Question: Forward Strand Or Reverse Strand
5
gravatar for siyu
6.3 years ago by
siyu130
China
siyu130 wrote:

How can I know that a read is mapping to forward strand or reverse strand ?

strand read • 15k views
ADD COMMENTlink modified 15 months ago by h.mon24k • written 6.3 years ago by siyu130

Hello everyone,

This post is really interesting, I actually need to analyse (ht-seq counts) .bam files (sorted and indexed) that have been mapped by someone else before me. My files are as this form: - sample.merged.mdup.forward.bam - sample.merged.mdup.reverse.bam.

When looking for ht-seq counts tutorials, there was only one .bam (and the .gff) as input... ->Does someone, working on paired-end mode, have any idea of how can we manipulate the forward.bam and the reverse.bam for counting step

Many thanks in advance

ADD REPLYlink written 18 months ago by lucilepain0
14
gravatar for Pierre Lindenbaum
6.3 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum118k wrote:

A sam record contains a binary flag (aka SAM-FLAG). The bit for 16 is set if the read is mapped on the reverse strand

see the SAM specification http://samtools.sourceforge.net/SAM1.pdf ( section "FLAG: bitwise FLAG")

see also "explain SAM flag " http://picard.sourceforge.net/explain-flags.html

ADD COMMENTlink written 6.3 years ago by Pierre Lindenbaum118k
1

"explain SAM flag" has moved on https://broadinstitute.github.io/picard/explain-flags.html

ADD REPLYlink written 2.6 years ago by Carlo Yague4.4k
8
gravatar for Alex Reynolds
6.3 years ago by
Alex Reynolds27k
Seattle, WA USA
Alex Reynolds27k wrote:

One way to extract stranded reads from a BAM file:

samtools view Reads.bam | gawk '(and(16, $2))' > forwardStrandReads.sam
samtools view Reads.bam | gawk '(! and(16, $2))' > reverseStrandReads.sam
ADD COMMENTlink modified 6.3 years ago • written 6.3 years ago by Alex Reynolds27k
3

faster:

samtools view -F 16 Reads.bam 
samtools view -f 16 Reads.bam
ADD REPLYlink modified 6.3 years ago • written 6.3 years ago by Pierre Lindenbaum118k
2

Won't this and Alex's solution will only work if you already have filtered away the unmapped reads using the

-F 4

option? Otherwise for with the '-F 16' option in SAMtools view or 'not containing 16' filtering by gawk you just get everything that isn't mapping to the reverse strand which could be more than those that map to the forward strand.

This post Samtools View: Only Forward Or Reverse Strand  seems to work for single paired reads to sort forward reads without that initial step because it excludes anything that is mapping to reverse strand or unmapped.

ADD REPLYlink modified 3.9 years ago • written 3.9 years ago by Wayne320

Nice, cheers for the suggestion.

ADD REPLYlink written 6.3 years ago by Alex Reynolds27k

Could you give me a detail description about this ? I don't understand why the Flag 16 stands for reverse strand.Thank you !

ADD REPLYlink written 6.3 years ago by siyu130

Read the PDF link in Pierre's answer for an explanation of the flags. Hexadecimal 0x10 is decimal 16.

ADD REPLYlink modified 6.3 years ago • written 6.3 years ago by Alex Reynolds27k

Is this still a SAM file? can it be converted to BAM? because it is missing the headers
 

ADD REPLYlink written 4.6 years ago by ilyco50

SAM is text-based, so perhaps use samtools to excise the headers to a separate file, then cat the header file and a filtered per-strand SAM file. Then convert the per-strand, headered SAM back to BAM.

ADD REPLYlink written 4.6 years ago by Alex Reynolds27k
0
gravatar for liangqinsi
23 months ago by
liangqinsi20
China
liangqinsi20 wrote:

Pierre and Alex's answers are great.

But flag :

16 or 0x10 means SEQ being reverse complemented

Which means it should be:

samtools view Reads.bam | gawk '(and(16, $2))' > reverseStrandReads.sam

samtools view Reads.bam | gawk '(! and(16,$2))' > forwardStrandReads.sam

------

I found it's strange when I grep primer in forward and reverse sam ( most forward primer in reverse sam, and vice versa ), it comes out to be misunderstanding of the flag.

wc -l *sam.out

2551 ATACTGGTATGTATTTAACCATGCAGATCC.forwardStrandReads.sam.out

1001 TTTCAGTTGATTTGCTTGAGATCAAGATTG.reverseStrandReads.sam.out

ADD COMMENTlink modified 23 months ago • written 23 months ago by liangqinsi20
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 752 users visited in the last hour