strange problem for cufflinks

sterding · 12-02-2010, 09:56 PM

HI, all

I am facing a strange problem for cufflinks: Reads that are aligned to exons obviously can NOT be counted by cufflinks. Please see if you have the same problem when running the following test.

A simple test SAM file (htt.sam)

Code:

HWUSI-EAS623:1:91:16574:8737#0	16	chr4	3129018	255	40M	*	0	0	AGCTGGCTGCTTCTTCAGGGGTTTCCACTCCAGGGTCAGC	@BBBCE-EEEDGGGDDFGEFGFGAGGGGGEGGGGGGGGDF	NM:i:0	NH:i:1	XS:A:-
HWUSI-EAS623:1:5:14642:7499#0	0	chr4	3129040	255	40M	*	0	0	TTCCACTCCAGGGTCAGCAGGTCATGACATCATCACAGAA	GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFF	NM:i:0	NH:i:1	XS:A:-
HWUSI-EAS623:1:27:15057:20736#0	16	chr4	3129048	255	40M	*	0	0	CAGGGTCAGCAGGTCATGACATCATCACAGAACAGCCACG	EAEEEEEEBEEEGEGEEEBEEEEEEFEGBGGEGFGGGGGG	NM:i:0	NH:i:1	XS:A:-
HWUSI-EAS623:1:117:13132:3189#0	16	chr4	3129072	255	40M	*	0	0	TCACATAACAGCCACGGTCACAGCACACACTGCACGCGGA	@B@CB7EEEEGGGEGGGGGGGGGGGECECECDDD*BBBBB	NM:i:2	NH:i:1	XS:A:-

Here is the GTF file (for test purpose, only transcript ENST00000355072 of gene HTT):
htt_ENST00000355072.gtf

The reads in SAM file are obviously resided in the exonic region of GTF file (see UCSC screenshot)

However, running cufflinks command

Code:

cufflinks --GTF htt_ENST00000355072.gtf htt.sam

will only give 0 FPKM value in the output file (genes.expr):

Code:

$ cat genes.expr 
gene_id	bundle_id	chr	left	right	FPKM	FPKM_conf_lo	FPKM_conf_hi	status
ENSG00000197386	3	chr4	3076406	3245676	0	0	0	OK

I found this problem by double-checking the cufflinks results of a genome-wide RNA-seq analysis; HTT gene has well-aligned reads, but its FPKM is zero. This is not the only case. I am simplifying the problem here by only showing the part of SAM file and GTF file.

Thanks for helps ahead.

plabaj · 12-03-2010, 02:23 AM

I had somehow similar problem.

I had few thousands reads which fall on specific splice junction belonging to just one transcript.
However, as mentioned above, Cufflinks was blind to them. The FPKM value was suspiciously close to zero for this transcript.

I reported this to Cole Trapnell few weeks ago. So far no response

If somebody know what is the problem I would appreciate such knowledge.

Thomas Doktor · 12-03-2010, 02:37 AM

I'm not certain, but I think this is because Cufflinks considers 4 reads in total within a gene as below detection threshold. In cuffdiff, you can adjust this by the "-c" options, but cufflinks doesn't have this, probably because it needs a certain amount of reads to assemble transcripts and doesn't simply count reads within a gene.

Try running Cuffdiff with "-c 2", it should be able to produce an FPKM value then (which will be very close to zero in any case if there is only 4 reads in total). The default -c value is 1000, btw.

sterding · 12-03-2010, 09:44 AM

Thanks a lot, Tom

I am re-running cuffdiff using -c option as you suggested, see what FPKM it will get.

meanwhile, I checked the problematic genes; there are reasonable numbers of reads mapped with its exons (see below). Actually the signal is quite clear for me. I wish Tole can open an option for minimal number of reads in cufflinks.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Strange maq 7.1 problem	homeriq5	Bioinformatics	7	11-26-2013 07:19 AM
Cufflinks / Cuffdiff problem	Morten	Bioinformatics	42	07-13-2013 10:51 PM
Problem with Cufflinks output	admiral	RNA Sequencing	0	06-01-2011 01:19 PM
strange bowtie index building and mapping problem	Gangcai	Bioinformatics	0	08-04-2010 05:02 PM
Problem with Tophat and Cufflinks	gongchang	RNA Sequencing	2	06-08-2010 04:16 PM

12-02-2010, 09:56 PM	#1
sterding Member Location: Boston Join Date: Sep 2010 Posts: 36	strange problem for cufflinks HI, all I am facing a strange problem for cufflinks: Reads that are aligned to exons obviously can NOT be counted by cufflinks. Please see if you have the same problem when running the following test. A simple test SAM file (htt.sam) Code: HWUSI-EAS623:1:91:16574:8737#0 16 chr4 3129018 255 40M * 0 0 AGCTGGCTGCTTCTTCAGGGGTTTCCACTCCAGGGTCAGC @BBBCE-EEEDGGGDDFGEFGFGAGGGGGEGGGGGGGGDF NM:i:0 NH:i:1 XS:A:- HWUSI-EAS623:1:5:14642:7499#0 0 chr4 3129040 255 40M * 0 0 TTCCACTCCAGGGTCAGCAGGTCATGACATCATCACAGAA GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFF NM:i:0 NH:i:1 XS:A:- HWUSI-EAS623:1:27:15057:20736#0 16 chr4 3129048 255 40M * 0 0 CAGGGTCAGCAGGTCATGACATCATCACAGAACAGCCACG EAEEEEEEBEEEGEGEEEBEEEEEEFEGBGGEGFGGGGGG NM:i:0 NH:i:1 XS:A:- HWUSI-EAS623:1:117:13132:3189#0 16 chr4 3129072 255 40M * 0 0 TCACATAACAGCCACGGTCACAGCACACACTGCACGCGGA @B@CB7EEEEGGGEGGGGGGGGGGGECECECDDDBBBBB NM:i:2 NH:i:1 XS:A:- Here is the GTF file (for test purpose, only transcript ENST00000355072 of gene HTT): htt_ENST00000355072.gtf The reads in SAM file are obviously resided in the exonic region of GTF file (see UCSC screenshot) However, running cufflinks command Code: cufflinks --GTF htt_ENST00000355072.gtf htt.sam will only give 0 FPKM value in the output file (genes.expr): Code: $ cat genes.expr gene_id bundle_id chr left right FPKM FPKM_conf_lo FPKM_conf_hi status ENSG00000197386 3 chr4 3076406 3245676 0 0 0 OK I found this problem by double-checking the cufflinks results of a genome-wide RNA-seq analysis; HTT gene has well-aligned reads, but its FPKM is zero. This is not the only case. I am simplifying the problem here by only showing the part of SAM file and GTF file. Thanks for helps ahead. Last edited by sterding; 12-02-2010 at 10:03 PM.*

12-03-2010, 02:23 AM	#2
plabaj Member Location: Vienna Join Date: Oct 2010 Posts: 86	I had somehow similar problem. I had few thousands reads which fall on specific splice junction belonging to just one transcript. However, as mentioned above, Cufflinks was blind to them. The FPKM value was suspiciously close to zero for this transcript. I reported this to Cole Trapnell few weeks ago. So far no response If somebody know what is the problem I would appreciate such knowledge. __________________ Pawel Labaj

12-03-2010, 02:37 AM	#3
Thomas Doktor Senior Member Location: University of Southern Denmark (SDU), Denmark Join Date: Apr 2009 Posts: 105	I'm not certain, but I think this is because Cufflinks considers 4 reads in total within a gene as below detection threshold. In cuffdiff, you can adjust this by the "-c" options, but cufflinks doesn't have this, probably because it needs a certain amount of reads to assemble transcripts and doesn't simply count reads within a gene. Try running Cuffdiff with "-c 2", it should be able to produce an FPKM value then (which will be very close to zero in any case if there is only 4 reads in total). The default -c value is 1000, btw.

12-03-2010, 09:44 AM	#4
sterding Member Location: Boston Join Date: Sep 2010 Posts: 36	Thanks a lot, Tom I am re-running cuffdiff using -c option as you suggested, see what FPKM it will get. meanwhile, I checked the problematic genes; there are reasonable numbers of reads mapped with its exons (see below). Actually the signal is quite clear for me. I wish Tole can open an option for minimal number of reads in cufflinks.