strange problem for cufflinks

[sterding]

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)

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

cufflinks --GTF htt_ENST00000355072.gtf htt.sam

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

$ 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]

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]

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]

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.