Calculate forward and reverse strand reads from bam

Question

[bioguy24]

I am trying to calculate the forward and reverse mapped reads from bam. I am eventually going to use a bed of targeted regions but started by using a script from a thread from @Alex Reynolds.

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

output for forwardStrandReads

 1QXCL:01920:06690  16  chr1    10177   0   34S41M  *   0   0   TCGTTACCTACCTACCTACCTACCTAACCTACCTACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA +6+(+.*44456775:44-:44-74/5/:773855/9965*78737*88848+9994<4?@>8<6;<<7808987 ZP:B:f,0.00718598,0.00485627,0.000543075    ZM:B:s,286,-24,286,12,14,256,-6,274,10,206,226,0,488,514,18,2,250,0,24,494,16,268,222,238,498,0,2,234,4,-26,658,-18,468,218,-28,700,476,260,6,734,456,6,-6,-14,242,682,28,-30,454,-2,-24,250,8,16,658,18,440,2,-12,236,34,-4,648,0,472,252,-12,564,346,210,48,574,332,10,-24,-12,244,530,56,-10,398,-34,30,224,-26,84,544,20,344,42,-38,228,82,-22,524,-26,340,222,-30,484,282,192,56,526,236,46,-22,0,320,482,76,12,324,14,48,382,80,180,258,168,124,102,-70,96,310,54,314,10,178,226,54,224,110,126,160   ZF:i:24 MD:Z:41 NM:i:0  AS:i:41 XM:i:41 XA:Z:map4-1 XS:i:41 RG:Z:1QXCL.IonXpress_001-7EBD8685   PG:Z:tmap-1F1907BA

I guess my question is what am I looking for as I am nor sure how to interpret this. Is there documentation that might help? Thank you :).

[GenoMax]

A BBMap suite solution: A: How to Separate a SAM File Containing Forward and Reverse Reads into Two Separat

[Carlo Yague]

Perhaps this can help too : Tutorial: How To Separate Illumina Based Strand Specific Rna-Seq Alignments By Strand

[lucilepain]

Hello everyone,

I am trying to process to counting from data "forward.bam" and "reverse.bam" that have been previously mapped (paired-end mode). I would like to use HT-seq counts , that takes only one .bam as input (forward+reverse). Does anyone know about some tools or documentations that allow to reunite a forward.bam and a reverse.bam in a single one?

Thank by advance

[GenoMax]

Can you clarify as to what your BAM files are? If the data was paired-end there should be just one BAM.

[Carlo Yague]

samtools cat can be used to concatenate two bam files. But as genomax pointed it out, your situation is a little bit weird and you might be doing something wrong.

Answer 1

[michael.ante]

Hi cmccabe,

just go to http://broadinstitute.github.io/picard/explain-flags.html. There you enter the number seen in the second column. In the given case, that's 16. It'll give you a tick at "read reverse strand".

If you are dealing with RNA-Seq data, have a look at RSeQC's infer_experiment tool, since it copes with paired-end information and the read-orientation with regard to the transcripts' orientation. Otherwise, have a look at bam_stat, it gives you the number of reads mapping to each strand.

Cheers,

Michael