c9941d5c3c0198fea692d72db7fa44a711858598
kent
  Mon Jun 20 13:51:30 2022 -0700
Fixing off by one bug in traceback

diff --git src/hg/oneShot/testLiftRechainer/testLiftRechainer.c src/hg/oneShot/testLiftRechainer/testLiftRechainer.c
index 890052d..4bc2c2d 100644
--- src/hg/oneShot/testLiftRechainer/testLiftRechainer.c
+++ src/hg/oneShot/testLiftRechainer/testLiftRechainer.c
@@ -216,30 +216,46 @@
 struct chain *chain = rechain->chain;
 if (chain == NULL)
     return 0;
 else
     return chain->id;
 }
 
 struct chainPart
 /* A part of a chain */
     {
     struct chainPart *next;
     struct chain *chain;
     int tStart, tEnd;
     };
 
+#ifdef DEBUG
+void showBacks(struct dlList *list, FILE *f)
+/* Show backwards pointers of all chains on list */
+{
+struct dlNode *el;
+for (el = list->head; !dlEnd(el); el = el->next)
+    {
+    struct rechain *rechain = el->val;
+    fprintf(f, "backs for %d:\n", rechainId(rechain));
+    struct crossover *cross;
+    for (cross = rechain->crossList; cross != NULL; cross = cross->next)
+       fprintf(f, "\t%d to %d\n", cross->tPos, rechainId(cross->rechain));
+    }
+}
+#endif /* DEBUG */
+
 void rechainOneTarget(struct chainTarget *target, struct dnaSeq *tSeq, struct hash *qSeqHash,
     struct hash *qRcSeqHash, FILE *f)
 /* Do rechaining on a single target */
 {
 int tSize = tSeq->size;
 
 /* Make rechain lists and put in the "no chain" on on the active list. */
 struct dlList *archiveList = dlListNew();   /* List of rechains no longer active */
 struct dlList *activeList = dlListNew();    /* List of active rechains. */
 rechainZero(tSize, activeList);
 
 struct chain *chain = target->chainList;
 int tPos;
 for (tPos=0; tPos < tSize; ++tPos)
     {
@@ -357,52 +373,54 @@
 	bestChain = rechain;
 	bestScore = rechain->score;
     verbose(1, "active %d: score %lld cross %d last from %d\n", rechainId(rechain), rechain->score, 
 	slCount(rechain->crossList), rechainId(rechain->crossList->rechain));
 	}
     }
 
 /* Do traceback */
 struct chainPart *partList = NULL, *part;
 struct rechain *rechain;
 int tEnd = tSeq->size;
 for (rechain = bestChain; rechain != NULL; )
     {
     // uglyf("rechain id %d tEnd %d\n", rechainId(rechain), tEnd);
     struct rechain *prev = NULL;
-    struct crossover *cross;
+    struct crossover *cross, *nextCross = NULL;
     int tStart = 0;
-    for (cross = rechain->crossList; cross != NULL; cross = cross->next)
+    for (cross = rechain->crossList; cross != NULL; cross = nextCross)
         {
+	nextCross = cross->next;
 	// uglyf("  crossing back to chain %d at %d\n", rechainId(cross->rechain), cross->tPos);
-	if (cross->tPos < tEnd && cross->rechain != rechain)
+	if (cross->tPos <= tEnd)
 	    {
 	    tStart = cross->tPos;
-	    prev = cross->rechain;
+	    if (nextCross != NULL)
+		prev = nextCross->rechain;
 	    rechain->crossList = cross;
 	    break;
 	    }
 	}
     chain = rechain->chain;
     if (chain != NULL)
         {
 	AllocVar(part);
 	part->chain = chain;
 	part->tStart = tStart;
 	part->tEnd = tEnd;
 	slAddHead(&partList, part);
-	uglyf("chainId %d (%d-%d) qSeq (%s:%d-%d) [%s:%d-%d] %d\n", chain->id, chain->tStart, chain->tEnd, chain->qName, chain->qStart, chain->qEnd, chain->tName, tStart, tEnd, tEnd-tStart);
+	// uglyf("chainId %d (%d-%d) qSeq (%s:%d-%d) [%s:%d-%d] %d\n", chain->id, chain->tStart, chain->tEnd, chain->qName, chain->qStart, chain->qEnd, chain->tName, tStart, tEnd, tEnd-tStart);
 	}
     rechain = prev;
     tEnd = tStart;
     }
 for (part = partList; part != NULL; part = part->next)
     {
     struct chain *chainToFree = NULL, *subchain;
     chainSubsetOnT(part->chain, part->tStart, part->tEnd, &subchain,  &chainToFree);
     chainWrite(subchain, f);
     chainFree(&chainToFree);
     }
 slFreeList(&partList);
 }
 
 void testLiftRechainer(char *chainFile, char *tTwoBit, char *qTwoBit, char *out)
@@ -429,27 +447,28 @@
 verbose(1, "Read %lld query bases in %d sequences\n", qTotal, qSeqHash->elCount);
 
 /* Open output now that all input has been at least checked */
 FILE *f = mustOpen(out, "w");
 
 struct twoBitFile *targetTwoBit = twoBitOpen(tTwoBit);
 struct chainTarget *target;
 for (target = targetList; target != NULL; target = target->next)
     {
     int tSize = hashIntVal(tSizeHash, target->name);
     struct dnaSeq *tSeq = twoBitReadSeqFragLower(targetTwoBit, target->name, 0, tSize);
     verbose(1, "Read %d bases in %s\n", tSeq->size, target->name);
     rechainOneTarget(target, tSeq, qSeqHash, qRcSeqHash, f);
     dnaSeqFree(&tSeq);
     }
+
 carefulClose(&f);
 }
 
 int main(int argc, char *argv[])
 /* Process command line. */
 {
 optionInit(&argc, argv, options);
 if (argc != 5)
     usage();
 testLiftRechainer(argv[1], argv[2], argv[3], argv[4]);
 return 0;
 }