/* align - blatz alignment module. The strategy is to scan the index * for msps in the gapless module, chain the msps, band extend, and * then fill in gaps by doing it again on a more detailed index of * the gap region. */ /* Copyright 2005 Jim Kent. All rights reserved. */ #include "common.h" #include "linefile.h" #include "hash.h" #include "options.h" #include "spacedSeed.h" #include "axt.h" #include "dnautil.h" #include "localmem.h" #include "gapCalc.h" #include "chain.h" #include "chainBlock.h" #include "chainConnect.h" #include "fuzzyFind.h" #include "bandExt.h" #include "boxClump.h" #include "boxLump.h" #include "bzp.h" #include "blatz.h" #include "dynamic.h" // LX static struct chain *chainsCreate(struct gapCalc *gapCalc, struct axtScoreScheme *ss, struct dnaSeq *query, char strand, struct dnaSeq *target, struct cBlock **pBlockList) /* Call chainer on blocks and clean things up a bit. */ { struct chain *chain, *chainList; struct cBlock *block; struct chainConnect cc; /* Set up variables for gap/connect costs. */ ZeroVar(&cc); cc.target = target; cc.query = query; cc.ss = ss; cc.gapCalc = gapCalc; /* We might not have scored the block list properly yet. */ for (block = *pBlockList; block != NULL; block = block->next) { block->score = chainScoreBlock(query->dna + block->qStart, target->dna + block->tStart, block->qEnd - block->qStart, ss->matrix); } /* Call the fancy KD-tree based chainer. */ chainList = chainBlocks(query->name, query->size, strand, target->name, target->size, pBlockList, (ConnectCost)chainConnectCost, (GapCost)chainConnectGapCost, &cc, NULL); /* Clean up overlaps */ for (chain = chainList; chain != NULL; chain = chain->next) { chainRemovePartialOverlaps(chain, query, target, ss->matrix); chainMergeAbutting(chain); } return chainList; } static void insertNewBlocks(struct cBlock *blockList, struct cBlock **pInsertPoint) /* Insert new block list at insertion point. */ { if (blockList != NULL) { struct cBlock *lastInList = slLastEl(blockList); lastInList->next = *pInsertPoint; *pInsertPoint = blockList; } } static struct cBlock *extendInRegion(struct bzp *bzp, struct dnaSeq *query, int qStart, int qEnd, struct dnaSeq *target, int tStart, int tEnd, int dir, char *qSym, char *tSym, int symAlloc) /* Do banded extension in region and return as a list of blocks. */ { int maxExtend = bzp->maxExtend; int symCount; int qSize = qEnd - qStart; int tSize = tEnd - tStart; int qs, ts, qExtStart, tExtStart; if (qSize > 0 && tSize > 0) { if (qSize > maxExtend) qSize = maxExtend; if (tSize > maxExtend) tSize = maxExtend; if (dir < 0) { qStart = qEnd - qSize; tStart = tEnd - tSize; } bandExt(FALSE, bzp->ss, bzp->maxBandGap, query->dna + qStart, qSize, target->dna + tStart, tSize, dir, symAlloc, &symCount, qSym, tSym, &qs, &ts); if (dir > 0) { qExtStart = qStart; tExtStart = tStart; } else { qExtStart = qEnd - qs - 1; tExtStart = tEnd - ts - 1; } return cBlocksFromAliSym(symCount, qSym, tSym, qExtStart, tExtStart); } else return NULL; } static void chainBandExtend(struct bzp *bzp, struct chain *chain, struct dnaSeq *query, struct dnaSeq *target) /* Call banded extender on each block of chain. */ { struct cBlock *oldFirst, *mid, *block, *nextBlock, *newList, *endNewList; int symAlloc = 2 * bzp->maxExtend; char *qSym, *tSym; qSym = needMem(symAlloc); tSym = needMem(symAlloc); #ifdef DEBUG uglyf("chainBandExtend t:%d-%d, q:%d-%d\n", chain->tStart, chain->tEnd, chain->qStart, chain->qEnd); chainWrite(chain, uglyOut); #endif /* DEBUG */ /* Do extension before first block in chain. */ oldFirst = chain->blockList; newList = extendInRegion(bzp, query, 0, chain->qStart, target, 0, chain->tStart, -1, qSym, tSym, symAlloc); #ifdef DEBUG uglyf("got %d new blocks before first\n", slCount(newList)); #endif /* DEBUG */ if (newList != NULL) insertNewBlocks(newList, &chain->blockList); /* Do extensions between blocks. */ for (block = oldFirst; ; block = nextBlock) { nextBlock = block->next; if (nextBlock == NULL) break; /* Extend into beginning of block. */ newList = extendInRegion(bzp, query, block->qEnd, nextBlock->qStart, target, block->tEnd, nextBlock->tStart, 1, qSym, tSym, symAlloc); if (newList != NULL) { endNewList = slLastEl(newList); endNewList->next = block->next; block->next = newList; mid = endNewList; } else mid = block; newList = extendInRegion(bzp, query, block->qEnd, nextBlock->qStart, target, block->tEnd, nextBlock->tStart, -1, qSym, tSym, symAlloc); #ifdef DEBUG uglyf("got %d new blocks in middle\n", slCount(newList)); #endif /* DEBUG */ if (newList != NULL) insertNewBlocks(newList, &mid->next); } /* Do extension after last block in chain. */ newList = extendInRegion(bzp, query, chain->qEnd, query->size, target, chain->tEnd, target->size, 1, qSym, tSym, symAlloc); #ifdef DEBUG if (bzpTimeOn) uglyf("got %d new blocks after last\n", slCount(newList)); #endif /* DEBUG */ if (newList != NULL) insertNewBlocks(newList, &block->next); chainMergeAbutting(chain); chainCalcBounds(chain); /* Clean up. */ freez(&qSym); freez(&tSym); } static void shrinkBlocks(struct chain *chain, struct dnaSeq *query, struct dnaSeq *target, int matrix[256][256]) /* Trim away up to 50 bp from each end of each block in * so that band extension can correct problems introduced * by over-aggressive unbanded extension. */ { struct cBlock *block, *lastBlock = NULL; for (block = chain->blockList; block != NULL; block = block->next) { int size = block->qEnd - block->qStart; int newSize = size - 200; int sizeDiff; int addStart, subEnd; if (newSize < 1) newSize = 1; sizeDiff = size - newSize; addStart = (sizeDiff>>1); subEnd = sizeDiff - addStart; block->qStart += addStart; block->tStart += addStart; block->qEnd -= subEnd; block->tEnd -= subEnd; lastBlock = block; } block = chain->blockList; chain->qStart = block->qStart; chain->tStart = block->tStart; chain->qEnd = lastBlock->qEnd; chain->tEnd = lastBlock->tEnd; } static void bandExtendChains(struct bzp *bzp, struct chain *chainList, struct dnaSeq *query, struct dnaSeq *target) /* Loop through chain list doing banded extensions. */ { struct chain *chain; for (chain = chainList; chain != NULL; chain = chain->next) { shrinkBlocks(chain, query, target, bzp->ss->matrix); chainBandExtend(bzp, chain, query, target); } } static void addBlocksToList(struct cBlock **pList, struct cBlock **toAdd) /* Add blocks in 'toAdd' to pList. This will end up scrambling * orders. */ { struct cBlock *block, *next; for (block = *toAdd; block != NULL; block = next) { next = block->next; slAddHead(pList, block); } *toAdd = NULL; } static void chainsToBlocks(struct chain **pChainList, struct cBlock **pBlockList) /* Convert chains to just a list of blocks. Free chains. */ { struct chain *chain; for (chain = *pChainList; chain != NULL; chain = chain->next) addBlocksToList(pBlockList, &chain->blockList); chainFreeList(pChainList); } static void thresholdChains(struct chain **pChainList, int minScore) /* Remove chains under the threshold from chainList. */ { struct chain *chain, *chainList = NULL, *next; for (chain = *pChainList; chain != NULL; chain = next) { next = chain->next; if (chain->score >= minScore) { slAddHead(&chainList, chain); } else chainFree(&chain); } slReverse(&chainList); *pChainList = chainList; } static void blatzChainAndExtend(struct bzp *bzp, struct dnaSeq *query, char strand, struct dnaSeq *target, struct cBlock **pBlockList, struct chain **retChainList) /* Process msp block list into chains. */ { struct cBlock *blockList = NULL; struct chain *chain, *chainList = NULL, *next; #ifdef DEBUG struct cBlock *block; { for (block = *pBlockList; block != NULL; block = block->next) fprintf(uglyOut, "block: %s %d-%d, %s %d-%d\n", query->name, block->qStart, block->qEnd, target->name, block->tStart, block->tEnd); } #endif /* DEBUG */ chainList = chainsCreate(bzp->cheapGap, bzp->ss, query, strand, target, pBlockList); bzpTime("chainsCreate %d initial chains", slCount(chainList)); if (bzp->bestChainOnly) { if (chainList != NULL) chainFreeList(&chainList->next); } else if (bzp->maxChainsToExplore > 0) { struct chain *lastChain = slElementFromIx(chainList, bzp->maxChainsToExplore-1); if (lastChain != NULL) { verbose(2, "Trimming chain list from %d", slCount(chainList)); chainFreeList(&lastChain->next); verbose(2, " to %d\n", slCount(chainList)); } } thresholdChains(&chainList, bzp->minChain); bzpTime("%d chains passed threshold of %d", slCount(chainList), bzp->minChain); if (bzp->maxExtend > 0) { bandExtendChains(bzp, chainList, query, target); bzpTime("banded extension"); } #ifdef DEBUG { fprintf(uglyOut, "---- Chains after banded extension ----\n"); for (chain = chainList; chain != NULL; chain = chain->next) chainWrite(chain, uglyOut); } #endif /* DEBUG */ chainsToBlocks(&chainList, &blockList); chainList = chainsCreate(bzp->gapCalc, bzp->ss, query, '+', target, &blockList); bzpTime("rechaining with normal gap costs"); #ifdef OLD /* Move to cumulative chain list. */ for (chain = chainList; chain != NULL; chain = next) { next = chain->next; if (chain->score >= bzp->minScore) { slAddHead(retChainList, chain); } else chainFree(&chain); } #endif /* OLD */ /* Move to cumulative chain list. */ for (chain = chainList; chain != NULL; chain = next) { next = chain->next; slAddHead(retChainList, chain); } chainList = NULL; } static struct chain *blatzChainAgainstIndex(struct bzp *bzp, struct blatzIndex *index, struct dnaSeq *query, char strand) /* Align query sequences against indexed target. */ { struct cBlock *blockList = NULL; struct chain *chainList = NULL; /* Chains here. */ /* Collect MSPs. */ bzpTime("Before index"); blatzGaplessScan(bzp, index, query, &blockList); bzpTime("Scanned index"); if (blockList != NULL) blatzChainAndExtend(bzp, query, strand, index->target, &blockList, &chainList); slReverse(&chainList); return chainList; } int uglyCount, uglyTotalQ, uglyTotalT; double uglyArea; static void expandInRegion(struct bzp *bzp, struct dnaSeq *query, int qStart, int qEnd, char strand, struct dnaSeq *target, int tStart, int tEnd, struct cBlock **pBlockList) /* Do alignment on a region inside of query/target governed by bzp parameters. */ { int qSize = qEnd - qStart; int tSize = tEnd - tStart; struct dnaSeq qNew, tNew; struct blatzIndex *index; struct chain *chainList, *chain; #ifdef DEBUG uglyf("----- expandInRegion t:%d-%d, q:%d-%d weight %d -----\n", tStart, tEnd, qStart, qEnd, bzp->weight); #endif /* DEBUG */ { int qSize = qEnd - qStart; int tSize = tEnd - tStart; uglyCount += 1; uglyTotalQ += qSize; uglyTotalT += tSize; uglyArea += (double)qSize * tSize; } /* Create subsequences of query and target. */ ZeroVar(&qNew); ZeroVar(&tNew); qNew.dna = query->dna + qStart; qNew.size = qSize; qNew.name = query->name; tNew.dna = target->dna + tStart; tNew.size = tSize; tNew.name = target->name; index = blatzIndexOne(&tNew, tStart, tEnd, target->size, bzp->weight); chainList = blatzChainAgainstIndex(bzp, index, &qNew, strand); blatzIndexFree(&index); for (chain = chainList; chain != NULL; chain = chain->next) { cBlocksAddOffset(chain->blockList, qStart, tStart); addBlocksToList(pBlockList, &chain->blockList); } chainFreeList(&chainList); } static struct cBlock *focusOnClumps(struct bzp *bzp, struct dnaSeq *query, char strand, struct dnaSeq *target, struct boxClump *clumpList) /* Do alignment process on each clump at parameters more sensitive than * the passed in bzp. Return list of blocks. */ { struct boxClump *clump; struct bzp bzpNew = *bzp; struct cBlock *blockList = NULL; /* Create more sensitive alignment seed parameters (but less aggressive extension). */ if (bzpNew.weight >= 9) bzpNew.weight = bzpNew.weight - 3; else bzpNew.weight = bzpNew.weight - 2; bzpNew.minScore -= 500; bzpNew.minGapless -= 500; bzpNew.minChain -= 500; bzpNew.maxExtend /= 6; bzpNew.maxBandGap /= 6; bzpNew.bestChainOnly = TRUE; for (clump = clumpList; clump != NULL; clump = clump->next) { #ifdef DEBUG uglyf("---- focusing on clump t:%d-%d (%d) q:%d-%d (%d) ----\n", clump->tStart, clump->tEnd, clump->tEnd - clump->tStart, clump->qStart, clump->qEnd, clump->qEnd - clump->qStart); #endif /* DEBUG */ expandInRegion(&bzpNew, query, clump->qStart, clump->qEnd, strand, target, clump->tStart, clump->tEnd, &blockList); } return blockList; } static void addBox(struct boxIn **pList, int qs, int qe, int ts, int te) /* Allocate box and add to list. */ { struct boxIn *box; AllocVar(box); box->qStart = qs; box->qEnd = qe; box->tStart = ts; box->tEnd = te; slAddHead(pList, box); } static void addClippedBox(struct boxIn **pList, int minSize, int winSize, int qStart, int qEnd, int tStart, int tEnd) /* Add box if it is big enough. Split and add if too big. */ { int qSize = qEnd - qStart, tSize = tEnd - tStart; if (qSize >= minSize && tSize >= minSize) { if (qSize > winSize || tSize > winSize) { /* Must split box in two. */ int qe,qs,te,ts; /* Add box at beginning. */ qs = qStart; qe = qStart + winSize; if (qe > qEnd) qe = qEnd; ts = tStart; te = tStart + winSize; if (te > tEnd) te = tEnd; addBox(pList, qs, qe, ts, te); /* Add box at end. */ qe = qEnd; qs = qEnd - winSize; if (qs < qStart) qs = qStart; te = tEnd; ts = tEnd - winSize; if (ts < tStart) ts = tStart; addBox(pList, qs, qe, ts, te); } else addBox(pList, qStart, qEnd, tStart, tEnd); } } struct boxClump *boxFindLumps(struct boxIn **pBoxList); // uglyf static struct boxClump *clumpGapsAndEnds(struct bzp *bzp, struct chain *chainList) /* Given a list of chains (all on same target and query) create a list of * regions where we'll want redo with smaller seeds. Clump these regions * together. */ { struct chain *chain; struct boxIn *boxList = NULL; struct cBlock *b1, *b2; struct boxClump *clumpList = NULL; int winSize = bzp->expandWindow; int minSize = bzp->weight * 2; for (chain = chainList; chain != NULL; chain = chain->next) { if (chain->score >= bzp->minExpand) { int qs, qe, ts, te; qe = chain->qStart; qs = qe - winSize; if (qs < 0) qs = 0; te = chain->tStart; ts = te - winSize; if (ts < 0) ts = 0; addClippedBox(&boxList, minSize, winSize, qs, qe, ts, te); b1 = chain->blockList; for (b2 = b1->next; b2 != NULL; b2 = b2->next) { addClippedBox(&boxList, minSize, winSize, b1->qEnd, b2->qStart, b1->tEnd, b2->tStart); b1 = b2; } qs = chain->qEnd; qe = qs + winSize; if (qe > chain->qSize) qe = chain->qSize; ts = chain->tEnd; te = ts + winSize; if (te > chain->tSize) te = chain->tSize; addClippedBox(&boxList, minSize, winSize, qs, qe, ts, te); } } clumpList = boxLump(&boxList); { /* uglyf */ struct boxClump *clump; int boxCount = 0; int clumpCount = 0; long clumpArea = 0; for (clump = clumpList; clump != NULL; clump = clump->next) { boxCount += clump->boxCount; clumpArea += (clump->qEnd - clump->qStart) * (clump->tEnd - clump->tStart); clumpCount += 1; } bzpTime("%d boxes, %d clumps, %ld clump area", boxCount, clumpCount, clumpArea); } /* uglyf */ return clumpList; } static void removeSimpleOverlaps(struct cBlock **pBlockList) /* Remove from list blocks that overlap and are on same diagonal. */ { struct cBlock *newList = NULL, *b, *next, *last = NULL; slSort(pBlockList, cBlockCmpDiagQuery); for (b = *pBlockList; b != NULL; b = next) { next = b->next; if (last != NULL) { /* Check same diagonal. */ if ((b->tStart - b->qStart) == (last->tStart - last->qStart)) { /* Check for overlap. */ if (last->qEnd >= b->qStart) { /* Fold this block into previous one. */ if (last->qEnd < b->qEnd) last->qEnd = b->qEnd; if (last->tEnd < b->tEnd) last->tEnd = b->tEnd; freez(&b); } } } if (b != NULL) { slAddHead(&newList, b); last = b; } } slReverse(&newList); *pBlockList = newList; } static void findBestGapPos(struct axtScoreScheme *ss, DNA *qStart, int qGap, int qSize, DNA *tStart, int tGap, int tSize, int *retBestPos, int *retBestScore) /* Figure out where to position gap for maximal score. */ { int maxSize = max(qSize, tSize); int gapSize = qGap + tGap; /* One of qGap or tGap is zero */ int matchSize = maxSize - gapSize; int pos, bestPos = 0, rightPos; int bestScore = -BIGNUM; int score; for (pos=0; pos<=matchSize; ++pos) { rightPos = matchSize-pos; score = axtScoreUngapped(ss, qStart, tStart, pos); score += axtScoreUngapped(ss, qStart+qSize-rightPos, tStart+tSize-rightPos, rightPos); if (score > bestScore) { bestScore = score; bestPos = pos; } } *retBestPos = bestPos; *retBestScore = bestScore; } static void reduceGaps(struct bzp *bzp, struct dnaSeq *query, struct dnaSeq *target, struct chain *chain) /* Go through chain looking for small double-sided gaps that * might be turned into single-sided gaps (or even eliminated) * in a way that would improve the score. */ { struct cBlock *b1, *b2; #ifdef DEBUG chain->score = chainCalcScore(chain, bzp->ss, bzp->gapCalc, query, target); // uglyf uglyf("score %f vs %f before reduce gaps\n", chain->score, chainCalcScore(chain, bzp->ss, bzp->gapCalc, query, target) ); #endif /* DEBUG */ b1 = chain->blockList; for (b2 = b1->next; b2 != NULL; b2 = b2->next) { int dq = b2->qStart - b1->qEnd; int dt = b2->tStart - b1->tEnd; if (dq != 0 && dt != 0) { if (dt == dq && dt < 200) { int curScore = -gapCalcCost(bzp->gapCalc, dq, dt); int score = axtScoreUngapped(bzp->ss, query->dna + b1->qEnd, target->dna + b1->tEnd, dt); if (score >= curScore) { /* Merge b2 into b1, pop b2 out of list, and free it. */ b1->qEnd = b2->qEnd; b1->tEnd = b2->tEnd; b1->next = b2->next; freeMem(b2); b2 = b1; chain->score += score - curScore; } } else if (dq < 30 || dt < 30) { int curScore = -gapCalcCost(bzp->gapCalc, dq, dt); int minGap = min(dq,dt); int newDq = dq - minGap; int newDt = dt - minGap; int newGapScore = -gapCalcCost(bzp->gapCalc, newDq, newDt); int matchScore, bestPos; int scoreImprovement; findBestGapPos(bzp->ss, query->dna + b1->qEnd, newDq, dq, target->dna + b1->tEnd , newDt, dt, &bestPos, &matchScore); scoreImprovement = (newGapScore + matchScore) - curScore; #ifdef DEBUG uglyf(" score improvement %d\n", scoreImprovement); #endif /* DEBUG */ if (scoreImprovement >= 0) { b1->qEnd += bestPos; b1->tEnd += bestPos; b2->qStart = b1->qEnd + newDq; b2->tStart = b1->tEnd + newDt; chain->score += scoreImprovement; } } } b1 = b2; } #ifdef DEBUG uglyf("score %f vs %f after reduce gaps\n", chain->score, chainCalcScore(chain, bzp->ss, bzp->gapCalc, query, target) ); #endif /* DEBUG */ } struct ffAli *chainToFfAli(struct chain *chain, struct dnaSeq *query, struct dnaSeq *target) /* Construct ffAli's from chain's blocks. */ { struct ffAli *ff, *ffList = NULL; DNA *q = query->dna, *t = target->dna; struct cBlock *block; for (block = chain->blockList; block != NULL; block = block->next) { AllocVar(ff); ff->nStart = q + block->qStart; ff->nEnd = q + block->qEnd; ff->hStart = t + block->tStart; ff->hEnd = t + block->tEnd; ff->left = ffList; ffList = ff; } ffList = ffMakeRightLinks(ffList); return ffList; } static void slideIntronsInChain(struct chain *chain, struct dnaSeq *query, struct dnaSeq *target) /* Slide around gaps to maximize gt/ag-ness of ends. */ { if (chain->blockList->next != NULL) /* Don't waste time on single-exons */ { /* Call existing code that works on ffAli's. */ struct ffAli *ff, *ffList = chainToFfAli(chain, query, target); struct cBlock *block; DNA *q = query->dna, *t = target->dna; ffSlideIntrons(ffList); for (ff = ffList, block = chain->blockList; ff != NULL; ff = ff->right, block = block->next) { block->qStart = ff->nStart - q; block->qEnd = ff->nEnd - q; block->tStart = ff->hStart - t; block->tEnd = ff->hEnd - t; } ffFreeAli(&ffList); } } static struct chain *blatzAlignOne(struct bzp *bzp, struct blatzIndex *index, struct dnaSeq *query, char strand) /* Align query sequence against indexed target sequence. */ { struct chain *chain, *chainList; struct dnaSeq *target = index->target; struct cBlock *blockList = NULL; /* Do main alignment and chaining. */ chainList = blatzChainAgainstIndex(bzp, index, query, strand); /* Build up list of regions between blocks and before and after * chains and run aligner at more sensitive settings on these regions. */ if (bzp->expandWindow > 0) { struct boxClump *missingClumps = clumpGapsAndEnds(bzp, chainList); bzpTimeOn = FALSE; blockList = focusOnClumps(bzp, query, strand, target, missingClumps); bzpTimeOn = TRUE; boxClumpFreeList(&missingClumps); #ifdef DEBUG uglyf("focus stats: count %d, area %f, tTotal %d, qTotal %d\n", uglyCount, uglyArea, uglyTotalT, uglyTotalQ); #endif /* DEBUG */ bzpTime("expand chains from finer seeds - got %d new blocks", slCount(blockList)); } /* Add blocks from original chains to the mix. */ chainsToBlocks(&chainList, &blockList); /* Create final list of chains. */ removeSimpleOverlaps(&blockList); chainList = chainsCreate(bzp->gapCalc, bzp->ss, query, strand, target, &blockList); thresholdChains(&chainList, bzp->minScore); bzpTime("final chaining"); for (chain = chainList; chain != NULL; chain = chain->next) { reduceGaps(bzp, query, target, chain); if (bzp->rna) slideIntronsInChain(chain, query, target); } #ifdef DEBUG for (chain = chainList; chain != NULL; chain = chain->next) { double oldScore = chain->score; chain->score = chainCalcScore(chain, bzp->ss, bzp->gapCalc, query, target); if (oldScore != chain->score) warn("score inconsistency %f vs %f", oldScore, chain->score); } #endif /* DEBUG */ bzpTime("reduced double to single gaps"); return chainList; } static struct chain *blatzAlignStrand(struct bzp *bzp, struct blatzIndex *indexList, struct dnaSeq *query, char strand) /* Align query sequence against list of indexed target sequences. * The order of chains returned is close to random. */ { struct blatzIndex *index; struct chain *chainList = NULL; int j; // LX /* Loop through indexes, aligning against each, and collecting results * into chainList. */ for (index = indexList; index != NULL; index = index->next) { struct chain *oneList = NULL, *chain, *nextChain; // LX BEG if (bzp->dynaLimitTdynaLimitT; // allocate zeroed memory for hit counters if necessary if (index->counter == NULL) { // index->counter could be set later? AllocArray(index->counter, index->target->size); globalCounter = index->counter; } } if (bzp->dynaLimitQsize); } // LX END oneList = blatzAlignOne(bzp, index, query, strand); // LX BEG // Add the contents of dynaCountQtemp to dynaCountQ if (bzp->dynaLimitQsize;j++) { dynaCountQ[j] = dynaCountQ[j]+dynaCountQtemp[j]; } } freez(&dynaCountQtemp); // LX END for (chain = oneList; chain != NULL; chain = nextChain) { nextChain = chain->next; slAddHead(&chainList, chain); } } return chainList; } struct chain *blatzAlign(struct bzp *bzp, struct blatzIndex *indexList, struct dnaSeq *query) /* Align both strands of query against index using parameters in bzp. * Return chains sorted by score (highest scoring first) */ { struct chain *chainList, *plusChains, *minusChains; plusChains = blatzAlignStrand(bzp, indexList, query, '+'); reverseComplement(query->dna, query->size); minusChains = blatzAlignStrand(bzp, indexList, query, '-'); reverseComplement(query->dna, query->size); chainList = slCat(plusChains, minusChains); slSort(&chainList, chainCmpScore); return chainList; }