1
$\begingroup$

I'm using pysam (and also samtools) to find all reads in BAM files that are within a specific region py_bamfile.fetch('4',42266768,42268410) and the following read was returned:

NB501288:541:HCWKVBGXB:2:21107:25304:16421  0   #16 41910661    255 27M167903N19M266888N9M  *   0   0   TTTCTTTTTTTTTTTTTTTTTTTTTTGAAAGGGGGTTTTTTTGTTAAACCAAGGG array('B', [32, 32, 32, 32, 32, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 14, 14, 14, 14, 14, 14, 14, 36, 32, 14, 36, 36, 36, 36, 14, 32, 32, 14, 14, 14, 14, 14, 14, 14, 27, 14, 14, 14, 14, 14])    [('NH', 1), ('HI', 1), ('AS', 40), ('nM', 1), ('RE', 'I'), ('xf', 0), ('li', 0), ('CR', 'CATTTCACAGGGACTA'), ('CY', 'AAAAAEEEEEEEEEEE'), ('CB', 'CATTTCACAGGGACTA-1'), ('UR', 'AATTTACAGC'), ('UY', 'EEEEAEEEEE'), ('UB', 'AATTTACAGC'), ('RG', 'Batch1:0:1:HCWKVBGXB:2')]

This read has start position is outside of the region and its insert size (isize) or fragment size (template_length) are 0. I'm wondering why it is considered to be within the region?

Thanks for any support.

$\endgroup$

1 Answer 1

1
$\begingroup$

Looking at the CIGAR alignment, there is a 167,903bp gap followed by a 266,888bp gap after the starting match point. That second gap spans the region that you have specified.

> 41910661 + 27 + 167903 + 19
[1] 42078610
> 41910661 + 27 + 167903 + 19 + 266888
[1] 42345498

In other words, it's mapping to the region you have identified because a portion of the alignment (between the alignment start point and alignment end point) intersects with the region you have identified.

As a general rule, it's a good idea in bioinformatics to verify that results make sense, as you are doing in asking this question.

This particular match to your target region seems unhelpful, especially given that the part of the alignment that "mapped" is a large skip. If I were to come across such an alignment, I would go back to the mapper (i.e. the thing that created this BAM file) and modify it to be less eager to map over large spans.

$\endgroup$
1
  • $\begingroup$ Thanks for your reply. $\endgroup$
    – Tien
    Commented Apr 25 at 6:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.