Question: I'm having a hard time implementing a calculation of the edit distance from a sam record.
For a bwa mem alignment, I would resort to taking the NM tag (explicitly defined as the edit distance), but this tag is not set by all aligners, so that's where the trouble starts... Alternatively, I would parse the MD tag to get the number of matches and subtract that from the aligned length. But conceptually that doesn't seem equivalent to the MD tag. I scale both to the **aligned** read length.

I should clarify this is for analysis of Oxford Nanopore data, which is dominated by indel errors.

I'll add my python code for calculating both "definitions" of the edit distances. I use pysam, but it's more about the logic than the actual code.

Using the NM tag:

read.get_tag("NM")/read.query_alignment_length
So I:

1. take the NM tag which should be the edit distance
2. divide it by the length of the alignment. This excludes clipped sequences from the calculation.

Using the MD tag:

(read.query_alignment_length - sum([int(item) for item in re.split('[ACTG^]', read.get_tag("MD")) if not item == '']))/read.query_alignment_length
So I:

1. Take the MD tag and sum all sequence 'matches' by removing all non-matching `ACTG^` characters from the string (is this conceptually correct?)
2. I subtract the number of matches from the aligned read length
3. I divide this by the length of the alignment. This excludes clipped sequences from the calculation.

To me, both calculations sound logical, but as a sanity check, I decided to plot them against each other for a sample dataset aligned with bwa mem to see how these metrics correlate. And this result seems okay (correlated) for many reads but deviates quite a lot for a subset. I think my mistake is conceptual or a misunderstanding of those tags. I assume something about how indels are taken into account is tripping me up.

Below I add a (fancy) scatter plot:
(note the pearsonr is quite okay)

![Scatter plot NM vs MD calculation of edit distance][1]

Since this plot is overcrowded by dots here is also the kernel density estimation of the same data:
![kde plot NM vs MD calculation of edit distance][2]

I think ideally I would take query and reference sequence and calculate the Levenshtein distance, but I'm afraid that's not very efficient for processing thousands of reads. Do you have better solutions?

[1]: http://i64.tinypic.com/156eao1.png
[2]: http://i65.tinypic.com/339mu8k.png

Answer: Hi Wouter, first the defintion of NM from [bowtie][1] manual

**NM:i:** The edit distance; that is, the minimal number of **one-nucleotide edits** (substitutions, insertions and deletions) needed to transform the read string into the reference string. Only present if SAM record is for an aligned read.

Note that edit distances are **one-nucleotide edits**.

Now, let's take a little complicated case; one of my reads aligned by BWA

**MD:Z:24^A24C0A0G10^AA39**

If we parse (number sepated from symbols) it will become

24 ^A 24 C 0 A 0 G 10 ^AA 39

If I understood your code correctly (I don't use python, sorry!), then you will calculte MD as

100 - (24 + 24 + 0 + 0 + 10 + 39) = 3

But this will give only the INDEL's edit distance. There are 3 SNPs, (C, A, G) missing from this calculation. This is also evident from the plots that your MD is lesser than NM. In fact, the NM for this is 6, as reported by BWA.

IMO, the correct way to calculate Edit Distance is to calculate the number of transitions while reading the MD from left to right. I'll explain it better:

24 -> ^A : First transition from 24 match to **one** deletion. NM = 1

24 -> C : 2nd transition. A SNP here. NM = 2

0 -> A : Third transition. Another (consecutive) SNP here. NM = 3

0 -> G : Fourth transition. Yet another (consecutive) SNP. NM = 4

10 ^ AA : Deletions of **two** = two transitions => NM = 6

You may also like to see this page regarding the INSERTIONS in MD tag

https://github.com/vsbuffalo/devnotes/wiki/The-MD-Tag-in-BAM-Files

[1]: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml