c709a48585f153d6aab69081cbc5ff3fee7b6b3f angie Mon Sep 26 15:21:34 2011 -0700 Samtools CIGAR strings have two new opcodes: '=' for exactly matchingbases, 'X' for mismatching bases. These are more explicit versions of the 'M' opcode for matching and/or mismatching bases. diff --git src/hg/lib/hgBam.c src/hg/lib/hgBam.c index 36b2a2a..6a8b39b 100644 --- src/hg/lib/hgBam.c +++ src/hg/lib/hgBam.c @@ -58,30 +58,32 @@ if (isRc) reverseComplement(target->dna, target->size); DNA *haystack = target->dna; unsigned int *cigarPacked = bam1_cigar(bam); int tStart = targetOffset, qStart = 0, i; // If isRc, need to go through the CIGAR ops backwards, but sequence offsets still count up. int iStart = isRc ? (core->n_cigar - 1) : 0; int iIncr = isRc ? -1 : 1; for (i = iStart; isRc ? (i >= 0) : (i < core->n_cigar); i += iIncr) { char op; int size = bamUnpackCigarElement(cigarPacked[i], &op); switch (op) { case 'M': // match or mismatch (gapless aligned block) + case '=': // match + case 'X': // mismatch AllocVar(ff); ff->left = ffList; ffList = ff; ff->nStart = needle + qStart; ff->nEnd = ff->nStart + size; ff->hStart = haystack + tStart - targetOffset; ff->hEnd = ff->hStart + size; tStart += size; qStart += size; break; case 'I': // inserted in query case 'S': // skipped query bases at beginning or end ("soft clipping") qStart += size; break; case 'D': // deleted from query