| Subject: | behaviour of qpos in deletion regions |
the current behaviour of qpos when entering a deletion in a read is to
keep its last value.
Using the SNP caller example code from the POD of Bio::DB::Sam, take
this example of a read with one deletion at base 10 (Perl debugger output):
x $b->padded_alignment
0
'GTTTTTGAAGACGACCTTCGATTTTTTTAATTTTATTTTCATTAACAGCAAAATCTCTGAAATCGACTATATATACTTTTAATTTATTCAATAAATGCAGATATGAATTGCATTCAAAAATATTTTCATAAAGGGATCTATATTTATT'
1 '|||||||||
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||'
2
'GTTTTTGAA-ACGACCTTCGATTTTTTTAATTTTATTTTCATTAACAGCAAAATCTCTGAAATCGACTATATATACTTTTAATTTATTCAATAAATGCAGATATGAATTGCATTCAAAAATATTTTCATAAAGGGATCTATATTTATT'
traversing the pileups from the position before the deletion (continuing
with the example SNP caller code), $pileup->qpos returns the following
values on each call of the callback sub for the above read:
8, 9, 9, 10, 11 ....
So it stays at 9 where the read has the deletion.
A problem arises when one uses the code in your example to get to the
base of the read
substr( $b->qseq, $pileup->qpos, 1 )
which would return an A inside the deletion where the read has no base.
I know that I could use is_del to check for this case and skip it but I
think it would be safer for qpos to return undef where the read has no
base, which would prevent incorrect bases to be returned and would
correctly reflect the fact that the query has no position that matches
the current reference position.
Would also be great if the example callback code could show how to
properly work with indels as this is not entirely intuitive I think.
Thanks for looking into this!