Question: Calculate forward and reverse strand reads from bam
1
3.8 years ago by
bioguy24190
Chicago
bioguy24190 wrote:

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 :).

ngs • 2.8k views
ADD COMMENTlink
modified 3.8 years ago by michael.ante3.6k • written 3.8 years ago by bioguy24190
1

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

ADD REPLYlink modified 3.8 years ago • written 3.8 years ago by genomax91k
1

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

ADD REPLYlink written 3.8 years ago by Carlo Yague5.0k

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

ADD REPLYlink written 3.1 years ago by lucilepain0

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

ADD REPLYlink written 3.1 years ago by genomax91k

`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.

ADD REPLYlink written 3.1 years ago by Carlo Yague5.0k
3
3.8 years ago by
michael.ante3.6k
Austria/Vienna
michael.ante3.6k wrote:

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

ADD COMMENTlink written 3.8 years ago by michael.ante3.6k
Please log in to add an answer.

Content
Help
Access

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