cufflinks - having difficulty reading reference GTF

tguo · 09-30-2010, 10:56 PM

Hi all,

I'm using cufflinks 0.9.0 beta. But the problem persists when I use v0.8.4.

Everything seems to work fine if I don't use a reference GTF.
But when I do use a reference (downloaded from the Ensembl FTP "current GTF" folder, Mouse release 59), cuffcompare doesn't seem to read from it.
(I did manually delete the giant transcript ENSMUST00000127664 from the GTF file to preclude the error message.)

When I run:
cuffcompare -r ../Mus_musculus.NCBIM37.56.gtf A.gtf B.gtf
where A.gtf and B.gtf are the two "transcripts.gtf" generated by cufflinks

the resulting "combined.gtf" file from cuffcompare does not contain any gene ID or gene name from the mouse reference GTF, and the refmap file is completely empty, too, besides the header row.

Does anyone have any idea what is going on?

Another weird thing is, if I run cufflinks with
"cufflinks -G ../Mus_musculus.NCBIM37.56.gtf A.sam",
then cufflinks does read the reference GTF, but in the output "transcripts.gtf", the estimated transcript levels are 0 for all transcripts...
However, if I run cufflinks without "-G", then the level estimation is correct.

Any ideas or help are highly appreciated!

TG

ldong · 10-01-2010, 07:27 AM

Hi, tguo,
It looks like you don't have enough reads as input. Please check how many reads are in your fastq file? To test cufflinks, I first used 1000 reads and got the same problem as yours. Later on, I used 10000 reads, then cuffcompare produces transcript expression level. But the problem is that I don't know how to get the gene expression level. Do you have any idea on this?
Best,
ldong

tguo · 10-01-2010, 08:21 AM

Hi ldong,

I guess gene expression levels are in the "genes.expr" output file.

I didn't think too few reads were the reason causing the problem. I was using a 10M read RNA-seq result as an input.. the SAM input into cufflinks were several gigabytes..

Thanks,

Ting

ldong · 10-01-2010, 09:27 AM

New version of tophat works!

Pejman · 10-07-2010, 11:33 AM

Quote:

Originally Posted by tguo

Hi all,

Everything seems to work fine if I don't use a reference GTF.

I have the same problem. I'm using cufflinks-0.9.1.Linux_x86_64 and trying to get it to produce expression levels based on a GTF file. If I run it without a reference GTF file, it works fine, gives some expressions. But with reference:

Code:

$cufflinks -G hg18_annotation.gtf MyFile.sam

it gives:

Code:

[20:24:17] Inspecting reads and determining fragment length distribution.
> Processing Locus chr6:27032750-27099732      [*****************        ]  69%
Error: this SAM file doesn't appear to be correctly sorted!
        current hit is at chr10:272501, last one was at chr9:139953334
Cufflinks requires that if your file has SQ records in
the SAM header that they appear in the same order as the chromosomes names 
in the alignments.
If there are no SQ records in the header, or if the header is missing,
the alignments must be sorted lexicographically by chromsome
name and by position.

The GTF file I downloaded from UCSC genomebrowser portal and .sam file is produced by Bowtie, and then sorted by Samtools. Does anybody have any clue?

ldong · 10-07-2010, 12:35 PM

Hi, Pejman,
I got the same error. To walk around it I use the gtf file when I call cuffcompare. See our wiki site for detail: https://sites.google.com/a/brown.edu...nks-and-tophat
Best,
ldong

adarob · 10-07-2010, 02:13 PM

Pejman,

Try remapping with the newest version of TopHat and your problem should disappear. Alternatively, you can add a proper header to your SAM file.

-Adam

Pejman · 10-08-2010, 03:44 AM

@ ldong Thanks for the nice link!

I managed to get around the problem, by manually adding the SQ headers to the sorted same file! They were dumped out during the SAM->BAM->sort->SAM process

But about TopHAt, I have tried to use it, so far have not been so successful, it gives errors. I'm using SOLiD data.

Pejman · 10-08-2010, 05:33 AM

@Adam:
I'm not interested in new splice junctions, just looking for differential expression of the known transcripts using Cufflinks. However, I kinda think that it would be still better to use tophat (without new splice junction prediction) for the alignment part than bowtie, in order to get more accurate results on the genes that have isoforms. What do you think as a Cufflinker?

ldong · 10-08-2010, 08:29 AM

Hi adarob,
Thank you for pointing to the right direction. I downloaded new version tophat binary. Now it can read gtf file. I assume in this way, tophat will be more accurate. But when I compare the output file accepted_hit.sam between using and not using gtf on a small data set, there is no difference. By the way, the new tophat does not output sam, I have to use samtools to convert the file:
samtools view -h -o accepted_hits.sam accepted_hits.bam
Any commnents on this? Best, ldong

adarob · 10-08-2010, 08:33 AM

Idong,

I don't work on TopHat so I can't answer your question about the GTF. You can contact Daehwan from the TopHat page. As for the bam output, conversion is not necessary since Cufflinks 0.9.x takes bam as input.

-Adam

ldong · 10-08-2010, 09:52 AM

Hi, Adam,
Thanks. I didn't realize new cufflnks can read bam. Now I am trying to see why UCSC genome browser can not read the bam file. Best, ldong

tguo · 10-09-2010, 07:46 AM

I think I've figured out what the problem was...

I have used a UCSC mm9 reference assembly against an EMSEMBL GTF, and that apparently was the problem. Some attributes of the two apparently doesn't match (didn't look in to see exactly why).

If instead, I use a UCSC GTF for mm9 reference, or conversely, use the EMSEMBL m_musculus_ncbi37 reference for EMSEMBL GTF, it seems to work.

ldong · 10-09-2010, 08:26 AM

Hi, tguo,
I think the reason is that reference sequence ID in mm9 does not match the chromosome Id in the GTF file.
By the way, where did you download the UCSC GTF file? Could you remind me? I could not find it anywhere. I use a GFF file from AceView.
Thanks,
ldong

tguo · 10-09-2010, 02:27 PM

Sure, from the UCSC table browser:

http://genome.ucsc.edu/cgi-bin/hgTab...a_doMainPage=1

adarob · 10-10-2010, 07:44 AM

Idong,

Did you use the --no-novel-juncs command?

tguo · 10-10-2010, 09:18 AM

No, I didn't. Would that help?

adarob · 10-10-2010, 10:23 AM

I'm sorry I actually should have addressed that to Pejman.

The idea is that if you're supplying a GTF in TopHat and you don't know want it find any splice junctions that are not in the GTF, you should use the --no-novel-juncs option.

Pejman · 10-11-2010, 05:13 AM

Hi, thanks for the hint. Yeah I know that, but apparently TopHat+colorspace has just been released in too much rush. It crashes on SOLiD data:
http://seqanswers.com/forums/showthread.php?t=7118
http://seqanswers.com/forums/showthread.php?p=26692

nkwuji · 10-28-2010, 04:51 AM

Hi,

In the tophat wiki page "https://sites.google.com/a/brown.edu/bioinformatics-in-biomed/cufflinks-and-tophat#read", I found some very helpful perl scripts, for example, find_known_gene_and_alternative_splicing_expression_from_cuffcompare_result.pl. Do you konw where I can download them?

Cheers,
Jun

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Difficulty aligning to reference with AVA	cog	Bioinformatics	8	07-07-2019 07:16 AM
cufflinks : analysis comparison with and without a gtf reference file	sohnic	Bioinformatics	3	07-07-2019 06:40 AM
GTF reference files that work with TopHat/Cufflinks	marcora	Bioinformatics	23	01-15-2014 12:10 AM
Cufflinks' computation of FPKM for --GTF and --GTF-guide estimation	burt	Bioinformatics	0	08-24-2011 12:59 AM
Difficulty mapping reads with non-reference allele?	krobison	Genomic Resequencing	3	10-09-2009 11:48 AM

09-30-2010, 10:56 PM	#1
tguo Junior Member Location: maryland Join Date: Sep 2010 Posts: 5	cufflinks - having difficulty reading reference GTF Hi all, I'm using cufflinks 0.9.0 beta. But the problem persists when I use v0.8.4. Everything seems to work fine if I don't use a reference GTF. But when I do use a reference (downloaded from the Ensembl FTP "current GTF" folder, Mouse release 59), cuffcompare doesn't seem to read from it. (I did manually delete the giant transcript ENSMUST00000127664 from the GTF file to preclude the error message.) When I run: cuffcompare -r ../Mus_musculus.NCBIM37.56.gtf A.gtf B.gtf where A.gtf and B.gtf are the two "transcripts.gtf" generated by cufflinks the resulting "combined.gtf" file from cuffcompare does not contain any gene ID or gene name from the mouse reference GTF, and the refmap file is completely empty, too, besides the header row. Does anyone have any idea what is going on? Another weird thing is, if I run cufflinks with "cufflinks -G ../Mus_musculus.NCBIM37.56.gtf A.sam", then cufflinks does read the reference GTF, but in the output "transcripts.gtf", the estimated transcript levels are 0 for all transcripts... However, if I run cufflinks without "-G", then the level estimation is correct. Any ideas or help are highly appreciated! TG Last edited by tguo; 09-30-2010 at 11:03 PM.

10-01-2010, 09:27 AM	#4
ldong Member Location: USA Join Date: May 2010 Posts: 15	New version of tophat works! Last edited by ldong; 10-08-2010 at 08:49 AM.

10-01-2010, 07:27 AM	#2
ldong Member Location: USA Join Date: May 2010 Posts: 15	Hi, tguo, It looks like you don't have enough reads as input. Please check how many reads are in your fastq file? To test cufflinks, I first used 1000 reads and got the same problem as yours. Later on, I used 10000 reads, then cuffcompare produces transcript expression level. But the problem is that I don't know how to get the gene expression level. Do you have any idea on this? Best, ldong

10-01-2010, 08:21 AM	#3
tguo Junior Member Location: maryland Join Date: Sep 2010 Posts: 5	Hi ldong, I guess gene expression levels are in the "genes.expr" output file. I didn't think too few reads were the reason causing the problem. I was using a 10M read RNA-seq result as an input.. the SAM input into cufflinks were several gigabytes.. Thanks, Ting

10-07-2010, 12:35 PM	#6
ldong Member Location: USA Join Date: May 2010 Posts: 15	Hi, Pejman, I got the same error. To walk around it I use the gtf file when I call cuffcompare. See our wiki site for detail: https://sites.google.com/a/brown.edu...nks-and-tophat Best, ldong

10-07-2010, 02:13 PM	#7
adarob Member Location: Berkeley, CA Join Date: Jul 2010 Posts: 71	Pejman, Try remapping with the newest version of TopHat and your problem should disappear. Alternatively, you can add a proper header to your SAM file. -Adam

10-08-2010, 03:44 AM	#8
Pejman Member Location: Switzerland Join Date: Jul 2010 Posts: 23	@ ldong Thanks for the nice link! I managed to get around the problem, by manually adding the SQ headers to the sorted same file! They were dumped out during the SAM->BAM->sort->SAM process But about TopHAt, I have tried to use it, so far have not been so successful, it gives errors. I'm using SOLiD data.

10-08-2010, 05:33 AM	#9
Pejman Member Location: Switzerland Join Date: Jul 2010 Posts: 23	@Adam: I'm not interested in new splice junctions, just looking for differential expression of the known transcripts using Cufflinks. However, I kinda think that it would be still better to use tophat (without new splice junction prediction) for the alignment part than bowtie, in order to get more accurate results on the genes that have isoforms. What do you think as a Cufflinker?

10-08-2010, 08:29 AM	#10
ldong Member Location: USA Join Date: May 2010 Posts: 15	Hi adarob, Thank you for pointing to the right direction. I downloaded new version tophat binary. Now it can read gtf file. I assume in this way, tophat will be more accurate. But when I compare the output file accepted_hit.sam between using and not using gtf on a small data set, there is no difference. By the way, the new tophat does not output sam, I have to use samtools to convert the file: samtools view -h -o accepted_hits.sam accepted_hits.bam Any commnents on this? Best, ldong

10-08-2010, 08:33 AM	#11
adarob Member Location: Berkeley, CA Join Date: Jul 2010 Posts: 71	Idong, I don't work on TopHat so I can't answer your question about the GTF. You can contact Daehwan from the TopHat page. As for the bam output, conversion is not necessary since Cufflinks 0.9.x takes bam as input. -Adam

10-08-2010, 09:52 AM	#12
ldong Member Location: USA Join Date: May 2010 Posts: 15	Hi, Adam, Thanks. I didn't realize new cufflnks can read bam. Now I am trying to see why UCSC genome browser can not read the bam file. Best, ldong

10-09-2010, 07:46 AM	#13
tguo Junior Member Location: maryland Join Date: Sep 2010 Posts: 5	I think I've figured out what the problem was... I have used a UCSC mm9 reference assembly against an EMSEMBL GTF, and that apparently was the problem. Some attributes of the two apparently doesn't match (didn't look in to see exactly why). If instead, I use a UCSC GTF for mm9 reference, or conversely, use the EMSEMBL m_musculus_ncbi37 reference for EMSEMBL GTF, it seems to work.

10-09-2010, 08:26 AM	#14
ldong Member Location: USA Join Date: May 2010 Posts: 15	Hi, tguo, I think the reason is that reference sequence ID in mm9 does not match the chromosome Id in the GTF file. By the way, where did you download the UCSC GTF file? Could you remind me? I could not find it anywhere. I use a GFF file from AceView. Thanks, ldong

10-09-2010, 02:27 PM	#15
tguo Junior Member Location: maryland Join Date: Sep 2010 Posts: 5	Sure, from the UCSC table browser: http://genome.ucsc.edu/cgi-bin/hgTab...a_doMainPage=1

10-10-2010, 07:44 AM	#16
adarob Member Location: Berkeley, CA Join Date: Jul 2010 Posts: 71	Idong, Did you use the --no-novel-juncs command?

10-10-2010, 09:18 AM	#17
tguo Junior Member Location: maryland Join Date: Sep 2010 Posts: 5	No, I didn't. Would that help?

10-10-2010, 10:23 AM	#18
adarob Member Location: Berkeley, CA Join Date: Jul 2010 Posts: 71	I'm sorry I actually should have addressed that to Pejman. The idea is that if you're supplying a GTF in TopHat and you don't know want it find any splice junctions that are not in the GTF, you should use the --no-novel-juncs option.

10-11-2010, 05:13 AM	#19
Pejman Member Location: Switzerland Join Date: Jul 2010 Posts: 23	Hi, thanks for the hint. Yeah I know that, but apparently TopHat+colorspace has just been released in too much rush. It crashes on SOLiD data: http://seqanswers.com/forums/showthread.php?t=7118 http://seqanswers.com/forums/showthread.php?p=26692

10-28-2010, 04:51 AM	#20
nkwuji Member Location: Dublin Join Date: Mar 2010 Posts: 19	Hi, In the tophat wiki page "https://sites.google.com/a/brown.edu/bioinformatics-in-biomed/cufflinks-and-tophat#read", I found some very helpful perl scripts, for example, find_known_gene_and_alternative_splicing_expression_from_cuffcompare_result.pl. Do you konw where I can download them? Cheers, Jun