Speed up bar's base aligners: pecan, minipoa and abPOA

.gitmodules

@@ -21,7 +21,7 @@
 	url = https://github.com/ComparativeGenomicsToolkit/kyoto.git
 [submodule "submodules/abPOA"]
 	path = submodules/abPOA
-	url = https://github.com/yangao07/abPOA.git
+	url = https://github.com/ComparativeGenomicsToolkit/abPOA.git
 [submodule "submodules/lastz"]
 	path = submodules/lastz
 	url = https://github.com/ComparativeGenomicsToolkit/lastz.git

Makefile

@@ -221,6 +221,14 @@ evolver_test_minipoa_local: all ${CWD}/test/primates-truth.maf
 evolver_test_minipoa_mammals_local: all ${CWD}/test/mammals-truth.maf
 	PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverMinipoaMammalsLocal
 
+# pecan, same dataset and tolerance as evolver_test_poa_local, so the two are comparable
+evolver_test_pecan_local: all ${CWD}/test/primates-truth.maf
+	PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPecanLocal
+
+# the mammals head-to-head against evolver_test_local, which runs the same data through abpoa
+evolver_test_pecan_mammals_local: all ${CWD}/test/mammals-truth.maf
+	PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPecanMammalsLocal
+
 evolver_test_refmap_local: all ${CWD}/test/primates-truth.maf
 	PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverRefmapLocal
 
@@ -258,7 +266,7 @@ evolver_test_primates_pangenome_resume_local: all ${CWD}/test/primates-truth.maf
 evolver_test_primates_pangenome_extend_local: all ${CWD}/test/primates-truth.maf
 	PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeExtendLocal
 
-evolver_test_all_local: evolver_test_local evolver_test_prepare_toil evolver_test_decomposed_local evolver_test_prepare_no_outgroup_local evolver_test_poa_local evolver_test_minipoa_local evolver_test_minipoa_mammals_local evolver_test_refmap_local evolver_test_minigraph_local
+evolver_test_all_local: evolver_test_local evolver_test_prepare_toil evolver_test_decomposed_local evolver_test_prepare_no_outgroup_local evolver_test_poa_local evolver_test_minipoa_local evolver_test_minipoa_mammals_local evolver_test_pecan_local evolver_test_pecan_mammals_local evolver_test_refmap_local evolver_test_minigraph_local
 
 yeast_test_local:
 	PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testYeastPangenomeLocal

bar/impl/bar.c

@@ -49,6 +49,28 @@ PairwiseAlignmentParameters *pairwiseAlignmentParameters_constructFromCactusPara
     p->alignAmbiguityCharacters = cactusParams_get_int(params, 3, "bar", "pecan", "alignAmbiguityCharacters");
     p->useMumAnchors = cactusParams_get_int(params, 3, "bar", "pecan", "useMumAnchors");
     p->recursiveMums = cactusParams_get_int(params, 3, "bar", "pecan", "recursiveMums");
+    // anchorMethod supersedes useMumAnchors when it is given; without it, useMumAnchors decides as it always has
+    if (cactusParams_has(params, 3, "bar", "pecan", "anchorMethod")) {
+        char *method = cactusParams_get_string(params, 3, "bar", "pecan", "anchorMethod");
+        if (strcmp(method, "seed") == 0) {
+            p->anchorMethod = PAIRWISE_ANCHOR_SEED;
+        } else if (strcmp(method, "mum") == 0 || strcmp(method, "lastz") == 0) {
+            p->anchorMethod = PAIRWISE_ANCHOR_LEGACY;
+            p->useMumAnchors = strcmp(method, "mum") == 0;
+        } else {
+            st_errAbort("Unknown <pecan anchorMethod=\"%s\">; expected seed, mum or lastz", method);
+        }
+        free(method);
+    }
+    if (cactusParams_has(params, 3, "bar", "pecan", "seedHspThreshold")) {
+        p->seedHspThreshold = cactusParams_get_int(params, 3, "bar", "pecan", "seedHspThreshold");
+    }
+    if (cactusParams_has(params, 3, "bar", "pecan", "seedHspThresholdMin")) {
+        p->seedHspThresholdMin = cactusParams_get_int(params, 3, "bar", "pecan", "seedHspThresholdMin");
+    }
+    if (cactusParams_has(params, 3, "bar", "pecan", "seedRecursionDepth")) {
+        p->seedRecursionDepth = cactusParams_get_int(params, 3, "bar", "pecan", "seedRecursionDepth");
+    }
     return p;
 }
 
@@ -110,7 +132,9 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList *
         st_logInfo("bar: base aligner %s, window %" PRIi64 ", maskFilter %" PRIi64 ", bandingLimit %" PRIi64 "\n",
                    baseAligner_toString(engine), poaWindow, maskFilter, maximumLength);
     } else {
-        st_logInfo("bar: base aligner pecan, bandingLimit %" PRIi64 "\n", maximumLength);
+        st_logInfo("bar: base aligner pecan, anchors %s, bandingLimit %" PRIi64 "\n",
+                   pairwiseAlignmentParameters->anchorMethod == PAIRWISE_ANCHOR_SEED ? "seed" :
+                   (pairwiseAlignmentParameters->useMumAnchors ? "mum" : "lastz"), maximumLength);
     }
 
     //////////////////////////////////////////////

bar/impl/flowerAligner.c

@@ -5,6 +5,7 @@
  */
 
 #include "endAligner.h"
+#include "flowerAligner.h" // so the compiler holds the definitions here to the declarations
 #include "cactus.h"
 #include "sonLib.h"
 #include "adjacencySequences.h"
@@ -505,7 +506,11 @@ static stSortedSet *getEndsToAlign(Flower *flower, int64_t maxSequenceLength) {
      */
     stSortedSet *endsToAlign = stSortedSet_construct();
     End *dominantEnd = getDominantEnd(flower); //an end to which all adjacencies are incident
-    if (dominantEnd != NULL && getMaxAdjacencyLength(flower) <= 2 * maxSequenceLength) {
+    // Aligning only the dominant end is a saving only if it loses nothing.  Its alignment reaches
+    // maxSequenceLength into each adjacency, from its own side, where aligning every end reaches
+    // that far from both sides -- so an adjacency longer than maxSequenceLength (this used to allow
+    // twice that) would have its far part left unaligned.  The poa path makes the same test.
+    if (dominantEnd != NULL && getMaxAdjacencyLength(flower) <= maxSequenceLength) {
         stSortedSet_insert(endsToAlign, dominantEnd);
     } else {
         End *end;
@@ -555,7 +560,7 @@ static void computeMissingEndAlignments(StateMachine *sM, Flower *flower, stHash
 
 stSortedSet *makeFlowerAlignment(StateMachine *sM, Flower *flower, int64_t spanningTrees, int64_t maxSequenceLength,
         bool useProgressiveMerging, float gapGamma,
-        PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments, int64_t poaWindow) {
+        PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments) {
     stHash *endAlignments = stHash_construct2(NULL, (void(*)(void *)) stSortedSet_destruct);
     computeMissingEndAlignments(sM, flower, endAlignments, spanningTrees, maxSequenceLength,
             useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters);

bar/impl/poaBarAligner.c

@@ -661,15 +661,14 @@ static void run_abpoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_param
  * last_scoring.py's learned matrix reaches both (it writes into <poa>, which local_alignment.py
  * also reads to score FastGA PAFs).
  *
- * The gaps have to be separate.  abPOA's gap model is convex -- min(open1 + L*ext1, open2 +
- * L*ext2) -- so with the shipped 400/30 and 1200/1 its effective extension beyond L~28 is 1.
- * minipoa has a single affine piece, so handing it abPOA's first-piece extension of 30 prices
- * long gaps about 30x above what abPOA charges, and it responds by packing bases into shared
- * columns instead of opening a gap.  That cost about six points of mafComparator accuracy on
- * evolver mammals.
+ * The gaps are separate so each engine can be priced for itself.  Both can be convex --
+ * min(open1 + L*ext1, open2 + L*ext2) -- and <minipoa> ships abPOA's 400/30 and 1200/1.  With
+ * the second piece off, minipoa has a single affine piece, and abPOA's first piece alone prices
+ * long gaps about 30x above what abPOA charges: it packs bases into shared columns instead of
+ * opening a gap, which cost about six points of mafComparator accuracy on evolver mammals.
  *
- * last-train reaches minipoa too: it fits a single affine model, which is minipoa's model
- * exactly, so last_scoring.py writes the learned open/extend straight into <minipoa>.  The
+ * last-train reaches minipoa too: it fits a single affine model, so last_scoring.py writes the
+ * learned open/extend into <minipoa>'s first piece and turns the second off.  The
  * GapOpen2/GapExtend2 pair it synthesises for abPOA is a stability workaround for that aligner
  * and is deliberately not passed on.
  */
@@ -704,14 +703,16 @@ static minipoa_para_t *minipoaParameters_constructFromCactusParams(CactusParams
     }
     free(submat_string);
 
-    /*
-     * Only the first gap piece.  abPOA takes min(open1 + L*ext1, open2 + L*ext2) and minipoa has
-     * no second piece at all, so gaps longer than where the two cross (~28bp with the shipped
-     * 400/30 and 1200/1) are penalised more heavily here than abPOA would.
-     */
     out->gapOpen = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaGapOpenPenalty");
     out->gapExt = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaGapExtensionPenalty");
     minipoa_set_gap(mpt, out->gapOpen, out->gapExt);
+    // The optional second piece, which makes the gap model convex, as abPOA's is.  Absent or 0 is
+    // the single piece above.
+    out->gapOpen2 = cactusParams_has(params, 3, "bar", "minipoa", "minipoaGapOpenPenalty2") ?
+        cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaGapOpenPenalty2") : 0;
+    out->gapExt2 = cactusParams_has(params, 3, "bar", "minipoa", "minipoaGapExtensionPenalty2") ?
+        cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaGapExtensionPenalty2") : 0;
+    minipoa_set_gap2(mpt, out->gapOpen2, out->gapExt2);
 
     out->bandConstant = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaBandConstant");
     out->bandFraction = cactusParams_get_float(params, 3, "bar", "minipoa", "minipoaBandFraction");
@@ -726,10 +727,9 @@ static minipoa_para_t *minipoaParameters_constructFromCactusParams(CactusParams
     minipoa_set_seeding(mpt, out->seeding, out->minimizerK, out->minimizerW, out->anchorWindow);
 
     /*
-     * Progressive ordering, on by default, because abPOA runs with it on
-     * (<poa partialOrderAlignmentProgressiveMode>) and the order sequences are added to a POA
-     * graph changes the alignment -- markedly so on diverged input.  Leaving it off here was worth
-     * about six points of mafComparator accuracy on the evolver mammals set.
+     * Progressive ordering, as abPOA does it (<poa partialOrderAlignmentProgressiveMode>): the
+     * order sequences are added to a POA graph changes the alignment, markedly so on diverged
+     * input, and adding the closest pair first is what a guide tree buys.
      */
     out->progressive = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaProgressiveMode");
     out->progressiveMaxRows = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaProgressiveMaxRows");
@@ -762,6 +762,11 @@ static char *dump_minipoa_input(Msa *msa, PoaParameters *pp, uint8_t **bseqs, ch
     char *command = st_malloc(4096 * sizeof(char));
     sprintf(command, "minipoa %s -m %s -O -%d -E -%d -b %d -f %d -r 1 -t 1",
             input_path, matrix_path, pp->gapOpen, pp->gapExt, pp->bandConstant, f);
+    if (pp->gapOpen2 != 0 || pp->gapExt2 != 0) {
+        char gap2_opts[128];
+        sprintf(gap2_opts, " --gap_open2 -%d --gap_ext2 -%d", pp->gapOpen2, pp->gapExt2);
+        strcat(command, gap2_opts);
+    }
     if (pp->seeding) {
         char kw_opts[128];
         sprintf(kw_opts, " -S -k %d -w %d", pp->minimizerK, pp->minimizerW);

bar/inc/poaBarAligner.h

@@ -78,6 +78,7 @@ typedef struct _PoaParameters {
      */
     int mat[25];
     int gapOpen, gapExt;
+    int gapOpen2, gapExt2; /* minipoa's optional convex second piece; 0 = off */
     int bandConstant;
     double bandFraction;
     bool seeding;

doc/progressive.md

@@ -427,7 +427,7 @@ These are the most relevant options for running on a cluster
 * `--consMemory`: Override the memory for each `cactus_consolidated` job. Can be useful if Cactus's estimates are wrong, but a job that runs out of memory is already retried with double the memory by default, so this is rarely needed.
 * `--consRetainPages`: Whether `cactus_consolidated` keeps the memory pages it frees rather than returning them to the system (`auto`, `1` or `0`). Keeping them is much faster but takes 2-3x the peak memory. The default, `auto`, keeps them when twice the memory estimate, plus one further estimate of headroom, fits what the job can be given (the system memory on a single machine, or `--maxMemory`), and then requests twice the estimate; otherwise the pages are returned and the request is the estimate itself. Where that ceiling is not known, the pages are returned. Corresponds to `<consolidated retain_pages>` in the configuration.
 * **Large, repeat-rich genomes**: turn on `partialOrderAlignmentMaskFilter` in the `<bar><poa>` section of the config (e.g. `1000`; the default `-1` is off). On 20-30 Gb salamander genomes, leaving it off cost ~38x the BAR time and ran out of memory on a 2 TB node; with it on, BAR took under an hour. `partialOrderAlignmentWindow` is a far cheaper lever: halving it to `5000` saved ~25% memory and ~30% BAR time.
-* **Which base aligner**: `<bar baseAligner>` selects it -- `abpoa` (the default), `minipoa`, or the deprecated `pecan`. The two levers above live in `<bar><poa>` and apply to **abpoa only**; under `minipoa` the equivalents are `minipoaMaskFilter` and `minipoaWindow` in `<bar><minipoa>` (and `--barMaskFilter` writes to whichever engine is selected). minipoa uses roughly half abpoa's memory at the *same* window, and ships the same `10000`: its memory is quadratic-ish in the window as well, and on human/chimp chr10 a 50 kb window is already OOM-killed while BAR time is flat from 10 kb to 30 kb. Raising `minipoaWindow` buys nothing measured and costs a lot; if you do raise it, note the `cactus_consolidated` memory estimate uses abpoa's fitted exponent, which under-requests badly for minipoa at larger windows. minipoa shares abpoa's substitution matrix (from `<bar><poa>`) but has its own gap penalties, `minipoaGapOpenPenalty`/`minipoaGapExtensionPenalty` in `<bar><minipoa>`. They are deliberately much lower than abpoa's: abpoa's gap model is convex, so past roughly 28bp its effective extension is 1 rather than 30, and giving minipoa's single affine piece abpoa's first-piece value of 30 prices long gaps about 30x too high. `--lastTrain` reaches both engines: it fits a single affine gap model, which is minipoa's model exactly, so minipoa gets the learned open/extend verbatim while abPOA additionally gets the synthesised second piece it needs to stay stable. minipoa inherits the learned substitution matrix too, via the empty `minipoaSubMatrix`.
+* **Which base aligner**: `<bar baseAligner>` selects it -- `abpoa` (the default), `minipoa`, or the deprecated `pecan`. The two levers above live in `<bar><poa>` and apply to **abpoa only**; under `minipoa` the equivalents are `minipoaMaskFilter` and `minipoaWindow` in `<bar><minipoa>` (and `--barMaskFilter` writes to whichever engine is selected). minipoa uses roughly half abpoa's memory at the *same* window, and ships the same `10000`: its memory is quadratic-ish in the window as well, and on human/chimp chr10 a 50 kb window is already OOM-killed while BAR time is flat from 10 kb to 30 kb. Raising `minipoaWindow` buys nothing measured and costs a lot; if you do raise it, note the `cactus_consolidated` memory estimate uses abpoa's fitted exponent, which under-requests badly for minipoa at larger windows. minipoa shares abpoa's substitution matrix (from `<bar><poa>`) but has its own gap penalties in `<bar><minipoa>`: `minipoaGapOpenPenalty`/`minipoaGapExtensionPenalty`, and a second piece, `minipoaGapOpenPenalty2`/`minipoaGapExtensionPenalty2`, that makes the gap cost convex exactly as abpoa's is. They ship with abpoa's values (400/30 and 1200/1), and together with the guide tree (`minipoaProgressiveMode`) that makes minipoa as accurate as abpoa on the evolver tests (mammals 0.854/0.652 against 0.853/0.652). Setting the second piece to 0 leaves a single affine piece, about a quarter cheaper in time and memory; it then needs a much lower extension than abpoa's first piece (the best found was 600/4), since past roughly 28bp abpoa's effective extension is 1 rather than 30. `--lastTrain` reaches both engines: it fits a single affine gap model, so minipoa gets the learned open/extend as its only piece, while abPOA additionally gets the synthesised second piece it needs to stay stable. minipoa inherits the learned substitution matrix too, via the empty `minipoaSubMatrix`.
 
 On a cluster with partitions and/or time limits, make sure to use