Bar scoring by divergence: per-row abPOA and per-pair pecan models, trained per pair of genomes (prototype)

bar/impl/bar.c

@@ -108,6 +108,8 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList *
     float matchGamma = cactusParams_get_float(params, 3, "bar", "pecan", "matchGamma");
     PairwiseAlignmentParameters *pairwiseAlignmentParameters = pairwiseAlignmentParameters_constructFromCactusParams(params);
     StateMachine *sM = stateMachine5_construct(fiveState);
+    // per-pair pecan models, if the config has any (see endAligner.h)
+    PecanPairModels *pecanPairModels = usePoa ? NULL : pecanPairModels_constructFromCactusParams(params);
     bool pruneOutStubAlignments = cactusParams_get_int(params, 3, "bar", "pecan", "pruneOutStubAlignments");
 
     // Poa params.  The window and the mask filter come from whichever engine is selected:
@@ -233,7 +235,7 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList *
                                                    poaMaxProgRows, poaMaxLenDiff, poaParameters);
             st_logDebug("Created the poa alignments: %" PRIi64 " poa alignment blocks for flower\n", stList_length(alignments));
         } else {
-            alignments = makeFlowerAlignment3(sM, flower, listOfEndAlignmentFiles, spanningTrees, maximumLength,
+            alignments = makeFlowerAlignment3(sM, pecanPairModels, flower, listOfEndAlignmentFiles, spanningTrees, maximumLength,
                                               useProgressiveMerging, matchGamma, pairwiseAlignmentParameters,
                                               pruneOutStubAlignments);
             st_logDebug("Created the alignment: %" PRIi64 " pairs for flower\n", stSortedSet_size(alignments));
@@ -356,6 +358,7 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList *
 
     pairwiseAlignmentBandingParameters_destruct(pairwiseAlignmentParameters);
     stateMachine_destruct(sM);
+    pecanPairModels_destruct(pecanPairModels);
 
     poaParameters_destruct(poaParameters);
 }

bar/impl/endAligner.c

@@ -5,6 +5,7 @@
  */
 
 #include "endAligner.h"
+#include <math.h>
 #include "multipleAligner.h"
 #include "adjacencySequences.h"
 #include "pairwiseAligner.h"
@@ -52,16 +53,209 @@ int alignedPair_cmpFn(const AlignedPair *alignedPair1, const AlignedPair *aligne
     return i;
 }
 
+struct _barDistances {
+    stTree *tree;   // the reconstruction tree, or NULL
+    stHash *nodes;  // its nodes by label
+};
+
+static void barDistances_addNodes(stHash *nodes, stTree *tree) {
+    if (stTree_getLabel(tree) != NULL) {
+        stHash_insert(nodes, (void *)stTree_getLabel(tree), tree);
+    }
+    for (int64_t i = 0; i < stTree_getChildNumber(tree); i++) {
+        barDistances_addNodes(nodes, stTree_getChild(tree, i));
+    }
+}
+
+BarDistances *barDistances_constructFromCactusParams(CactusParams *params) {
+    BarDistances *distances = st_calloc(1, sizeof(BarDistances));
+    if (cactusParams_has(params, 2, "reference", "reconstructionTree")) {
+        char *newick = cactusParams_get_string(params, 2, "reference", "reconstructionTree");
+        distances->tree = stTree_parseNewickString(newick);
+        distances->nodes = stHash_construct3(stHash_stringKey, stHash_stringEqualKey, NULL, NULL);
+        barDistances_addNodes(distances->nodes, distances->tree);
+        free(newick);
+    }
+    return distances;
+}
+
+bool barDistances_hasReconstructionTree(BarDistances *distances) {
+    return distances->tree != NULL;
+}
+
+double barDistances_get(BarDistances *distances, Event *event1, Event *event2) {
+    stTree *node1 = distances->nodes == NULL ? NULL : stHash_search(distances->nodes, (void *)event_getHeader(event1));
+    stTree *node2 = distances->nodes == NULL ? NULL : stHash_search(distances->nodes, (void *)event_getHeader(event2));
+    if (node1 != NULL && node2 != NULL) {
+        double d1 = 0.0;
+        for (stTree *a = node1; a != NULL; a = stTree_getParent(a)) {
+            double d2 = 0.0;
+            for (stTree *b = node2; b != NULL; b = stTree_getParent(b)) {
+                if (a == b) {
+                    return d1 + d2;
+                }
+                d2 += stTree_getBranchLength(b);
+            }
+            d1 += stTree_getBranchLength(a);
+        }
+    }
+    Event *ancestor = eventTree_getCommonAncestor(event1, event2);
+    double d = 0.0;
+    for (Event *e = event1; e != ancestor; e = event_getParent(e)) {
+        d += event_getBranchLength(e);
+    }
+    for (Event *e = event2; e != ancestor; e = event_getParent(e)) {
+        d += event_getBranchLength(e);
+    }
+    return d;
+}
+
+void barDistances_destruct(BarDistances *distances) {
+    if (distances->tree != NULL) {
+        stHash_destruct(distances->nodes);
+        stTree_destruct(distances->tree);
+    }
+    free(distances);
+}
+
+void barModel_readGenomes(CactusParams *params, const char *section, const char *name, char **genomes) {
+    genomes[0] = genomes[1] = NULL;
+    if (!cactusParams_has(params, 4, "bar", section, name, "genomes")) {
+        return;
+    }
+    char *value = cactusParams_get_string(params, 4, "bar", section, name, "genomes");
+    stList *tokens = stString_split(value);
+    if (stList_length(tokens) != 2) {
+        st_errAbort("bar model %s/%s has genomes=\"%s\", which is not two genome names", section, name, value);
+    }
+    genomes[0] = stString_copy(stList_get(tokens, 0));
+    genomes[1] = stString_copy(stList_get(tokens, 1));
+    stList_destruct(tokens);
+    free(value);
+}
+
+bool barModel_isForGenomes(char **genomes, const char *a, const char *b) {
+    return genomes[0] != NULL && ((strcmp(genomes[0], a) == 0 && strcmp(genomes[1], b) == 0) ||
+                                  (strcmp(genomes[0], b) == 0 && strcmp(genomes[1], a) == 0));
+}
+
+struct _pecanPairModels {
+    int64_t modelNo;
+    double *maxDistance; // NAN for a model with none
+    char **genomes; // model m's two genomes at 2m and 2m+1, NULL for a model with none
+    StateMachine **stateMachines;
+    BarDistances *distances;
+};
+
+PecanPairModels *pecanPairModels_constructFromCactusParams(CactusParams *params) {
+    int64_t modelNo = cactusParams_has(params, 3, "bar", "pecan", "pecanPairModels") ?
+        cactusParams_get_int(params, 3, "bar", "pecan", "pecanPairModels") : 0;
+    if (modelNo <= 0) {
+        return NULL;
+    }
+    PecanPairModels *pairModels = st_calloc(1, sizeof(PecanPairModels));
+    pairModels->modelNo = modelNo;
+    pairModels->maxDistance = st_malloc(sizeof(double) * modelNo);
+    pairModels->genomes = st_malloc(sizeof(char *) * 2 * modelNo);
+    pairModels->stateMachines = st_malloc(sizeof(StateMachine *) * modelNo);
+    int64_t keyed = 0;
+    double lastMaxDistance = -INFINITY;
+    for (int64_t m = 0; m < modelNo; m++) {
+        char *name = stString_print("pairModel%" PRIi64, m);
+        barModel_readGenomes(params, "pecan", name, pairModels->genomes + 2 * m);
+        keyed += pairModels->genomes[2 * m] != NULL;
+        pairModels->maxDistance[m] = NAN;
+        if (cactusParams_has(params, 4, "bar", "pecan", name, "maxDistance")) {
+            pairModels->maxDistance[m] = cactusParams_get_float(params, 4, "bar", "pecan", name, "maxDistance");
+            if (pairModels->maxDistance[m] < lastMaxDistance) {
+                st_errAbort("pecan pair models must be in ascending maxDistance, but %s has %f after %f", name,
+                            pairModels->maxDistance[m], lastMaxDistance);
+            }
+            lastMaxDistance = pairModels->maxDistance[m];
+        } else if (pairModels->genomes[2 * m] == NULL) {
+            st_errAbort("pecan pair model %s has neither genomes nor a maxDistance", name);
+        }
+        char *json = cactusParams_get_string(params, 4, "bar", "pecan", name, "hmm");
+        Hmm *hmm = hmm_jsonParse(json, strlen(json));
+        if (hmm->type != fiveState) {
+            st_errAbort("pecan pair model %s is a type %i HMM; bar's pecan uses the five-state model (type %i)",
+                        name, (int)hmm->type, (int)fiveState);
+        }
+        pairModels->stateMachines[m] = hmm_getStateMachine(hmm);
+        hmm_destruct(hmm);
+        free(json);
+        free(name);
+    }
+    pairModels->distances = barDistances_constructFromCactusParams(params);
+    st_logInfo("bar: %" PRIi64 " pecan pair models (%" PRIi64 " for pairs of genomes), distances on the %s tree\n",
+               modelNo, keyed, barDistances_hasReconstructionTree(pairModels->distances) ? "reconstruction" : "event");
+    return pairModels;
+}
+
+void pecanPairModels_destruct(PecanPairModels *pairModels) {
+    if (pairModels == NULL) {
+        return;
+    }
+    for (int64_t m = 0; m < pairModels->modelNo; m++) {
+        stateMachine_destruct(pairModels->stateMachines[m]);
+    }
+    for (int64_t i = 0; i < 2 * pairModels->modelNo; i++) {
+        free(pairModels->genomes[i]);
+    }
+    free(pairModels->genomes);
+    free(pairModels->stateMachines);
+    free(pairModels->maxDistance);
+    barDistances_destruct(pairModels->distances);
+    free(pairModels);
+}
+
+/*
+ * cPecan's per-pair callback: the model for the two sequences' genomes, else for the distance between them.
+ */
+typedef struct {
+    PecanPairModels *pairModels;
+    stList *events; // the genome of each sequence, by its index in the alignment
+} PairModelChoice;
+
+static StateMachine *choosePairModel(int64_t seqX, int64_t seqY, void *extraArgs) {
+    PairModelChoice *choice = extraArgs;
+    Event *eventX = stList_get(choice->events, seqX), *eventY = stList_get(choice->events, seqY);
+    if (eventX == eventY) {
+        return NULL;
+    }
+    PecanPairModels *pairModels = choice->pairModels;
+    for (int64_t m = 0; m < pairModels->modelNo; m++) {
+        if (barModel_isForGenomes(pairModels->genomes + 2 * m, event_getHeader(eventX), event_getHeader(eventY))) {
+            return pairModels->stateMachines[m];
+        }
+    }
+    double distance = barDistances_get(pairModels->distances, eventX, eventY);
+    for (int64_t m = 0; m < pairModels->modelNo; m++) {
+        if (distance <= pairModels->maxDistance[m]) { // false for a NAN maxDistance
+            return pairModels->stateMachines[m];
+        }
+    }
+    return NULL;
+}
+
 stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees, int64_t maxSequenceLength,
         bool useProgressiveMerging, float gapGamma,
         PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
+    return makeEndAlignmentWithPairModels(sM, NULL, end, spanningTrees, maxSequenceLength, useProgressiveMerging, gapGamma,
+                                          pairwiseAlignmentBandingParameters);
+}
+
+stSortedSet *makeEndAlignmentWithPairModels(StateMachine *sM, PecanPairModels *pairModels, End *end, int64_t spanningTrees,
+        int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma,
+        PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
     //Make an alignment of the sequences in the ends
 
     //Get the adjacency sequences to be aligned.
     Cap *cap;
     End_InstanceIterator *it = end_getInstanceIterator(end);
     stList *sequences = stList_construct3(0, (void (*)(void *))adjacencySequence_destruct);
     stList *seqFrags = stList_construct3(0, (void (*)(void *))seqFrag_destruct);
+    stList *events = stList_construct();
     stHash *endInstanceNumbers = stHash_construct2(NULL, free);
     while((cap = end_getNext(it)) != NULL) {
         if(cap_getSide(cap)) {
@@ -72,6 +266,7 @@ stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees,
         assert(cap_getAdjacency(cap) != NULL);
         End *otherEnd = end_getPositiveOrientation(cap_getEnd(cap_getAdjacency(cap)));
         stList_append(seqFrags, seqFrag_construct(adjacencySequence->string, 0, end_getName(otherEnd)));
+        stList_append(events, cap_getEvent(cap));
         //Increase count of seqfrags with a given end.
         int64_t *c = stHash_search(endInstanceNumbers, otherEnd);
         if(c == NULL) {
@@ -84,7 +279,10 @@ stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees,
     end_destructInstanceIterator(it);
 
     //Get the alignment.
-    MultipleAlignment *mA  = makeAlignment(sM, seqFrags, spanningTrees, 100000000, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters);
+    PairModelChoice choice = { pairModels, events };
+    MultipleAlignment *mA  = makeAlignmentWithPairStateMachines(sM, pairModels == NULL ? NULL : choosePairModel, &choice,
+            seqFrags, spanningTrees, 100000000, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters);
+    stList_destruct(events);
 
     //Build an array of weights to reweight pairs in the alignment.
     int64_t *pairwiseAlignmentsPerSequenceNonCommonEnds = st_calloc(stList_length(seqFrags), sizeof(int64_t));

bar/impl/flowerAligner.c

@@ -528,7 +528,7 @@ static stSortedSet *getEndsToAlign(Flower *flower, int64_t maxSequenceLength) {
  * then call the makeFlowerAlignment2 consistency generating function.
  */
 
-static void computeMissingEndAlignments(StateMachine *sM, Flower *flower, stHash *endAlignments, int64_t spanningTrees,
+static void computeMissingEndAlignments(StateMachine *sM, PecanPairModels *pairModels, Flower *flower, stHash *endAlignments, int64_t spanningTrees,
         int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma,
         PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) {
     /*
@@ -546,7 +546,7 @@ static void computeMissingEndAlignments(StateMachine *sM, Flower *flower, stHash
                 stHash_insert(
                         endAlignments,
                         end,
-                        makeEndAlignment(sM, end, spanningTrees, maxSequenceLength,
+                        makeEndAlignmentWithPairModels(sM, pairModels, end, spanningTrees, maxSequenceLength,
                                 useProgressiveMerging, gapGamma,
                                 pairwiseAlignmentBandingParameters));
             } else {
@@ -562,7 +562,7 @@ stSortedSet *makeFlowerAlignment(StateMachine *sM, Flower *flower, int64_t spann
         bool useProgressiveMerging, float gapGamma,
         PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments) {
     stHash *endAlignments = stHash_construct2(NULL, (void(*)(void *)) stSortedSet_destruct);
-    computeMissingEndAlignments(sM, flower, endAlignments, spanningTrees, maxSequenceLength,
+    computeMissingEndAlignments(sM, NULL, flower, endAlignments, spanningTrees, maxSequenceLength,
             useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters);
     return makeFlowerAlignment2(flower, endAlignments, pruneOutStubAlignments);
 }
@@ -583,14 +583,14 @@ static void loadEndAlignments(Flower *flower, stHash *endAlignments, stList *lis
     }
 }
 
-stSortedSet *makeFlowerAlignment3(StateMachine *sM, Flower *flower, stList *listOfEndAlignmentFiles, int64_t spanningTrees,
+stSortedSet *makeFlowerAlignment3(StateMachine *sM, PecanPairModels *pairModels, Flower *flower, stList *listOfEndAlignmentFiles, int64_t spanningTrees,
         int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma,
         PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments) {
     stHash *endAlignments = stHash_construct2(NULL, (void(*)(void *)) stSortedSet_destruct);
     if(listOfEndAlignmentFiles != NULL) {
         loadEndAlignments(flower, endAlignments, listOfEndAlignmentFiles);
     }
-    computeMissingEndAlignments(sM, flower, endAlignments, spanningTrees, maxSequenceLength,
+    computeMissingEndAlignments(sM, pairModels, flower, endAlignments, spanningTrees, maxSequenceLength,
             useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters);
     return makeFlowerAlignment2(flower, endAlignments, pruneOutStubAlignments);
 }