How to filter out cross alignments from a BED file?

Question

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?

Answer 1

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'