From 12ce97891423c2b80a148d2a57452d4d44e0fc77 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 27 Sep 2026 19:52:06 -0400 Subject: [PATCH] Faster, exact forward-backward; in-process seed anchoring; array-based progressive MSA Forward-backward. The posterior probabilities (getPosteriorProbsWithBanding with diagonalCalculationPosteriorMatchProbs or ...PosteriorProbs) now go through a probability-space version of the same banded, chunked recursion. Each diagonal keeps a log scale that is moved only when its values drift out of [1e-100,1e100], the state machine is flattened into arrays once and cached on the StateMachine, and each state of a diagonal is one vectorised loop, with AVX2 gathers for the emissions. The band, the traceback points and which diagonals get posteriors in which chunk are exactly those of the log-space code, which remains for EM. It matches an exact reference to 1e-7; the old logAdd, a piecewise cubic that drops terms under e^-7.5, was off by up to ~0.01. The chunks are planned from the band before any dp, so buffers are sized once, and a chunk whose forward values would pass 16MB keeps only a pair of diagonals before each block and recomputes the block on the way back: an unanchored 10kb x 10kb pair drops from 3.9GB to 116MB. 5-11x faster per pairwise alignment. Seed anchoring (anchorMethod = PAIRWISE_ANCHOR_SEED): lastz's recipe -- 12of19 spaced seeds allowing a transition, ungapped HOXD70 x-drop extension, the heaviest colinear chain of the HSPs -- in process, recursing into the gaps with a threshold that falls with their area. MUM anchors accept "unique" matches of a dozen bases; at 0.3 substitutions per site per branch half of them were off the true alignment, and the banded dp followed them. In simulation at that divergence the MSA goes from 0.68/0.32 precision/recall with MUMs to 0.79/0.48. The default is unchanged (PAIRWISE_ANCHOR_LEGACY); cactus selects seed. Progressive MSA. getMultipleSequenceAlignmentProgressive and filterMultipleAlignedPairs work over arrays indexed by (sequence, position) instead of hashes and AVL trees, with the weights in one block, and no longer maintain the weight-ordered set only the greedy method reads. Same alignment, same order; 2.3x faster and half the memory on 6 x 100kb. column_hashFn no longer sums sequence and position, which made every bucket hold one of each sequence. Also: mum_destruct is a loop rather than a recursion as deep as the MUM chain; a zero-length HSP cannot enter the chain; threshold 0 reports every cell as the log path does. Tests: the fast path against the log-space one (random bands, machines, ragged ends, and unanchored pairs big enough to be recomputed in blocks) and seed anchor validity. test_filterToRemoveOverlap compared each pair with itself, so the set it expected was always empty; it now skips that. Co-Authored-By: Claude Opus 5.5 --- impl/multipleAligner.c | 439 ++++++++++++- impl/pairwiseAligner.c | 1213 ++++++++++++++++++++++++++++++++++- impl/stateMachine.c | 3 + inc/pairwiseAligner.h | 15 + inc/stateMachine.h | 4 + tests/pairwiseAlignerTest.c | 121 ++++ 6 files changed, 1758 insertions(+), 37 deletions(-) diff --git a/impl/multipleAligner.c b/impl/multipleAligner.c index 5829a86..7baacb7 100644 --- a/impl/multipleAligner.c +++ b/impl/multipleAligner.c @@ -63,7 +63,13 @@ int column_cmp(Column *c, Column *c2) { } uint64_t column_hashFn(Column *c) { - return c->seqName + c->position; + // Mixed, not summed: seqName + position put every sequence's positions on top of one another, + // so each bucket held a position from every sequence, and lookups here were a fifth of the + // runtime of a multiple alignment. + uint64_t h = (uint64_t) c->seqName * 0x9E3779B97F4A7C15ULL ^ (uint64_t) c->position; + h ^= h >> 31; + h *= 0xBF58476D1CE4E5B9ULL; + return h ^ (h >> 29); } int column_equalsFn(Column *c, Column *c2) { @@ -190,16 +196,20 @@ stSortedSet *makeOrderedSetOfAlignmentWeights(stHash *alignmentWeightAdjLists) { static void removeAlignmentFromSortedAlignmentWeights(AlignmentWeight *aW, stSortedSet *alignmentWeightsOrderedByWeight) { /* - * Removes the weight from the ordered set. + * Removes the weight from the ordered set, if one is being kept (the progressive merge does not need it). */ - stSortedSet_remove(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight); + if (alignmentWeightsOrderedByWeight != NULL) { + stSortedSet_remove(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight); + } } static void insertAlignmentIntoSortedAlignmentWeights(AlignmentWeight *aW, stSortedSet *alignmentWeightsOrderedByWeight) { /* - * Adds weight to the ordered set. + * Adds weight to the ordered set, if one is being kept. */ - stSortedSet_insert(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight); + if (alignmentWeightsOrderedByWeight != NULL) { + stSortedSet_insert(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight); + } } static void alignmentWeight_destruct(AlignmentWeight *aW) { @@ -509,18 +519,359 @@ stList *makeColumnSequences(stList *seqFrags, stSet *columns) { return columnSequences; } +/* + * The progressive merge, over arrays. + * + * The functions above keep a column's alignment weights in an AVL tree, find a column's tree + * through a hash, and find a column from its sequence position through another; the merge + * below spent most of its time in those lookups, and most of its memory -- over 200 bytes for + * each aligned pair -- on the trees' nodes and the weights' own allocations. Here every column + * position has an id, (sequence, position) numbered in order, and the columns, their weights + * and the scratch the merge needs are arrays indexed by it. The weights, two per aligned pair, + * sit in one block. The merge makes the same choices in the same order with the same arithmetic + * as the functions above, which remain for the greedy method and the tests, and gives the same + * alignment. + */ + +typedef struct _weightList { // a column's weights, in no particular order + AlignmentWeight **items; + int64_t length, capacity; + bool owned; // items was allocated for this list, rather than being its slice of the shared block +} WeightList; + +typedef struct _progressiveMsa { + int64_t *seqOffsets; // the id of (seq, pos) is seqOffsets[seq] + pos + int64_t columnNumber; + Column **columns; // by id + bool *isHead; // by id: is the column's first position, and so the column, in the alignment + WeightList *weights; // by id, the weights of the column headed there + AlignmentWeight *weightPool; // the weights, two per aligned pair + AlignmentWeight **itemBlock; // the weight lists' first slices + AlignmentWeight **mark; // by id, scratch for merging: the weight to the column headed there, if any; else NULL + int64_t *yIndex, *yTag, tag; // by id, scratch for aligning two column-sequences: the index of the column in Y, if yTag is tag +} ProgressiveMsa; + +static inline int64_t progressiveMsa_id(ProgressiveMsa *m, Column *c) { + return m->seqOffsets[c->seqName] + c->position; +} + +static inline WeightList *progressiveMsa_weights(ProgressiveMsa *m, Column *c) { + return &m->weights[progressiveMsa_id(m, c)]; +} + +static void weightList_append(WeightList *l, AlignmentWeight *aW) { + if (l->length == l->capacity) { + int64_t capacity = l->capacity * 2 + 4; + AlignmentWeight **items = st_malloc(sizeof(AlignmentWeight *) * capacity); + memcpy(items, l->items, sizeof(AlignmentWeight *) * l->length); + if (l->owned) { + free(l->items); + } + l->items = items; + l->capacity = capacity; + l->owned = 1; + } + l->items[l->length++] = aW; +} + +static void weightList_remove(WeightList *l, AlignmentWeight *aW) { + for (int64_t i = 0; i < l->length; i++) { + if (l->items[i] == aW) { + l->items[i] = l->items[--l->length]; + return; + } + } + assert(0); +} + +static ProgressiveMsa *progressiveMsa_construct(stList *seqFrags, stList *multipleAlignedPairs) { + ProgressiveMsa *m = st_calloc(1, sizeof(ProgressiveMsa)); + int64_t seqNumber = stList_length(seqFrags); + m->seqOffsets = st_malloc(sizeof(int64_t) * (seqNumber + 1)); + m->columnNumber = 0; + for (int64_t seq = 0; seq < seqNumber; seq++) { + m->seqOffsets[seq] = m->columnNumber; + m->columnNumber += ((SeqFrag *) stList_get(seqFrags, seq))->length; + } + m->seqOffsets[seqNumber] = m->columnNumber; + int64_t n = m->columnNumber; + // The columns, one per position, as makeColumns makes them + m->columns = st_malloc(sizeof(Column *) * (n + 1)); + m->isHead = st_malloc(sizeof(bool) * (n + 1)); + for (int64_t seq = 0; seq < seqNumber; seq++) { + for (int64_t pos = 0; pos < m->seqOffsets[seq + 1] - m->seqOffsets[seq]; pos++) { + Column *c = st_malloc(sizeof(Column)); + c->seqName = seq; + c->position = pos; + c->nColumn = NULL; + m->columns[m->seqOffsets[seq] + pos] = c; + m->isHead[m->seqOffsets[seq] + pos] = 1; + } + } + // The weights, made in the order, and with the random tie breaking, that makeAlignmentWeightAdjacencyLists uses + int64_t pairNumber = stList_length(multipleAlignedPairs); + m->weights = st_calloc(n + 1, sizeof(WeightList)); + for (int64_t i = 0; i < pairNumber; i++) { // each list's size, so each can have its own slice of one block + stIntTuple *aP = stList_get(multipleAlignedPairs, i); + m->weights[m->seqOffsets[stIntTuple_get(aP, 1)] + stIntTuple_get(aP, 2)].capacity++; + m->weights[m->seqOffsets[stIntTuple_get(aP, 3)] + stIntTuple_get(aP, 4)].capacity++; + } + m->itemBlock = st_malloc(sizeof(AlignmentWeight *) * (2 * pairNumber + 1)); + for (int64_t i = 0, offset = 0; i < n; i++) { + m->weights[i].items = m->itemBlock + offset; + offset += m->weights[i].capacity; + } + m->weightPool = st_malloc(sizeof(AlignmentWeight) * (2 * pairNumber + 1)); + for (int64_t i = 0; i < pairNumber; i++) { + /*Tuple of score, seq1, pos1, seq2, pos2 */ + stIntTuple *aP = stList_get(multipleAlignedPairs, i); + assert(stIntTuple_length(aP) == 5); + int64_t score = stIntTuple_get(aP, 0); + int64_t id1 = m->seqOffsets[stIntTuple_get(aP, 1)] + stIntTuple_get(aP, 2); + int64_t id2 = m->seqOffsets[stIntTuple_get(aP, 3)] + stIntTuple_get(aP, 4); + AlignmentWeight *aW = &m->weightPool[2 * i], *rW = &m->weightPool[2 * i + 1]; + aW->column = m->columns[id1]; + aW->numberOfWeights = 1; + aW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; // as makeAlignmentWeight + rW->column = m->columns[id2]; + rW->numberOfWeights = 1; + rW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; + aW->rWeight = rW; + rW->rWeight = aW; + // A weight lives in the list of the column at the other end of it, as insertWeight puts it + m->weights[id2].items[m->weights[id2].length++] = aW; + m->weights[id1].items[m->weights[id1].length++] = rW; + } + m->mark = st_calloc(n + 1, sizeof(AlignmentWeight *)); + m->yIndex = st_malloc(sizeof(int64_t) * (n + 1)); + m->yTag = st_calloc(n + 1, sizeof(int64_t)); + m->tag = 0; + return m; +} + +static void progressiveMsa_destruct(ProgressiveMsa *m) { + for (int64_t i = 0; i < m->columnNumber; i++) { + if (m->weights[i].owned) { + free(m->weights[i].items); + } + } + free(m->weights); + free(m->itemBlock); + free(m->weightPool); + free(m->seqOffsets); + free(m->columns); // the Columns themselves belong to whoever has the set of them + free(m->isHead); + free(m->mark); + free(m->yIndex); + free(m->yTag); + free(m); +} + +/* + * As mergeColumns. + */ +static Column *progressiveMsa_merge(ProgressiveMsa *m, AlignmentWeight *aW) { + if (progressiveMsa_weights(m, aW->column)->length < progressiveMsa_weights(m, aW->rWeight->column)->length) { + aW = aW->rWeight; + } + Column *c1 = aW->column, *c2 = aW->rWeight->column; + assert(c1 != c2); + WeightList *l1 = progressiveMsa_weights(m, c1), *l2 = progressiveMsa_weights(m, c2); + assert(l1->length >= l2->length); + // Merge the columns + Column *c = c1; + while (c->nColumn != NULL) { + c = c->nColumn; + } + c->nColumn = c2; + m->isHead[progressiveMsa_id(m, c2)] = 0; + // Drop the merging weight + weightList_remove(l1, aW->rWeight); + weightList_remove(l2, aW); + // Merge each of c2's remaining weights into c1's to the same column, or move it across if c1 has none + for (int64_t i = 0; i < l1->length; i++) { + m->mark[progressiveMsa_id(m, l1->items[i]->column)] = l1->items[i]; + } + for (int64_t i = 0; i < l2->length; i++) { + AlignmentWeight *aW2 = l2->items[i]; + AlignmentWeight *aW1 = m->mark[progressiveMsa_id(m, aW2->column)]; + if (aW1 != NULL) { + weightList_remove(progressiveMsa_weights(m, aW2->column), aW2->rWeight); + aW1->avgWeight = ((aW1->avgWeight * aW1->numberOfWeights) + (aW2->avgWeight * aW2->numberOfWeights)) / + (aW1->numberOfWeights + aW2->numberOfWeights); + aW1->numberOfWeights += aW2->numberOfWeights; + aW1->rWeight->avgWeight = aW1->avgWeight; + aW1->rWeight->numberOfWeights = aW1->numberOfWeights; + } else { + aW2->rWeight->column = c1; + weightList_append(l1, aW2); + } + } + for (int64_t i = 0; i < l1->length; i++) { + m->mark[progressiveMsa_id(m, l1->items[i]->column)] = NULL; + } + if (l2->owned) { + free(l2->items); + } + l2->items = NULL; + l2->length = 0; + l2->capacity = 0; + l2->owned = 0; + return c1; +} + +static int64_t progressiveMsa_totalWeights(ProgressiveMsa *m, stList *seqColumns) { + int64_t totalWeights = 0; + for (int64_t i = 0; i < stList_length(seqColumns); i++) { + totalWeights += progressiveMsa_weights(m, stList_get(seqColumns, i))->length; + } + return totalWeights; +} + +/* + * As pairwiseAlignColumns. + */ +static stList *progressiveMsa_alignColumns(ProgressiveMsa *m, stList *seqXColumns, stList *seqYColumns, double matchGamma) { + //Switch seqX and seqY if seqX has more alignment weights associated with it. This is critical to ensure linear scaling, + //else worse case performance is quadratic + if (progressiveMsa_totalWeights(m, seqXColumns) > progressiveMsa_totalWeights(m, seqYColumns)) { + stList *l = seqYColumns; + seqYColumns = seqXColumns; + seqXColumns = l; + } + + //Index the columns of Y + m->tag++; + for (int64_t i = 0; i < stList_length(seqYColumns); i++) { + int64_t id = progressiveMsa_id(m, stList_get(seqYColumns, i)); + m->yIndex[id] = i; + m->yTag[id] = m->tag; + } + + //Best scoring pairs + stSortedSet *bestScoringAlignments = stSortedSet_construct3(columnPair_cmpByYIndex, (void(*)(void *)) columnPair_destruct); + //Add in buffering first and last pairs + ColumnPair *minPair = columnPair_construct(-1, -1, 0, NULL, NULL); + stSortedSet_insert(bestScoringAlignments, minPair); + stSortedSet_insert(bestScoringAlignments, columnPair_construct(stList_length(seqXColumns), stList_length(seqYColumns), INT64_MAX, minPair, NULL)); + + //For each column in X. + stList *l = stList_construct(); + for (int64_t i = 0; i < stList_length(seqXColumns); i++) { + WeightList *aWsX = progressiveMsa_weights(m, stList_get(seqXColumns, i)); + //We first get all the valid new column pairs. + for (int64_t k = 0; k < aWsX->length; k++) { + AlignmentWeight *aWX = aWsX->items[k]; + //Add pair if exceeds the gap gamma. + if (aWX->avgWeight >= matchGamma && aWX->avgWeight > 0.0) { //Must be greater than zero else screws up dynamic programming assumptions + //The column weight may point to a column not in the Y column sequence, if so ignore. + int64_t id = progressiveMsa_id(m, aWX->column); + if (m->yTag[id] == m->tag) { + ColumnPair cP; + cP.yIndex = m->yIndex[id]; + //Search for highest scoring point up to but less than that index. + ColumnPair *cPP = stSortedSet_searchLessThan(bestScoringAlignments, &cP); + assert(cPP != NULL); + assert(i - cPP->xIndex > 0); + assert(cP.yIndex - cPP->yIndex > 0); + assert(cPP->score + aWX->avgWeight * aWX->numberOfWeights > cPP->score); + stList_append(l, columnPair_construct(i, cP.yIndex, /*new score */ cPP->score + aWX->avgWeight * aWX->numberOfWeights, cPP, aWX)); //Make first to increase ref-count of previous position. + } + } + } + //We now work through the new column pairs, from right-to-left along Y. + stList_sort(l, columnPair_cmpByYIndex); + while(stList_length(l) > 0) { + ColumnPair *cP = stList_pop(l); + //Find point that is equal or to the right of cP->yIndex + ColumnPair *cPN = stSortedSet_searchGreaterThanOrEqual(bestScoringAlignments, cP); + assert(cPN != NULL); + if (cP->score >= cPN->score || cPN->yIndex > cP->yIndex) { + //Remove points that overlap or are to the right that score more poorly and clean them up. + while (cP->score >= cPN->score) { + ColumnPair *cPNN = stSortedSet_searchGreaterThan(bestScoringAlignments, cPN); + assert(cPNN != NULL); + stSortedSet_remove(bestScoringAlignments, cPN); + columnPair_destruct(cPN); + cPN = cPNN; + } + //Insert new point. + assert(stSortedSet_search(bestScoringAlignments, cP) == NULL); + stSortedSet_insert(bestScoringAlignments, cP); + } + else { //The new cP is redundant. + columnPair_destruct(cP); + } + } + } + stList_destruct(l); + + //Link the right-most Y pair to the next-right most. + ColumnPair *maxPair = stSortedSet_getLast(bestScoringAlignments); + assert(maxPair != NULL); + stSortedSet_remove(bestScoringAlignments, maxPair); + assert(stSortedSet_getLast(bestScoringAlignments) != NULL); + maxPair->pPair = stSortedSet_getLast(bestScoringAlignments); + maxPair->pPair->refCount++; + + //Now traceback from highest scoring point to generate the alignment + ColumnPair *cP = maxPair; + stList *alignment = stList_construct(); + int64_t merges = 0; + while (1) { + assert(cP->pPair != NULL); + //Add any unaligned Y columns + assert(cP->yIndex > cP->pPair->yIndex); + while (--cP->yIndex > cP->pPair->yIndex) { + stList_append(alignment, stList_get(seqYColumns, cP->yIndex)); + } + //Add any unaligned X columns + assert(cP->xIndex > cP->pPair->xIndex); + while (--cP->xIndex > cP->pPair->xIndex) { + stList_append(alignment, stList_get(seqXColumns, cP->xIndex)); + } + //Now move to previous pair + cP = cP->pPair; + //If this is the final pair we're done + if (cP == minPair) { + break; + } + //Merge two columns. + stList_append(alignment, progressiveMsa_merge(m, cP->aW)); + merges++; + } + assert(stList_length(alignment) + merges == stList_length(seqXColumns) + stList_length(seqYColumns)); + //Make the list of columns left-to-right + stList_reverse(alignment); + + //Cleanup + assert(maxPair->refCount == 1); + columnPair_destruct(maxPair); + stSortedSet_destruct(bestScoringAlignments); + assert(minPair->refCount == 1); + columnPair_destruct(minPair); + stList_destruct(seqXColumns); + stList_destruct(seqYColumns); + + return alignment; +} + stSet *getMultipleSequenceAlignmentProgressive(stList *seqFrags, stList *multipleAlignedPairs, double matchGamma, stList *seqPairSimilarityScores) { - //Get the data-structures needed for the pairwise alignments - stSet *columns = makeColumns(seqFrags); - stHash *alignmentWeightAdjLists = makeAlignmentWeightAdjacencyLists(columns, multipleAlignedPairs); - stSortedSet *alignmentWeightsOrderedByWeight = makeOrderedSetOfAlignmentWeights(alignmentWeightAdjLists); + ProgressiveMsa *m = progressiveMsa_construct(seqFrags, multipleAlignedPairs); //sort list of pairwise distances seqPairSimilarityScores = stList_copy(seqPairSimilarityScores, NULL); stList_sort(seqPairSimilarityScores, (int(*)(const void *, const void *)) stIntTuple_cmpFn); //get list of column-sequences - stList *columnSequences = makeColumnSequences(seqFrags, columns); + stList *columnSequences = stList_construct(); + for (int64_t seq = 0; seq < stList_length(seqFrags); seq++) { + stList *columnSequence = stList_construct(); + for (int64_t id = m->seqOffsets[seq]; id < m->seqOffsets[seq + 1]; id++) { + stList_append(columnSequence, m->columns[id]); + } + stList_append(columnSequences, columnSequence); + } //do n-1 merges while (stList_length(seqPairSimilarityScores) > 0) { @@ -533,7 +884,7 @@ stSet *getMultipleSequenceAlignmentProgressive(stList *seqFrags, stList *multipl stList *seqXColumns = stList_get(columnSequences, seqX); stList *seqYColumns = stList_get(columnSequences, seqY); if (seqXColumns != seqYColumns) { - stList *seqColumns = pairwiseAlignColumns(seqXColumns, seqYColumns, alignmentWeightAdjLists, columns, alignmentWeightsOrderedByWeight, matchGamma); + stList *seqColumns = progressiveMsa_alignColumns(m, seqXColumns, seqYColumns, matchGamma); for(int64_t i=0; icolumnNumber; id++) { + if (m->isHead[id]) { + stSet_insert(columns, m->columns[id]); + } + } + //Clean up - stSortedSet_destruct(alignmentWeightsOrderedByWeight); - stHash_destruct(alignmentWeightAdjLists); if(stList_length(columnSequences) > 0) { stList_destruct(stList_peek(columnSequences)); //This is because we have repeated copies of the same column-sequence in the list } stList_destruct(columnSequences); stList_destruct(seqPairSimilarityScores); - //Return final set of columns. + progressiveMsa_destruct(m); return columns; } @@ -559,44 +917,61 @@ stSet *getMultipleSequenceAlignmentProgressive(stList *seqFrags, stList *multipl * Methods to extract consistent pairs. */ -static Column *getColumn2(stHash *columns, int64_t seq, int64_t pos) { - Column c; - c.seqName = seq; - c.position = pos; - return stHash_search(columns, &c); -} - stList *filterMultipleAlignedPairs(stSet *columns, stList *multipleAlignedPairs) { /* * Processes the list of multipleAlignedPairs and places those that align pairs within the same column in a list which is * returned. Pairs that do not make the list are cleaned up, as is the input list. + * + * Positions are numbered as in the progressive merge, from the columns themselves, so the column each is in is an array lookup. */ - //Build hash of positions to columns + int64_t seqNumber = 0; stSetIterator *it = stSet_getIterator(columns); Column *c; - stHash *positionsToColumns = stHash_construct3((uint64_t(*)(const void *)) column_hashFn, - (int(*)(const void *, const void *)) column_equalsFn, NULL, NULL); while ((c = stSet_getNext(it)) != NULL) { - Column *c2 = c; - do { - stHash_insert(positionsToColumns, c2, c); - c2 = c2->nColumn; - } while (c2 != NULL); + for (Column *c2 = c; c2 != NULL; c2 = c2->nColumn) { + seqNumber = c2->seqName + 1 > seqNumber ? c2->seqName + 1 : seqNumber; + } + } + stSet_destructIterator(it); + int64_t *seqOffsets = st_calloc(seqNumber + 1, sizeof(int64_t)); // first the length of each sequence, then the offsets + it = stSet_getIterator(columns); + while ((c = stSet_getNext(it)) != NULL) { + for (Column *c2 = c; c2 != NULL; c2 = c2->nColumn) { + seqOffsets[c2->seqName] = c2->position + 1 > seqOffsets[c2->seqName] ? c2->position + 1 : seqOffsets[c2->seqName]; + } + } + stSet_destructIterator(it); + int64_t total = 0; + for (int64_t seq = 0; seq <= seqNumber; seq++) { + int64_t length = seqOffsets[seq]; + seqOffsets[seq] = total; + total += length; + } + Column **positionsToColumns = st_calloc(total + 1, sizeof(Column *)); + it = stSet_getIterator(columns); + while ((c = stSet_getNext(it)) != NULL) { + for (Column *c2 = c; c2 != NULL; c2 = c2->nColumn) { + positionsToColumns[seqOffsets[c2->seqName] + c2->position] = c; + } } stSet_destructIterator(it); //Now walk through pairs stList *filteredMultipleAlignedPairs = stList_construct3(0, (void(*)(void *)) stIntTuple_destruct); while (stList_length(multipleAlignedPairs) > 0) { stIntTuple *mAP = stList_pop(multipleAlignedPairs); - if (getColumn2(positionsToColumns, stIntTuple_get(mAP, 1), stIntTuple_get(mAP, 2)) == getColumn2(positionsToColumns, - stIntTuple_get(mAP, 3), stIntTuple_get(mAP, 4))) { + int64_t seq1 = stIntTuple_get(mAP, 1), pos1 = stIntTuple_get(mAP, 2); + int64_t seq2 = stIntTuple_get(mAP, 3), pos2 = stIntTuple_get(mAP, 4); + Column *c1 = seq1 < seqNumber && pos1 < seqOffsets[seq1 + 1] - seqOffsets[seq1] ? positionsToColumns[seqOffsets[seq1] + pos1] : NULL; + Column *c2 = seq2 < seqNumber && pos2 < seqOffsets[seq2 + 1] - seqOffsets[seq2] ? positionsToColumns[seqOffsets[seq2] + pos2] : NULL; + if (c1 == c2) { stList_append(filteredMultipleAlignedPairs, mAP); } else { stIntTuple_destruct(mAP); } } //Cleanup - stHash_destruct(positionsToColumns); + free(positionsToColumns); + free(seqOffsets); stList_destruct(multipleAlignedPairs); return filteredMultipleAlignedPairs; } diff --git a/impl/pairwiseAligner.c b/impl/pairwiseAligner.c index db8f509..e6093a3 100644 --- a/impl/pairwiseAligner.c +++ b/impl/pairwiseAligner.c @@ -755,6 +755,793 @@ static void diagonalCalculationExpectations(StateMachine *sM, int64_t xay, DpMat dpMatrix_getDiagonal(forwardDpMatrix, xay - 2), sX, sY, cell_calculateExpectation, extraArgs2); } +/////////////////////////////////// +/////////////////////////////////// +//Fast banded forward-backward +// +//The log space recursion below goes through three levels of function pointer per cell (the +//diagonal, the cell and the transition) and an approximate log-add per transition. For the +//posterior probabilities -- which is all bar asks for -- the same banded, chunked recursion is +//computed here in probability space instead: each diagonal is rescaled so its largest value is +//one and the scale is kept as a log, the state machine is flattened into arrays once per call, +//and each state of a diagonal is a plain loop over its cells. The band, the traceback points +//and which diagonals get posteriors in which chunk are exactly those of the log space version; +//only the arithmetic differs, and this is the more accurate of the two, since logAdd is a +//piecewise cubic fit that also drops any term under e^-7.5 of the larger. +/////////////////////////////////// +/////////////////////////////////// + +#if defined(__SSE2__) +#include +#endif +#if defined(__AVX2__) +#include +#endif + +#define PAIR_HMM_MAX_STATES 5 + +enum { PAIR_HMM_MATCH = 0, PAIR_HMM_GAP_X = 1, PAIR_HMM_GAP_Y = 2 }; + +typedef struct _pairHmmTerms { // the nonzero transitions into (or out of) a state + int64_t number; + int64_t state[PAIR_HMM_MAX_STATES]; + double t[PAIR_HMM_MAX_STATES]; +} PairHmmTerms; + +typedef struct _pairHmm { + int64_t stateNumber; + int64_t type[PAIR_HMM_MAX_STATES]; // what a state emits: a match, a gap in y (x only) or a gap in x (y only) + double t[PAIR_HMM_MAX_STATES][PAIR_HMM_MAX_STATES]; // transition probabilities, [from][to] + double eMatch[SYMBOL_NUMBER * SYMBOL_NUMBER], eGapX[SYMBOL_NUMBER], eGapY[SYMBOL_NUMBER]; + double start[PAIR_HMM_MAX_STATES], raggedStart[PAIR_HMM_MAX_STATES]; + double end[PAIR_HMM_MAX_STATES], raggedEnd[PAIR_HMM_MAX_STATES]; + PairHmmTerms into[PAIR_HMM_MAX_STATES], outOf[PAIR_HMM_MAX_STATES]; // t, sparse, by destination and by source + bool ok; // false if the state machine does not fit this form, in which case the log space code is used +} PairHmm; + +typedef struct { + PairHmm *hmm; + double *lower, *middle, *upper; + int64_t x, y; +} PairHmmExtraction; + +static void pairHmm_setOnce(double *slot, double value, PairHmm *hmm) { + if (*slot < 0.0) { + *slot = value; + } else if (*slot != value) { + hmm->ok = 0; + } +} + +/* + * A doTransition that records the transition instead of computing with it. Which neighbour + * the transition reads from says what the destination state emits. + */ +static void pairHmm_recordTransition(double *fromCells, double *toCells, int64_t from, int64_t to, double eP, double tP, + void *extraArgs) { + PairHmmExtraction *e = extraArgs; + PairHmm *hmm = e->hmm; + if (from < 0 || from >= hmm->stateNumber || to < 0 || to >= hmm->stateNumber) { + hmm->ok = 0; + return; + } + int64_t type = fromCells == e->lower ? PAIR_HMM_GAP_X : (fromCells == e->middle ? PAIR_HMM_MATCH : PAIR_HMM_GAP_Y); + if (hmm->type[to] == -1) { + hmm->type[to] = type; + } else if (hmm->type[to] != type) { // the recursion needs each state entered from one neighbour only + hmm->ok = 0; + } + pairHmm_setOnce(&hmm->t[from][to], exp(tP), hmm); // must not depend on the symbols + double *emission = type == PAIR_HMM_MATCH ? &hmm->eMatch[e->x * SYMBOL_NUMBER + e->y] : + (type == PAIR_HMM_GAP_X ? &hmm->eGapX[e->x] : &hmm->eGapY[e->y]); + pairHmm_setOnce(emission, exp(eP), hmm); // must not depend on the state transitioned from +} + +/* + * Flattens the state machine by running its cell calculation once for every pair of symbols + * and recording the transitions, rather than by reading its fields, so it works for any state + * machine whose states each emit in one way and whose transitions do not depend on the symbols. + */ +static void pairHmm_construct(PairHmm *hmm, StateMachine *sM) { + hmm->ok = sM->stateNumber > 0 && sM->stateNumber <= PAIR_HMM_MAX_STATES; + if (!hmm->ok) { + return; + } + hmm->stateNumber = sM->stateNumber; + for (int64_t s = 0; s < PAIR_HMM_MAX_STATES; s++) { + hmm->type[s] = -1; + for (int64_t s2 = 0; s2 < PAIR_HMM_MAX_STATES; s2++) { + hmm->t[s][s2] = -1.0; + } + } + for (int64_t i = 0; i < SYMBOL_NUMBER * SYMBOL_NUMBER; i++) { + hmm->eMatch[i] = -1.0; + } + for (int64_t i = 0; i < SYMBOL_NUMBER; i++) { + hmm->eGapX[i] = -1.0; + hmm->eGapY[i] = -1.0; + } + double current[PAIR_HMM_MAX_STATES], lower[PAIR_HMM_MAX_STATES], middle[PAIR_HMM_MAX_STATES], upper[PAIR_HMM_MAX_STATES]; + PairHmmExtraction e = { hmm, lower, middle, upper, 0, 0 }; + for (e.x = 0; e.x < SYMBOL_NUMBER; e.x++) { + for (e.y = 0; e.y < SYMBOL_NUMBER; e.y++) { + sM->cellCalculate(sM, current, lower, middle, upper, (Symbol) e.x, (Symbol) e.y, pairHmm_recordTransition, &e); + } + } + for (int64_t s = 0; s < hmm->stateNumber; s++) { + if (hmm->type[s] == -1) { + hmm->ok = 0; + } + for (int64_t s2 = 0; s2 < hmm->stateNumber; s2++) { + if (hmm->t[s][s2] < 0.0) { + hmm->t[s][s2] = 0.0; // no such transition + } + } + hmm->start[s] = exp(sM->startStateProb(sM, s)); + hmm->raggedStart[s] = exp(sM->raggedStartStateProb(sM, s)); + hmm->end[s] = exp(sM->endStateProb(sM, s)); + hmm->raggedEnd[s] = exp(sM->raggedEndStateProb(sM, s)); + } + for (int64_t i = 0; i < SYMBOL_NUMBER * SYMBOL_NUMBER; i++) { + hmm->ok = hmm->ok && hmm->eMatch[i] >= 0.0; + } + for (int64_t i = 0; i < SYMBOL_NUMBER; i++) { + hmm->ok = hmm->ok && hmm->eGapX[i] >= 0.0 && hmm->eGapY[i] >= 0.0; + } + for (int64_t s = 0; s < hmm->stateNumber; s++) { + hmm->into[s].number = 0; + hmm->outOf[s].number = 0; + } + for (int64_t from = 0; from < hmm->stateNumber; from++) { // in increasing order of state, both ways + for (int64_t to = 0; to < hmm->stateNumber; to++) { + double t = hmm->t[from][to]; + if (t > 0.0) { + PairHmmTerms *i = &hmm->into[to], *o = &hmm->outOf[from]; + i->state[i->number] = from; + i->t[i->number++] = t; + o->state[o->number] = to; + o->t[o->number++] = t; + } + } + }} + +/* + * A diagonal of the scaled dp matrix. Cells are stored state-major, so the cell for state s + * at xmy = xmyL + 2i is cells[s * width + i], and the true value is that times exp(logScale). + */ +enum { FAST_STORED, FAST_RING, FAST_BLOCK }; // where a forward diagonal's cells are, see FastForward + +typedef struct _fastDiagonal { + int64_t xmyL, width; + int64_t where, offset; // which buffer holds the cells, and where in it (for the ring, which of its three) + double logScale; +} FastDiagonal; + +/* + * The dp between one traceback and the next. A chunk computes forward values for diagonals + * (start, end], walks the backward values from end down to tracedBackTo, takes posteriors up + * to tracedBackFrom, and hands the forward values from tracedBackFrom on to the next chunk. + * + * A large chunk does not keep all of its forward diagonals: only the ones it hands on and a + * pair before every blockSize'th, from which the walk back recomputes each block of them when + * it gets there. That is a second forward pass over the chunk, but it holds a few hundred + * diagonals rather than thousands, and what the walk reads is still in cache. Keeping them + * all, an unanchored 3000x3000 matrix holds 360MB per thread and runs at the speed of memory. + */ +typedef struct _fastChunk { + int64_t start, end, tracedBackTo, tracedBackFrom; + int64_t blockSize; // 0 if every forward diagonal is kept +} FastChunk; + +#define FAST_CHECKPOINT_BYTES (16 * 1024 * 1024) // chunks whose forward values would take more are recomputed in blocks + +/* + * Whether the chunk keeps diagonal xay's forward values, rather than computing them only to feed + * the next two. Blocks start at start + 1 + j * blockSize. + */ +static bool fastChunk_stores(FastChunk *c, int64_t xay) { + if (c->blockSize == 0 || xay >= c->tracedBackFrom || xay <= c->start) { + return 1; + } + for (int64_t blockStart = xay + 1; blockStart <= xay + 2; blockStart++) { // the pair before a block + if (blockStart > c->start + 1 && blockStart < c->tracedBackFrom && (blockStart - c->start - 1) % c->blockSize == 0) { + return 1; + } + } + return 0; +} + +/* + * The forward diagonals the current chunk holds, and the buffers they are held in. + */ +typedef struct _fastForward { + FastDiagonal *diagonals; // diagonals[xay - base], for every diagonal of the chunk however it is held + int64_t base, length, capacity; + double *cells; // FAST_STORED: kept diagonals, contiguous and in order + int64_t cellsLength, cellsCapacity; + double *ring[3]; // FAST_RING: diagonals computed only to feed the next two, round robin + double *block; // FAST_BLOCK: the block the walk back last recomputed + int64_t blockStart, blockCellsLength, blockCellsCapacity; +} FastForward; + +static FastDiagonal *fastForward_get(FastForward *f, int64_t xay) { + return xay >= f->base && xay < f->base + f->length ? &f->diagonals[xay - f->base] : NULL; +} + +static double *fastForward_cells(FastForward *f, FastDiagonal *d) { + return d->where == FAST_STORED ? f->cells + d->offset : (d->where == FAST_RING ? f->ring[d->offset] : f->block + d->offset); +} + +static FastDiagonal *fastForward_add(FastForward *f, Diagonal diagonal, int64_t stateNumber, bool store) { + assert(f->length < f->capacity); // all the buffers are sized up front, from the band + int64_t xay = diagonal_getXay(diagonal); + assert(xay == f->base + f->length); + FastDiagonal *d = &f->diagonals[f->length++]; + d->xmyL = diagonal_getMinXmy(diagonal); + d->width = diagonal_getWidth(diagonal); + d->logScale = 0.0; + if (store) { + d->where = FAST_STORED; + d->offset = f->cellsLength; + f->cellsLength += d->width * stateNumber; + assert(f->cellsLength <= f->cellsCapacity); + } else { + d->where = FAST_RING; + d->offset = xay % 3; + } + return d; +} + +/* + * Plans the chunks. Where the tracebacks fall depends only on the band, so the chunks, and so + * the buffers they need, are known before any dp is done: nothing is grown as the dp goes. + */ +static FastChunk *fastForward_construct(FastForward *f, Band *band, int64_t diagonalNumber, int64_t stateNumber, + PairwiseAlignmentParameters *p, int64_t *chunkNumber) { + int64_t chunkCapacity = 16; + FastChunk *chunks = st_malloc(sizeof(FastChunk) * chunkCapacity); + *chunkNumber = 0; + int64_t maxDiagonals = 1, maxCells = 1, maxBlockCells = 1, maxWidth = 1; + int64_t tracedBackTo = 0, previousEnd = 0; + for (int64_t xay = 0; xay <= diagonalNumber; xay++) { + int64_t width = diagonal_getWidth(band->diagonals[xay]); + maxWidth = width > maxWidth ? width : maxWidth; + bool atEnd = xay == diagonalNumber; + bool tracebackPoint = xay >= tracedBackTo + p->minDiagsBetweenTraceBack && width <= p->diagonalExpansion * 2 + 1; + if (xay == 0 || (!atEnd && !tracebackPoint)) { + continue; + } + FastChunk c = { previousEnd, xay, tracedBackTo, xay - (atEnd ? 0 : p->traceBackDiagonals + 1), 0 }; + int64_t cells = 0; + for (int64_t d = c.tracedBackTo; d <= c.end; d++) { + cells += diagonal_getWidth(band->diagonals[d]); + } + if (cells * stateNumber * (int64_t) sizeof(double) > FAST_CHECKPOINT_BYTES) { + c.blockSize = (int64_t) sqrt((double) (c.tracedBackFrom - c.start)); + c.blockSize = c.blockSize < 8 ? 8 : c.blockSize; + cells = 0; + int64_t blockCells = 0; + for (int64_t d = c.tracedBackTo; d <= c.end; d++) { + int64_t w = diagonal_getWidth(band->diagonals[d]); + cells += fastChunk_stores(&c, d) ? w : 0; + if (d > c.start && d < c.tracedBackFrom) { + blockCells = (d - c.start - 1) % c.blockSize == 0 ? w : blockCells + w; + maxBlockCells = blockCells > maxBlockCells ? blockCells : maxBlockCells; + } + } + } + maxCells = cells > maxCells ? cells : maxCells; + maxDiagonals = c.end - c.tracedBackTo + 1 > maxDiagonals ? c.end - c.tracedBackTo + 1 : maxDiagonals; + if (*chunkNumber == chunkCapacity) { + chunkCapacity *= 2; + chunks = st_realloc(chunks, sizeof(FastChunk) * chunkCapacity); + } + chunks[(*chunkNumber)++] = c; + tracedBackTo = c.tracedBackFrom; + previousEnd = xay; + } + f->capacity = maxDiagonals + 1; + f->diagonals = st_malloc(sizeof(FastDiagonal) * f->capacity); + f->cellsCapacity = maxCells * stateNumber; + f->cells = st_malloc(sizeof(double) * f->cellsCapacity); + for (int64_t i = 0; i < 3; i++) { + f->ring[i] = st_malloc(sizeof(double) * maxWidth * stateNumber); + } + f->blockCellsCapacity = maxBlockCells * stateNumber; + f->block = st_malloc(sizeof(double) * f->blockCellsCapacity); + f->blockStart = -1; + f->blockCellsLength = 0; + f->base = 0; + f->length = 0; + f->cellsLength = 0; + return chunks; +} + +static void fastForward_destruct(FastForward *f) { + free(f->diagonals); + free(f->cells); + for (int64_t i = 0; i < 3; i++) { + free(f->ring[i]); + } + free(f->block); +} + +/* + * Drops the diagonals before xay, moving the rest -- all kept -- to the front of the buffer. + */ +static void fastForward_keepFrom(FastForward *f, int64_t xay) { + assert(xay >= f->base && xay < f->base + f->length); + int64_t first = xay - f->base; + assert(f->diagonals[first].where == FAST_STORED); + int64_t cellStart = f->diagonals[first].offset; + memmove(f->cells, f->cells + cellStart, sizeof(double) * (f->cellsLength - cellStart)); + f->cellsLength -= cellStart; + memmove(f->diagonals, f->diagonals + first, sizeof(FastDiagonal) * (f->length - first)); + f->length -= first; + f->base = xay; + for (int64_t i = 0; i < f->length; i++) { + assert(f->diagonals[i].where == FAST_STORED); + f->diagonals[i].offset -= cellStart; + } + f->blockStart = -1; +} + +/* + * Where cell i of a diagonal finds its neighbour in another diagonal, whose xmy is the cell's + * plus xmyDelta: the neighbour is cell i + shift, and exists for i in [lo, hi). + */ +static void fastDiagonal_neighbours(FastDiagonal *cur, FastDiagonal *other, int64_t xmyDelta, + int64_t *shift, int64_t *lo, int64_t *hi) { + if (other == NULL) { + *shift = 0; + *lo = 0; + *hi = 0; + return; + } + assert((cur->xmyL + xmyDelta - other->xmyL) % 2 == 0); + *shift = (cur->xmyL + xmyDelta - other->xmyL) / 2; + *lo = *shift < 0 ? -*shift : 0; + *hi = other->width - *shift < cur->width ? other->width - *shift : cur->width; + if (*hi < *lo) { + *hi = *lo; + } +} + +/* + * The values of a diagonal are kept relative to a scale, exp(logScale), chosen so they neither + * underflow nor overflow. Rescaling every diagonal cost a pass over it for nothing: probabilities + * shrink by a bounded factor per diagonal, so the scale only needs moving every few dozen. Moves + * it, if the diagonal's largest value has drifted out of range, and returns the log of the factor + * the cells were divided by. + */ +#define FAST_SCALE_LOW 1e-100 +#define FAST_SCALE_HIGH 1e100 + +static double fastDiagonal_rescale(double *restrict cells, int64_t n, double max) { + if (max <= 0.0 || (max >= FAST_SCALE_LOW && max <= FAST_SCALE_HIGH)) { + return 0.0; + } + double r = 1.0 / max; + for (int64_t i = 0; i < n; i++) { + cells[i] *= r; + } + return log(max); +} + +/* + * out[i] = em[i] * sum_k t[k] * in[k][i] for i in [0, n), or without the em factor if em is NULL, + * returning the largest out[i]. One loop per number of terms, so each case vectorises and the + * transitions a state lacks cost nothing. The terms are summed in order, as the separate passes + * this replaced did. + */ +#define FAST_SUM_LOOP(EXPR) \ + _Pragma("omp simd reduction(max:max)") \ + for (int64_t i = 0; i < n; i++) { \ + double v = (EXPR); \ + out[i] = v; \ + max = v > max ? v : max; \ + } + +static double fastSum(double *restrict out, const double *restrict em, const double *const *in, const double *t, + int64_t terms, int64_t n) { + double max = 0.0; + const double *restrict i0 = terms > 0 ? in[0] : NULL, *restrict i1 = terms > 1 ? in[1] : NULL; + const double *restrict i2 = terms > 2 ? in[2] : NULL, *restrict i3 = terms > 3 ? in[3] : NULL; + const double *restrict i4 = terms > 4 ? in[4] : NULL; + double t0 = terms > 0 ? t[0] : 0, t1 = terms > 1 ? t[1] : 0, t2 = terms > 2 ? t[2] : 0; + double t3 = terms > 3 ? t[3] : 0, t4 = terms > 4 ? t[4] : 0; + if (em != NULL) { + switch (terms) { + case 0: FAST_SUM_LOOP(0.0); break; + case 1: FAST_SUM_LOOP(em[i] * (t0 * i0[i])); break; + case 2: FAST_SUM_LOOP(em[i] * (t0 * i0[i] + t1 * i1[i])); break; + case 3: FAST_SUM_LOOP(em[i] * (t0 * i0[i] + t1 * i1[i] + t2 * i2[i])); break; + case 4: FAST_SUM_LOOP(em[i] * (t0 * i0[i] + t1 * i1[i] + t2 * i2[i] + t3 * i3[i])); break; + default: FAST_SUM_LOOP(em[i] * (t0 * i0[i] + t1 * i1[i] + t2 * i2[i] + t3 * i3[i] + t4 * i4[i])); + } + } else { + switch (terms) { + case 0: FAST_SUM_LOOP(0.0); break; + case 1: FAST_SUM_LOOP(t0 * i0[i]); break; + case 2: FAST_SUM_LOOP(t0 * i0[i] + t1 * i1[i]); break; + case 3: FAST_SUM_LOOP(t0 * i0[i] + t1 * i1[i] + t2 * i2[i]); break; + case 4: FAST_SUM_LOOP(t0 * i0[i] + t1 * i1[i] + t2 * i2[i] + t3 * i3[i]); break; + default: FAST_SUM_LOOP(t0 * i0[i] + t1 * i1[i] + t2 * i2[i] + t3 * i3[i] + t4 * i4[i]); + } + } + return max; +} + +static void fastZero(double *restrict out, int64_t n) { + for (int64_t i = 0; i < n; i++) { + out[i] = 0.0; + } +} + +/* + * The two sequences as bytes, with Y reversed, so that walking along a diagonal -- x up, y + * down -- reads both forwards, and the emission lookups below vectorise. + */ +typedef struct _fastSymbols { + uint8_t *x, *yReversed; + int64_t lX, lY; +} FastSymbols; + +static FastSymbols fastSymbols_construct(const SymbolString sX, const SymbolString sY) { + FastSymbols f; + f.lX = sX.length; + f.lY = sY.length; + f.x = st_malloc(sX.length + 1); + f.yReversed = st_malloc(sY.length + 1); + for (int64_t i = 0; i < sX.length; i++) { + f.x[i] = (uint8_t) sX.sequence[i]; + } + for (int64_t i = 0; i < sY.length; i++) { + f.yReversed[i] = (uint8_t) sY.sequence[sY.length - 1 - i]; + } + return f; +} + +static void fastSymbols_destruct(FastSymbols f) { + free(f.x); + free(f.yReversed); +} + +/* + * The emission of cells [lo, hi) of a diagonal into a state of the given type, into em[lo, hi). + * x0/y0 are the coordinates of cell 0, and cell i is at (x0 + i, y0 - i). For the forward + * recursion a cell emits its own symbols (the ones before it, at x-1 and y-1); for the backward + * it is the successor's (at x and y) -- pass x0 and y0 one larger in that case. + */ +static void fastDiagonal_emissions(PairHmm *hmm, int64_t type, FastSymbols *symbols, + int64_t x0, int64_t y0, int64_t lo, int64_t hi, double scale, double *restrict em) { + double *restrict e = em + lo; + int64_t n = hi - lo, i = 0; + // Only the sequences the type emits are read: at the edge of the matrix the other one's index + // is off its end. cX[i] is X[x0 - 1 + lo + i] and cY[i] is Y[y0 - 1 - lo - i]. + const uint8_t *restrict cX = type != PAIR_HMM_GAP_Y ? symbols->x + (x0 - 1 + lo) : NULL; + const uint8_t *restrict cY = type != PAIR_HMM_GAP_X ? symbols->yReversed + (symbols->lY - y0 + lo) : NULL; +#if defined(__AVX2__) + // Table lookups the compiler will not vectorise on its own + __m256d vScale = _mm256_set1_pd(scale); + __m128i vSymbols = _mm_set1_epi32(SYMBOL_NUMBER); + for (; i + 4 <= n; i += 4) { + int32_t bytes; + __m256d v; + if (type == PAIR_HMM_GAP_X) { + memcpy(&bytes, cX + i, sizeof(int32_t)); + v = _mm256_i32gather_pd(hmm->eGapX, _mm_cvtepu8_epi32(_mm_cvtsi32_si128(bytes)), 8); + } else if (type == PAIR_HMM_GAP_Y) { + memcpy(&bytes, cY + i, sizeof(int32_t)); + v = _mm256_i32gather_pd(hmm->eGapY, _mm_cvtepu8_epi32(_mm_cvtsi32_si128(bytes)), 8); + } else { + memcpy(&bytes, cX + i, sizeof(int32_t)); + __m128i vx = _mm_cvtepu8_epi32(_mm_cvtsi32_si128(bytes)); + memcpy(&bytes, cY + i, sizeof(int32_t)); + __m128i vy = _mm_cvtepu8_epi32(_mm_cvtsi32_si128(bytes)); + v = _mm256_mul_pd(_mm256_i32gather_pd(hmm->eMatch, _mm_add_epi32(_mm_mullo_epi32(vx, vSymbols), vy), 8), vScale); + } + _mm256_storeu_pd(e + i, v); + } +#endif + if (type == PAIR_HMM_GAP_X) { + for (; i < n; i++) { + e[i] = hmm->eGapX[cX[i]]; + } + } else if (type == PAIR_HMM_GAP_Y) { + for (; i < n; i++) { + e[i] = hmm->eGapY[cY[i]]; + } + } else { + for (; i < n; i++) { + e[i] = hmm->eMatch[cX[i] * SYMBOL_NUMBER + cY[i]] * scale; + } + } +} + +/* + * Forward values for diagonal cur, from the diagonals one (m1) and two (m2) before it. + */ +static void fastForwardDiagonal(PairHmm *hmm, int64_t xay, FastDiagonal *cur, double *restrict curCells, + FastDiagonal *m1, const double *m1Cells, FastDiagonal *m2, const double *m2Cells, + FastSymbols *symbols, double *restrict em) { + int64_t w = cur->width; + int64_t x0 = (xay + cur->xmyL) / 2, y0 = (xay - cur->xmyL) / 2; + // Matches read from two diagonals back, whose scale can differ from the previous diagonal's, + // which this one's is taken relative to. Fold the difference into the match emissions. + double matchScale = m2 != NULL ? (m2->logScale == m1->logScale ? 1.0 : exp(m2->logScale - m1->logScale)) : 0.0; + double max = 0.0; + for (int64_t type = 0; type < 3; type++) { + FastDiagonal *src = type == PAIR_HMM_MATCH ? m2 : m1; + const double *srcCells = type == PAIR_HMM_MATCH ? m2Cells : m1Cells; + int64_t shift, lo, hi; + fastDiagonal_neighbours(cur, src, type == PAIR_HMM_GAP_X ? -1 : (type == PAIR_HMM_GAP_Y ? 1 : 0), &shift, &lo, &hi); + if (hi > lo) { + fastDiagonal_emissions(hmm, type, symbols, x0, y0, lo, hi, matchScale, em); + } + for (int64_t s = 0; s < hmm->stateNumber; s++) { + if (hmm->type[s] != type) { + continue; + } + double *restrict out = curCells + s * w; + fastZero(out, lo); + fastZero(out + hi, w - hi); + if (hi > lo) { + const double *in[PAIR_HMM_MAX_STATES]; + for (int64_t k = 0; k < hmm->into[s].number; k++) { + in[k] = srcCells + hmm->into[s].state[k] * src->width + shift + lo; + } + double m = fastSum(out + lo, em + lo, in, hmm->into[s].t, hmm->into[s].number, hi - lo); + max = m > max ? m : max; + } + } + } + cur->logScale = m1->logScale + fastDiagonal_rescale(curCells, w * hmm->stateNumber, max); +} + +/* + * Backward values for diagonal cur, from the diagonals one (p1) and two (p2) after it. w is + * scratch for stateNumber * width values. + */ +static void fastBackwardDiagonal(PairHmm *hmm, int64_t xay, FastDiagonal *cur, double *restrict curCells, + FastDiagonal *p1, const double *p1Cells, FastDiagonal *p2, const double *p2Cells, + FastSymbols *symbols, double *restrict w, double *restrict em) { + int64_t width = cur->width; + int64_t x0 = (xay + cur->xmyL) / 2, y0 = (xay - cur->xmyL) / 2; + double matchScale = p2 != NULL ? (p2->logScale == p1->logScale ? 1.0 : exp(p2->logScale - p1->logScale)) : 0.0; + // w[s] = what entering state s from this cell is worth: the successor's emission and backward value + for (int64_t type = 0; type < 3; type++) { + FastDiagonal *succ = type == PAIR_HMM_MATCH ? p2 : p1; + const double *succCells = type == PAIR_HMM_MATCH ? p2Cells : p1Cells; + int64_t shift, lo, hi; + fastDiagonal_neighbours(cur, succ, type == PAIR_HMM_GAP_X ? 1 : (type == PAIR_HMM_GAP_Y ? -1 : 0), &shift, &lo, &hi); + if (hi > lo) { + fastDiagonal_emissions(hmm, type, symbols, x0 + 1, y0 + 1, lo, hi, matchScale, em); + } + for (int64_t s = 0; s < hmm->stateNumber; s++) { + if (hmm->type[s] != type) { + continue; + } + double *restrict ws = w + s * width; + fastZero(ws, lo); + fastZero(ws + hi, width - hi); + if (hi > lo) { + const double *restrict in = succCells + s * succ->width + shift + lo; + double *restrict out = ws + lo; + const double *restrict e = em + lo; + for (int64_t i = 0; i < hi - lo; i++) { + out[i] = e[i] * in[i]; + } + } + } + } + double max = 0.0; + for (int64_t f = 0; f < hmm->stateNumber; f++) { + const double *in[PAIR_HMM_MAX_STATES]; + for (int64_t k = 0; k < hmm->outOf[f].number; k++) { + in[k] = w + hmm->outOf[f].state[k] * width; + } + double m = fastSum(curCells + f * width, NULL, in, hmm->outOf[f].t, hmm->outOf[f].number, width); + max = m > max ? m : max; + } + cur->logScale = p1->logScale + fastDiagonal_rescale(curCells, width * hmm->stateNumber, max); +} + +/* + * Posteriors of the states the caller wants lists for, over one diagonal. + */ +static void fastPosteriors(StateMachine *sM, int64_t xay, FastDiagonal *fd, const double *fCells, + FastDiagonal *bd, const double *bCells, double logTotal, PairwiseAlignmentParameters *p, + stList *alignedPairs, stList *gapXPairs, stList *gapYPairs) { + assert(fd->xmyL == bd->xmyL && fd->width == bd->width); + int64_t w = fd->width; + int64_t x0 = (xay + fd->xmyL) / 2, y0 = (xay - fd->xmyL) / 2; + double logNorm = fd->logScale + bd->logScale - logTotal; + bool direct = logNorm < 700.0; // else exp(logNorm) would overflow, and it is done cell by cell in log space + double norm = direct ? exp(logNorm) : 0.0; + // Most cells fall short of the threshold, so test them against it in their own scale -- one + // multiply and one compare -- a hair low, so addPosteriorProb makes the exact decision. + double fbThreshold = direct ? p->threshold * exp(-logNorm) * (1.0 - 1e-9) : 0.0; + bool everyCell = p->threshold <= 0.0; // then the log space path reports every cell, zeros included + int64_t states[3] = { sM->matchState, sM->gapXState, sM->gapYState }; + stList *lists[3] = { alignedPairs, gapXPairs, gapYPairs }; + for (int64_t k = 0; k < 3; k++) { + if (lists[k] == NULL) { + continue; + } + // Matches and gaps in y need x > 0, i.e. i >= 1 - x0; matches and gaps in x need y > 0, i.e. i < y0 + int64_t lo = k != 2 && 1 - x0 > 0 ? 1 - x0 : 0; + int64_t hi = k != 1 && y0 < w ? y0 : w; + const double *f = fCells + states[k] * w, *b = bCells + states[k] * w; + for (int64_t i = lo; i < hi; i++) { + double fb = f[i] * b[i]; + if (everyCell || (fb >= fbThreshold && fb > 0.0)) { + addPosteriorProb(x0 + i, y0 - i, direct ? fb * norm : (fb > 0.0 ? exp(log(fb) + logNorm) : 0.0), lists[k], p); + } + } + } +} + +/* + * The fast equivalent of getPosteriorProbsWithBanding with diagonalCalculationPosteriorMatchProbs + * (gapXPairs and gapYPairs NULL) or diagonalCalculationPosteriorProbs. Returns false, having done + * nothing, if the state machine is not one it can use. + */ +/* + * The machine's flattened form, made the first time it is asked for and kept on the machine: + * flattening runs the cell calculation for every pair of symbols and takes hundreds of exps, which + * cost more than the dp itself for the thousands of few-base alignments a bar run makes. Threads + * share machines, so a thread that loses the race to store its copy frees it and uses the winner's. + */ +static PairHmm *pairHmm_get(StateMachine *sM) { + PairHmm *hmm = __atomic_load_n((PairHmm **) &sM->flat, __ATOMIC_ACQUIRE); + if (hmm == NULL) { + PairHmm *made = st_malloc(sizeof(PairHmm)); + pairHmm_construct(made, sM); + PairHmm *expected = NULL; + if (__atomic_compare_exchange_n((PairHmm **) &sM->flat, &expected, made, 0, __ATOMIC_ACQ_REL, __ATOMIC_ACQUIRE)) { + hmm = made; + } else { + free(made); + hmm = expected; + } + } + return hmm; +} + +static bool getPosteriorProbsWithBandingFast(StateMachine *sM, stList *anchorPairs, const SymbolString sX, const SymbolString sY, + PairwiseAlignmentParameters *p, bool alignmentHasRaggedLeftEnd, bool alignmentHasRaggedRightEnd, + stList *alignedPairs, stList *gapXPairs, stList *gapYPairs) { + PairHmm hmm = *pairHmm_get(sM); + if (!hmm.ok) { + return 0; + } + int64_t diagonalNumber = sX.length + sY.length; + if (diagonalNumber == 0) { + return 1; + } + int64_t ns = hmm.stateNumber; + +#if defined(__SSE2__) + // Cells far off the alignment decay towards zero; flush them rather than crawl through denormals + unsigned int oldCsr = _mm_getcsr(); + _mm_setcsr(oldCsr | 0x8040); // flush-to-zero and denormals-are-zero +#endif + + Band *band = p->dynamicAnchorExpansion ? band_constructDynamic(anchorPairs, sX.length, sY.length) : + band_construct(anchorPairs, sX.length, sY.length, p->diagonalExpansion); + FastForward fwd; + int64_t chunkNumber; + FastChunk *chunks = fastForward_construct(&fwd, band, diagonalNumber, ns, p, &chunkNumber); + FastSymbols symbols = fastSymbols_construct(sX, sY); + + // scratch, sized to the widest diagonal + int64_t maxWidth = 1; + for (int64_t d = 0; d <= diagonalNumber; d++) { + int64_t w = diagonal_getWidth(band->diagonals[d]); + maxWidth = w > maxWidth ? w : maxWidth; + } + double *em = st_malloc(sizeof(double) * maxWidth); + double *scratch = st_malloc(sizeof(double) * maxWidth * ns); + FastDiagonal back[3]; // the backward diagonals xay, xay + 1 and xay + 2, round robin + double *backCells[3]; + for (int64_t i = 0; i < 3; i++) { + backCells[i] = st_malloc(sizeof(double) * maxWidth * ns); + } + + FastDiagonal *d0 = fastForward_add(&fwd, band->diagonals[0], ns, 1); + for (int64_t s = 0; s < ns; s++) { + double v = alignmentHasRaggedLeftEnd ? hmm.raggedStart[s] : hmm.start[s]; + for (int64_t i = 0; i < d0->width; i++) { + fwd.cells[d0->offset + s * d0->width + i] = v; + } + } + + for (int64_t c = 0; c < chunkNumber; c++) { + FastChunk *chunk = &chunks[c]; + for (int64_t xay = chunk->start + 1; xay <= chunk->end; xay++) { + FastDiagonal *cur = fastForward_add(&fwd, band->diagonals[xay], ns, fastChunk_stores(chunk, xay)); + FastDiagonal *m1 = fastForward_get(&fwd, xay - 1), *m2 = fastForward_get(&fwd, xay - 2); + fastForwardDiagonal(&hmm, xay, cur, fastForward_cells(&fwd, cur), m1, fastForward_cells(&fwd, m1), + m2, m2 != NULL ? fastForward_cells(&fwd, m2) : NULL, &symbols, em); + } + + // Treat every cell of the chunk's last diagonal as an end point, and walk back to where the + // last traceback stopped. Posteriors are taken only up to traceBackDiagonals before the + // end, where that pretence has worn off; the rest are redone by the next traceback. + bool atEnd = chunk->end == diagonalNumber; + const double *endProbs = atEnd && alignmentHasRaggedRightEnd ? hmm.raggedEnd : hmm.end; + FastDiagonal *last = fastForward_get(&fwd, chunk->end); + const double *lastCells = fastForward_cells(&fwd, last); + double total = 0.0; + for (int64_t s = 0; s < ns; s++) { + for (int64_t i = 0; i < last->width; i++) { + total += lastCells[s * last->width + i] * endProbs[s]; + } + } + double logTotal = last->logScale + log(total); + bool havePosteriors = total > 0.0; + + for (int64_t xay2 = chunk->end; xay2 > chunk->tracedBackTo; xay2--) { + FastDiagonal *b = &back[xay2 % 3]; + double *bCells = backCells[xay2 % 3]; + FastDiagonal *f = fastForward_get(&fwd, xay2); + b->xmyL = f->xmyL; + b->width = f->width; + if (xay2 == chunk->end) { + for (int64_t s = 0; s < ns; s++) { + for (int64_t i = 0; i < b->width; i++) { + bCells[s * b->width + i] = endProbs[s]; + } + } + b->logScale = 0.0; + } else { + FastDiagonal *p1 = &back[(xay2 + 1) % 3], *p2 = xay2 + 2 <= chunk->end ? &back[(xay2 + 2) % 3] : NULL; + fastBackwardDiagonal(&hmm, xay2, b, bCells, p1, backCells[(xay2 + 1) % 3], + p2, p2 != NULL ? backCells[(xay2 + 2) % 3] : NULL, &symbols, scratch, em); + } + if (xay2 > chunk->tracedBackFrom || !havePosteriors) { + continue; + } + if (f->where != FAST_STORED && !(f->where == FAST_BLOCK && xay2 >= fwd.blockStart)) { + // Recompute the block this diagonal is in, from the kept pair before it. The same + // code on the same inputs, so the same values, scales and all, as the first time. + int64_t blockStart = chunk->start + 1 + (xay2 - chunk->start - 1) / chunk->blockSize * chunk->blockSize; + fwd.blockStart = blockStart; + fwd.blockCellsLength = 0; + for (int64_t xay3 = blockStart; xay3 <= xay2; xay3++) { + FastDiagonal *cur = fastForward_get(&fwd, xay3); + cur->where = FAST_BLOCK; + cur->offset = fwd.blockCellsLength; + fwd.blockCellsLength += cur->width * ns; + assert(fwd.blockCellsLength <= fwd.blockCellsCapacity); + FastDiagonal *m1 = fastForward_get(&fwd, xay3 - 1), *m2 = fastForward_get(&fwd, xay3 - 2); + fastForwardDiagonal(&hmm, xay3, cur, fastForward_cells(&fwd, cur), m1, fastForward_cells(&fwd, m1), + m2, m2 != NULL ? fastForward_cells(&fwd, m2) : NULL, &symbols, em); + } + } + fastPosteriors(sM, xay2, f, fastForward_cells(&fwd, f), b, bCells, logTotal, p, alignedPairs, gapXPairs, gapYPairs); + } + if (!atEnd) { + fastForward_keepFrom(&fwd, chunk->tracedBackFrom); + } + } + + free(chunks); + fastForward_destruct(&fwd); + fastSymbols_destruct(symbols); + free(em); + free(scratch); + for (int64_t i = 0; i < 3; i++) { + free(backCells[i]); + } + band_destruct(band); +#if defined(__SSE2__) + _mm_setcsr(oldCsr); +#endif + return 1; +} + /////////////////////////////////// /////////////////////////////////// //Banded alignment routine to calculate posterior match probs @@ -779,6 +1566,18 @@ void getPosteriorProbsWithBanding(StateMachine *sM, stList *anchorPairs, const S return; } + //The posterior probabilities have a fast path. The expectations for EM training do not. + if (diagonalPosteriorProbFn == diagonalCalculationPosteriorMatchProbs && + getPosteriorProbsWithBandingFast(sM, anchorPairs, sX, sY, p, alignmentHasRaggedLeftEnd, alignmentHasRaggedRightEnd, + ((void **) extraArgs)[0], NULL, NULL)) { + return; + } + if (diagonalPosteriorProbFn == diagonalCalculationPosteriorProbs && + getPosteriorProbsWithBandingFast(sM, anchorPairs, sX, sY, p, alignmentHasRaggedLeftEnd, alignmentHasRaggedRightEnd, + ((void **) extraArgs)[0], ((void **) extraArgs)[2], ((void **) extraArgs)[4])) { + return; + } + //Primitives for the forward matrix recursion Band *band = p->dynamicAnchorExpansion ? band_constructDynamic(anchorPairs, sX.length, sY.length) : band_construct(anchorPairs, sX.length, sY.length, p->diagonalExpansion); BandIterator *forwardBandIterator = bandIterator_construct(band); @@ -1224,6 +2023,9 @@ stList *getAnchorPairsForPairwiseAlignmentParameters(const char *sX, const char if ((int64_t) lX * lY <= p->anchorMatrixBiggerThanThis) { return stList_construct(); } + if (p->anchorMethod == PAIRWISE_ANCHOR_SEED) { + return getSeedAnchors(sX, sY, lX, lY, p); + } if(p->useMumAnchors) { return getAlignedMums(sX, sY, lX, lY, p, 0, 0); } @@ -1386,6 +2188,11 @@ PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters_construct() { p->recursiveMums = 1; p->k = 50; p->u = 1; + p->anchorMethod = PAIRWISE_ANCHOR_LEGACY; + p->seedHspThreshold = 2000; + p->seedHspThresholdMin = 1000; + p->seedRecursionDepth = 2; + p->seedXDrop = 910; return p; } @@ -1442,6 +2249,21 @@ PairwiseAlignmentParameters *pairwiseAlignmentParameters_jsonParse(char *buf, si else if (strcmp(keyString, "dynamicAnchorExpansion") == 0) { params->dynamicAnchorExpansion = stJson_parseBool(js, tokens, ++tokenIndex); } + else if (strcmp(keyString, "anchorMethod") == 0) { + params->anchorMethod = stJson_parseInt(js, tokens, ++tokenIndex); + } + else if (strcmp(keyString, "seedHspThreshold") == 0) { + params->seedHspThreshold = stJson_parseInt(js, tokens, ++tokenIndex); + } + else if (strcmp(keyString, "seedHspThresholdMin") == 0) { + params->seedHspThresholdMin = stJson_parseInt(js, tokens, ++tokenIndex); + } + else if (strcmp(keyString, "seedRecursionDepth") == 0) { + params->seedRecursionDepth = stJson_parseInt(js, tokens, ++tokenIndex); + } + else if (strcmp(keyString, "seedXDrop") == 0) { + params->seedXDrop = stJson_parseInt(js, tokens, ++tokenIndex); + } else { st_errAbort("ERROR: Unrecognised key in pairwise alignment parameters json: %s\n", keyString); } @@ -1948,13 +2770,15 @@ static Mum *mum_construct(int64_t x, int64_t y, int64_t length) { } void mum_destruct(Mum *mum) { - assert(mum->refCount > 0); - mum->refCount--; - if(mum->refCount == 0) { - if(mum->pMum != NULL) { - mum_destruct(mum->pMum); + // A loop, not recursion: a chain of mums along a megabase is tens of thousands long + while (mum != NULL) { + assert(mum->refCount > 0); + if (--mum->refCount > 0) { + return; } + Mum *pMum = mum->pMum; free(mum); + mum = pMum; } } @@ -2121,3 +2945,382 @@ stList *getAlignedMums(const char *sX, const char *sY, int64_t lX, int64_t lY, P } +/////////////////////////////////// +/////////////////////////////////// +//Anchoring by spaced seeds, ungapped extension and chaining +// +//What the lastz call above does -- 12of19 spaced seeds allowing one transition, HOXD70 scored +//ungapped extension, and the heaviest colinear chain of the resulting HSPs -- but in process. +//lastz costs a process, two temporary files and ~50 ms of startup for every pairwise alignment +//it anchors. The MUM anchors, which were the way around that, accept "unique" matches as short +//as a dozen bases, and at the divergence of a deep branch most of those are random: simulated +//at 0.3 substitutions per site per branch, half the MUM anchors were off the true alignment, +//and the banded dp, which cannot leave the band they define, followed them. +/////////////////////////////////// +/////////////////////////////////// + +#define SEED_SPAN 19 +#define SEED_WEIGHT 12 +#define SEED_MAX_WORD_COUNT 64 // Y positions a seed word may have before it is taken to be a repeat and skipped + +static const int64_t seedCare[SEED_WEIGHT] = { 0, 1, 2, 4, 7, 8, 11, 13, 15, 16, 17, 18 }; // 1110100110010101111 + +static const int64_t seedHoxd70[5][5] = { // A C G T, then N against anything + { 91, -114, -31, -123, -100 }, + { -114, 100, -125, -31, -100 }, + { -31, -125, 100, -114, -100 }, + { -123, -31, -114, 91, -100 }, + { -100, -100, -100, -100, -100 } }; + +typedef struct _seedSequence { + uint8_t *codes; // 0-3 for ACGT, 4 for anything else + int32_t *masked; // masked[i] is the number of lower case bases before i +} SeedSequence; + +static SeedSequence seedSequence_construct(const char *s, int64_t l) { + SeedSequence seq; + seq.codes = st_malloc(l + 1); + seq.masked = st_malloc(sizeof(int32_t) * (l + 1)); + seq.masked[0] = 0; + for (int64_t i = 0; i < l; i++) { + switch (s[i]) { + case 'A': case 'a': seq.codes[i] = 0; break; + case 'C': case 'c': seq.codes[i] = 1; break; + case 'G': case 'g': seq.codes[i] = 2; break; + case 'T': case 't': seq.codes[i] = 3; break; + default: seq.codes[i] = 4; + } + seq.masked[i + 1] = seq.masked[i] + (islower((unsigned char) s[i]) ? 1 : 0); + } + return seq; +} + +static void seedSequence_destruct(SeedSequence seq) { + free(seq.codes); + free(seq.masked); +} + +/* + * The seed word at i, over the seed's care positions. False if the span has an N, or, when + * repeat masking, a lower case base. + */ +static inline bool seed_word(const SeedSequence *seq, int64_t i, bool repeatMask, uint32_t *word) { + if (repeatMask && seq->masked[i + SEED_SPAN] != seq->masked[i]) { + return 0; + } + uint32_t w = 0; + for (int64_t j = 0; j < SEED_WEIGHT; j++) { + uint8_t c = seq->codes[i + seedCare[j]]; + if (c > 3) { + return 0; + } + w |= ((uint32_t) c) << (2 * j); + } + *word = w; + return 1; +} + +typedef struct _seedEntry { + uint32_t word; + int32_t pos; +} SeedEntry; + +static int seedEntry_cmp(const void *a, const void *b) { + const SeedEntry *e = a, *f = b; + if (e->word != f->word) { + return e->word < f->word ? -1 : 1; + } + return e->pos < f->pos ? -1 : (e->pos > f->pos ? 1 : 0); +} + +typedef struct _seedHsp { + int64_t x, y, length, score; +} SeedHsp; + +/* + * Extends a seed hit at (x, y) both ways along its diagonal, stopping each way once the score + * falls xDrop below the best seen, and returns the best scoring segment. + */ +static SeedHsp seed_extend(const uint8_t *cX, int64_t x0, int64_t x1, const uint8_t *cY, int64_t y0, int64_t y1, + int64_t x, int64_t y, int64_t xDrop) { + int64_t s = 0, best = 0, right = 0; + for (int64_t k = 0; x + k < x1 && y + k < y1; k++) { + s += seedHoxd70[cX[x + k]][cY[y + k]]; + if (s > best) { + best = s; + right = k + 1; + } else if (s < best - xDrop) { + break; + } + } + int64_t s2 = 0, best2 = 0, left = 0; + for (int64_t k = 1; x - k >= x0 && y - k >= y0; k++) { + s2 += seedHoxd70[cX[x - k]][cY[y - k]]; + if (s2 > best2) { + best2 = s2; + left = k; + } else if (s2 < best2 - xDrop) { + break; + } + } + SeedHsp hsp = { x - left, y - left, left + right, best + best2 }; + return hsp; +} + +/* + * The HSPs scoring at least threshold between X[x0, x1) and Y[y0, y1). + */ +static SeedHsp *seed_getHsps(SeedSequence *sX, int64_t x0, int64_t x1, SeedSequence *sY, int64_t y0, int64_t y1, + bool repeatMask, int64_t threshold, int64_t xDrop, int64_t *hspNumber) { + *hspNumber = 0; + int64_t lX = x1 - x0, lY = y1 - y0; + if (lX < SEED_SPAN || lY < SEED_SPAN) { + return NULL; + } + // Index the words of Y, sorted, with a table of where each bucket of leading word bits starts + SeedEntry *entries = st_malloc(sizeof(SeedEntry) * (lY - SEED_SPAN + 1)); + int64_t entryNumber = 0; + for (int64_t y = y0; y <= y1 - SEED_SPAN; y++) { + uint32_t w; + if (seed_word(sY, y, repeatMask, &w)) { + entries[entryNumber].word = w; + entries[entryNumber++].pos = (int32_t) y; + } + } + qsort(entries, entryNumber, sizeof(SeedEntry), seedEntry_cmp); + int64_t bucketBits = 4; + while (bucketBits < 20 && ((int64_t) 1 << bucketBits) < entryNumber) { + bucketBits++; + } + int64_t shift = 2 * SEED_WEIGHT - bucketBits; + int64_t *bucketStart = st_malloc(sizeof(int64_t) * (((int64_t) 1 << bucketBits) + 1)); + for (int64_t b = 0, e = 0; b <= ((int64_t) 1 << bucketBits); b++) { + while (e < entryNumber && (int64_t) (entries[e].word >> shift) < b) { + e++; + } + bucketStart[b] = e; + } + + // Walk X, looking up each word and the words one transition away from it. A hit on a + // diagonal already extended past it is part of an HSP that has been found. + int32_t *diagonalEnd = st_malloc(sizeof(int32_t) * (lX + lY)); + for (int64_t i = 0; i < lX + lY; i++) { + diagonalEnd[i] = INT32_MIN; + } + int64_t hspCapacity = 64; + SeedHsp *hsps = st_malloc(sizeof(SeedHsp) * hspCapacity); + for (int64_t x = x0; x <= x1 - SEED_SPAN; x++) { + uint32_t w; + if (!seed_word(sX, x, repeatMask, &w)) { + continue; + } + for (int64_t v = 0; v <= SEED_WEIGHT; v++) { + uint32_t w2 = v == 0 ? w : w ^ (2u << (2 * (v - 1))); // A<->G and C<->T are the same bit flip + int64_t b = w2 >> shift, lo = bucketStart[b], hi = bucketStart[b + 1]; + while (lo < hi && entries[lo].word < w2) { + lo++; + } + int64_t e = lo; + while (e < hi && entries[e].word == w2) { + e++; + } + if (e - lo > SEED_MAX_WORD_COUNT) { + continue; + } + for (int64_t k = lo; k < e; k++) { + int64_t y = entries[k].pos; + int64_t d = (x - x0) - (y - y0) + lY - 1; + if (diagonalEnd[d] > x) { + continue; + } + SeedHsp hsp = seed_extend(sX->codes, x0, x1, sY->codes, y0, y1, x, y, xDrop); + diagonalEnd[d] = hsp.x + hsp.length > x + 1 ? (int32_t) (hsp.x + hsp.length) : (int32_t) (x + 1); + if (hsp.length > 0 && hsp.score >= threshold) { // a zero-length "HSP", possible only with a threshold of 0, would break the chain's sweep + if (*hspNumber == hspCapacity) { + hspCapacity *= 2; + hsps = st_realloc(hsps, sizeof(SeedHsp) * hspCapacity); + } + hsps[(*hspNumber)++] = hsp; + } + } + } + } + free(entries); + free(bucketStart); + free(diagonalEnd); + return hsps; +} + +static int seedHsp_cmpByX(const void *a, const void *b) { + const SeedHsp *h = a, *i = b; + if (h->x != i->x) { + return h->x < i->x ? -1 : 1; + } + return h->y < i->y ? -1 : (h->y > i->y ? 1 : 0); +} + +static int seed_cmpInt64(const void *a, const void *b) { + int64_t k = *(int64_t *) a, l = *(int64_t *) b; + return k < l ? -1 : (k > l ? 1 : 0); +} + +/* + * Replaces hsps with the heaviest chain of them that is colinear and non-overlapping in both + * sequences, in order. A sweep along X, with a Fenwick tree of the best chain ending at or + * before each Y position. + */ +static void seed_chain(SeedHsp *hsps, int64_t *hspNumber) { + int64_t n = *hspNumber; + if (n <= 1) { + return; + } + qsort(hsps, n, sizeof(SeedHsp), seedHsp_cmpByX); + int64_t *ends = st_malloc(sizeof(int64_t) * n), *byEnd = st_malloc(sizeof(int64_t) * n); + int64_t *yEnds = st_malloc(sizeof(int64_t) * n); + for (int64_t i = 0; i < n; i++) { + ends[i] = hsps[i].x + hsps[i].length; + byEnd[i] = i; + yEnds[i] = hsps[i].y + hsps[i].length; + } + // The HSPs in order of where they end in X, by sorting (end, index) pairs + int64_t *pairs = st_malloc(sizeof(int64_t) * 2 * n); + for (int64_t i = 0; i < n; i++) { + pairs[2 * i] = ends[i]; + pairs[2 * i + 1] = i; + } + qsort(pairs, n, 2 * sizeof(int64_t), seed_cmpInt64); // compares the end, the first of each pair + for (int64_t i = 0; i < n; i++) { + byEnd[i] = pairs[2 * i + 1]; + } + free(pairs); + qsort(yEnds, n, sizeof(int64_t), seed_cmpInt64); + int64_t m = 0; + for (int64_t i = 0; i < n; i++) { + if (m == 0 || yEnds[m - 1] != yEnds[i]) { + yEnds[m++] = yEnds[i]; + } + } + int64_t *fenwickScore = st_calloc(m + 1, sizeof(int64_t)), *fenwickIndex = st_malloc(sizeof(int64_t) * (m + 1)); + int64_t *best = st_malloc(sizeof(int64_t) * n), *prev = st_malloc(sizeof(int64_t) * n); + int64_t e = 0, bestEnd = -1; + for (int64_t i = 0; i < n; i++) { + // Make the chains ending at or before this HSP's start in X available + while (e < n && ends[byEnd[e]] <= hsps[i].x) { + int64_t j = byEnd[e++]; + int64_t yEnd = hsps[j].y + hsps[j].length; + int64_t lo = 0, hi = m; // rank of yEnd, 1 based + while (lo < hi) { + int64_t mid = (lo + hi) / 2; + if (yEnds[mid] < yEnd) { + lo = mid + 1; + } else { + hi = mid; + } + } + for (int64_t r = lo + 1; r <= m; r += r & -r) { + if (best[j] > fenwickScore[r]) { + fenwickScore[r] = best[j]; + fenwickIndex[r] = j; + } + } + } + // The best of them ending at or before its start in Y + int64_t lo = 0, hi = m; // number of yEnds <= hsps[i].y + while (lo < hi) { + int64_t mid = (lo + hi) / 2; + if (yEnds[mid] <= hsps[i].y) { + lo = mid + 1; + } else { + hi = mid; + } + } + int64_t s = 0, p = -1; + for (int64_t r = lo; r > 0; r -= r & -r) { + if (fenwickScore[r] > s) { + s = fenwickScore[r]; + p = fenwickIndex[r]; + } + } + best[i] = hsps[i].score + s; + prev[i] = p; + if (bestEnd == -1 || best[i] > best[bestEnd]) { + bestEnd = i; + } + } + int64_t chainLength = 0; + for (int64_t i = bestEnd; i != -1; i = prev[i]) { + chainLength++; + } + SeedHsp *chain = st_malloc(sizeof(SeedHsp) * chainLength); + for (int64_t i = bestEnd, k = chainLength - 1; i != -1; i = prev[i], k--) { + chain[k] = hsps[i]; + } + memcpy(hsps, chain, sizeof(SeedHsp) * chainLength); + *hspNumber = chainLength; + free(chain); + free(ends); + free(byEnd); + free(yEnds); + free(fenwickScore); + free(fenwickIndex); + free(best); + free(prev); +} + +/* + * The HSP score needed in a box of the given area. How high a score chance alone reaches grows + * with the log of the area searched, so a whole megabase-square problem needs a high one -- on + * unrelated random sequence, 2000 still lets 16 HSPs into the chain of a 1 Mb x 1 Mb search, and + * 2500 none -- while the gap between two anchors, a few hundred bases square and already known to + * sit between homologous sequence, can take a much lower one: 1000 finds nothing by chance in a + * 1 kb square. Rises 250 per tenfold of area from minThreshold at 1e6, capped at maxThreshold. + */ +static int64_t seed_threshold(int64_t area, int64_t minThreshold, int64_t maxThreshold) { + double t = minThreshold + (area > 1000000 ? 250.0 * log10((double) area / 1000000.0) : 0.0); + return t > maxThreshold ? maxThreshold : (int64_t) t; +} + +/* + * Appends to anchors, in order, the anchors between X[x0, x1) and Y[y0, y1): the chained HSPs, + * less constraintDiagonalTrim at each end, and, down to seedRecursionDepth levels, the anchors + * of the gaps between them that are big enough to be worth anchoring themselves. The whole box + * takes seedHspThreshold; the gaps within it, flanked by anchors, the lower seed_threshold. + */ +static void seed_anchorBox(SeedSequence *sX, int64_t x0, int64_t x1, SeedSequence *sY, int64_t y0, int64_t y1, + PairwiseAlignmentParameters *p, bool repeatMask, int64_t depth, stList *anchors) { + int64_t hspNumber; + int64_t threshold = depth == 0 ? p->seedHspThreshold : + seed_threshold((x1 - x0) * (y1 - y0), p->seedHspThresholdMin, p->seedHspThreshold); + SeedHsp *hsps = seed_getHsps(sX, x0, x1, sY, y0, y1, repeatMask, threshold, p->seedXDrop, &hspNumber); + seed_chain(hsps, &hspNumber); + int64_t pX = x0, pY = y0; + for (int64_t i = 0; i <= hspNumber; i++) { + int64_t nX = i < hspNumber ? hsps[i].x : x1, nY = i < hspNumber ? hsps[i].y : y1; + int64_t gapSize = (nX - pX) * (nY - pY); + if (depth < p->seedRecursionDepth && gapSize > p->anchorMatrixBiggerThanThis && hspNumber > 0) { // not the whole box again + seed_anchorBox(sX, pX, nX, sY, pY, nY, p, gapSize > p->repeatMaskMatrixBiggerThanThis, depth + 1, anchors); + } + if (i < hspNumber) { + for (int64_t k = p->constraintDiagonalTrim; k < hsps[i].length - p->constraintDiagonalTrim; k++) { + stList_append(anchors, stIntTuple_construct3(hsps[i].x + k, hsps[i].y + k, p->diagonalExpansion)); + } + pX = hsps[i].x + hsps[i].length; + pY = hsps[i].y + hsps[i].length; + } + } + free(hsps); +} + +stList *getSeedAnchors(const char *sX, const char *sY, int64_t lX, int64_t lY, PairwiseAlignmentParameters *p) { + stList *anchors = stList_construct3(0, (void (*)(void *)) stIntTuple_destruct); + if (lX >= INT32_MAX || lY >= INT32_MAX) { // the index keeps positions in 32 bits; bar never asks for more than a megabase + st_logCritical("getSeedAnchors: sequences of %" PRIi64 " and %" PRIi64 " bases are too long to anchor\n", lX, lY); + return anchors; + } + SeedSequence seqX = seedSequence_construct(sX, lX), seqY = seedSequence_construct(sY, lY); + // Masked or not by the same rule as the gaps within: repeatMaskMatrixBiggerThanThis + seed_anchorBox(&seqX, 0, lX, &seqY, 0, lY, p, lX * lY > p->repeatMaskMatrixBiggerThanThis, 0, anchors); + seedSequence_destruct(seqX); + seedSequence_destruct(seqY); + return anchors; +} diff --git a/impl/stateMachine.c b/impl/stateMachine.c index dc0d548..7d7e4fd 100644 --- a/impl/stateMachine.c +++ b/impl/stateMachine.c @@ -516,6 +516,7 @@ StateMachine *stateMachine5_construct(StateMachineType type) { sM5->model.raggedStartStateProb = stateMachine5_raggedStartStateProb; sM5->model.raggedEndStateProb = stateMachine5_raggedEndStateProb; sM5->model.cellCalculate = stateMachine5_cellCalculate; + sM5->model.flat = NULL; return (StateMachine *) sM5; } @@ -740,6 +741,7 @@ StateMachine *stateMachine3_construct(StateMachineType type) { sM3->model.raggedStartStateProb = stateMachine3_raggedStartStateProb; sM3->model.raggedEndStateProb = stateMachine3_raggedEndStateProb; sM3->model.cellCalculate = stateMachine3_cellCalculate; + sM3->model.flat = NULL; return (StateMachine *) sM3; } @@ -819,5 +821,6 @@ StateMachine *hmm_getStateMachine(Hmm *hmm) { } void stateMachine_destruct(StateMachine *stateMachine) { + free(stateMachine->flat); free(stateMachine); } diff --git a/inc/pairwiseAligner.h b/inc/pairwiseAligner.h index efb7e9f..c89047e 100644 --- a/inc/pairwiseAligner.h +++ b/inc/pairwiseAligner.h @@ -51,8 +51,16 @@ typedef struct _pairwiseAlignmentBandingParameters { int64_t k; // k-mer length used by MUM anchoring approach int64_t u; // a MUM must be 1 + u longer than any other match to be considered unique int64_t recursiveMums; // anchor recursively + int64_t anchorMethod; // PAIRWISE_ANCHOR_SEED, or PAIRWISE_ANCHOR_LEGACY to use lastz or MUMs as useMumAnchors says + int64_t seedHspThreshold; // HOXD70 score an ungapped seed extension needs to be an anchor candidate (seed anchoring) + int64_t seedHspThresholdMin; // the least that threshold falls to in the small gaps between anchors (seed anchoring) + int64_t seedRecursionDepth; // how many levels of gaps between anchors are searched again (seed anchoring) + int64_t seedXDrop; // x-drop for the ungapped extensions (seed anchoring) } PairwiseAlignmentParameters; +#define PAIRWISE_ANCHOR_LEGACY 0 +#define PAIRWISE_ANCHOR_SEED 1 + PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters_construct(); void pairwiseAlignmentBandingParameters_destruct(PairwiseAlignmentParameters *p); @@ -267,6 +275,13 @@ int cmpKmers(const char *k1, const char *k2, int64_t k, int64_t *matchLength); stList *getAlignedMums(const char *sX, const char *sY, int64_t lX, int64_t lY, PairwiseAlignmentParameters *p, int64_t offsetX, int64_t offsetY); +/* + * Anchors from spaced seeds, ungapped extension and chaining, done in process: the same kind of + * anchors as the lastz call makes, without the process. Returns (x, y, expansion) tuples, + * strictly increasing in both x and y. + */ +stList *getSeedAnchors(const char *sX, const char *sY, int64_t lX, int64_t lY, PairwiseAlignmentParameters *p); + stList *getBlastPairs(const char *sX, const char *sY, int64_t lX, int64_t lY, PairwiseAlignmentParameters *p, bool repeatMask); stList *getBlastPairsForPairwiseAlignmentParameters(const char *sX, const char *sY, const int64_t lX, const int64_t lY, diff --git a/inc/stateMachine.h b/inc/stateMachine.h index 683fcaf..053a9eb 100644 --- a/inc/stateMachine.h +++ b/inc/stateMachine.h @@ -52,6 +52,10 @@ struct _stateMachine { //Cells (states at a given coordinate( void (*cellCalculate)(StateMachine *sM, double *current, double *lower, double *middle, double *upper, Symbol cX, Symbol cY, void(*doTransition)(double *, double *, int64_t, int64_t, double, double, void *), void *extraArgs); + + // The machine flattened into arrays for the probability space forward-backward, made by pairwiseAligner.c the + // first time it is needed. A machine's parameters are not to be changed once it has been used. + void *flat; }; /* diff --git a/tests/pairwiseAlignerTest.c b/tests/pairwiseAlignerTest.c index c473cd0..b71bc65 100644 --- a/tests/pairwiseAlignerTest.c +++ b/tests/pairwiseAlignerTest.c @@ -525,6 +525,9 @@ void test_filterToRemoveOverlap(CuTest *testCase) { int64_t y = stIntTuple_get(pair, 1); bool nonOverlapping = 1; for (int64_t j = 0; j < stList_length(pairs); j++) { + if (j == i) { // a pair does not overlap itself -- without this every pair did, and the expected set was always empty + continue; + } stIntTuple *pair2 = stList_get(pairs, j); int64_t x2 = stIntTuple_get(pair2, 0); int64_t y2 = stIntTuple_get(pair2, 1); @@ -1201,6 +1204,122 @@ void test_computeForwardProbability(CuTest *testCase) { } } +/* + * A posterior callback the fast path does not recognise, so getPosteriorProbsWithBanding takes the + * original log space route with it. + */ +static void logSpacePosteriorMatchProbs(StateMachine *sM, int64_t xay, DpMatrix *forwardDpMatrix, DpMatrix *backwardDpMatrix, + const SymbolString sX, const SymbolString sY, double totalProbability, + PairwiseAlignmentParameters *p, void *extraArgs) { + diagonalCalculationPosteriorMatchProbs(sM, xay, forwardDpMatrix, backwardDpMatrix, sX, sY, totalProbability, p, extraArgs); +} + +static void checkPosteriorsAgree(CuTest *testCase, stList *pairs, stList *otherPairs, double minProb, double tolerance) { + stHash *other = stHash_construct3((uint64_t (*)(const void *)) stIntTuple_hashKey, + (int (*)(const void *, const void *)) stIntTuple_equalsFn, + (void (*)(void *)) stIntTuple_destruct, NULL); + for (int64_t i = 0; i < stList_length(otherPairs); i++) { + stIntTuple *pair = stList_get(otherPairs, i); + stHash_insert(other, stIntTuple_construct2(stIntTuple_get(pair, 1), stIntTuple_get(pair, 2)), pair); + } + for (int64_t i = 0; i < stList_length(pairs); i++) { + stIntTuple *pair = stList_get(pairs, i); + double prob = (double) stIntTuple_get(pair, 0) / PAIR_ALIGNMENT_PROB_1; + stIntTuple *key = stIntTuple_construct2(stIntTuple_get(pair, 1), stIntTuple_get(pair, 2)); + stIntTuple *otherPair = stHash_search(other, key); + stIntTuple_destruct(key); + if (prob >= minProb) { + CuAssertTrue(testCase, otherPair != NULL); + } + if (otherPair != NULL) { + double otherProb = (double) stIntTuple_get(otherPair, 0) / PAIR_ALIGNMENT_PROB_1; + CuAssertTrue(testCase, fabs(prob - otherProb) <= tolerance); + } + } + stHash_destruct(other); +} + +void test_fastPosteriorsMatchLogSpace(CuTest *testCase) { + /* + * The probability space forward-backward against the log space one it replaced, over random + * bands, traceback settings, ragged ends and state machines. They differ only by the error of + * the log space version's approximate log-add, a few thousandths. The last few are long and + * unanchored, which is what makes the fast path recompute its forward values in blocks. + */ + StateMachineType types[4] = { fiveState, fiveStateAsymmetric, threeState, threeStateAsymmetric }; + for (int64_t test = 0; test < 60; test++) { + bool long_ = test >= 56; + char *sX = getRandomSequence(long_ ? st_randomInt(1000, 1500) : st_randomInt(0, 150)); + char *sY = evolveSequence(sX); + int64_t lX = strlen(sX), lY = strlen(sY); + SymbolString sX2 = symbolString_construct(sX, lX); + SymbolString sY2 = symbolString_construct(sY, lY); + PairwiseAlignmentParameters *p = pairwiseAlignmentBandingParameters_construct(); + p->traceBackDiagonals = st_randomInt(1, 10); + p->minDiagsBetweenTraceBack = p->traceBackDiagonals + st_randomInt(2, 10); + p->diagonalExpansion = st_randomInt(0, 10) * 2; + p->dynamicAnchorExpansion = st_random() > 0.5; + StateMachineType type = types[st_randomInt(0, 4)]; + StateMachine *sM = type == fiveState || type == fiveStateAsymmetric ? stateMachine5_construct(type) : stateMachine3_construct(type); + stList *anchorPairs = long_ ? stList_construct3(0, (void (*)(void *)) stIntTuple_destruct) : getRandomAnchorPairs(lX, lY); + bool raggedLeft = st_random() > 0.5, raggedRight = st_random() > 0.5; + + stList *fastPairs = stList_construct3(0, (void (*)(void *)) stIntTuple_destruct); + stList *logPairs = stList_construct3(0, (void (*)(void *)) stIntTuple_destruct); + void *fastArgs[1] = { fastPairs }, *logArgs[1] = { logPairs }; + getPosteriorProbsWithBanding(sM, anchorPairs, sX2, sY2, p, raggedLeft, raggedRight, diagonalCalculationPosteriorMatchProbs, fastArgs); + getPosteriorProbsWithBanding(sM, anchorPairs, sX2, sY2, p, raggedLeft, raggedRight, logSpacePosteriorMatchProbs, logArgs); + st_logInfo("Fast and log space posteriors for lengths %" PRIi64 " and %" PRIi64 ": %" PRIi64 " and %" PRIi64 " pairs\n", + lX, lY, stList_length(fastPairs), stList_length(logPairs)); + checkAlignedPairs(testCase, fastPairs, lX, lY, 0, 0); + checkPosteriorsAgree(testCase, fastPairs, logPairs, 0.05, 0.02); + checkPosteriorsAgree(testCase, logPairs, fastPairs, 0.05, 0.02); + + stList_destruct(fastPairs); + stList_destruct(logPairs); + stList_destruct(anchorPairs); + stateMachine_destruct(sM); + pairwiseAlignmentBandingParameters_destruct(p); + free(sX); + free(sY); + free(sX2.sequence); + free(sY2.sequence); + } +} + +void test_getSeedAnchors(CuTest *testCase) { + /* + * The anchors are in range, strictly increasing in both sequences, and for a sequence against + * itself all on the one diagonal. + */ + for (int64_t test = 0; test < 20; test++) { + // Half against an evolved copy, half against itself. getRandomSequence is half lower case, + // which seeds do not start in, so the copies against themselves are upper case throughout. + char *seqX = test % 2 == 0 ? getRandomSequence(st_randomInt(0, 20000)) : getRandomACGTSequence(st_randomInt(0, 20000)); + char *seqY = test % 2 == 0 ? evolveSequence(seqX) : stString_copy(seqX); + int64_t lX = strlen(seqX), lY = strlen(seqY); + PairwiseAlignmentParameters *p = pairwiseAlignmentBandingParameters_construct(); + p->constraintDiagonalTrim = st_randomInt(0, 15); + p->diagonalExpansion = st_randomInt(0, 5) * 2; + stList *anchors = getSeedAnchors(seqX, seqY, lX, lY, p); + st_logInfo("Got %" PRIi64 " seed anchors for lengths %" PRIi64 " and %" PRIi64 "\n", stList_length(anchors), lX, lY); + checkBlastPairs(testCase, anchors, lX, lY, p->diagonalExpansion, 1); + if (test % 2 == 1) { + for (int64_t i = 0; i < stList_length(anchors); i++) { + stIntTuple *anchor = stList_get(anchors, i); + CuAssertIntEquals(testCase, stIntTuple_get(anchor, 0), stIntTuple_get(anchor, 1)); + } + if (lX > 1000) { // a random sequence against itself has plenty to anchor on + CuAssertTrue(testCase, stList_length(anchors) > lX / 2); + } + } + stList_destruct(anchors); + free(seqX); + free(seqY); + pairwiseAlignmentBandingParameters_destruct(p); + } +} + CuSuite* pairwiseAlignmentTestSuite(void) { CuSuite* suite = CuSuiteNew(); @@ -1213,6 +1332,8 @@ CuSuite* pairwiseAlignmentTestSuite(void) { SUITE_ADD_TEST(suite, test_dpMatrix); SUITE_ADD_TEST(suite, test_diagonalDPCalculations); SUITE_ADD_TEST(suite, test_getAlignedPairsWithBanding); + SUITE_ADD_TEST(suite, test_fastPosteriorsMatchLogSpace); + SUITE_ADD_TEST(suite, test_getSeedAnchors); SUITE_ADD_TEST(suite, test_getBlastPairs); SUITE_ADD_TEST(suite, test_getBlastPairsWithRecursion); SUITE_ADD_TEST(suite, test_filterToRemoveOverlap);