<p>Hi,</p>

<p>I want to extract read-pairs that aligned concordantly exactly 1 time to the reference genome. Is there any way to parse the output <a href='http://samtools.sourceforge.net/SAM1.pdf'>SAM</a> file? I would really appreciate your suggestions.</p>

<p>Best regards,</p>

<p>Rahul </p>

<p>The simplest way is likely something along the lines of:</p>

<pre><code>samtools view -hf 0x2 alignments.bam | grep -v "XS:i:" &gt; filtered.alignments.sam
</code></pre>

<p>As I mentioned in the comment above, the <code>-f 0x2</code> part will get only "properly paired" alignments, which will effectively be the concordant alignments. For the option #2 definition of "map only once", you can take advantage of the fact that bowtie2 will add an <code>XS</code> auxiliary tag to reads that have another "valid" mapping. So, a quick inverse grep (<code>grep -v</code>) can get rid of those. One possible problem with this is if only one read of a pair can map to multiple places (e.g., the original fragment partly overlapped a simple tandem repeat) then you'd end up with orphans. The easiest fix to that (if it's a problem) would be to just check if the read names of pairs of alignments are the same. I'm sure someone will come up with a way to do that in awk, but the following python script is likely simpler:</p>

<pre><code>#!/usr/bin/env python
import csv
import sys

f = csv.reader(sys.stdin, dialect="excel-tab")
of = csv.writer(sys.stdout, dialect="excel-tab")
last_read = None
for line in f :
    #take care of the header
    if(line[0][0] == "@") :
        of.writerow(line)
        continue

    if(last_read == None) : 
        last_read = line
    else :
        if(last_read[0] == line[0]) :
            of.writerow(last_read)
            of.writerow(line)
            last_read = None
        else :
            last_read = line
</code></pre>

<p>If you saved that as <code>foo.py</code> and made it executable and in your PATH, then the following would solve the aforementioned issue:</p>

<pre><code>samtools view -hf 0x2 alignments.bam | grep -v "XS:i:" | foo.py &gt; filtered.alignments.sam
</code></pre>

<p>Note that this won't work with coordinate-sorted alignments, but bowtie2 doesn't produce those.</p>
