extract assembled transcripts

[zchou]

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

[john_mu]

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

[zchou]

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?

[john_mu]

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.

[zchou]

Exactly, that's I want. It is necessary to get consensus transcript sequence for each rebuild transcript for further analysis.

[dariober]

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):

## 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:

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]

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?

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):

## 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:

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]

(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:

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:

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

[zchou]

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

[k-gun12]

bump

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

Thanks!

[adarob]

I believe you can upload the GTF to the UCSC genome browser and then extract the sequences from the Table tab.

[lilin001]

Can you share those perl codes for the sequence extraction according to the GTF file? Thanks a lot!

[adarob]

This software will do what you want.

http://deweylab.biostat.wisc.edu/rsem/rsem-prepare-reference-bowtie.html

[edge]

Hi k-gun12,
mind do share the script that you have written to extract the rebuild transcripts sequence?
Thanks in advance.

[edge]

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