6.3 Transcriptome reconstruction and expression using cufflinks

There are a number of tools that perform reconstruction of the transcriptome and for this workshop you are going to use cufflinks, which can do transcriptome assembly with and without a reference annotation. It also quantifies the isoform expression in FPKMs. Let’s import the cufflinks module and check its (rather long) list of parameters:

$ module load cufflinks/2.1.1
$ cufflinks --help

You are going to use cufflinks with the UCSC annotation file for zebrafish limiting the transcriptome reconstruction strictly to the available annotations. You could also use the annotations as a guide but letting the algorithm look also for transcribed loci not included in the annotations. In the first case cufflinks will only report isoforms that are included in the annotation, while in the latter case it will report novel isoforms as well. First, copy the annotation file from UCSC for chromosome 12 of Danio rerio:

$ cp –r /ibers/repository/public/courses/Rna-Seq/annotations .  

Then let’s create a folder for the cufflinks output storage:

$ mkdir cuff_out

Now you are ready to run cufflinks. The general format of the cufflinks command is:

cufflinks [options]* <aligned_reads.(sam/bam)>

where the input is the aligned reads (either in SAM or BAM format). Prepare a script file to launch cufflinks, by opening it with a text editor (you can call it for example cufflinks_2cells.sh). Copy the same header for the scheduler that you used previously, use the commands to load the required modules, and point the script to your working space. Your script should look something like this:

#$ -S /bin/sh
#$ -cwd
#$ -q amd.q,large.q,intel.q
#$ -l h_vmem=24G
#$ -e cff_12.e
#$ -N cff
#$ -o cff_12.o
module load cufflinks/2.2.1
cd /ibers/ernie/home/vpl/zebrafish

Then you can write (all in one line) the command:

cufflinks -o cuff_out/2cells -G annotations/Zebrafish_refGene_chr12.gtf \
-b /ibers/ernie/home/vpl/zebrafish/bt2Index/chr12.fa \
-u --library-type fr-unstranded /ibers/ernie/scratch/vpl/zebra_fish/tophat_out_2cells/accepted_hits.bam

The options that you used above are: -o: output directory

-G: tells cufflinks to use the supplied annotation in order to estimate isoform annotation.

-b: instructs cufflinks to run a bias detection and correction algorithm which can significantly improve accuracy of transcript abundance estimates (the genome sequence is required to perform this).

-u: tells cufflinks to do an initial estimation procedure to more accurately weight reads mapping to multiple locations in the genome.

--library-type: specify the employed sequencing strategy

Save, close and run using qsub. With your small dataset, it should not take more than few minutes. While waiting, you can prepare and run a similar file for the 6-hours embryo sample (again, remember to change the output folder name, otherwise the second run will overwrite the first run files). You should know now how to do that, so try by yourself before looking at the command below:

cufflinks -o cuff_out/6hours -G annotations/Zebrafish_refGene_chr12.gtf \
-b /ibers/ernie/home/vpl/zebrafish/bt2Index/chr12.fa \
-u --library-type fr-unstranded /ibers/ernie/scratch/vpl/zebra_fish/tophat_out_6hours/accepted_hits.bam

When both have finished, you can take a look at the output folders that have been created. The results from cufflinks are stored in 4 different files: genes.fpkm_tracking: This file contains the estimated geneAlevel expression values in the generic FPKM Tracking Format.

isoforms.fpkm_tracking: This file contains the estimated isoform level expression values in the generic FPKM Tracking Format.

transcripts.gtf: This GTF file contains the assembled isoforms.

skipped.gtf: This GTF file contains the skipped loci (it’s often empty)

Let’s have a look at the file genes.fpkm_tracking, where gene expression is stored. The first column reports the UCSC identifier of the gene, the fifth its common name, the seventh its genomic coordinates, and the tenth its expression in FPKM. You would like to know which are the most expressed genes in the sample. This can be done in many ways, for example opening the file as a spreadsheet with Excel, but let’s do it using Unix commands. You can use the sort command, which can order a file based on the values in a specific column:

$ sort -t$'\t' -r -g -k 10 cuff_out/2cells/genes.fpkm_tracking > cuff_out/2cells/genes.sorted.fpkm_tracking

If you want to better understand the command that you just typed, check the manual page for the sort command using:

$ man sort

The first parameter of sort tells the command that the file is divided in columns separated by tabs; the –k parameter specify which column you are interested in; the –g tells sort that the column contains numbers (the default sort behavior is to order by lexicographic order); the –r tells sort to order from the largest to the smallest value. The sort output is redirected into a new file. If you open this new file, the genes are now ordered by decreasing expression. If you now want to know, for example, which is the most expressed gen, you can copy UCSC identifier in the first row of the sorted file, the UCSC genome browser to get more information about it . Do that for a few genes.

Now do the same also for the 6-hours embryo sample. Which are the most expressed genes in this other sample? Are they in common with the 2-cells embryo?

Now that you have all the strongly expressed genes in the two embryos, you might wonder in which biological processes they participate. Let’s see how you can use an automatic web-based functional annotation tool (DAVID) to help you in the rationalization of which are the most represented process in a biological sample. To use DAVID you need to submit a list of gene identifiers. In the genes.sorted.fpkm_tracking, the first column contains the UCSC gene ids of the genes in the annotation GTF file. You need to extract these identifiers from this file. Let’s take the top 100 of the genes sorted by FPKM. Let’s use the command head to select the top 100 genes, and the command cut to extract only the first column from these rows (if you want, use man to have more details on how these commands work):

$ head -100 cuff_out/2cells/genes.sorted.fpkm_tracking | cut -f 1 > cuff_out/2cells/top100byFpkm.txt

Remember that the > operator redirects the output of a command into a file. Now open the top100byFpkm.txt file, verify that it’s a single-column file with the UCSC gene identifiers, and copy its content. Open any browser and go to: [1]