How to filter out cross alignments from a BED file? 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?
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=="=" &amp;&amp; $3==$7 | 
    samtolls view -b &gt; 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'