diff --git a/src/alignment.cpp b/src/alignment.cpp index 49a39e0787..b9dce5a260 100644 --- a/src/alignment.cpp +++ b/src/alignment.cpp @@ -1292,14 +1292,20 @@ pair compute_template_lengths(const int64_t& pos1, const vecto int64_t here = pos; for (auto& item : cigar) { // Trace along the cigar - if (item.second == 'M') { - // Bases are matched. Count them in the bounds and execute the operation - low = min(low, here); - here += item.first; - high = max(high, here); - } else if (item.second == 'D') { - // Only other way to advance in the reference - here += item.first; + switch (item.second) { + case 'M': + case 'X': + case '=': + // Bases are matched or mismatched. Count them in the bounds and execute the operation + low = min(low, here); + high = max(high, here + item.first); + // Fall through. + case 'D': + case 'N': + // Even for deletions/gaps, we advance on the reference. + here += item.first; + break; + // Inserts and various clips/paddings don't do anything. } } @@ -1308,16 +1314,12 @@ pair compute_template_lengths(const int64_t& pos1, const vecto auto bounds1 = find_bounds(pos1, cigar1); auto bounds2 = find_bounds(pos2, cigar2); - - // Compute the separation - int32_t dist = 0; - if (bounds1.first < bounds2.second) { - // The reads are in order - dist = bounds2.second - bounds1.first; - } else if (bounds2.first < bounds1.second) { - // The reads are out of order so the other bounds apply - dist = bounds1.second - bounds2.first; - } + + // We wanted the distance between the outermost points. So find them. + auto min_start = std::min(bounds1.first, bounds2.first); + auto max_end = std::max(bounds1.second, bounds2.second); + // And find the distance + int32_t dist = max_end - min_start; if (pos1 < pos2) { // Count read 1 as the overall "leftmost", so its value will be positive