Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-26T10:18:26.377Z
Answers.Text: 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`.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-19T17:44:29.193Z
Answers.LastActivityDate: 2017-05-21T23:21:14.657Z
Answers.LastEditDate: 2017-05-21T23:21:14.657Z
Answers.UserDisplayName: terdon
Body: <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>
  
Tags: <sam><bed><genomics>
Title: How to filter out cross alignments from a BED file?
ThreadText: # 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](https://samtools.github.io/hts-specs/SAMv1.pdf), 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`.  
  
  ---
  
  