196d206a199ee47dc303b5ddf54170f89e4daa75 hiram Thu Nov 3 17:13:25 2016 -0700 correctly counting matches and mis matches refs #13673 diff --git src/lib/bamFile.c src/lib/bamFile.c index beeb94c..001705c 100644 --- src/lib/bamFile.c +++ src/lib/bamFile.c @@ -614,30 +614,34 @@ { char op; int n = bamUnpackCigarElement(cigar[i], &op); switch (op) { case 'X': // mismatch (gapless aligned block) case '=': // match (gapless aligned block) case 'M': // match or mismatch (gapless aligned block) blockSizes[blockCount] = n; qStarts[blockCount] = qPos; tStarts[blockCount] = tPos; blockCount++; tPos += n; qPos += n; qLength += n; + if (op == 'X') + psl->misMatch += n; + else + psl->match += n; break; case 'I': // inserted in query qPos += n; qLength += n; break; case 'D': // deleted from query case 'N': // long deletion from query (intron as opposed to small del) tPos += n; break; case 'S': // skipped query bases at beginning or end ("soft clipping") qPos += n; qLength += n; break; case 'H': // skipped query bases not stored in record's query sequence ("hard clipping") case 'P': // P="silent deletion from padded reference sequence" -- ignore these. @@ -653,18 +657,19 @@ return NULL; } psl->tSize = hdr->target_len[core->tid]; psl->tStart = tStarts[0]; psl->tEnd = tStarts[blockCount-1] + blockSizes[blockCount-1]; psl->qSize = qLength; psl->qStart = qStarts[0]; psl->qEnd = qStarts[blockCount-1] + blockSizes[blockCount-1]; if (isRc) reverseIntRange(&psl->qStart, &psl->qEnd, psl->qSize); psl->blockCount = blockCount; psl->blockSizes = blockSizes; psl->qStarts = qStarts; psl->tStarts = tStarts; +pslComputeInsertCounts(psl); return psl; }