SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Pileup / extract information from BAM/SAM files liu_xt005 Bioinformatics 4 01-19-2015 09:14 AM
Extract mapped information from bam files empyrean Bioinformatics 0 03-23-2014 10:45 AM
how to know whether a bam file is strand specific? jay2008 Bioinformatics 0 11-27-2012 11:12 PM
How to extract SNP information of a specific population from 1KG data? wwmm933 Bioinformatics 0 12-22-2011 10:00 AM
Re-write strand information in BAM file droog_22 Bioinformatics 1 11-07-2011 03:25 AM

Reply
 
Thread Tools
Old 07-15-2014, 05:41 AM   #1
paperblu
Junior Member
 
Location: Italy

Join Date: Jul 2014
Posts: 2
Default extract strand information from bam file in a specific position

Hi All
I am new on this Forum and I have an issue that I hope you can help me to solve.
I have strand specific paired-end RNA sequencing and I would like to extract the strand information of reads that fall in a particular genomic region (from bam file)
Obviously I have to consider the first in pair to determine the strand orientation (because in paired-end one of reads is always on the opposite strand).

For example, I would like to know if reads falling in chr1:122592431-123294331 are located in plus or minus strand (in bam file).

How would you do it?Any suggestion?

Thank you
paperblu is offline   Reply With Quote
Old 07-15-2014, 06:35 AM   #2
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 310
Default

Quote:
Originally Posted by paperblu View Post
Hi All
I am new on this Forum and I have an issue that I hope you can help me to solve.
I have strand specific paired-end RNA sequencing and I would like to extract the strand information of reads that fall in a particular genomic region (from bam file)
Obviously I have to consider the first in pair to determine the strand orientation (because in paired-end one of reads is always on the opposite strand).

For example, I would like to know if reads falling in chr1:122592431-123294331 are located in plus or minus strand (in bam file).

How would you do it?Any suggestion?

Thank you
It's a bit unclear to me what output format you want. If you want to know how many first-in-pair alignments in a region map to the forward and to the reverse strand, you could use something on these lines:

Code:
region='chr1:1-10000'
bam='myaln.bam'

samtools view -f 64 -F 4 $bam $region \
| gawk '{if (and($2, 16)) {cnt_reverse++} else {cnt_forward++}}END{print "Reads on forward: " cnt_forward "\nReads on reverse: " cnt_reverse}' 

## Example output:
Reads on forward: 2970
Reads on reverse: 3950
dariober is offline   Reply With Quote
Old 07-15-2014, 07:35 AM   #3
paperblu
Junior Member
 
Location: Italy

Join Date: Jul 2014
Posts: 2
Default

Hi dariober
thank you for your reply and I am sorry if I have been a little unclear.
But yes, my desired output should be something like your (I would like to obtain the number just of reads firs in pair or second in pair, that are, if I am correct, reads that mark genes in positive strand or negative).

I made this:

Code:
samtools view -b $mybam $myregion | samtools view -f 0x0040 - | wc -l | awk '{ print$0"\t+"}'
output

Code:
7    +
That means that in $myregion I have 7 reads first in pair.Isn't it?
The problem is that for each coordinate I have to perfom the command twice to know first in pair and second in pair.
Any other suggestion ?

Thank you
paperblu 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:43 AM.


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