Just copied over the bug report received to make sure I'll handle it later on:
. An example using
columns 6,10,17,18 of a SAM record created by MALT with the custom ZI
tag showing the % identity integer as calculated by MALT:
37M4D8M CTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA ZI:i:85
MD:Z:G36^CCTTTC6
My manual reconstruction of the reference sequence from
CIGAR and MD:
CTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA----AAAAAAAA read
|||||||||||||||||||||||||||||||||||||----|||||||| CIGAR
G||||||||||||||||||||||||||||||||||||CCTTTC|||||| MD
GTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACCTTTCAAAAAA reference
The problem with the MD tag is the ambiguous ^CCTTTC string. It could be
interpreted as deletion of 6 bases although the CIGAR string clearly
shows only 4 deletions. SAMtools circumvent this ambiguity by including
another 0 (zero) between the 4 deleted characters and the two SNPs that
directly follow it:
before: ^CCTTTC
after: ^CCTT0TC
as this clearly means "deletion of CCTT, then 0 matches, then mismatch
T, then mismatch C". However this is not specified in the SAM standard
so not all programs will include the 0 divider. MALT does not but rather
expects you to interpret the CIGAR string first before checking the MD
tag for further details.
An example by someone else having the same problem:
http://seqanswers.com/forums/showthread.php?t=8978
From your code, line 315:
MDlist=re.findall('(\d+|\D+)',MD)
splits the MD tag into contiguous chunks of digits or non-digits. This
would consider the entire ^CCTTTC as one element, ignoring the fact that
according to the CIGAR only the first 4 are deleted.
From your code, lines 328-331:
elif '^' in e:
ef=e.lstrip('^')
alignment += ef
continue
will then add the six characters to the alignment, thus shifting the
entire alignment by two positions?
CTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA------AAAAAAAA read
GTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACCTTTCAAAAAA reference
I think this may also explain why my modified divergence filter raised
an exception on different lengths of read and reference. However I don't
fully understand your method to recreate the reference from CIGAR and MD
so I might be wrong.
This has implications for all processed SAM files produced by tools
which do not use a 0 divider in their MD tag. Could you please check if
my assumptions are correct and let me know how to proceed? I'll try to
fix the bug by myself, probably by inverting the parsing order of CIGAR
and MD, if I don't hear from you.
Just copied over the bug report received to make sure I'll handle it later on:
. An example using
columns 6,10,17,18 of a SAM record created by MALT with the custom ZI
tag showing the % identity integer as calculated by MALT:
37M4D8M CTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA ZI:i:85
MD:Z:G36^CCTTTC6
My manual reconstruction of the reference sequence from
CIGAR and MD:
CTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA----AAAAAAAA read
|||||||||||||||||||||||||||||||||||||----|||||||| CIGAR
G||||||||||||||||||||||||||||||||||||CCTTTC|||||| MD
GTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACCTTTCAAAAAA reference
The problem with the MD tag is the ambiguous ^CCTTTC string. It could be
interpreted as deletion of 6 bases although the CIGAR string clearly
shows only 4 deletions. SAMtools circumvent this ambiguity by including
another 0 (zero) between the 4 deleted characters and the two SNPs that
directly follow it:
before: ^CCTTTC
after: ^CCTT0TC
as this clearly means "deletion of CCTT, then 0 matches, then mismatch
T, then mismatch C". However this is not specified in the SAM standard
so not all programs will include the 0 divider. MALT does not but rather
expects you to interpret the CIGAR string first before checking the MD
tag for further details.
An example by someone else having the same problem:
http://seqanswers.com/forums/showthread.php?t=8978
From your code, line 315:
MDlist=re.findall('(\d+|\D+)',MD)
splits the MD tag into contiguous chunks of digits or non-digits. This
would consider the entire ^CCTTTC as one element, ignoring the fact that
according to the CIGAR only the first 4 are deleted.
From your code, lines 328-331:
elif '^' in e:
ef=e.lstrip('^')
alignment += ef
continue
will then add the six characters to the alignment, thus shifting the
entire alignment by two positions?
CTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA------AAAAAAAA read
GTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACCTTTCAAAAAA reference
I think this may also explain why my modified divergence filter raised
an exception on different lengths of read and reference. However I don't
fully understand your method to recreate the reference from CIGAR and MD
so I might be wrong.
This has implications for all processed SAM files produced by tools
which do not use a 0 divider in their MD tag. Could you please check if
my assumptions are correct and let me know how to proceed? I'll try to
fix the bug by myself, probably by inverting the parsing order of CIGAR
and MD, if I don't hear from you.