SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Re-write strand information in BAM file droog_22 Bioinformatics 1 11-07-2011 03:25 AM
cufflinks transcript on wrong strand ornam Bioinformatics 2 12-13-2010 12:33 PM
Cufflinks library-type for SOLiD strand specific single data gengen Bioinformatics 1 10-19-2010 08:41 AM
tophat islands information syslm01 Bioinformatics 0 08-11-2010 03:48 AM
Bowtie, BWA, strand information Zimbobo Bioinformatics 6 04-24-2010 08:30 AM

Reply
 
Thread Tools
Old 05-10-2010, 05:54 PM   #21
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

Haneko,

Does BioScope report a strand for the spliced alignments which correspond to the strand of the "parent" transcript? If so you could insert the XS tag only in the spliced reads. I would not insert XS tag in the reads contained within exons since it would probably only confuse Cufflinks or give you faulty estimates of the expression of transcripts.

Alternatively you could ask Cole (the principal author of TopHat and Cufflinks) if there is any way to run TopHat on colorspace sequences. I'm pretty sure Bowtie supports it in any case, so TopHat might eventually as well (or it's already there and I just haven't seen the option for it).
Thomas Doktor is offline   Reply With Quote
Old 05-11-2010, 02:21 AM   #22
damiankao
Member
 
Location: UK

Join Date: Jan 2010
Posts: 49
Default

Bioscope will only do split reads according to the .gtf you've given it. It will not to do de novo split reads.
damiankao is offline   Reply With Quote
Old 05-12-2010, 11:12 AM   #23
xhuister
Member
 
Location: China

Join Date: Apr 2010
Posts: 41
Default

Hi all,

I tested the strand-specific RNA-seq using TopHat and cufflinks in order to get enriched regions or transcript boundaries.

1) TopHat output accepted_hits.sam file, where the splice reads have "XS:A:+-" tag while other reads don't.

2) I seperated accepted_hits.sam into two files accepted_hits.plus.sam and accepted_hits.minus.sam. The XS:A:+ reads and reads with flag=0 are put into accepted_hits.plus.sam and XS:A:-... to minus.sam.
Here, I didn't add XS:A tag to exon-reads mannually.

3) Cufflinks to get regions:
cufflinks accepted_hits.minus.sam
cufflinks accepted_hits.plus.sam

In one output transcript.expr of cufflinks:
trans_id bundle_id chr left right FPKM FMI frac FPKM_conf_lo FPKM_conf_hi coverage length
CUFF.1.1 21679 chloroplast 165 1522 67673.5 1 1 66397.6 68949.5 414.607 1357
CUFF.3.1 21683 chloroplast 1788 1889 2424.11 1 1 1538.95 3309.27 14.8515 101

In this post http://seqanswers.com/forums/showthr...links+coverage, Cole said that "Multiplying the average depth of coverage by the transcript length will give you the *estimated* number of reads assigned to each transcript. "

But, when I multiplied the above coverage and length: 414.607*1357=562621.699, it just didn't meet the actual number of reads in that transcript which is 11341.
$ awk -F'\t' '$3=="chloroplast" && $4>=165 && $4<=1522' accepted_hits.minus.sam | wc -l
11341

Am I wrong in calculating the number of reads in one transcript?

Another question:
In the output accepted_hits.sam of Tophat, there are 408,360 rows. But, when I used Bowtie to do the read mapping with at most 2 mismatches and only one location, the number is much larger 2,288,016.
Does Tophat filter the mapping results and how? I thought that the number of mapped reads in tophat should larger than Bowtie's result, because Tophat contains the junction reads.

Thank you in advance.

Last edited by xhuister; 05-12-2010 at 11:23 AM.
xhuister is offline   Reply With Quote
Old 05-13-2010, 01:42 PM   #24
xguo
Member
 
Location: Maryland

Join Date: Jul 2008
Posts: 48
Default

Quote:
Originally Posted by xhuister View Post
Hi all,

I tested the strand-specific RNA-seq using TopHat and cufflinks in order to get enriched regions or transcript boundaries.

1) TopHat output accepted_hits.sam file, where the splice reads have "XS:A:+-" tag while other reads don't.

2) I seperated accepted_hits.sam into two files accepted_hits.plus.sam and accepted_hits.minus.sam. The XS:A:+ reads and reads with flag=0 are put into accepted_hits.plus.sam and XS:A:-... to minus.sam.
Here, I didn't add XS:A tag to exon-reads mannually.

3) Cufflinks to get regions:
cufflinks accepted_hits.minus.sam
cufflinks accepted_hits.plus.sam

In one output transcript.expr of cufflinks:
trans_id bundle_id chr left right FPKM FMI frac FPKM_conf_lo FPKM_conf_hi coverage length
CUFF.1.1 21679 chloroplast 165 1522 67673.5 1 1 66397.6 68949.5 414.607 1357
CUFF.3.1 21683 chloroplast 1788 1889 2424.11 1 1 1538.95 3309.27 14.8515 101

In this post http://seqanswers.com/forums/showthr...links+coverage, Cole said that "Multiplying the average depth of coverage by the transcript length will give you the *estimated* number of reads assigned to each transcript. "

But, when I multiplied the above coverage and length: 414.607*1357=562621.699, it just didn't meet the actual number of reads in that transcript which is 11341.
$ awk -F'\t' '$3=="chloroplast" && $4>=165 && $4<=1522' accepted_hits.minus.sam | wc -l
11341

Am I wrong in calculating the number of reads in one transcript?

Another question:
In the output accepted_hits.sam of Tophat, there are 408,360 rows. But, when I used Bowtie to do the read mapping with at most 2 mismatches and only one location, the number is much larger 2,288,016.
Does Tophat filter the mapping results and how? I thought that the number of mapped reads in tophat should larger than Bowtie's result, because Tophat contains the junction reads.

Thank you in advance.
The coverage in cufflink output seems to be base-level coverage. If the read length is 50, 11341*50 = 567050 is roughly the same as 562621 derived by the multiplification of coverage and transcript length. Does it make more sense to use FPKM to get expected reads for each transcript?
xguo is offline   Reply With Quote
Old 05-17-2010, 02:08 PM   #25
joung
Junior Member
 
Location: USA

Join Date: May 2010
Posts: 1
Default

Hi everyone,

I am still not sure about XS tag.
How does tophat generate XS:A:+/- ?
joung is offline   Reply With Quote
Old 05-17-2010, 05:03 PM   #26
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

It looks at the genomic sequence of the splice site spanning reads to identify the strand of the transcript. XS tells you the strand of the transcript that produced the read, not the read itself.
Thomas Doktor is offline   Reply With Quote
Old 07-30-2011, 08:42 PM   #27
johnwu
Junior Member
 
Location: Taiwan

Join Date: Jun 2011
Posts: 5
Question

Hi! I'm new to the RNA-seq field, and now I'm trying to use tophat-cufflinks pipeline to analyze strand-specific libraries (illumina).

After reading above posts in this thread, I'm still confused by the existence of such alignment (in SAM format):

ERR030871.1444950 16 chr10 83984 255 71M689N29M * 0 0 CACTGACTCCATCAGCTCCGCGCCTTCGGTGTAGTGTCCCTTAGCCCAGTTGTTTCCGGCCCCACACTGACCGCTGGCCTCGTTGTAGTACACGTTGATG G;EHHEGDGCFGHAFGBHIHHHFHHDHHHHHHHHHHIHHHHHHHHHHHHHFHHEHHHHHHHHHAHHHHHHHHHHHHFHHHHHHHHHHHHHHHHHHHHHHH NM:i:1 XS:A:- NH:i:1

According to the read's FLAG (0x10 SEQ being reverse complemented) and the XS tag (XS:A:-), the complement sequence of the read (according to its FLAG: 0x10 SEQ being reverse complemented) could be aligned to negative strand of chromosome 10.

What I don't get is that since the reverse complement sequence of the read could be aligned to negative strand of chromosome, then the sequence of the read should be directly aligned to positive strand as well, why doesn't tophat just report so? like:

ERR030871.1444950 0 chr10 83984 255 71M689N29M * 0 0 CACTGACTCCATCAGCTCCGCGCCTTCGGTGTAGTGTCCCTTAGCCCAGTTGTTTCCGGCCCCACACTGACCGCTGGCCTCGTTGTAGTACACGTTGATG G;EHHEGDGCFGHAFGBHIHHHFHHDHHHHHHHHHHIHHHHHHHHHHHHHFHHEHHHHHHHHHAHHHHHHHHHHHHFHHHHHHHHHHHHHHHHHHHHHHH NM:i:1 XS:A:+ NH:i:1

By the way, I'm using Tophat version 1.30.

Thanks for your help
johnwu is offline   Reply With Quote
Old 04-03-2012, 07:46 PM   #28
sterding
Member
 
Location: Boston

Join Date: Sep 2010
Posts: 36
Default

Very useful post!

I followed this post due to the same problem: For strand-specific RNAseq data, the accepted_hits.sam (converted from BAM) from Tophat does not contain XS:A:+/- tag for unspliced alignment (yes, for spliced one, there are XS:A tag). This will result in the loss of strand information if I run cufflinks directly using the Tophat output file.

To fix the issue, I found it's useful if I manually add XS:A tag for those non-spliced alignments. For example:

Code:
samtools view accepted_hits.bam | awk '{if($0 ~ /XS:A:/) print $0; else {if($2==16) print $0"\tXS:A:-"; else print $0"\tXS:A:+";}}' > accepted_hits.sam
Again, as Thomas Doktor said, it's safe to do so only if the lib is from strand-specific RNAseq. For non-strand-specific one, it makes no sense to do so.

What I want to comment here is, there is still slightly problem for the solution -- For multiple mapped reads, the FLAG (e.g. 2nd column of SAM) is 256. In that case, we can't really tell which strand the reads map to. This is especially the case for output from aligner allowing a number of mismatches, or reads mapped in repeat region. To guarantee the assembly of de novo transcript, a dirty way to do that might be to duplicate the reads and assign both XS:A:- and XS:A:+. I've not tested this solution. If anyone else has better solution, I would like to know.
sterding is offline   Reply With Quote
Old 04-04-2012, 06:43 AM   #29
sterding
Member
 
Location: Boston

Join Date: Sep 2010
Posts: 36
Default

I just noticed that there are FLAG values of 256, 272 (not just 0, 16, 4 for single-strand RNAseq). Since 256=256+0, 272=256+16, so we may assign 256 to plus (+) strand, and minus (-) for 272 reads.

So, the above code for adding XS:A tag should be modified as:

samtools view accepted_hits.bam | awk '{if($0 ~ /XS:A:/) print $0; else {if($2%256==16) print $0"\tXS:A:-"; if($2%256==0) print $0"\tXS:A:+";}}' > accepted_hits.sam


UPDATE: This only works for single-end lib. For pair-end lib, the FLAG should be odd number, but in any case, reads on minus strand always have 1 on the 5th bit of binary code (e.g. 0x10 =10000). Thanks to Wei's suggestion on this. Here is the updated code:

samtools view -h accepted_hits.bam | awk '{if($0 ~ /XS:A:/ || $1 ~ /^@/) print $0; else {if(and($2,0x10)) print $0"\tXS:A:-"; else print $0"\tXS:A:+";}}' accepted_hits.sam

Last edited by sterding; 05-09-2012 at 06:11 AM.
sterding is offline   Reply With Quote
Old 04-20-2012, 12:35 AM   #30
blanco
Member
 
Location: Iceland

Join Date: Apr 2012
Posts: 28
Default

I agree, a very useful thread indeed.
Still it has been some time since it started and now tophat has a -library--type option.

I have sequences from a dUTP strand specific library and I ran it through tophat with the -library--type option. By using this option nearly all (99.9%) lines in the .sam file have the XS:A tag.

By applying this -library--type option to a strand specific library can one safely use the XS:A tag to say from which strand the read came from?
blanco 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 08:32 PM.


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