SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
find ORFs for cufflinks assembled transcripts jiexiong Bioinformatics 1 05-08-2014 12:26 PM
fetch transcripts assembled by cufflinks asling Bioinformatics 6 09-27-2012 09:46 PM
Improvements to assembled genome NPalopoli De novo discovery 6 01-31-2012 09:00 AM
Assembled sequence submission to Genbank? Melissa General 0 04-26-2011 12:54 AM
454 assembled reads katussa10 Bioinformatics 7 12-02-2010 01:19 PM

Reply
 
Thread Tools
Old 05-13-2010, 11:55 AM   #1
zchou
Member
 
Location: detroit, US

Join Date: May 2009
Posts: 13
Default extract assembled transcripts

Hi All,

I use Bowtie/TopHat/cufflink to analyze the RNA-Seq. I want to extract the rebuild transcripts by these softwares. However, I can only the BED file. Can anyone give idea to extract these assembled transcript sequence?

Thanks,
ZC
zchou is offline   Reply With Quote
Old 05-20-2010, 11:11 PM   #2
john_mu
Member
 
Location: Stanford, CA

Join Date: May 2010
Posts: 88
Default

If you want to assemble the transcript with cufflinks you need to use the SAM file in the output of Tophat. If Tophat is not outputting a SAM file.. then maybe you need to email the authors for help.

Cufflinks takes SAM files as input.

Alternatively, you could try SpliceMap as another option, followed by Cufflinks
__________________
SpliceMap: De novo detection of splice junctions from RNA-seq
Download SpliceMap Comment here
john_mu is offline   Reply With Quote
Old 05-21-2010, 06:19 AM   #3
zchou
Member
 
Location: detroit, US

Join Date: May 2009
Posts: 13
Default

Thanks. I used tophat/cufflink to rebuild the transcripts. tophat generates sam format results and then, cufflink rebuilds transcripts. However, cufflink can only give the rebuild transcript name. If we want to analyze infurther based on these rebuild transcripts, we have to get rebuilt transcript sequence. Is there tool to get seq from gtf file, except write script by ourself?
zchou is offline   Reply With Quote
Old 05-21-2010, 09:53 AM   #4
john_mu
Member
 
Location: Stanford, CA

Join Date: May 2010
Posts: 88
Default

Quote:
Originally Posted by zchou View Post
Thanks. I used tophat/cufflink to rebuild the transcripts. tophat generates sam format results and then, cufflink rebuilds transcripts. However, cufflink can only give the rebuild transcript name. If we want to analyze infurther based on these rebuild transcripts, we have to get rebuilt transcript sequence. Is there tool to get seq from gtf file, except write script by ourself?
So for example, if cufflinks outputs the following assembled transcript

chr21 Cufflinks transcript 39473730 39474150 1000 - . gene_id "CUFF.171"; transcript_id "CUFF.171.1"; FPKM "15397.8812515398"; frac "1.000000"; conf_lo "0.000000"; conf_hi "33177.823023"; cov "3.252033";

chr21 Cufflinks exon 39473730 39473782 1000 - . gene_id "CUFF.171"; transcript_id "CUFF.171.1"; exon_number "1"; FPKM "15397.8812515398"; frac "1.000000"; conf_lo "0.000000"; conf_hi "33177.823023"; cov "3.252033";

chr21 Cufflinks exon 39474081 39474150 1000 - . gene_id "CUFF.171"; transcript_id "CUFF.171.1"; exon_number "2"; FPKM "15397.8812515398"; frac "1.000000"; conf_lo "0.000000"; conf_hi "33177.823023"; cov "3.252033";

You would like a tool to output some kind of consensus sequence of this transcript with two exons based on all reads that mapped there?

I'm not aware of such a tool, but I see how it could be useful.
__________________
SpliceMap: De novo detection of splice junctions from RNA-seq
Download SpliceMap Comment here
john_mu is offline   Reply With Quote
Old 05-21-2010, 10:22 AM   #5
zchou
Member
 
Location: detroit, US

Join Date: May 2009
Posts: 13
Default

Exactly, that's I want. It is necessary to get consensus transcript sequence for each rebuild transcript for further analysis.
zchou is offline   Reply With Quote
Old 05-23-2010, 08:45 AM   #6
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 311
Default

Hello,

One option could be to use samtools to generate a pileup file of the mapped reads (output of tophat). Then from the pileup format you could extract the column of the consensus base (the 4th column I think) at the genomic positions corresponding to your exons/transcripts.

Something like (pseudocode):

Code:
## Use -c option to output consensus sequence, not just the reference sequence
samtools pilup -c tophat_output_accepted_hits.sam/bam -f myreference_genome.fa
Pileup file should like:
Code:
MT      1       C       C       36      0       60      3       ^~.^~.^~.       ABA     ~~~
MT      2       A       A       12      0       60      6       ...^~T^~.^~.    =BCBAB  ~~~~~~
etc..
It would take a little bit of coding to reconstruct the sequences from the pileup file but it's not too bad.

Hope this helps

Dario
dariober is offline   Reply With Quote
Old 05-23-2010, 05:06 PM   #7
zchou
Member
 
Location: detroit, US

Join Date: May 2009
Posts: 13
Default

samtools can only give consensus sequence information, but where is the reconstructed sequence information for each transcript? Moreover, cufflink or other RNA-Seq software can get different isoforms for each gene. How can we extract these information?



Quote:
Originally Posted by dariober View Post
Hello,

One option could be to use samtools to generate a pileup file of the mapped reads (output of tophat). Then from the pileup format you could extract the column of the consensus base (the 4th column I think) at the genomic positions corresponding to your exons/transcripts.

Something like (pseudocode):

Code:
## Use -c option to output consensus sequence, not just the reference sequence
samtools pilup -c tophat_output_accepted_hits.sam/bam -f myreference_genome.fa
Pileup file should like:
Code:
MT      1       C       C       36      0       60      3       ^~.^~.^~.       ABA     ~~~
MT      2       A       A       12      0       60      6       ...^~T^~.^~.    =BCBAB  ~~~~~~
etc..
It would take a little bit of coding to reconstruct the sequences from the pileup file but it's not too bad.

Hope this helps

Dario
zchou is offline   Reply With Quote
Old 05-24-2010, 12:27 AM   #8
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 311
Default

(Just rephrasing my previous post...)

"samtools pileup" with -c option gives at each genomic position the base of the reference *and* the base reconstructed from the reads. This latter is what you need, isn't it?

As I said, to reconstruct the sequences of transcripts and isoforms you would need to write some script that extracts the regions from the pileup file that correspond to the transcripts in the cufflinks gtf file.

So, if the GTF file looks (vaguely) like this:

Code:
chr    start    end    feature          exon_number           transcript_id
-------------------------------------------------------------------------
chr1    1         9        transcript                                cuff_1
chr1    1         3         exon             1                       cuff_1
chr1    6         9         exon             2                       cuff_1
etc.
And the pileup file looks (vaguely) like this:

Code:
chr    position        ref_base    consensus
-------------------------------------------
chr1     1                 A             T
chr1     2                 C             T
chr1     3                 G             T
chr1     4                 A             A
chr1     5                 C             A
chr1     6                 G             G
chr1     7                 A             G
chr1     8                 C             G
chr1     9                 G             G
etc.
Than you need a script that extracts the bases in the consensus column corresponding to chr1:1-3 (TTT) and to chr1:6-9 (GGGG) to have the reconstructed sequence of the transcript cuff_1.

Apologies if I'm misunderstanding your problem

Dario
dariober is offline   Reply With Quote
Old 05-24-2010, 05:28 AM   #9
zchou
Member
 
Location: detroit, US

Join Date: May 2009
Posts: 13
Default

Thanks dario,

It is very clear to illustrate how to get transcript sequence from sam and gtf file. I want to write scripts to extract sequence from gtf and fasta file before. Thanks a lot,

ZC
zchou is offline   Reply With Quote
Old 10-21-2010, 07:31 AM   #10
k-gun12
Member
 
Location: NJ

Join Date: Feb 2010
Posts: 54
Default

bump

I'm looking to do this myself.. anyone have a script handy? If not, I'll get started.

Thanks!
k-gun12 is offline   Reply With Quote
Old 10-21-2010, 09:06 AM   #11
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

I believe you can upload the GTF to the UCSC genome browser and then extract the sequences from the Table tab.
adarob is offline   Reply With Quote
Old 12-09-2010, 09:09 AM   #12
lilin001
Junior Member
 
Location: minnesota

Join Date: Dec 2010
Posts: 3
Default

Can you share those perl codes for the sequence extraction according to the GTF file? Thanks a lot!
lilin001 is offline   Reply With Quote
Old 12-09-2010, 09:15 AM   #13
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

This software will do what you want.

http://deweylab.biostat.wisc.edu/rse...ce-bowtie.html
adarob is offline   Reply With Quote
Old 05-11-2011, 10:31 PM   #14
edge
Senior Member
 
Location: China

Join Date: Sep 2009
Posts: 199
Default

Hi k-gun12,
mind do share the script that you have written to extract the rebuild transcripts sequence?
Thanks in advance.
edge is offline   Reply With Quote
Old 05-11-2011, 11:40 PM   #15
edge
Senior Member
 
Location: China

Join Date: Sep 2009
Posts: 199
Default

Hi zchou,
Do you already figure out the problem to extract assembled transcript?
Can I have a reference on how you do that?
I'm facing the same problems as well now
edge is offline   Reply With Quote
Old 05-17-2011, 09:24 AM   #16
k-gun12
Member
 
Location: NJ

Join Date: Feb 2010
Posts: 54
Default

It turns out that a rather nice one is shipped with EVidence Modeler:
http://evidencemodeler.sourceforge.net/

.. in the utilities directory.
k-gun12 is offline   Reply With Quote
Old 05-17-2011, 02:16 PM   #17
adarob
Member
 
Location: Berkeley, CA

Join Date: Jul 2010
Posts: 71
Default

There is a new utility included with cufflinks now, called "gffread", which should be able to extract spliced transcript sequences from a GFF/GTF. This is now documented on this page (scroll down to the "Extracting transcript sequences" section at the end):

http://cufflinks.cbcb.umd.edu/gff.html
adarob is offline   Reply With Quote
Old 10-25-2011, 01:56 AM   #18
Anelda
Member
 
Location: Cape Town, South Africa

Join Date: May 2010
Posts: 30
Default

Quote:
Originally Posted by adarob View Post
There is a new utility included with cufflinks now, called "gffread", which should be able to extract spliced transcript sequences from a GFF/GTF. This is now documented on this page (scroll down to the "Extracting transcript sequences" section at the end):

http://cufflinks.cbcb.umd.edu/gff.html
The problem with this script is that it returns the reference DNA sequence (from the genome fasta file used for mapping) for the transcripts and not the assembly from the mapped reads.

Does anyone know of a tool that can give the assembled transcript sequence from the RNA-seq data? Or can someone help me if I'm misunderstanding GFFRead?

Thanks!
Anelda is offline   Reply With Quote
Old 11-23-2011, 08:15 AM   #19
pkerrwall
Junior Member
 
Location: Cary, NC

Join Date: Dec 2009
Posts: 2
Default

Quote:
Originally Posted by Anelda View Post
The problem with this script is that it returns the reference DNA sequence (from the genome fasta file used for mapping) for the transcripts and not the assembly from the mapped reads.

Does anyone know of a tool that can give the assembled transcript sequence from the RNA-seq data? Or can someone help me if I'm misunderstanding GFFRead?
Here is the process that I am using:

samtools mpileup -uf ref.fa accepted_hits.bam | bcftools view -cg - | vcfutils.pl vcf2fq | fq2fa.pl > new_ref.fa
gffread -w transcripts.fa -g new_ref.fa transcripts.gtf

where fq2fa.pl is a bioperl script to convert from fq to fasta

I also have an email into the cufflinks developers to see if there is a way that the gffread utility can be enhanced to get the consensus sequence from the bam file and not the reference genome.
pkerrwall is offline   Reply With Quote
Old 12-11-2011, 04:14 PM   #20
huangjun
Member
 
Location: Wuhan China

Join Date: Dec 2011
Posts: 13
Default

the cufflinks only give us the rebuild assembled sequence based on the reference genome, but i want to know how to get unmapped contigs also , only that we could get some new transcripts.
huangjun 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 02:08 AM.


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