SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 09-30-2010, 09:56 PM   #1
tguo
Junior Member
 
Location: maryland

Join Date: Sep 2010
Posts: 5
Default 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 10:03 PM.
tguo is offline   Reply With Quote
Old 10-01-2010, 06:27 AM   #2
ldong
Member
 
Location: USA

Join Date: May 2010
Posts: 15
Default

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
ldong is offline   Reply With Quote
Old 10-01-2010, 07:21 AM   #3
tguo
Junior Member
 
Location: maryland

Join Date: Sep 2010
Posts: 5
Default

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
tguo is offline   Reply With Quote
Old 10-01-2010, 08:27 AM   #4
ldong
Member
 
Location: USA

Join Date: May 2010
Posts: 15
Default

New version of tophat works!

Last edited by ldong; 10-08-2010 at 07:49 AM.
ldong is offline   Reply With Quote
Old 10-07-2010, 10:33 AM   #5
Pejman
Member
 
Location: Switzerland

Join Date: Jul 2010
Posts: 23
Default

Quote:
Originally Posted by tguo View Post
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?
Pejman is offline   Reply With Quote
Old 10-07-2010, 11:35 AM   #6
ldong
Member
 
Location: USA

Join Date: May 2010
Posts: 15
Default

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
ldong is offline   Reply With Quote
Old 10-07-2010, 01:13 PM   #7
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

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
adarob is offline   Reply With Quote
Old 10-08-2010, 02:44 AM   #8
Pejman
Member
 
Location: Switzerland

Join Date: Jul 2010
Posts: 23
Default

@ 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 is offline   Reply With Quote
Old 10-08-2010, 04:33 AM   #9
Pejman
Member
 
Location: Switzerland

Join Date: Jul 2010
Posts: 23
Default

@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?
Pejman is offline   Reply With Quote
Old 10-08-2010, 07:29 AM   #10
ldong
Member
 
Location: USA

Join Date: May 2010
Posts: 15
Default

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
ldong is offline   Reply With Quote
Old 10-08-2010, 07:33 AM   #11
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

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
adarob is offline   Reply With Quote
Old 10-08-2010, 08:52 AM   #12
ldong
Member
 
Location: USA

Join Date: May 2010
Posts: 15
Default

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
ldong is offline   Reply With Quote
Old 10-09-2010, 06:46 AM   #13
tguo
Junior Member
 
Location: maryland

Join Date: Sep 2010
Posts: 5
Default

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.
tguo is offline   Reply With Quote
Old 10-09-2010, 07:26 AM   #14
ldong
Member
 
Location: USA

Join Date: May 2010
Posts: 15
Default

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
ldong is offline   Reply With Quote
Old 10-09-2010, 01:27 PM   #15
tguo
Junior Member
 
Location: maryland

Join Date: Sep 2010
Posts: 5
Default

Sure, from the UCSC table browser:

http://genome.ucsc.edu/cgi-bin/hgTab...a_doMainPage=1
tguo is offline   Reply With Quote
Old 10-10-2010, 06:44 AM   #16
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

Idong,

Did you use the --no-novel-juncs command?
adarob is offline   Reply With Quote
Old 10-10-2010, 08:18 AM   #17
tguo
Junior Member
 
Location: maryland

Join Date: Sep 2010
Posts: 5
Default

No, I didn't. Would that help?
tguo is offline   Reply With Quote
Old 10-10-2010, 09:23 AM   #18
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

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.
adarob is offline   Reply With Quote
Old 10-11-2010, 04:13 AM   #19
Pejman
Member
 
Location: Switzerland

Join Date: Jul 2010
Posts: 23
Default

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
Pejman is offline   Reply With Quote
Old 10-28-2010, 03:51 AM   #20
nkwuji
Member
 
Location: Dublin

Join Date: Mar 2010
Posts: 19
Default

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
nkwuji is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 10:06 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO