Cufflinks/Cuffdiff different FPKM values for multiple genes at one location
 

I am doing my first RNA-seq for Drosophila melanogaster (I normally deal with human or mouse data). It turns out there are a lot of fly genes that have identical coordinates as other genes. In other words, the same location has multiple genes assigned to it.

If there are multiple genes that have the same exact coordinates, they should have the same FPKM values. However, running Cuffdiff using a GTF like that does not yield the same values for all genes.

Is there a way to force Cuffdiff to assign the same values to all overlapping genes? I could not find any arguments that may do that. Is there a proper way of dealing with such situations? Do I need to optimize the GTF file? I use the one from iGenomes, which is endorsed by Cufflinks, so it seems like it should be fine.

Hi,

The FPKM values for the genes are a sum of the FPKM values found for transcripts of that gene... so I guess you could expect differing values come from the difference in transcripts and the reads/fragments covering them.

To test I would look at the transcripts FPKM values for each gene that have the same loci.

The FPKM values for multiple transcripts with identical coordinates are very different. Additonally, they can even be called significantly different between two samples and in different directions.

I stopped using cufflinks/cuffdiff 3 months ago as the latest version was producing implausible results. I would recommend using tophat2 + htseq-count + edgeR (or DESeq). I based my workflow on this nice tutorial: http://www-huber.embl.de/pub/pdf/nprot.2013.099.pdf

I agree with feralBiologist, but would switch to featureCounts instead of HTSeq-count for performance reasons (can be run multithreaded and does not require resorted SAM file).

Thanks for this suggestion. Have you checked whether HTSeq-count and featureCounts produce the same results? HTSeq-count supports certain counting models. Basically, the main issue is how to count reads that hit overlapping genes. I chose HTSeq-count because it is written by the author of DESeq and in the above paper has been "approved" also by the authors of edgeR. If featureCounts produces the same results than the switch would be painless but if there are differences than you need to look into the details.

EDIT: I realised that featureCounts is written by the authors of edgeR so it shall be straight-forward to substitute HTSeq-count. Thanks again to rboettcher.

I suggest to have a look a their manuscript on arXiv where the authors made such as comparison. From my experience both tools produce similar results, see http://arxiv.org/abs/1305.3347

EDIT: another nice feature is that fC outputs gene length, so computation of RPKM is straight forward.

I am going to try these new methods out, thanks.