Solution found: For Cuffdiff/links 2.0.2 Make cufflinks FPKM match Cuffdiff FPKM

[NGSfan]

This is in reference to cuffdiff v2.0.2 (3522)

Hi, I have been banging my head on why cuffdiff and cufflinks have different FPKM calculations. I finally stumbled upon the settings to make cufflinks match the FPKM values created with default settings of cuffdiff:

cuffdiff (defaults) mygenes.gtf mysample1.bam mysample1.bam
cufflinks --compatible-hits-norm --max-frag-multihits 1 -G mygenes.gtf mysample1.bam

I think there is a mistake in the cuffdiff options list, where --max-frag-multihits is listed as unlimited in the default state, but it is probably actually set to 1.

also, I add a pseudocount RPKM of 0.01 to both, otherwise cufflinks makes very small RPKMs, unlike cuffdiff, where it sets it to zero.



At least now it is possible for me to reproduce the FPKMs that cuffdiff makes on default settings. I hope the authors fix this.

Note the plot is on log2 scale for the RPKMs, but the match is exact when on the linear scale as well (with exception of a few - those genes with extremely small FPKMs ~ e-200 as calculated by cufflinks).

btw, I carefully compared all the options and saw that the options for --compatible-hits-norm and --total-hits-norm are flipped between the two programs:

CUFFLINKS: --compatible-hits-norm count hits compatible with reference RNAs only [ default: FALSE ]
CUFFDIFF: --compatible-hits-norm count hits compatible with reference RNAs only [ default: TRUE ]

CUFFLINKS: --total-hits-norm count all hits for normalization [ default: TRUE ]
CUFFDIFF: --total-hits-norm count all hits for normalization [ default: FALSE ]

Interestingly, there is another normalization only offered by Cuffdiff and set to true:

CUFFDIFF only: --geometric-norm use geometric mean normalization [ default: TRUE ]


Both Cufflinks and Cuffdiff say the default for max-frag-multihits is unlimited, but I have found setting it to 1 in Cufflinks brings the two into agreement:

CUFFLINKS --max-frag-multihits Maximum number of alignments allowed per fragment [ default: unlim ]
CUFFDIFF --max-frag-multihits Maximum number of alignments allowed per fragment [ default: unlim ]

[syfo]

Thanks for sharing

[NGSfan]

Alternatively, instead of pseudocounts to deal with some of the very small RPKMs generated by Cufflinks, you can also replace any RPKM < 0.0001 with 0 and this will give you the same values as reported by Cuffdiff.

I don't know why on some genes, the RPKMs get so small (1e-300 for example) in Cufflinks, whereas in Cuffdiff they are 0 RPKM.

Has anybody compared Cufflinks to eXpress? I notice that there are subsets of genes where cufflinks reports 1000 RPKM while eXpress reports 0 and vice versa... many look to be small RNAs.. just wondering what others are seeing out there. This happens with or without effective length correction.

I will try RSEM next.

I am a bit concerned though with aligning only to transcripts, wouldn't this result in some reads being forcibly aligned to a transcript reference? A lot of transcription goes on that is not annotated (as we know from the ENCODE project) and this could result in skewing RPKMs for some genes. Perhaps the genomic reference approach is better? With accurate aligners like STAR the argument for aligning solely to the transcriptome seems weaker.

[fuad193]

Very Useful I will try it now. Just wondering why you make multi-hit only 1, this seems odd because the aligner will stop for looking of other possibility in mapping when it encounter the first hit.

[NGSfan]

Very Useful I will try it now. Just wondering why you make multi-hit only 1, this seems odd because the aligner will stop for looking of other possibility in mapping when it encounter the first hit.

If my interpretation or understanding of the --max-frag-multihits 1 option is correct, this tells cuffdiff to limit itself to using reads that align uniquely to one locus, and to ignore reads that map to more than 1 locus. In the case of tophat/STAR, reads that are uniquely mapped get a score of 255, so it should be counting those only.

I have not played around with changing this value. I would assume that unique reads are better to use for differential gene expression.