<p>According to the <a href="https://samtools.github.io/hts-specs/SAMv1.pdf" rel="nofollow noreferrer">SAM specification</a>, the 3rd field of a SAM line (<code>RNAME</code>) is:</p>

<blockquote>
  <p>RNAME: Reference sequence NAME of the alignment. If @SQ header lines
  are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
  tag. An unmapped segment without coordinate has a ‘*’ at this field.
  However, an unmapped segment may also have an ordinary coordinate such
  that it can be placed at a desired position after sorting. If RNAME is
  ‘*’, no assumptions can be made about POS and CIGAR.</p>
</blockquote>

<p>And the 7th field is (emphasis mine, missing "to" theirs):</p>

<blockquote>
  <p>RNEXT: Reference sequence name of the primary alignment of the NEXT
  read in the template. For the last read, the next read is the first
  read in the template. If @SQ header lines are present, RNEXT (if not
  ‘*’ or ‘=’) must be present in one of the SQ-SN tag. This field is set
  as ‘*’ when the information is unavailable, <strong>and set as ‘=’ if RNEXT is
  identical RNAME</strong>. If not ‘=’ and the next read in the template has one
  primary mapping (see also bit 0x100 in FLAG), this field is identical
  to RNAME at the primary line of the next read. If RNEXT is ‘*’, no
  assumptions can be made on PNEXT and bit 0x20</p>
</blockquote>

<p>So, you want to remove those lines whose 7th field isn't <code>=</code> and, just in case, those lines whose 7th field isn't <code>=</code> <em>and</em> isn't the same as the 3rd field. You can therefore use something like this:</p>

<pre><code>samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7
</code></pre>

<p>And, to save as a bam file again:</p>

<pre><code>samtools view -L test.bed test.bam | awk '$7=="=" &amp;&amp; $3==$7 | 
    samtolls view -b &gt; fixed.bam
</code></pre>

<hr>

<p>On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use <code>\|</code> (or <code>|</code> with the <code>-E</code> or <code>-P</code> options) to separate them. Something like:</p>

<pre><code>samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'
</code></pre>

<p>Or</p>

<pre><code>samtools view -L test.bed test.bam | grep -Ev 'chrx|chr2|chr10|chrN'
</code></pre>

CC BY-SA 3.0
2017-05-26T10:18:26.377Z
If you do it like this, the `fixed.bam` file is missing the header, which in my experience generates a lot of problems. I recommend always adding the header back in; either by specifying `-h` when reading the original BAM, or by adding it separately: `(samtools view -H infile.bam; samtools view …) > samtools view -b > outfile.bam`.
Konrad Rudolph
CC BY-SA 3.0
2017-05-19T17:44:29.193Z
2017-05-21T23:21:14.657Z
2017-05-21T23:21:14.657Z
terdon
<p>I have a BAM file:</p>

<pre><code>@SQ SN:chr1 LN:248956422
@SQ SN:chrx LN:248956423
ST-E00110:348:HGVKKALXX:1:1201:5822:48670   323 chr1    9999    0   67H66M16H   chrx    1000    0   GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC  JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ  NM:i:0  MD:Z:66 AS:i:66 XS:i:65 SA:Z:chr5,18606834,-,73S76M,34,0;   RG:Z:g1
</code></pre>

<p>There is a read aligned to <code>chr1</code>, and it's mate aligned to <code>chrx</code>. </p>

<p>I have a BED file:</p>

<pre><code>chr1    0   100000  TestOnly
</code></pre>

<p>I would like to filter out everything that falls outside my BED region, that includes cross-alignments. In my example, although my read aligned to <code>chr1</code> but it's mate is not. I don't want this read.</p>

<p>When I do:</p>

<blockquote>
  <p>samtools view -L test.bed test.bam</p>
</blockquote>

<p>the command gives me the read because it doesn't check cross-alignments.</p>

<p>My solution:</p>

<blockquote>
  <p>samtools view -L test.bed test.bam | grep -v chrx</p>
</blockquote>

<p>but this is very slow and clumsy. In my production pipeline I would have to do something like:</p>

<blockquote>
  <p>samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...</p>
</blockquote>

<p><strong>Q:</strong> Is there a better solution?</p>

<sam><bed><genomics>
How to filter out cross alignments from a BED file?
# How to filter out cross alignments from a BED file?

(Asked by: username_0 on 2017-05-19)

I have a BAM file:
@SQ SN:chr1 LN:248956422
@SQ SN:chrx LN:248956423
ST-E00110:348:HGVKKALXX:1:1201:5822:48670   323 chr1    9999    0   67H66M16H   chrx    1000    0   GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC  JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ  NM:i:0  MD:Z:66 AS:i:66 XS:i:65 SA:Z:chr5,18606834,-,73S76M,34,0;   RG:Z:g1

There is a read aligned to chr1, and it's mate aligned to chrx. 
I have a BED file:
chr1    0   100000  TestOnly

I would like to filter out everything that falls outside my BED region, that includes cross-alignments. In my example, although my read aligned to chr1 but it's mate is not. I don't want this read.
When I do:

samtools view -L test.bed test.bam

the command gives me the read because it doesn't check cross-alignments.
My solution:

samtools view -L test.bed test.bam | grep -v chrx

but this is very slow and clumsy. In my production pipeline I would have to do something like:

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

Q: Is there a better solution?

**Tags:** sam, bed, genomics

---

## 1 answers

### Answer by username_1 on 2017-05-19. Score: 6 | Accepted answer

According to the SAM specification, the 3rd field of a SAM line (RNAME) is:

RNAME: Reference sequence NAME of the alignment. If @SQ header lines
  are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
  tag. An unmapped segment without coordinate has a ‘*’ at this field.
  However, an unmapped segment may also have an ordinary coordinate such
  that it can be placed at a desired position after sorting. If RNAME is
  ‘*’, no assumptions can be made about POS and CIGAR.

And the 7th field is (emphasis mine, missing "to" theirs):

RNEXT: Reference sequence name of the primary alignment of the NEXT
  read in the template. For the last read, the next read is the first
  read in the template. If @SQ header lines are present, RNEXT (if not
  ‘*’ or ‘=’) must be present in one of the SQ-SN tag. This field is set
  as ‘*’ when the information is unavailable, and set as ‘=’ if RNEXT is
  identical RNAME. If not ‘=’ and the next read in the template has one
  primary mapping (see also bit 0x100 in FLAG), this field is identical
  to RNAME at the primary line of the next read. If RNEXT is ‘*’, no
  assumptions can be made on PNEXT and bit 0x20

So, you want to remove those lines whose 7th field isn't = and, just in case, those lines whose 7th field isn't = and isn't the same as the 3rd field. You can therefore use something like this:
samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7

And, to save as a bam file again:
samtools view -L test.bed test.bam | awk '$7=="=" && $3==$7 | 
    samtolls view -b > fixed.bam


On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use \| (or | with the -E or -P options) to separate them. Something like:
samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'

Or
samtools view -L test.bed test.bam | grep -Ev 'chrx|chr2|chr10|chrN'

> **Comments:**
> **username_2**: If you do it like this, the `fixed.bam` file is missing the header, which in my experience generates a lot of problems. I recommend always adding the header back in; either by specifying `-h` when reading the original BAM, or by adding it separately: `(samtools view -H infile.bam; samtools view …) > samtools view -b > outfile.bam`.  

---

