extract strand information from bam file in a specific position

[paperblu]

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

[dariober]

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:

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

[paperblu]

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:

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

output

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