From 5cac8f8c41a29f070e20724ae53f28969a02bbf6 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 27 Sep 2026 19:53:30 -0400 Subject: [PATCH 1/5] Make pecan a practical base aligner for progressive bar Bumps cPecan to faster-pecan (12ce978): an exact, probability-space forward-backward 5-11x faster per pairwise alignment that recomputes large unanchored chunks in blocks instead of holding them (3.9GB -> 116MB for a 10kb pair), an array-based progressive MSA (same output, 2x faster, half the memory), and in-process seed anchoring. The anchoring is what decides pecan's accuracy on deep branches. The MUM anchors bar used accept "unique" matches of a dozen bases, and at 0.3 substitutions per site per branch half of them are off the true alignment and drag the banded dp with them. is now shipped: lastz's recipe (spaced seeds, ungapped extension, chaining) without a lastz process per alignment. anchorMethod="mum" or "lastz" select the old ones; without the attribute useMumAnchors decides as before. Evolver mammals, mafComparator: abpoa 0.853/0.652, pecan before 0.859/0.654, pecan now 0.869/0.666, with bar taking 1-3s per ancestor against abpoa's 0-3s (stock pecan 1-10s). Evolver primates: abpoa 0.9973/0.9844, pecan 0.9977/0.9846. 8Mb of human/mouse: bar 26s/0.5GB against abpoa's 59s/6.7GB, aligning slightly more of human to mouse. On evolver-like inputs pecan still spends about twice abpoa's CPU in bar, so baseAligner stays abpoa. use deepAligner for ancestors whose two farthest children are at least that far apart in the scaled tree the workflow already uses for its divergence settings. Off by default (0); with branchScale 2, 0.5 picks the deepest 15% of the VGP 577-way's ancestors. Fixes in bar: aligning only the dominant end is now taken only when no adjacency is longer than bandingLimit -- it allowed twice that, leaving the far part of such adjacencies unaligned; flowerAligner.c includes its own header, which shows makeFlowerAlignment's definition had a parameter its declaration lacks. The bar log line names the anchor method. New evolver tests run pecan on primates and mammals. Co-Authored-By: Claude Opus 5.5 --- Makefile | 10 ++++++- bar/impl/bar.c | 26 ++++++++++++++++- bar/impl/flowerAligner.c | 9 ++++-- src/cactus/cactus_progressive_config.xml | 34 ++++++++++++++++++++-- src/cactus/pipeline/cactus_workflow.py | 37 ++++++++++++++++++++---- submodules/cPecan | 2 +- test/evolverTest.py | 19 ++++++++++++ 7 files changed, 125 insertions(+), 12 deletions(-) diff --git a/Makefile b/Makefile index da6cb8126..6cdaecdcf 100644 --- a/Makefile +++ b/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 diff --git a/bar/impl/bar.c b/bar/impl/bar.c index e342ea9f0..59a54e268 100644 --- a/bar/impl/bar.c +++ b/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 ; 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); } ////////////////////////////////////////////// diff --git a/bar/impl/flowerAligner.c b/bar/impl/flowerAligner.c index f1221c282..390d9bc3e 100644 --- a/bar/impl/flowerAligner.c +++ b/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); diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index f765540b7..966feb7c0 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -249,13 +249,25 @@ + @@ -264,6 +276,8 @@ runBar="1" bandingLimit="1000000" baseAligner="abpoa" + deepAligner="pecan" + deepAlignerDivergence="0" minimumBlockDegree="2" minimumIngroupDegree="1" minimumOutgroupDegree="0" @@ -286,7 +300,21 @@ - + + + diff --git a/src/cactus/pipeline/cactus_workflow.py b/src/cactus/pipeline/cactus_workflow.py index 6492ad058..8fa71ee53 100644 --- a/src/cactus/pipeline/cactus_workflow.py +++ b/src/cactus/pipeline/cactus_workflow.py @@ -46,6 +46,31 @@ def cons_core_scale(cores, baseline=CONS_CORE_BASELINE, parallel=CONS_PARALLEL_F return (1.0 - parallel) + parallel * (float(baseline) / cores) return 1.0 +def ingroup_divergence(tree, ancestor_event): + """ How far apart the genomes bar aligns into an ancestor are: the distance in the tree between its + two farthest-apart children, i.e. its two longest child branches. The tree is the one the workflow + hands consolidated, which is scaled by --branchScale and, with upweightAncestorDistances, has each + ancestral child's height added to its branch -- so this is measured as cactus's other + divergence-dependent settings are. """ + node = tree.getNodeId(ancestor_event) + lengths = sorted([tree.getWeight(node, child) for child in tree.getChildren(node)], reverse=True) + return sum(lengths[:2]) if len(lengths) >= 2 else 0.0 + +def bar_base_aligner(tree, ancestor_event, config_node, name=None): + """ The base aligner bar uses for this ancestor: when its ingroups are at least + apart, else . """ + bar_node = findRequiredNode(config_node, 'bar') + base_aligner = getOptionalAttrib(bar_node, 'baseAligner', typeFn=str, default='abpoa') + deep_aligner = getOptionalAttrib(bar_node, 'deepAligner', typeFn=str, default='') + deep_divergence = getOptionalAttrib(bar_node, 'deepAlignerDivergence', typeFn=float, default=0.0) + if deep_aligner and deep_aligner != base_aligner and deep_divergence > 0: + divergence = ingroup_divergence(tree, ancestor_event) + if divergence >= deep_divergence: + RealtimeLogger.info('cactus_consolidated({}): ingroups {:g} apart, at or above deepAlignerDivergence {:g}, so bar uses {} rather than {}'.format( + name if name else ancestor_event, divergence, deep_divergence, deep_aligner, base_aligner)) + return deep_aligner + return base_aligner + def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_map, og_map, paf_id, cons_cores = None, cons_memory = None, intermediate_results_url = None, chrom_name = None, cons_retain_pages = None): @@ -92,8 +117,7 @@ def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_ma ramp = min(1.0, input_gb / 0.02) estimate = int(mem_coef * (input_gb ** mem_exp) * ramp * 2**30) if input_gb > 0 else 0 # A POA aligner (abPOA or minipoa) needs a table even for tiny alignments - bar_node = findRequiredNode(config_node, 'bar') - base_aligner = getOptionalAttrib(bar_node, 'baseAligner', typeFn=str, default='abpoa') + base_aligner = bar_base_aligner(tree, ancestor_event, config_node, name) if base_aligner != 'pecan': estimate = max(estimate, int(4e9)) @@ -224,17 +248,17 @@ def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_ma cons_job = job.addChildJobFn(cactus_cons, tree, ancestor_event, config_node, seq_id_map, og_map, paf_id, intermediate_results_url=intermediate_results_url, chrom_name=chrom_name, cores = cons_cores, memory=cactus_clamp_memory(mem), disk=disk, retain_pages=retain_pages, - poa_window=poa_window, walltime=cactus_walltime(walltime_secs)) + poa_window=poa_window, base_aligner=base_aligner, walltime=cactus_walltime(walltime_secs)) return cons_job.rv() def cactus_cons(job, tree, ancestor_event, config_node, seq_id_map, og_map, paf_id, intermediate_results_url = None, chrom_name = None, retain_pages = None, - poa_window = None): + poa_window = None, base_aligner = None): ''' run cactus_consolidated ''' # cactus_consolidated reads its settings from the config, so the resolved page retention # goes into the copy it is given (this job's copy of the node, so nothing else sees it) - if retain_pages is not None or poa_window is not None: + if retain_pages is not None or poa_window is not None or base_aligner is not None: config_node = copy.deepcopy(config_node) if retain_pages is not None: findRequiredNode(config_node, 'consolidated').set('retain_pages', str(retain_pages)) @@ -243,6 +267,9 @@ def cactus_cons(job, tree, ancestor_event, config_node, seq_id_map, og_map, paf_ poa_node = findRequiredNode(config_node, 'bar').find('poa') if poa_window is not None and poa_node is not None: poa_node.set('partialOrderAlignmentWindow', str(poa_window)) + # likewise the base aligner, if chose it for this ancestor + if base_aligner is not None: + findRequiredNode(config_node, 'bar').set('baseAligner', base_aligner) # Build up a genome -> fasta map. work_dir = job.fileStore.getLocalTempDir() diff --git a/submodules/cPecan b/submodules/cPecan index f7e7d6320..12ce97891 160000 --- a/submodules/cPecan +++ b/submodules/cPecan @@ -1 +1 @@ -Subproject commit f7e7d6320878b4dbcfaa16cbf49d03c2d8a038ac +Subproject commit 12ce97891423c2b80a148d2a57452d4d44e0fc77 diff --git a/test/evolverTest.py b/test/evolverTest.py index f339896ca..0f15ee9b9 100644 --- a/test/evolverTest.py +++ b/test/evolverTest.py @@ -1913,6 +1913,25 @@ def testEvolverMinipoaMammalsLocal(self): chromInfoDict = {'simChow' : 'X,Y', 'simDog' : 'X', 'simRat' : 'Y', 'simHuman' : 'X,Y,Z'}) self._check_maf_accuracy(self._out_hal(name), delta=(0.05,0.13)) + def testEvolverPecanLocal(self): + """ testEvolverPOALocal with baseAligner="pecan": the same data and tolerance, so the two + engines are directly comparable. + """ + pecan_config_path = self._write_base_aligner_config("pecan", star=True) + name = "local" + self._run_evolver_primates_star(name, configFile = pecan_config_path) + self._check_maf_accuracy(self._out_hal("local"), delta=(0.0025,0.0075), dataset='primates') + + def testEvolverPecanMammalsLocal(self): + """ The mammals head-to-head with baseAligner="pecan", as testEvolverMinipoaMammalsLocal is + for minipoa. + """ + pecan_config_path = self._write_base_aligner_config("pecan") + name = "local" + self._run_evolver(name, configFile = pecan_config_path, + chromInfoDict = {'simChow' : 'X,Y', 'simDog' : 'X', 'simRat' : 'Y', 'simHuman' : 'X,Y,Z'}) + self._check_maf_accuracy(self._out_hal(name), delta=(0.05,0.13)) + def testEvolverRefmapLocal(self): """ Use the new minimap pangenome pipeline to create an alignment of the primates, then compare with the baseline """ From c28e345acd56483606d94dcd3fa340ba06225b9b Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Mon, 28 Sep 2026 07:24:00 -0400 Subject: [PATCH 2/5] Make minipoa as accurate as abpoa in bar, and several times faster Bumps minipoa to 0416daa (ai-optimization), which vectorises and streams its DP, keeps its scratch between windows, centres its band better, adds an optional convex gap, and fixes its guide tree below 8 sequences. now ships the convex gap abpoa uses (400/30 and 1200/1, via the new minipoaGapOpenPenalty2/minipoaGapExtensionPenalty2) and progressive mode on. Evolver mammals: 0.846/0.646 -> 0.854/0.652, against abpoa's 0.853/0.652; primates unchanged. On 8Mb of human/mouse at 8 threads, bar takes 14s and 5.6GB against abpoa's 59s and 6.7GB (10s and 3.5GB with the second gap piece off, where the old minipoa took 29s and 3.1GB for identical output). last-train now turns the second piece off, so minipoa runs the trained single affine model as trained. Co-Authored-By: Claude Opus 5.5 --- bar/impl/poaBarAligner.c | 39 ++++++++------- bar/inc/poaBarAligner.h | 1 + doc/progressive.md | 2 +- src/cactus/cactus_progressive_config.xml | 64 +++++++++++++----------- src/cactus/paf/last_scoring.py | 14 ++++-- submodules/minipoa | 2 +- 6 files changed, 68 insertions(+), 54 deletions(-) diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index 68b91a91a..f936ed4ba 100644 --- a/bar/impl/poaBarAligner.c +++ b/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 , 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 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 . The + * last-train reaches minipoa too: it fits a single affine model, so last_scoring.py writes the + * learned open/extend into '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 - * () 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 (): 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); diff --git a/bar/inc/poaBarAligner.h b/bar/inc/poaBarAligner.h index 912a3d680..f28c59b22 100644 --- a/bar/inc/poaBarAligner.h +++ b/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; diff --git a/doc/progressive.md b/doc/progressive.md index 080bdeec0..3cf39b30e 100644 --- a/doc/progressive.md +++ b/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 `` in the configuration. * **Large, repeat-rich genomes**: turn on `partialOrderAlignmentMaskFilter` in the `` 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**: `` selects it -- `abpoa` (the default), `minipoa`, or the deprecated `pecan`. The two levers above live in `` and apply to **abpoa only**; under `minipoa` the equivalents are `minipoaMaskFilter` and `minipoaWindow` in `` (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 ``) but has its own gap penalties, `minipoaGapOpenPenalty`/`minipoaGapExtensionPenalty` in ``. 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**: `` selects it -- `abpoa` (the default), `minipoa`, or the deprecated `pecan`. The two levers above live in `` and apply to **abpoa only**; under `minipoa` the equivalents are `minipoaMaskFilter` and `minipoaWindow` in `` (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 ``) but has its own gap penalties in ``: `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 third 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 diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index 966feb7c0..0ef8ac866 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -255,12 +255,13 @@ 0.869/0.666 against abpoa's 0.853/0.652, evolver primates 0.9977/0.9846 against 0.9973/0.9844, and on 8Mb of human/mouse bar took 26s and 0.5GB against abpoa's 59s and 6.7GB while aligning slightly more of human to mouse. On evolver mammals it takes about twice abpoa's CPU time in bar. minipoa is a - minimizer-based POA that uses roughly half abpoa's memory and is somewhat faster at the same - window size. It does NOT afford a larger window: measured on human/chimp chr10, its memory - grows about as window x (bandConstant + bandFraction x window), so 10kb costs 3.5GB and 50kb - is already over 23GB and gets OOM-killed. minipoa shares abpoa's substitution matrix - (from ) but has its own gap penalties in , which - are deliberately different - see the note there for why. --> + POA like abpoa, with a static band and a vectorised dp that streams its matrices, and with the + shipped settings (the same convex gaps as abpoa, and a guide tree) it is as accurate: + evolver mammals 0.854/0.652 against abpoa's 0.853/0.652, primates equal. On 8Mb of human/mouse + at 8 threads its bar took 14s and 5.8GB against abpoa's 59s and 6.7GB (9s and 3.5GB with a + single affine gap). Its memory still grows about as window x (bandConstant + bandFraction x + window), so it does not afford a much larger window. minipoa shares abpoa's substitution matrix + (from ) but has its own gap penalties in . --> - + + + + On a tree-shaped window - two sister ingroups and outgroups, as progressive cactus aligns - adding + the closest sequences first raised recall of simulated truth by 1-2 points at moderate divergence + and 7-13 at high divergence, and was worth 0.2/0.2 on evolver mammals on top of the convex gap. + (Earlier measurements that found it a no-op predate a fix: below 8 sequences the tree was + built from no minimizers at all and came out as the input order.) --> diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index 534db52fa..f23e6de5a 100644 --- a/src/cactus/paf/last_scoring.py +++ b/src/cactus/paf/last_scoring.py @@ -133,15 +133,19 @@ def apply_scores_to_config(score_dict, config_xml): poa_node.attrib['partialOrderAlignmentSubMatrix'] = score_string - # minipoa gets the learned gaps verbatim. last-train fits a single affine gap model, which is - # exactly minipoa's model -- the GAP-OPEN-2/GAP-EXTEND-2 pair above is a synthesised second - # piece that exists only because abPOA is unstable without one, so it would be wrong to hand it - # on. minipoaSubMatrix is left alone: empty means inherit 's matrix, which is the learned - # one we just wrote, so the trained substitution scores reach minipoa too. + # minipoa gets the learned gaps verbatim, as a single affine piece: last-train fits a single + # affine gap model, and the GAP-OPEN-2/GAP-EXTEND-2 pair above is a synthesised second piece + # that exists only because abPOA is unstable without one, so it would be wrong to hand it on. + # minipoa's own second piece is turned off too: its shipped values are priced against the + # shipped first piece, not the learned one. minipoaSubMatrix is left alone: empty + # means inherit 's matrix, which is the learned one we just wrote, so the trained + # substitution scores reach minipoa too. minipoa_node = bar_node.find("minipoa") if minipoa_node is not None: minipoa_node.attrib['minipoaGapOpenPenalty'] = str(score_dict['GAP-OPEN']) minipoa_node.attrib['minipoaGapExtensionPenalty'] = str(score_dict['GAP-EXTEND']) + minipoa_node.attrib['minipoaGapOpenPenalty2'] = '0' + minipoa_node.attrib['minipoaGapExtensionPenalty2'] = '0' RealtimeLogger.info("Overriding minipoa scores with trained values: GapOpen {}; GapExtend {} (single affine, as trained)".format( minipoa_node.attrib['minipoaGapOpenPenalty'], minipoa_node.attrib['minipoaGapExtensionPenalty'])) diff --git a/submodules/minipoa b/submodules/minipoa index fd134725e..0416daa63 160000 --- a/submodules/minipoa +++ b/submodules/minipoa @@ -1 +1 @@ -Subproject commit fd134725ee16acb5de152f096a0a84e7664d51cc +Subproject commit 0416daa63e69021b008b37ab19f2e29ae6786225 From 6d15d7c6caed762cc09fe86f7028ea9e6b9f846d Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Mon, 28 Sep 2026 07:57:53 -0400 Subject: [PATCH 3/5] Bump minipoa to keep its horizontal gap states as bits minipoa's DP now keeps two bits a cell for each horizontal gap state instead of a 16-bit value, with alignments unchanged. On 8Mb of human/mouse at 8 threads, bar with the shipped settings goes from 14s and 5.6GB to 12s and 4GB (abpoa: 59s and 6.7GB), and with a single affine gap from 10s and 3.5GB to 9s and 2.8GB. Co-Authored-By: Claude Opus 5.5 --- doc/progressive.md | 2 +- src/cactus/cactus_progressive_config.xml | 6 +++--- submodules/minipoa | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/progressive.md b/doc/progressive.md index 3cf39b30e..ef8a4e8d5 100644 --- a/doc/progressive.md +++ b/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 `` in the configuration. * **Large, repeat-rich genomes**: turn on `partialOrderAlignmentMaskFilter` in the `` 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**: `` selects it -- `abpoa` (the default), `minipoa`, or the deprecated `pecan`. The two levers above live in `` and apply to **abpoa only**; under `minipoa` the equivalents are `minipoaMaskFilter` and `minipoaWindow` in `` (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 ``) but has its own gap penalties in ``: `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 third 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`. +* **Which base aligner**: `` selects it -- `abpoa` (the default), `minipoa`, or the deprecated `pecan`. The two levers above live in `` and apply to **abpoa only**; under `minipoa` the equivalents are `minipoaMaskFilter` and `minipoaWindow` in `` (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 ``) but has its own gap penalties in ``: `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 diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index 0ef8ac866..9473c1177 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -258,7 +258,7 @@ POA like abpoa, with a static band and a vectorised dp that streams its matrices, and with the shipped settings (the same convex gaps as abpoa, and a guide tree) it is as accurate: evolver mammals 0.854/0.652 against abpoa's 0.853/0.652, primates equal. On 8Mb of human/mouse - at 8 threads its bar took 14s and 5.8GB against abpoa's 59s and 6.7GB (9s and 3.5GB with a + at 8 threads its bar took 12s and 4GB against abpoa's 59s and 6.7GB (9s and 2.8GB with a single affine gap). Its memory still grows about as window x (bandConstant + bandFraction x window), so it does not afford a much larger window. minipoa shares abpoa's substitution matrix (from ) but has its own gap penalties in . --> @@ -397,8 +397,8 @@