Wednesday, January 30, 2013

What your mother never told you about cufflinks

In this post, I am trying to nail down the "best" way to run cufflinks, part of the so-called "tuxedo suite" of RNA-Seq tools (tophat, cufflinks, etc.)

I have a dataset composed of 3 technical replicates -- exactly the same library, run in 3 separate hi-seq lanes (multiplexed along with other samples, but I am ignoring those for now.) They were run on an Illumina HiSeq 2000, with 101 bp paired-end reads.

Number of reads, before and after quality trimming and preprocessing:

Replicate   Before      After
1           19,939,455  19,074,425
2           19,929,909  19,120,424
3           19,950,587  19,101,188
You can see, the lanes were very uniform in the number of reads.

For all of the following commands, I am using the iGenomes files (for the bowtie2 index, the GTF file, and the genome.fa file), available here: 
http://ccb.jhu.edu/software/tophat/igenomes.shtml

I ran tophat like this:
tophat -G genes.gtf --transcriptome-index=transcriptome -p 8 -o out genome rep1_1.fq.gz rep1_2.fq.gz

Cufflinks was run like this:
cufflinks -p 24 --no-effective-length-correction -b genome.fa -G genes.gtf -o out out/accepted_hits.bam

The two variables are highlighted.  

The cufflinks online manual says (note, "Cuffdiff" is a typo and should read "Cufflinks")


--no-effective-length-correctionCuffdiff will not employ its "effective" length normalization to transcript FPKM.

and:

-b/--frag-bias-correct Providing Cufflinks with the multifasta file your reads were mapped to via this option instructs it to run our bias detection and correction algorithm which can significantly improve accuracy of transcript abundance estimates. See How Cufflinks Works for more details.

For each run, I sorted the resulting genes.fpkm_tracking files, pulled out the FPKM column, and calculated the mean of the pairwise Pearson correlation for the 3 replicates. Here are the results:


effective length corr
fragment bias corr
Time
Correlation
Yes
Yes
43 min
0.9790755
Yes
No
16 min
0.9749819
No
Yes
43 min
0.9999465
No
No
16 min
0.9999626


Note: the effective length correction is on by default, whereas the frag bias correction is off by default.

The best way to run cufflinks is with --no-effective-length-correction !!!

frag-bias-correction makes the analysis much slower, without any discernable benefit in this trial run.

2 comments: