[billstevens]

Thank you Simon.

To echo turnersd, it would be very, very nice to have a comparison to the new Cufflinks/Cuffdiff, but the new version did JUST come out. Have you had a chance to take a peak at it at all, Simon?

This article was very useful in that I can now safely ignore the results that I got with the old version of cufflinks. Cole has stated the new version is much, much more accurate. However, I just used Cufflinks and Cuffdiff 2.0.2, and I have found the results more than little confusing in regards to FPKM and values (and yes, I'm aware of the common scale transformation). I'm going to give the Tophat results to EasyRNASeq and then run it on DESeq, and then compare what I've got.

[turnersd]

Well, for our take on it, see the Discussion section in our DEXSeq paper (Anders, Reyes, Huber; Genome Res., 2012).

Annoyingly, Genome Research still has not managed to put our supplement online, which might be relevant; hence, please look at it in the in the PDF available here. (Same PDF as the preprint, at the end.)

Thanks Simon.

Based on many threads here regarding huge differences folks are seeing in the new version of cufflinks/cuffdiff, it would be nice to see the comparisons in tables S1 and S2 supplemented with results from Cufflinks 2.0.2.

[Simon Anders]

Sorry Simon, would you mind commenting on the new Cufflinks and Cuffdiff? I don't mean to put you on the spot, since there is clearly a competition aspect of this, but in this nascent field, there are so few experts that can comment authoritatively on this.

Well, for our take on it, see the Discussion section in our DEXSeq paper (Anders, Reyes, Huber; Genome Res., 2012).

Annoyingly, Genome Research still has not managed to put our supplement online, which might be relevant; hence, please look at it in the in the PDF available here. (Same PDF as the preprint, at the end.)

[turnersd]

Sorry Simon, would you mind commenting on the new Cufflinks and Cuffdiff? I don't mean to put you on the spot, since there is clearly a competition aspect of this, but in this nascent field, there are so few experts that can comment authoritatively on this.

+1.

My current strategy is to use cufflinks/cuffdiff for isoform level differences, DEXSeq for differential exon usage, and DESeq for differential expression at the gene level, then presenting all the results to the core's users while explaining the limitations, caveats, gotchas, and known reproducibility issues with each software/method.

Further, I've spoken with others recently about alternative workflows. E.g., tophat for alignment, cufflinks or trinity for assembly, RSEM or eXpress for getting counts for isoforms, then edgeR or DESeq for differential expression testing. There are several permutations of this workflow, and I'm hesitant to use any of them in my core until I see something published, but the ideas are worth entertaining.

[billstevens]

Sorry Simon, would you mind commenting on the new Cufflinks and Cuffdiff? I don't mean to put you on the spot, since there is clearly a competition aspect of this, but in this nascent field, there are so few experts that can comment authoritatively on this.

Apart from increased knowledge, I ask for my own selfish purposes because I've fallen victim to what turnersd has said: "I've heard far too many examples of people using tophat-cufflinks-cuffmerge-cuffdiff for their RNA-seq workflow, and getting extremely different results from one version of cufflinks to the next."

[billstevens]

Wow, DESeq must be very, very conservative because the new version of Cufflinks is incredibly conservative. The number of differentially expressed genes went from the high hundreds to just 5 or so when I upgraded to the new version of Cufflinks.

[DJParker]

Cufflinks2 and DeSeq

Hello

I have been working with the Cufflinks2 Vs. DESeq issue for a while now and thus decided to test how both packages performed with my own data.

Essentially DESeq appears to be more conservative producing 73 DE genes to 173 from cufflinks at a cutoff of p=0.05

However of the 73 genes detected by DESeq 55 of them also appeared in the list produced by Cufflinks.

Thus, it looks to me that the 2 methods are producing broadly the same result.

I am still unsure which of the results to 'trust' more, however, I think areyes point about combining the approaches is a good move as I now have a fair bit of confidence that in the genes that at DE in both lists really are DE.


ALSO:

lshen and natstreet

The Binary release of cufflinks 2.0.0 has an issue with cuffdiff where it assigns lots of zeros and other strange counts. This has been fixed for the binary release in 2.0.1 and later see: http://cufflinks.cbcb.umd.edu/

[natstreet]

Cuffdiff zero FPKM with reference GFF

I have the exact same problem of having zero FPKM values for all genes that have a reference transcript. In my case, the GFF reference file comes from phytozome. I have run it through gffread to be sure it is being interpreted correctly.

If I run the samples through as individual samples rather than biological replicates in groups then things seem OK. Using previous versions of cufflinks/diff this worked fine.

My process was tophat2-->cufflinks2-->cuffmerge-->cuffdiff

I'm using version v2.0.0 of cufflinks and tophat-2.0.0.Linux_x86_64

[lshen]

I got another issue with new CUFFLINK 2:

When I directly quantify against ensembl gtf, the cufflinks returned 0 expression for most of them. This only occurred when I used replicates. single sample group is fine. And seems only when transcripts matched to known gene's annotation.


#command:
cufflinks-2.0.0.Linux_x86_64/cuffdiff -p 8 -L P1,P2 -c 1 -b anFam2.fa -o cuffdiff.P1.P2.ensembl canFam2.67.gtf TOPHAT2.C1.bam,TOPHAT2.C2.bam,TOPHAT2.C3.bam,TOPHAT2.C4.bam TOPHAT2.C5.bam,TOPHAT2.C6.bam,TOPHAT2.C7.bam,TOPHAT2.C8.bam,TOPHAT2.C9.bam

Here are the number of genes returned FPKM 0 in cuffdiff:

$8 is for treatment P1, $9 is for treatment P2 in output.

awk ' $8 ==0 { i++}; END {print i " of " NR " = " i/NR*100 "%"} ' cuffdiff.P1.P2.ensembl/gene_exp.diff

cufflinks-2.0.0:

P1: 24649 of 24661 = 99.9513%
P2: 24645 of 24661 = 99.9351%

cufflinks 1.3.0 seems right:

P1: 6345 of 24661 = 25.7289%
P2: 6564 of 24661 = 26.6169%

[epi]

The thread you linked to doesn't show a bug, just an incorrect use of Cufflinks, as I just noted in that thread.

Hi Cole, as reflected in the thread, it seems this not an incorrect use.

From the cufflinks manual: -g option

Output will include all reference transcripts as well as any novel genes and isoforms that are assembled.

I was wondering if there is any solution or resolution, as many others have reported the same issue from this thread

http://seqanswers.com/forums/showthread.php?t=18070

.

I also reported another issue in the same thread, can you please look and comment on it.

Hi Cole, Thanks again for replying. While your comment clarifies some things, it also raises some more question. First here is my command
I am still unclear why there are new predictions (with CUFF. IDs) that are identical to existing annotations. I have manually looked at one of the test runs and all new predictions (with CUFF. IDs) were original genes/transcripts in the source GTF annotation (I must have looked only about 20-25 ). So far I have not encountered any novel predictions which do not correspond to existing annotation. In some cases, original annotation still exists and in some cases not. Similar to the example about Trib3

[Xi Wang]

Hi Simon, Cole & Lior,

Glad to find you'er discussing the two main strategies in analysing RNA-seq data for differential expression. I think either exon-centric or isoform-centric strategy has its own pros and cons, and both perform reasonably well in real-world RNA-seq studies. What I concern is, however, which strategy is more preferred in paired-end (PE) RNA-sequencing. PE design will become the unique in RNA-seq, as PE is widely recognised better for identifying splice variants in transcriptomes. The advantage of PE can be easily and fully modelled in isoform-centric methods including Cufflinks and RSEM, but how exon-centric methods can take PE into account? Any idea? Thanks.

[gringer]

A good genomic reference sequence also helps. Isoforms are much easier to find and identify when you know what the genomic sequence looks like.

[areyes]

I think the combination of several methods its the way to go, depending on the question one may ask. I think a single method won't answer all questions. At the end all the methods have advantages and limitations. (e.g. edgeR and DESeq have very good control over in the false discovery rate aspect compared other methods)

Regarding the last point of longer/shorter isoforms as possible confounders of differential expression, maybe one could also look at differences in exon usage. e.g. DEXSeq would detect the exons that are differentially used in between the shorter (condition A) and longer isoform (condition B), independently if the gene is differentially expressed between A and B.

[lpachter]

I'd like to point out that Cole's example is not only showing that DESeq and edgeR's assumption of features being the same length is problematic because it might lead to a canceling out that leads to non-DE calls when they should be made. In fact, the reverse can happen as well. In the example, suppose that in condition B you actually have 200 reads (on the longer isoform), instead of 100 on the shorter. In this case, the 2 fold change in read count is certain to look like differential expression, but in reality, there has been no change in expression at all. So its not only that you might get a few splicing changes among your expression changes with DESeq or edgeR. You will also miss significant expression changes.

The numbers of counts in Cole's example are made up (to make the point), but the differences in transcript length are not. Its common for transcripts in a gene to have variable lengths, even much more than the factor of 2 in the example. In DESeq/edgeR every time that happens in the structure of a gene, the computed fold change in expression, which is the basis for calling differential expression, is incorrect.

In summary, I agree with Simon that DESeq and edgeR are useful tools for analysis of differential expression in single isoform genes where all reads map uniquely.

[Simon Anders]

Well, Cole is certainly write with both that this is an issue and that it is a contrived example. In practice, it will be rare that a change in isoform usage and a change in overall expression cancel out that well, and hence I would not worry about it. It is true, however, that DESeq and edgeR cannot distinguish changes in number of transcripts from changes in the average transcript length, and hence, if the tools report significant evidence for differential expression, this should be read as "significant evidence for differential expression and/or differential isoform usage". In many cases, you will be able to live with that -- either because you intend to have a closer look at the top hits anyway and would then discover this by manual inspection, or because you just want an overview and don't mind if a few splicing changes are among your expression changes.

And if you are interested in splicing changes, have a look at DEXSeq, too.