Test status in Cufflinks and Cuffdiff

[qiongyi]

Hi, I ran Cufflinks and Cuffdiff on some RNA-Seq datasets, and I got the result files "isoform_exp.diff" and "gene_exp.diff". When I looked at these two list, I found something interesting. Here I take gene "Phactr1" as an example:

gene level estimation:

test_id	gene_id	gene	locus	sample_1	sample_2	status	value_1	value_2	log2(fold_change)	test_stat	p_value	q_value	significant
Phactr1	Phactr1	Phactr1	chr13:42680622-43138512	KB	NB	OK	1.0429	1.49812	0.522548	-1.36759	0.171441	0.677408	no

isoform level estimation:

test_id	gene_id		locus	sample_1	sample_2	status	value_1	value_2	log2(fold_change)	test_stat	p_value	q_value	significant
NM_001005740	Phactr1	Phactr1	chr13:42680622-43138512	KB	NB	NOTEST	1.0429	1.49812	0.522548	-1.36759	0.171441	1	no
NM_001005748	Phactr1	Phactr1	chr13:42680622-43138512	KB	NB	NOTEST	0	0	0	0	1	1	no
NM_198419	Phactr1	Phactr1	chr13:42680622-43138512	KB	NB	NOTEST	0	0	0	0	1	1	no

Gene "Phactr1" has three isoforms (NM_001005740, NM_001005748, NM_198419).
Why the gene is given status "OK" while the isoform "NM_001005740" is given "NOTEST". It seems strange that cufflinks/cuffdiff extends different treatment to isoform and gene with the same expression value.
I wonder what the criteria to set the status tags for genes and isoforms? Different criteria?

[qiongyi]

Hope someone here could help to clarify this...

[drosoform]

I was wondering the exact same thing in my RNA-seq data from Cuffdiff. If anyone could explain why this happens, it would be very helpful!

[bhootnaath]

Here is a somewhat related behavior of Cuffdiff / cummerbund that gives erroneous results.

A) Ran Tophat v2.0.4 for 9 samples (3 conditions c1, c2, and c3 each with 3 biological replicates _1, _2 and _3), for example for c1_1:

tophat -p 8 -G genes.gtf -o tophat.output.c1_1 --no-novel-juncs --solexa1.3-quals bowtie_index c1_1.fq

Repeated the above for c1_2, c1_3, c2_1, c2_2, c2_3, c3_1, c3_2, and c3_3.

B) Ran Cuffdiff v2.0.2 (3522) for differential expression analysis for known list of genes/junctions:

cuffdiff -p 8 -b genome.fa -u genes.gtf -o cuffdiff.out -L c1,c2,c3 tophat.output.c1_1/accepted_hits.bam,tophat.output.c1_2/accepted_hits.bam,tophat.output.c1_3/accepted_hits.bam tophat.output.c2_1/accepted_hits.bam,tophat.output.c2_2/accepted_hits.bam,tophat.output.c2_3/accepted_hits.bam tophat.output.c3_1/accepted_hits.bam,tophat.output.c3_2/accepted_hits.bam,tophat.output.c3_3/accepted_hits.bam

C) Then issued the following commands in R/cummerbund (v. 1.99.2) to write an output file containing differentially expressed genes for c1 vs. c2 with alpha = 0.05:

> cuff <- readCufflinks()
> C1vsC2.sigGeneIds<-getSig(cuff,"c1","c2",alpha=0.05,level="genes")
> C1vsC2.sigGenes <- getGenes(cuff,C1vsC2.sigGeneIds)
> write.table(diffData(C1vsC2.sigGenes),"C1vsC2sigGenes.diff")

Now here is one of many such erroneously reported significant genes in the C1vsC2sigGenes.diff file. The following gene *fails* the alpha = 0.05 threshold for c1 vs c2 and yet gets reported!

gene_id	sample_1	sample_2	status	value_1	value_2	log2_fold_change	test_stat	p_value	q_value	significant
Gene_X	c1	c2	OK	1.12696	0.0921556	-3.61222	3.44126	0.000579016	[COLOR="Red"]0.0751255[/COLOR]	[COLOR="Red"]no[/COLOR]
Gene_X	c1	c3	OK	1.12696	0.264807	-2.08942	1.30123	0.19318	0.999998	no
Gene_X	c2	c3	OK	0.0921556	0.264807	1.5228	-0.990653	0.321855	0.999998	no

Added later:

This phenomenon has possibly been explained here and here.

[sdriscoll]

I've seen similar inconsistencies. I've had a gene knocked out and it was not found significant at the gene level but at the isoform level it was. The thing is at the isoform level there was only one isoform with expression and its expression was the same as the gene level expression. Seems like a flaw in gene's estimated variance verses the isoform variance. Whatever the case it makes no biological sense which is always gonna be the point.