a44421a79fb36cc2036fe116b97ea3bc9590cd0c braney Fri Dec 2 09:34:39 2011 -0800 removed rcsid (#295) diff --git src/lib/chainConnect.c src/lib/chainConnect.c index ce166df..3e513dc 100644 --- src/lib/chainConnect.c +++ src/lib/chainConnect.c @@ -1,367 +1,366 @@ /* chainConnect - Stuff to cope with connect costs and overlaps when * making chains. This works closely with the chainBlock module. */ #include "common.h" #include "chain.h" #include "axt.h" #include "gapCalc.h" #include "chainConnect.h" -static char const rcsid[] = "$Id: chainConnect.c,v 1.3 2005/08/07 23:25:26 baertsch Exp $"; double chainScoreBlock(char *q, char *t, int size, int matrix[256][256]) /* Score block through matrix. */ { double score = 0; int i; for (i=0; i<size; ++i) score += matrix[(int)q[i]][(int)t[i]]; return score; } double chainCalcScore(struct chain *chain, struct axtScoreScheme *ss, struct gapCalc *gapCalc, struct dnaSeq *query, struct dnaSeq *target) /* Calculate chain score freshly. */ { struct cBlock *b1, *b2; double score = 0; for (b1 = chain->blockList; b1 != NULL; b1 = b2) { score += chainScoreBlock(query->dna + b1->qStart, target->dna + b1->tStart, b1->tEnd - b1->tStart, ss->matrix); b2 = b1->next; if (b2 != NULL) score -= gapCalcCost(gapCalc, b2->qStart - b1->qEnd, b2->tStart - b1->tEnd); } return score; } double chainCalcScoreSubChain(struct chain *chain, struct axtScoreScheme *ss, struct gapCalc *gapCalc, struct dnaSeq *query, struct dnaSeq *target) /* Calculate chain score assuming query and target span the chained region rather than entire chrom. */ { struct cBlock *b1, *b2; double score = 0; for (b1 = chain->blockList; b1 != NULL; b1 = b2) { score += chainScoreBlock(query->dna + (b1->qStart) - chain->qStart, target->dna + (b1->tStart) - chain->tStart, b1->tEnd - b1->tStart, ss->matrix); b2 = b1->next; if (b2 != NULL) score -= gapCalcCost(gapCalc, b2->qStart - b1->qEnd, b2->tStart - b1->tEnd); } return score; } void cBlockFindCrossover(struct cBlock *left, struct cBlock *right, struct dnaSeq *qSeq, struct dnaSeq *tSeq, int overlap, int matrix[256][256], int *retPos, int *retScoreAdjustment) /* Find ideal crossover point of overlapping blocks. That is * the point where we should start using the right block rather * than the left block. This point is an offset from the start * of the overlapping region (which is the same as the start of the * right block). */ { int bestPos = 0; char *rqStart = qSeq->dna + right->qStart; char *lqStart = qSeq->dna + left->qEnd - overlap; char *rtStart = tSeq->dna + right->tStart; char *ltStart = tSeq->dna + left->tEnd - overlap; int i; double score, bestScore, rScore, lScore; /* Make sure overlap is not larger than either block size: */ if (overlap > (left->tEnd - left->tStart) || overlap > (right->tEnd - right->tStart)) errAbort("overlap is %d -- too large for one of these:\n" "qSize=%d tSize=%d\n" "left: qStart=%d qEnd=%d (%d) tStart=%d tEnd=%d (%d)\n" "right: qStart=%d qEnd=%d (%d) tStart=%d tEnd=%d (%d)\n", overlap, qSeq->size, tSeq->size, left->qStart, left->qEnd, left->qEnd - left->qStart, left->tStart, left->tEnd, left->tEnd - left->tStart, right->qStart, right->qEnd, right->qEnd - right->qStart, right->tStart, right->tEnd, right->tEnd - right->tStart); score = bestScore = rScore = chainScoreBlock(rqStart, rtStart, overlap, matrix); lScore = chainScoreBlock(lqStart, ltStart, overlap, matrix); for (i=0; i<overlap; ++i) { score += matrix[(int)lqStart[i]][(int)ltStart[i]]; score -= matrix[(int)rqStart[i]][(int)rtStart[i]]; if (score > bestScore) { bestScore = score; bestPos = i+1; } } *retPos = bestPos; *retScoreAdjustment = rScore + lScore - bestScore; } int chainConnectGapCost(int dq, int dt, struct chainConnect *cc) /* Calculate cost of non-overlapping gap. */ { return gapCalcCost(cc->gapCalc, dq, dt); } int chainConnectCost(struct cBlock *a, struct cBlock *b, struct chainConnect *cc) /* Calculate connection cost - including gap score * and overlap adjustments if any. */ { int dq = b->qStart - a->qEnd; int dt = b->tStart - a->tEnd; int overlapAdjustment = 0; if (a->qStart >= b->qStart || a->tStart >= b->tStart) { errAbort("a (%d %d) not strictly before b (%d %d)", a->qStart, a->tStart, b->qStart, b->tStart); } if (dq < 0 || dt < 0) { int bSize = b->qEnd - b->qStart; int aSize = a->qEnd - a->qStart; int overlap = -min(dq, dt); int crossover; if (overlap >= bSize || overlap >= aSize) { /* One of the blocks is enclosed completely on one dimension * or other by the other. Discourage this overlap. */ overlapAdjustment = 100000000; } else { cBlockFindCrossover(a, b, cc->query, cc->target, overlap, cc->ss->matrix, &crossover, &overlapAdjustment); dq += overlap; dt += overlap; } } return overlapAdjustment + gapCalcCost(cc->gapCalc, dq, dt); } static void checkStartBeforeEnd(struct chain *chain, char *message) /* Check qStart < qEnd, tStart < tEnd for each block. */ { struct cBlock *b; for (b = chain->blockList; b != NULL; b = b->next) { if (b->qStart >= b->qEnd || b->tStart >= b->tEnd) errAbort("Start after end in (%d %d) to (%d %d) %s", b->qStart, b->tStart, b->qEnd, b->tEnd, message); } } static void checkChainGaps(struct chain *chain, char *message) /* Check that gaps between blocks are non-negative. */ { struct cBlock *a, *b; a = chain->blockList; if (a == NULL) return; for (b = a->next; b != NULL; b = b->next) { if (a->qEnd > b->qStart || a->tEnd > b->tStart) { errAbort("Negative gap between (%d %d - %d %d) and (%d %d - %d %d) %s", a->qStart, a->tStart, a->qEnd, a->tEnd, b->qStart, b->tStart, b->qEnd, b->tEnd, message); } a = b; } } static void checkChainIncreases(struct chain *chain, char *message) /* Check that qStart and tStart both strictly increase * in chain->blockList. */ { struct cBlock *a, *b; a = chain->blockList; if (a == NULL) return; for (b = a->next; b != NULL; b = b->next) { if (a->qStart >= b->qStart || a->tStart >= b->tStart) { errAbort("a (%d %d) not before b (%d %d) %s", a->qStart, a->tStart, b->qStart, b->tStart, message); } a = b; } } void chainCalcBounds(struct chain *chain) /* Recalculate chain boundaries - setting qStart/qEnd/tStart/tEnd from * a scan of blockList. */ { struct cBlock *b = chain->blockList; if (chain->blockList == NULL) return; chain->qStart = b->qStart; chain->tStart = b->tStart; b = slLastEl(chain->blockList); chain->qEnd = b->qEnd; chain->tEnd = b->tEnd; } static boolean removeNegativeBlocks(struct chain *chain) /* Removing the partial overlaps occassional results * in all of a block being removed. This routine * removes the dried up husks of these blocks * and returns TRUE if it finds any. */ { struct cBlock *newList = NULL, *b, *next; boolean gotNeg = FALSE; for (b = chain->blockList; b != NULL; b = next) { next = b->next; if (b->qStart >= b->qEnd || b->tStart >= b->tEnd) { gotNeg = TRUE; freeMem(b); } else { slAddHead(&newList, b); } } slReverse(&newList); chain->blockList = newList; if (gotNeg) chainCalcBounds(chain); return gotNeg; } void setChainBounds(struct chain *chain) /* Set chain overall bounds to fit blocks. */ { struct cBlock *b = chain->blockList; chain->qStart = b->qStart; chain->tStart = b->tStart; while (b->next != NULL) b = b->next; chain->qEnd = b->qEnd; chain->tEnd = b->tEnd; } void chainRemovePartialOverlaps(struct chain *chain, struct dnaSeq *qSeq, struct dnaSeq *tSeq, int matrix[256][256]) /* If adjacent blocks overlap then find crossover points between them. */ { struct cBlock *a, *b; boolean totalTrimA, totalTrimB; assert(chain->blockList != NULL); /* Do an internal sanity check to make sure that things * really are sorted by both qStart and tStart. */ checkChainIncreases(chain, "before removePartialOverlaps"); /* Remove overlapping portion of blocks. In some * tricky repeating regions this can result in * complete blocks being removed. This complicates * the loop below in some cases, forcing us to essentially * start over when the first of the two blocks we * are examining gets trimmed out completely. */ for (;;) { totalTrimA = totalTrimB = FALSE; a = chain->blockList; b = a->next; for (;;) { int dq, dt; if (b == NULL) break; dq = b->qStart - a->qEnd; dt = b->tStart - a->tEnd; if (dq < 0 || dt < 0) { int overlap = -min(dq, dt); int aSize = a->qEnd - a->qStart; int bSize = b->qEnd - b->qStart; int crossover, invCross, overlapAdjustment; if (overlap >= aSize || overlap >= bSize) { totalTrimB = TRUE; } else { cBlockFindCrossover(a, b, qSeq, tSeq, overlap, matrix, &crossover, &overlapAdjustment); b->qStart += crossover; b->tStart += crossover; invCross = overlap - crossover; a->qEnd -= invCross; a->tEnd -= invCross; if (b->qEnd <= b->qStart) { totalTrimB = TRUE; } else if (a->qEnd <= a->qStart) { totalTrimA = TRUE; } } } if (totalTrimA == TRUE) { removeNegativeBlocks(chain); break; } else if (totalTrimB == TRUE) { b = b->next; freez(&a->next); a->next = b; totalTrimB = FALSE; } else { a = b; b = b->next; } } if (!totalTrimA) break; } /* Reset chain bounds - may have clipped them in this * process. */ setChainBounds(chain); /* Do internal sanity checks. */ checkChainGaps(chain, "after removePartialOverlaps"); checkStartBeforeEnd(chain, "after removePartialOverlaps"); } void chainMergeAbutting(struct chain *chain) /* Merge together blocks in a chain that abut each * other exactly. */ { struct cBlock *newList = NULL, *b, *last = NULL, *next; for (b = chain->blockList; b != NULL; b = next) { next = b->next; if (last == NULL || last->qEnd != b->qStart || last->tEnd != b->tStart) { slAddHead(&newList, b); last = b; } else { last->qEnd = b->qEnd; last->tEnd = b->tEnd; freeMem(b); } } slReverse(&newList); chain->blockList = newList; }