Skip to content

Per-pair state machines in the multiple aligner, and make it deterministic #10

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
glennhickey wants to merge 2 commits into from
Closed

Per-pair state machines in the multiple aligner, and make it deterministic #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ccdbc32
Per-pair state machines in the multiple aligner
glennhickey Oct 1, 2026
2ec01d1
Make the multiple aligner deterministic
glennhickey Oct 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 106 additions & 22 deletions impl/multipleAligner.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,38 @@ int column_cmp(Column *c, Column *c2) {
: (c->position < c2->position ? -1 : 0)));
}

/*
* Deterministic stand-ins for the random numbers this file used to draw from sonLib's global generator.
* Those made an alignment depend on that generator's state -- so on whatever had run before it, and in a
* multithreaded caller on thread timing -- and pointer comparisons made it depend on where things had been
* allocated. Now the same sequences give the same alignment.
*/
static uint64_t mix64(uint64_t h) {
h ^= h >> 30;
h *= 0xBF58476D1CE4E5B9ULL;
h ^= h >> 27;
h *= 0x94D049BB133111EBULL;
return h ^ (h >> 31);
}

static uint64_t pairHash(int64_t seq1, int64_t pos1, int64_t seq2, int64_t pos2) {
// the same either way round
if (seq1 > seq2 || (seq1 == seq2 && pos1 > pos2)) {
int64_t s = seq1, p = pos1;
seq1 = seq2, pos1 = pos2, seq2 = s, pos2 = p;
}
return mix64(mix64((uint64_t) seq1 * 0x9E3779B97F4A7C15ULL ^ (uint64_t) pos1) ^ ((uint64_t) seq2 * 0xC2B2AE3D27D4EB4FULL + (uint64_t) pos2));
}

/*
* The tie-breaking jitter added to an aligned pair's weight, in [0, 0.00001) like the random one it
* replaces (which kept the merge order from degenerating into badly unbalanced trees, without affecting
* accuracy). Both directions of the pair get the same value.
*/
static double weightJitter(int64_t seq1, int64_t pos1, int64_t seq2, int64_t pos2) {
return (double) (pairHash(seq1, pos1, seq2, pos2) >> 11) / 9007199254740992.0 * 0.00001;
}

uint64_t column_hashFn(Column *c) {
// 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
Expand Down Expand Up @@ -116,14 +148,23 @@ int alignmentWeight_cmpByPosition(AlignmentWeight *aW, AlignmentWeight *aW2) {

int alignmentWeight_cmpByWeight(AlignmentWeight *aW, AlignmentWeight *aW2) {
/*
* Cmp primarily by weight, then position.
* Cmp primarily by weight, then by the columns the weights join (not by address, which would make
* the order depend on where the weights were allocated).
*/
if (aW->avgWeight == aW2->avgWeight) {
return aW < aW2 ? -1 : (aW > aW2 ? 1 : 0);
int i = column_cmp(aW->column, aW2->column);
return i != 0 ? i : column_cmp(aW->rWeight->column, aW2->rWeight->column);
}
return aW->avgWeight > aW2->avgWeight ? 1 : -1;
}

/*
* Of a weight and its reverse, the one kept in the set ordered by weight: the one from the lesser column.
*/
static bool alignmentWeight_isCanonical(AlignmentWeight *aW) {
return column_cmp(aW->column, aW->rWeight->column) < 0;
}

void insertWeight(AlignmentWeight *aW, stHash *alignmentWeightAdjLists) {
/*
* Inserts an alignment weight to an adjacency list.
Expand All @@ -143,12 +184,12 @@ static Column *getColumn(stSet *columns, int64_t seqName, int64_t position) {
return stSet_search(columns, &c);
}

static AlignmentWeight *makeAlignmentWeight(stSet *columns, int64_t score, int64_t seqName, int64_t position) {
static AlignmentWeight *makeAlignmentWeight(stSet *columns, double weight, int64_t seqName, int64_t position) {
AlignmentWeight *aW = st_malloc(sizeof(AlignmentWeight));
aW->column = getColumn(columns, seqName, position);
assert(aW->column != NULL);
aW->numberOfWeights = 1;
aW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; //This randomness avoids nasty types of unbalanced trees and doesn't really affect accuracy
aW->avgWeight = weight;
return aW;
}

Expand All @@ -161,8 +202,10 @@ stHash *makeAlignmentWeightAdjacencyLists(stSet *columns, stList *multipleAligne
/*Tuple of score, seq1, pos1, seq2, pos2 */
stIntTuple *aP = stList_get(multipleAlignedPairs, i);
assert(stIntTuple_length(aP) == 5);
AlignmentWeight *aW = makeAlignmentWeight(columns, stIntTuple_get(aP, 0), stIntTuple_get(aP, 1), stIntTuple_get(aP, 2));
aW->rWeight = makeAlignmentWeight(columns, stIntTuple_get(aP, 0), stIntTuple_get(aP, 3), stIntTuple_get(aP, 4));
double weight = ((double) stIntTuple_get(aP, 0)) / PAIR_ALIGNMENT_PROB_1 +
weightJitter(stIntTuple_get(aP, 1), stIntTuple_get(aP, 2), stIntTuple_get(aP, 3), stIntTuple_get(aP, 4));
AlignmentWeight *aW = makeAlignmentWeight(columns, weight, stIntTuple_get(aP, 1), stIntTuple_get(aP, 2));
aW->rWeight = makeAlignmentWeight(columns, weight, stIntTuple_get(aP, 3), stIntTuple_get(aP, 4));
aW->rWeight->rWeight = aW;
insertWeight(aW, alignmentWeightAdjLists);
insertWeight(aW->rWeight, alignmentWeightAdjLists);
Expand All @@ -184,7 +227,7 @@ stSortedSet *makeOrderedSetOfAlignmentWeights(stHash *alignmentWeightAdjLists) {
AlignmentWeight *aW;
while ((aW = stSortedSet_getNext(aWIt)) != NULL) {
assert(aW->column != aW->rWeight->column);
if (aW->column < aW->rWeight->column) {
if (alignmentWeight_isCanonical(aW)) {
stSortedSet_insert(alignmentWeightsOrderedByWeight, aW);
}
}
Expand All @@ -199,7 +242,7 @@ static void removeAlignmentFromSortedAlignmentWeights(AlignmentWeight *aW, stSor
* Removes the weight from the ordered set, if one is being kept (the progressive merge does not need it).
*/
if (alignmentWeightsOrderedByWeight != NULL) {
stSortedSet_remove(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight);
stSortedSet_remove(alignmentWeightsOrderedByWeight, alignmentWeight_isCanonical(aW) ? aW : aW->rWeight);
}
}

Expand All @@ -208,7 +251,7 @@ static void insertAlignmentIntoSortedAlignmentWeights(AlignmentWeight *aW, stSor
* Adds weight to the ordered set, if one is being kept.
*/
if (alignmentWeightsOrderedByWeight != NULL) {
stSortedSet_insert(alignmentWeightsOrderedByWeight, aW->column < aW->rWeight->column ? aW : aW->rWeight);
stSortedSet_insert(alignmentWeightsOrderedByWeight, alignmentWeight_isCanonical(aW) ? aW : aW->rWeight);
}
}

Expand Down Expand Up @@ -608,7 +651,7 @@ static ProgressiveMsa *progressiveMsa_construct(stList *seqFrags, stList *multip
m->isHead[m->seqOffsets[seq] + pos] = 1;
}
}
// The weights, made in the order, and with the random tie breaking, that makeAlignmentWeightAdjacencyLists uses
// The weights, made in the order, and with the tie-breaking jitter, 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
Expand All @@ -630,12 +673,14 @@ static ProgressiveMsa *progressiveMsa_construct(stList *seqFrags, stList *multip
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];
double weight = ((double) score) / PAIR_ALIGNMENT_PROB_1 +
weightJitter(stIntTuple_get(aP, 1), stIntTuple_get(aP, 2), stIntTuple_get(aP, 3), stIntTuple_get(aP, 4));
aW->column = m->columns[id1];
aW->numberOfWeights = 1;
aW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001; // as makeAlignmentWeight
aW->avgWeight = weight;
rW->column = m->columns[id2];
rW->numberOfWeights = 1;
rW->avgWeight = ((double) score) / PAIR_ALIGNMENT_PROB_1 + st_random() * 0.00001;
rW->avgWeight = weight;
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
Expand Down Expand Up @@ -1025,22 +1070,37 @@ static stList *convertMultipleAlignedPairsToAlignedPairs(stList *multipleAligned
return alignedPairs;
}

static int64_t addMultipleAlignedPairs(StateMachine *sM, int64_t sequence1, int64_t sequence2, stList *seqFrags, stList *multipleAlignedPairs,
/*
* The machines a multiple alignment scores its pairs with: sM, or what pairStateMachine picks for the pair
* when it is given and does not return NULL.
*/
typedef struct _pairStateMachines {
StateMachine *sM;
PairStateMachineFn pairStateMachine;
void *extraArgs;
} PairStateMachines;

static StateMachine *pairStateMachines_get(PairStateMachines *sMs, int64_t sequence1, int64_t sequence2) {
StateMachine *sM = sMs->pairStateMachine == NULL ? NULL : sMs->pairStateMachine(sequence1, sequence2, sMs->extraArgs);
return sM == NULL ? sMs->sM : sM;
}

static int64_t addMultipleAlignedPairs(PairStateMachines *sMs, int64_t sequence1, int64_t sequence2, stList *seqFrags, stList *multipleAlignedPairs,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
/*
* Computes a pairwise alignment and returns the pairwise match probabilities as tuples of (score, seq1, pos1, seq2, pos2).
*/
SeqFrag *seqFrag1 = stList_get(seqFrags, sequence1);
SeqFrag *seqFrag2 = stList_get(seqFrags, sequence2);
stList *alignedPairs = getAlignedPairs(sM, seqFrag1->seq, seqFrag2->seq, pairwiseAlignmentBandingParameters,
stList *alignedPairs = getAlignedPairs(pairStateMachines_get(sMs, sequence1, sequence2), seqFrag1->seq, seqFrag2->seq, pairwiseAlignmentBandingParameters,
seqFrag1->leftEndId != seqFrag2->leftEndId, seqFrag1->rightEndId != seqFrag2->rightEndId);
alignedPairs = reweightAlignedPairs2(alignedPairs, seqFrag1->length, seqFrag2->length, pairwiseAlignmentBandingParameters->gapGamma);
int64_t distance = getAlignmentScore(alignedPairs, seqFrag1->length, seqFrag2->length);
convertAlignedPairsToMultipleAlignedPairs(alignedPairs, multipleAlignedPairs, sequence1, sequence2);
return distance;
}

stList *makeAllPairwiseAlignments(StateMachine *sM, stList *seqFrags, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, stList **seqPairSimilarityScores) {
static stList *makeAllPairwiseAlignments2(PairStateMachines *sMs, stList *seqFrags, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, stList **seqPairSimilarityScores) {
/*
* Generate the set of pairwise alignments between the sequences.
*/
Expand All @@ -1049,20 +1109,25 @@ stList *makeAllPairwiseAlignments(StateMachine *sM, stList *seqFrags, PairwiseAl
int64_t seqNo = stList_length(seqFrags);
for (int64_t seq1 = 0; seq1 < seqNo; seq1++) {
for (int64_t seq2 = seq1 + 1; seq2 < seqNo; seq2++) {
stList_append(*seqPairSimilarityScores, stIntTuple_construct3(addMultipleAlignedPairs(sM, seq1, seq2, seqFrags, multipleAlignedPairs, pairwiseAlignmentBandingParameters), seq1, seq2));
stList_append(*seqPairSimilarityScores, stIntTuple_construct3(addMultipleAlignedPairs(sMs, seq1, seq2, seqFrags, multipleAlignedPairs, pairwiseAlignmentBandingParameters), seq1, seq2));
}
}
return multipleAlignedPairs;
}

MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags,
stList *makeAllPairwiseAlignments(StateMachine *sM, stList *seqFrags, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, stList **seqPairSimilarityScores) {
PairStateMachines sMs = { sM, NULL, NULL };
return makeAllPairwiseAlignments2(&sMs, seqFrags, pairwiseAlignmentBandingParameters, seqPairSimilarityScores);
}

static MultipleAlignment *makeAlignmentUsingAllPairs2(PairStateMachines *sMs, stList *seqFrags,
bool useProgressiveMerging, float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
/*
* Generate a multiple alignment considering all pairs of sequences.
*/
MultipleAlignment *mA = st_calloc(1, sizeof(MultipleAlignment));
mA->alignedPairs = makeAllPairwiseAlignments(sM, seqFrags, pairwiseAlignmentBandingParameters, &mA->chosenPairwiseAlignments);
mA->alignedPairs = makeAllPairwiseAlignments2(sMs, seqFrags, pairwiseAlignmentBandingParameters, &mA->chosenPairwiseAlignments);
if(stList_length(seqFrags) == 2 || useProgressiveMerging) { //Compute an optimum exactly
mA->columns = getMultipleSequenceAlignmentProgressive(seqFrags, mA->alignedPairs, matchGamma, mA->chosenPairwiseAlignments);
}
Expand All @@ -1073,6 +1138,13 @@ MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags
return mA;
}

MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags,
bool useProgressiveMerging, float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
PairStateMachines sMs = { sM, NULL, NULL };
return makeAlignmentUsingAllPairs2(&sMs, seqFrags, useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters);
}

void multipleAlignment_destruct(MultipleAlignment *mA) {
stList_destruct(mA->alignedPairs);
stSet_destruct(mA->columns);
Expand Down Expand Up @@ -1244,7 +1316,9 @@ int64_t getNextBestPair(int64_t seq1, int64_t *distanceCounts, int64_t seqNo,
for (int64_t seq2 = 0; seq2 < seqNo; seq2++) {
if (seq1 != seq2) {
double gain = distances[seq2] - subsPerSite(seq1, seq2, distanceCounts, seqNo);
if (gain > maxGain || (gain == maxGain && st_random() > 0.5) /*lame attempt to reduce always picking the same seq*/) {
// ties go by a hash of the pair, which spreads the choice over the sequences (the point of the coin
// flip this replaces) without making it differ from one run to the next
if (gain > maxGain || (gain == maxGain && pairHash(seq1, 0, seq2, 0) > pairHash(seq1, 0, maxGainSeq, 0))) {
stIntTuple *pairToAlign = makePairToAlign(seq1, seq2);
if (stSortedSet_search(chosenPairsOfSequencesToAlign, pairToAlign) == NULL) { //So that any pair is unique
maxGain = gain;
Expand All @@ -1262,12 +1336,22 @@ int64_t getNextBestPair(int64_t seq1, int64_t *distanceCounts, int64_t seqNo,
MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, int64_t spanningTrees, int64_t maxPairsToConsider,
bool useProgressiveMerging, float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
return makeAlignmentWithPairStateMachines(sM, NULL, NULL, seqFrags, spanningTrees, maxPairsToConsider,
useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters);
}

MultipleAlignment *makeAlignmentWithPairStateMachines(StateMachine *sM, PairStateMachineFn pairStateMachine, void *extraArgs,
stList *seqFrags, int64_t spanningTrees, int64_t maxPairsToConsider,
bool useProgressiveMerging, float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
/*
* Computes an MSA, making up to "spanningTrees"*no of seqs pairwise alignments.
*/
PairStateMachines sMsValue = { sM, pairStateMachine, extraArgs };
PairStateMachines *sMs = &sMsValue;
int64_t seqNo = stList_length(seqFrags);
if (spanningTrees * (seqNo - 1) >= (seqNo * (seqNo - 1)) / 2) { //Do all pairs if we can
return makeAlignmentUsingAllPairs(sM, seqFrags, useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters);
return makeAlignmentUsingAllPairs2(sMs, seqFrags, useProgressiveMerging, matchGamma, pairwiseAlignmentBandingParameters);
}
MultipleAlignment *mA = st_calloc(1, sizeof(MultipleAlignment));
mA->alignedPairs = stList_construct3(0, (void(*)(void *)) stIntTuple_destruct); //pairwise alignment pairs, with sequence indices
Expand All @@ -1279,7 +1363,7 @@ MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, int64_t spa
int64_t seqX = stIntTuple_get(pairToAlign, 0);
int64_t seqY = stIntTuple_get(pairToAlign, 1);
//We get pairwise alignments, for this first alignment we filter the pairs greedily to make them consistent
stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sM, seqX, seqY, seqFrags,
stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sMs, seqX, seqY, seqFrags,
mA->alignedPairs, pairwiseAlignmentBandingParameters),
seqX, seqY));
}
Expand All @@ -1304,7 +1388,7 @@ MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags, int64_t spa
assert(seq != otherSeq);
stIntTuple *pairToAlign = makePairToAlign(seq, otherSeq);
assert(stSortedSet_search(chosenPairwiseAlignmentsSet, pairToAlign) == NULL);
stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sM, seq, otherSeq, seqFrags, mA->alignedPairs, pairwiseAlignmentBandingParameters), seq, otherSeq));
stList_append(mA->chosenPairwiseAlignments, stIntTuple_construct3(addMultipleAlignedPairs(sMs, seq, otherSeq, seqFrags, mA->alignedPairs, pairwiseAlignmentBandingParameters), seq, otherSeq));
stSortedSet_insert(chosenPairwiseAlignmentsSet, pairToAlign);
}
}
Expand Down
18 changes: 18 additions & 0 deletions inc/multipleAligner.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,24 @@ MultipleAlignment *makeAlignment(StateMachine *sM, stList *seqFrags,
float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters);

/*
* Picks the machine to align one pair of sequences with, by their indices in seqFrags: lets one multiple
* alignment score different pairs with different models, e.g. by how diverged the pair's genomes are.
* Called once for each pairwise alignment, possibly from several threads if the caller aligns several
* MSAs at once. The machine must outlive the alignment, and the same machine can be returned for any
* number of pairs. NULL means the alignment's default machine.
*/
typedef StateMachine *(*PairStateMachineFn)(int64_t seqX, int64_t seqY, void *extraArgs);

/*
* As makeAlignment, scoring each pair with the machine pairStateMachine picks for it (sM for any pair it
* returns NULL for, and for every pair when pairStateMachine is NULL).
*/
MultipleAlignment *makeAlignmentWithPairStateMachines(StateMachine *sM, PairStateMachineFn pairStateMachine, void *extraArgs,
stList *seqFrags, int64_t spanningTrees, int64_t maxPairsToConsider,
bool useProgressiveMerging, float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters);

MultipleAlignment *makeAlignmentUsingAllPairs(StateMachine *sM, stList *seqFrags,
bool useProgressiveMerging, float matchGamma,
PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters);
Expand Down
Loading