From 1c082a38364d1a9218f28e0506a40de57c8f1cae Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 16 Sep 2026 07:57:49 -0400 Subject: [PATCH 1/4] bar: minipoa as an optional third base aligner abPOA's memory use is what forces BAR's 10kb window. minipoa is a minimizer-based POA that partitions the graph on its consensus path and needs far less of it, so it can take a much larger window. On evolver mammals the two are equivalent on accuracy -- 0.8506 against abPOA's 0.8513, and 0.9974 against 0.9973 on primates -- at 1.0 GB peak RSS against 2.3 GB, with a 100kb window against abPOA's 10kb. abPOA stays the default. Selected with . There is no command line flag; the config attribute is the only selector. Configs written before the attribute existed keep working: both it and the older boolean are optional, and when only the boolean is present it decides as it always did. When both are present and disagree, baseAligner wins and says so. Scoring is the part that needed care. minipoa shares abPOA's substitution matrix, from , so substitutions score identically and last-train's learned matrix reaches both. The gap penalties are deliberately separate, because the same nominal values do not produce the same alignment: 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 past L~28 is 1, not 30. minipoa has a single affine piece, and charging it abPOA's first-piece 30 prices long gaps about 30x above what abPOA charges. It responds by packing bases into shared columns instead of opening a gap -- 9% fewer columns and 151k spurious aligned pairs measured across 229 real BAR windows -- which cost about six points of mafComparator accuracy. That matters here because 94% of the gap bases in the evolver mammals truth sit in runs of 28bp or longer; on data with mostly short gaps, such as human pangenomes, the two models differ far less. default to 600/4, which are fitted to that test and worth revisiting against other data. --lastTrain replaces them and is a better fit for minipoa than for abPOA: it trains a single affine model, which is minipoa's model exactly, whereas abPOA has to be handed a synthesised second piece to stay stable. Build: minipoa is a submodule built into lib/libminipoa.a, gated on a "minipoa" variable in include.mk that defaults on for x86 and off for ARM, where its simde/NEON path is unproven. With it off the submodule is not built, -lminipoa is not linked, and selecting the engine aborts with a message naming the switch, so a tree that cannot build minipoa still builds cactus. -lminipoa sits next to -labpoa, ahead of -lz: libminipoa.a has undefined references into zlib and ld resolves archives left to right. Verified: 20 unit tests over both engines; abPOA byte-for-byte identical to master across 229 windows, so the refactor preserves it; SSE2 and AVX2 produce identical alignments; full static link; docker image builds and contains a working minipoa; evolver primates and mammals accuracy both pass. Two things here are independent of minipoa and could be taken separately: - CACTUS_BAR_DUMP_DIR dumps every window handed to the base aligner, with a command line that replays it. This was a commented-out #define that also deleted its own output on success and named files by a pointer that got reused, so it was unusable for looking at a window that aligned badly rather than one that crashed. - last_scoring.py passed the extension factor to apply_long_gap in both the open and extension positions, so partialOrderAlignmentTrainedGapOpen2Factor was read from the config and discarded. Both default to 3, so the shipped config is unaffected; setting them differently silently did nothing, and setting the extension factor to 1 raised an assertion. Co-Authored-By: Claude Opus 5 (1M context) --- .gitmodules | 4 + Makefile | 23 +- README.md | 1 + ReleaseNotes.md | 8 + bar/impl/bar.c | 41 +- bar/impl/poaBarAligner.c | 534 +++++++++++++++++++---- bar/inc/poaBarAligner.h | 83 +++- bar/tests/poaBarTest.c | 179 ++++++-- doc/progressive.md | 1 + include.mk | 29 +- pipeline/cactus_consolidated.c | 4 +- src/cactus/cactus_progressive_config.xml | 79 +++- src/cactus/paf/last_scoring.py | 18 +- src/cactus/pipeline/cactus_workflow.py | 24 +- src/cactus/setup/cactus_align.py | 13 +- submodules/minipoa | 1 + test/evolverTest.py | 41 ++ 17 files changed, 952 insertions(+), 131 deletions(-) create mode 160000 submodules/minipoa diff --git a/.gitmodules b/.gitmodules index fc01b14c1..d2d8ac08b 100644 --- a/.gitmodules +++ b/.gitmodules @@ -46,3 +46,7 @@ [submodule "submodules/FASTAN"] path = submodules/FASTAN url = https://github.com/thegenemyers/FASTAN.git +[submodule "submodules/minipoa"] + path = submodules/minipoa + url = https://github.com/ComparativeGenomicsToolkit/minipoa.git + branch = cactus-integration diff --git a/Makefile b/Makefile index 58bfd58ce..cde0142dc 100644 --- a/Makefile +++ b/Makefile @@ -8,6 +8,9 @@ modules = api setup caf bar hal reference pipeline preprocessor # both cactus and sonLib # jemalloc is conditionally added based on include.mk settings submodules1 = sonLib cPecan hal matchingAndOrdering pinchesAndCacti abPOA lastz paffy red collapse-bubble FASTGA FASTAN alntools +ifeq ($(minipoa),on) +submodules1 += minipoa +endif submodules2 = cactus2hal submodules = ${submodules1} ${submodules2} @@ -208,6 +211,14 @@ evolver_test_update_branch_local: ${CWD}/test/mammals-truth.maf evolver_test_poa_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::testEvolverPOALocal +# minipoa, same dataset and tolerance as evolver_test_poa_local, so the two are comparable +evolver_test_minipoa_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::testEvolverMinipoaLocal + +# the mammals head-to-head against evolver_test_local, which runs the same data through abpoa +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 + 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 @@ -230,7 +241,7 @@ evolver_test_primates_pangenome_steps_mgsplit_docker: all ${CWD}/test/primates-t PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=docker CACTUS_DOCKER_MODE=1 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeStepByStepSplitDocker -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_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_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 @@ -315,6 +326,16 @@ suball.abPOA: ln -f submodules/abPOA/include/*.h ${INCLDIR} rm -fr ${INCLDIR}/simde && cp -r submodules/abPOA/include/simde ${INCLDIR} +# minipoa carries a plain Makefile for this, because cactus does not drive cmake and because +# upstream's CMakeLists picks its SIMD backend by probing the build host. It reads the same +# avx2/sse41/sse2/armv8 environment variables include.mk already exports for abPOA, and produces a +# fixed archive name, so there is no need for abPOA's cascade of if-file-exists links. +# Only minipoa_c.h is published: the rest of minipoa's headers are C++ and bar is C. +suball.minipoa: + cd submodules/minipoa && ${MAKE} + ln -f submodules/minipoa/lib/libminipoa.a ${LIBDIR}/libminipoa.a + ln -f submodules/minipoa/include/minipoa_c.h ${INCLDIR} + suball.lastz: suball.jemalloc # Inject ${LIBS} into lastz's link lines so jemalloc reaches it. This must not assume the # link line is still pristine: makeBinRelease seds 's/-lm/-lm -static/g' into this same file diff --git a/README.md b/README.md index 25d834b2a..8d60c9329 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,7 @@ Cactus uses many different algorithms and individual code contributions, princip - Melissa Jane Hubiz and Adam Siepel for halPhyloP and [Phast](http://compgen.cshl.edu/phast/). - B Gulhan, R Burhans, R Harris, M Kandemir, M Haeussler, A Nekrutenko for [KegAlign](https://github.com/galaxyproject/KegAlign), the GPU-accelerated version of LastZ. - Yan Gao et al. for [abPOA](https://github.com/yangao07/abPOA) +- Haodong Liu et al. for [minipoa](https://github.com/NCl3-lhd/minipoa), an optional alternative base aligner in BAR: "Minipoa: A minimizer-based method for fast and memory-efficient partial order alignment" ([doi:10.64898/2026.02.18.706716](https://doi.org/10.64898/2026.02.18.706716)) - Heng Li for [minigraph](https://github.com/lh3/minigraph), [minimap2](https://github.com/lh3/minimap2), [gfatools](https://github.com/lh3/gfatools) and [dna-brnn](https://github.com/lh3/dna-rnn) - Dany Doerr for [GFAffix](https://github.com/marschall-lab/GFAffix), used to optionally clean pangenome graphs. - The vg team for [vg](https://github.com/vgteam/vg), used to process pangenome graphs. diff --git a/ReleaseNotes.md b/ReleaseNotes.md index 8cac23f84..6f6cd4e59 100644 --- a/ReleaseNotes.md +++ b/ReleaseNotes.md @@ -1,3 +1,11 @@ +# Unreleased + +- `minipoa` added as an optional third base aligner in BAR, selected with `` in the config. abPOA remains the default. On the evolver tests the two are equivalent: primates 0.9974 against abPOA's 0.9973, mammals 0.8506 against 0.8513 -- at roughly half the peak memory (1.0 GB against 2.3 GB on evolver mammals). +- minipoa shares abPOA's substitution matrix (from ``) but has its own gap penalties in ``, because the same nominal values do not produce the same alignment. abPOA's gap model is convex, so its effective extension past ~28bp is 1, not 30; minipoa has a single affine piece and charging it 30 made long gaps ~30x dearer, costing six points of mammals accuracy. `--lastTrain` fits a *single* affine gap model, which is exactly minipoa's model, so the learned gaps are passed to minipoa verbatim; abPOA additionally gets a synthesised second piece it needs for stability, which minipoa does not. +- `` replaces the `` boolean as the way to choose between pecan, abpoa and minipoa. Configs with no `baseAligner` still work: the boolean is used as before. +- Fix `--lastTrain` ignoring `partialOrderAlignmentTrainedGapOpen2Factor`: `apply_long_gap` was being passed the extension factor in both the open and extension positions, so the open factor was read from the config and discarded. Both default to 3, so results with the shipped config are unchanged; setting them to different values previously did nothing (and setting the extension factor to 1 raised an assertion). +- `CACTUS_BAR_DUMP_DIR` dumps every window handed to the base aligner, with a command line that replays it. This was previously a commented-out `#define` that also deleted its own output. + # Release 3.3.0 2026-08-11 This release a few general patches and new pangenome functionality such as graph reference (gref), panacus and panpatch support. diff --git a/bar/impl/bar.c b/bar/impl/bar.c index c6dc446fa..c0860193e 100644 --- a/bar/impl/bar.c +++ b/bar/impl/bar.c @@ -60,8 +60,13 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * //Parse the many, many necessary parameters from the params file ////////////////////////////////////////////// + bar_dump_dir_init(); // before any thread starts, so the flower loop never calls getenv + int64_t maximumLength = cactusParams_get_int(params, 2, "bar", "bandingLimit"); - int64_t usePoa = cactusParams_get_int(params, 2, "bar", "partialOrderAlignment"); + BaseAligner engine = baseAligner_constructFromCactusParams(params); + // Every site that used to ask "poa or pecan?" still only needs that much. abpoa and minipoa + // both produce an MSA and so share the AlignmentBlock/stPinch path; only pecan differs. + bool usePoa = engine != BASE_ALIGNER_PECAN; // Pecan prams int64_t spanningTrees = cactusParams_get_int(params, 3, "bar", "pecan", "spanningTrees"); @@ -71,14 +76,30 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * StateMachine *sM = stateMachine5_construct(fiveState); bool pruneOutStubAlignments = cactusParams_get_int(params, 3, "bar", "pecan", "pruneOutStubAlignments"); - // Poa params - // toggle from pecan to abpoa for multiple alignment, by setting to non-zero - // Note that poa uses about N^2 memory, so maximum value is generally in 10s of kb - int64_t poaWindow = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentWindow"); - int64_t maskFilter = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMaskFilter"); + // Poa params. The window and the mask filter come from whichever engine is selected: + // abpoa's memory is quadratic in the window so it is held to 10s of kb, while minipoa is the + // reason to have a second engine at all and can take a much larger one. + int64_t poaWindow, maskFilter; + if (engine == BASE_ALIGNER_MINIPOA) { + poaWindow = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaWindow"); + maskFilter = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaMaskFilter"); + } else { + poaWindow = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentWindow"); + maskFilter = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMaskFilter"); + } + // abpoa-only progressive-mode guards; inert for minipoa, which has no such mode enabled. int64_t poaMaxProgRows = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentProgressiveMaxRows"); double poaMaxLenDiff = cactusParams_get_float(params, 3, "bar", "poa", "partialOrderAlignmentProgressiveMaxLengthDiff"); - abpoa_para_t *poaParameters = usePoa ? abpoaParamaters_constructFromCactusParams(params) : NULL; + PoaParameters *poaParameters = poaParameters_constructFromCactusParams(params, engine); + + // Say which aligner ran. Without this a report of a suspect alignment cannot be tied back to + // the engine or the settings that produced it. + if (usePoa) { + 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); + } ////////////////////////////////////////////// //Run the bar algorithm @@ -160,7 +181,7 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * void *alignments; if (usePoa) { /* - * This makes a consistent set of alignments using abPoa. + * This makes a consistent set of alignments using the selected POA engine. * * It does not use any precomputed alignments, if they are provided they will be ignored */ @@ -284,7 +305,5 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * pairwiseAlignmentBandingParameters_destruct(pairwiseAlignmentParameters); stateMachine_destruct(sM); - if (poaParameters) { - abpoa_free_para(poaParameters); - } + poaParameters_destruct(poaParameters); } diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index c1c4937d5..ea68835b6 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -6,15 +6,40 @@ #include #include "abpoa.h" +#ifdef HAVE_MINIPOA +#include "minipoa_c.h" +#endif #include "poaBarAligner.h" #include "flowerAligner.h" #include #include +#include + +/* + * Set CACTUS_BAR_DUMP_DIR in the environment and every window handed to the base aligner is + * written there as a FASTA, alongside the substitution matrix and a command line that replays it. + * That is how a suspect alignment gets out of a 64-thread run and onto a developer's laptop. + * + * This used to be a commented-out #define that also deleted its own output on success, so using + * it meant editing and rebuilding, and the files were gone by the time you looked. Read once in + * bar(), before the flower loop, so the OpenMP region never calls getenv. + */ +static const char *bar_dump_dir = NULL; +static int64_t bar_dump_counter = 0; + +void bar_dump_dir_init(void) { + bar_dump_dir = getenv("CACTUS_BAR_DUMP_DIR"); + if (bar_dump_dir != NULL && bar_dump_dir[0] == '\0') { + bar_dump_dir = NULL; // set-but-empty means off, not a dump into the current directory + } + if (bar_dump_dir != NULL) { + st_logInfo("bar: dumping base-aligner windows to %s\n", bar_dump_dir); + } +} -// FOR DEBUGGING ONLY: Specify directory where abpoa inputs get dumped -//#define CACTUS_ABPOA_MSA_DUMP_DIR "/home/hickey/dev/cactus/dump" -// FOR DEBUGGING ONLY: Run abpoa from command line instead of via API (only works with CACTUS_ABPOA_MSA_DUMP_DIR defined) +// FOR DEBUGGING ONLY: run abpoa from the command line instead of via the API. +// Requires CACTUS_BAR_DUMP_DIR to be set at run time. //#define CACTUS_ABPOA_FROM_COMMAND_LINE // OpenMP @@ -167,7 +192,79 @@ static inline uint8_t msa_to_rc(uint8_t n) { return rc_table[n]; } -#ifdef CACTUS_ABPOA_MSA_DUMP_DIR +/* + * Write the window as FASTA plus its 5x5 matrix. The matrix file format is the one both abpoa -t + * and minipoa -m read, so either aligner can be pointed straight at it. + */ +static void dump_window_fasta_and_matrix(Msa *msa, uint8_t **bseqs, const int *mat, + const char *input_path, const char *matrix_path) { + FILE *mat_file = fopen(matrix_path, "w"); + if (mat_file != NULL) { + fprintf(mat_file, "\tA\tC\tG\tT\tN\n"); + for (size_t i = 0; i < 5; ++i) { + fprintf(mat_file, "%c", "ACGTN"[i]); + for (size_t j = 0; j < 5; ++j) { + fprintf(mat_file, "\t%d", mat[i * 5 + j]); + } + fprintf(mat_file, "\n"); + } + fclose(mat_file); + } + FILE *fa_file = fopen(input_path, "w"); + if (fa_file == NULL) { + return; + } + for (int64_t i = 0; i < msa->seq_no; ++i) { + fprintf(fa_file, ">%" PRIi64 "\n", i); + for (int64_t j = 0; j < msa->seq_lens[i]; ++j) { + fputc(msa_to_base(bseqs[i][j]), fa_file); + } + fputc('\n', fa_file); + } + fclose(fa_file); +} + +/* + * A command line that replays this window through minipoa. Gap penalties are negated: minipoa + * maximises, so penalties are negative there, while cactus and abpoa carry them positive. + */ +static char *dump_minipoa_input(Msa *msa, PoaParameters *pp, uint8_t **bseqs, char *input_path, + char *matrix_path, char *command_path, char *output_path) { + dump_window_fasta_and_matrix(msa, bseqs, pp->mat, input_path, matrix_path); + + int f = pp->bandFraction > 0.0 ? (int)(1.0 / pp->bandFraction + 0.5) : 0; + 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->seeding) { + char kw_opts[128]; + sprintf(kw_opts, " -S -k %d -w %d", pp->minimizerK, pp->minimizerW); + strcat(command, kw_opts); + if (pp->anchorWindow > 0) { + sprintf(kw_opts, " -W %d", pp->anchorWindow); + strcat(command, kw_opts); + } + } + // -p and -B are on in the shipped config, and progressive ordering in particular changes the + // alignment, so a replay that omitted them would not reproduce the window it is meant to + // explain. + if (pp->progressive) { + strcat(command, " -p"); + } + if (pp->adaptiveBand) { + strcat(command, " -B"); + } + strcat(command, " > "); + strcat(command, output_path); + + FILE *cmd_file = fopen(command_path, "w"); + if (cmd_file != NULL) { + fprintf(cmd_file, "%s\n", command); + fclose(cmd_file); + } + return command; +} + // dump the abpoa input to files, and return a command line for running abpoa on them char* dump_abpoa_input(Msa* msa, abpoa_para_t* abpt, uint8_t **bseqs, char* abpoa_input_path, char* abpoa_matrix_path, char* abpoa_command_path, char* abpoa_output_path) { @@ -227,7 +324,6 @@ char* dump_abpoa_input(Msa* msa, abpoa_para_t* abpt, uint8_t **bseqs, char* abpo return abpoa_command; } -#endif #ifdef CACTUS_ABPOA_FROM_COMMAND_LINE void abpoa_msa_from_command_line(char* abpoa_command_line, char* abpoa_output_path, uint8_t*** msa_seq, int* col_no) { @@ -461,8 +557,364 @@ static void msa_fix_trimmed(Msa* msa) { msa->column_no -= empty_columns; } +/* + * picks the engine. It is read with cactusParams_has so that a config + * predating the attribute -- including any a user has saved -- still selects what it used to via + * the older boolean. + */ +BaseAligner baseAligner_constructFromCactusParams(CactusParams *params) { + /* + * Both attributes are optional, in both directions. An old config has only + * partialOrderAlignment; a config written from now on may reasonably have only baseAligner -- + * including one produced by following the warning below, which tells the user to delete the + * legacy attribute. Reading either unguarded would st_errAbort on the other's config. + */ + bool hasLegacy = cactusParams_has(params, 2, "bar", "partialOrderAlignment"); + int64_t usePoa = hasLegacy ? cactusParams_get_int(params, 2, "bar", "partialOrderAlignment") : 1; + if (!cactusParams_has(params, 2, "bar", "baseAligner")) { + return usePoa ? BASE_ALIGNER_ABPOA : BASE_ALIGNER_PECAN; + } + char *name = cactusParams_get_string(params, 2, "bar", "baseAligner"); + BaseAligner engine; + if (strcmp(name, "pecan") == 0) { + engine = BASE_ALIGNER_PECAN; + } else if (strcmp(name, "abpoa") == 0) { + engine = BASE_ALIGNER_ABPOA; + } else if (strcmp(name, "minipoa") == 0) { + engine = BASE_ALIGNER_MINIPOA; + } else { + st_errAbort("Unknown ; expected pecan, abpoa or minipoa", name); + engine = BASE_ALIGNER_ABPOA; /* not reached */ + } + free(name); + + /* + * Both attributes present and disagreeing is worth saying out loud. partialOrderAlignment="0" + * is how the config has always documented "use pecan", so someone who sets it and gets abpoa + * anyway should not have to discover that from the alignment. + */ + bool poaImplied = engine != BASE_ALIGNER_PECAN; + if (hasLegacy && (usePoa != 0) != poaImplied) { + st_logCritical("Warning: overrides , which asks for the opposite. baseAligner wins; remove the other to silence this.\n", + baseAligner_toString(engine), usePoa); + } + return engine; +} + +const char *baseAligner_toString(BaseAligner engine) { + switch (engine) { + case BASE_ALIGNER_PECAN: return "pecan"; + case BASE_ALIGNER_ABPOA: return "abpoa"; + case BASE_ALIGNER_MINIPOA: return "minipoa"; + } + return "unknown"; +} + +/* + * minipoa's substitution matrix comes from ; its gap penalties do not. + * + * Sharing the matrix keeps the two engines scoring substitutions identically, and means + * 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. + * + * 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 + * GapOpen2/GapExtend2 pair it synthesises for abPOA is a stability workaround for that aligner + * and is deliberately not passed on. + */ +#ifdef HAVE_MINIPOA +static minipoa_para_t *minipoaParameters_constructFromCactusParams(CactusParams *params, PoaParameters *out) { + minipoa_para_t *mpt = minipoa_init_para(); + if (mpt == NULL) { + st_errAbort("Failed to allocate minipoa parameters: %s", minipoa_last_error()); + } + + /* + * minipoa's own matrix, falling back to abPOA's when it is left empty -- that fallback is the + * only way to get last-train's learned scores into minipoa, since last_scoring.py writes into + * . + */ + char *submat_string = cactusParams_get_string(params, 3, "bar", "minipoa", "minipoaSubMatrix"); + if (submat_string == NULL || strlen(submat_string) == 0) { + free(submat_string); + submat_string = cactusParams_get_string(params, 3, "bar", "poa", "partialOrderAlignmentSubMatrix"); + } + if (submat_string != NULL && strlen(submat_string) > 0) { + int mat[25]; + int count = 0; + for (char *val = strtok(submat_string, " "); val != NULL && count < 25; val = strtok(NULL, " ")) { + mat[count++] = atoi(val); + } + if (count != 25) { + st_errAbort(" needs 25 values, got %d", count); + } + minipoa_set_score_matrix(mpt, mat); + memcpy(out->mat, mat, sizeof(mat)); + } + 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); + + out->bandConstant = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaBandConstant"); + out->bandFraction = cactusParams_get_float(params, 3, "bar", "minipoa", "minipoaBandFraction"); + minipoa_set_band(mpt, out->bandConstant, out->bandFraction); + out->adaptiveBand = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaAdaptiveBand"); + minipoa_set_adaptive_band(mpt, out->adaptiveBand); + + out->seeding = !cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaDisableSeeding"); + out->minimizerK = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaMinimizerK"); + out->minimizerW = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaMinimizerW"); + out->anchorWindow = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaAnchorWindow"); + 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. + */ + out->progressive = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaProgressiveMode"); + out->progressiveMaxRows = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaProgressiveMaxRows"); + minipoa_set_progressive(mpt, out->progressive); + return mpt; +} +#else +/* + * minipoa=off in include.mk. Selecting it is a configuration error rather than a crash, and the + * message has to name the build switch -- the config is portable between machines, the build is + * not. + */ +static void *minipoaParameters_constructFromCactusParams(CactusParams *params, PoaParameters *out) { + (void)params; (void)out; + st_errAbort(" was selected, but this cactus was built with " + "minipoa=off (see include.mk). Rebuild with minipoa=on, or choose abpoa or pecan."); + return NULL; +} +#endif + +PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, BaseAligner engine) { + if (engine == BASE_ALIGNER_PECAN) { + return NULL; + } + PoaParameters *poaParameters = st_calloc(1, sizeof(PoaParameters)); + poaParameters->engine = engine; + if (engine == BASE_ALIGNER_ABPOA) { + poaParameters->abpt = abpoaParamaters_constructFromCactusParams(params); + } else { + poaParameters->mpt = minipoaParameters_constructFromCactusParams(params, poaParameters); +#ifdef HAVE_MINIPOA + if (poaParameters->progressive) { + // Same settings with the guide tree off, for windows too wide to afford a dense NxN + // distance matrix. abPOA caps this the same way, via partialOrderAlignmentProgressiveMaxRows. + PoaParameters scratch = *poaParameters; + minipoa_para_t *plain = minipoaParameters_constructFromCactusParams(params, &scratch); + minipoa_set_progressive(plain, 0); + poaParameters->mptNoProgressive = plain; + } +#endif + } + return poaParameters; +} + +void poaParameters_destruct(PoaParameters *poaParameters) { + if (poaParameters == NULL) { + return; + } + if (poaParameters->abpt != NULL) { + abpoa_free_para(poaParameters->abpt); + } +#ifdef HAVE_MINIPOA + if (poaParameters->mpt != NULL) { + minipoa_free_para((minipoa_para_t *)poaParameters->mpt); + } + if (poaParameters->mptNoProgressive != NULL) { + minipoa_free_para((minipoa_para_t *)poaParameters->mptNoProgressive); + } +#endif + free(poaParameters); +} + +/* + * The minipoa half of run_poa_window(). + * + * minipoa_msa() returns a status rather than exiting: its internal band/backtrack traps used to + * call exit(0) -- a success code -- which from here is indistinguishable from a clean run with + * truncated output. A failure here is fatal for the job either way, but it says which window and + * leaves it on disk when dumping is on, which is the difference between a bug report and a shrug. + */ +static void run_minipoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_parameters) { +#ifndef HAVE_MINIPOA + (void)msa; (void)bseqs; (void)poa_parameters; + st_errAbort("minipoa was selected but this cactus was built with minipoa=off (see include.mk)"); +#else + char input_path[1024], matrix_path[1024], command_path[1024], output_path[1024]; + char *command_line = NULL; + if (bar_dump_dir != NULL) { + int64_t dump_id; +#if defined(_OPENMP) +#pragma omp atomic capture +#endif + dump_id = ++bar_dump_counter; + sprintf(input_path, "%s/bar_window_%d_%" PRIi64 ".fa", bar_dump_dir, (int)getpid(), dump_id); + sprintf(matrix_path, "%s.mat", input_path); + sprintf(command_path, "%s.cmd", input_path); + sprintf(output_path, "%s.out", input_path); + command_line = dump_minipoa_input(msa, poa_parameters, bseqs, input_path, matrix_path, + command_path, output_path); + } + + const minipoa_para_t *mpt = (const minipoa_para_t *)poa_parameters->mpt; + if (poa_parameters->mptNoProgressive != NULL && msa->seq_no > poa_parameters->progressiveMaxRows) { + mpt = (const minipoa_para_t *)poa_parameters->mptNoProgressive; + } + + int column_no = 0; + uint8_t **msa_seq = NULL; + int ret = minipoa_msa(mpt, (int)msa->seq_no, msa->seq_lens, bseqs, &msa_seq, &column_no); + if (ret != 0) { + st_errAbort("minipoa failed on a %" PRIi64 " x %d window: %s. %s", + msa->seq_no, msa->seq_lens[0], minipoa_last_error(), + command_line != NULL + ? command_path + : "Set CACTUS_BAR_DUMP_DIR to capture the window that did this."); + } + msa->msa_seq = msa_seq; + msa->column_no = column_no; + /* + * A global alignment cannot have fewer columns than its longest input row, so this catches a + * truncated or empty result -- which is the shape minipoa failures take when they do not + * return an error. O(rows), next to nothing against the DP, and worth it: an all-gap MSA + * produces zero alignment blocks and BAR would otherwise drop the flower in silence. + */ + int longest = 0; + for (int64_t i = 0; i < msa->seq_no; i++) { + if (msa->seq_lens[i] > longest) { + longest = msa->seq_lens[i]; + } + } + if (column_no < longest) { + st_errAbort("minipoa returned %d columns for a %" PRIi64 " x %d window, which cannot be a " + "global alignment of it. %s", column_no, msa->seq_no, longest, + command_line != NULL + ? command_path + : "Set CACTUS_BAR_DUMP_DIR to capture the window that did this."); + } + if (command_line != NULL) { + free(command_line); + } +#endif +} + +/* + * Align one window with whichever engine was selected, filling msa->msa_seq and msa->column_no. + * + * The contract both engines meet: seq_no rows of column_no bytes, one malloc per row plus one for + * the row array (msa_destruct frees them that way), rows in input order, values 0-4 for ACGTN and + * 5 for a gap. Everything around this -- the sliding window, the empty-sequence hack, trimming, + * stitching, block extraction -- is engine-neutral and shared. + * + * max_prog_rows and max_prog_length_diff are abPOA's progressive-mode guards; minipoa has no + * equivalent knob, so they are inert for it. + */ +static void run_poa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_parameters, + int64_t max_prog_rows, double max_prog_length_diff) { + if (poa_parameters->engine == BASE_ALIGNER_MINIPOA) { + run_minipoa_window(msa, bseqs, poa_parameters); + return; + } + + // init abpoa + abpoa_t *ab = abpoa_init(); + abpoa_para_t *abpt = copy_abpoa_params(poa_parameters->abpt); + if (msa->seq_no > max_prog_rows || + // note: these are sorted by length excep in unit tests + (1. - (double)msa->seq_lens[msa->seq_no-1] / (double)msa->seq_lens[0] > max_prog_length_diff)) { + abpt->progressive_poa = 0; + } + + // dump the input to file, if asked to at run time + char abpoa_input_path[1024], abpoa_matrix_path[1024], abpoa_command_path[1024], abpoa_output_path[1024]; + char *abpoa_command_line = NULL; + if (bar_dump_dir != NULL) { + // The old name keyed off the Msa pointer, so two windows that reused the same freed + // allocation silently overwrote each other. pid + counter is unique for the run. + int64_t dump_id; +#if defined(_OPENMP) +#pragma omp atomic capture +#endif + dump_id = ++bar_dump_counter; + sprintf(abpoa_input_path, "%s/bar_window_%d_%" PRIi64 ".fa", bar_dump_dir, (int)getpid(), dump_id); + sprintf(abpoa_matrix_path, "%s.mat", abpoa_input_path); + sprintf(abpoa_command_path, "%s.cmd", abpoa_input_path); + sprintf(abpoa_output_path, "%s.out", abpoa_input_path); + abpoa_command_line = dump_abpoa_input(msa, abpt, bseqs, + abpoa_input_path, abpoa_matrix_path, abpoa_command_path, abpoa_output_path); + } + +#ifdef CACTUS_ABPOA_FROM_COMMAND_LINE + // run abpoa from the command line + if (abpoa_command_line == NULL) { + st_errAbort("CACTUS_ABPOA_FROM_COMMAND_LINE needs CACTUS_BAR_DUMP_DIR set in the environment"); + } + abpoa_msa_from_command_line(abpoa_command_line, abpoa_output_path, &(msa->msa_seq), &(msa->column_no)); + + int test_cols = 0; + uint8_t** test_msa = NULL; + abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); + // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they + // end up in the ab->abc struct -- we extract them here + test_msa = ab->abc->msa_base; + ab->abc->msa_base = NULL; + test_cols = ab->abc->msa_len; + + // sanity check to make sure we get the same output + assert(msa->column_no == test_cols); + for (int i = 0; i < msa->seq_no; ++i) { + for (int j = 0; j < test_cols; ++j) { + //todo: not sure why this doesn't work anymore !!!! + //assert(test_msa[i][j] == msa->msa_seq[i][j]); + } + free(test_msa[i]); + } + free(test_msa); +#else + // perform abpoa-msa + abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); + // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they + // end up in the ab->abc struct -- we extract them here + msa->msa_seq = ab->abc->msa_base; + ab->abc->msa_base = NULL; + msa->column_no = ab->abc->msa_len; +#endif + + // The dumps are kept. Deleting them on success made the switch useless for its main + // job -- looking at a window that aligned badly rather than one that crashed. + if (abpoa_command_line != NULL) { + free(abpoa_command_line); + } + + // free abpoa + abpoa_free(ab); + abpoa_free_para(abpt); +} + Msa *msa_make_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no, int64_t window_size, - int64_t max_prog_rows, double max_prog_length_diff, abpoa_para_t *poa_parameters) { + int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { assert(seq_no > 0); @@ -562,71 +1014,7 @@ Msa *msa_make_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no } } - // init abpoa - abpoa_t *ab = abpoa_init(); - abpoa_para_t *abpt = copy_abpoa_params(poa_parameters); - if (msa->seq_no > max_prog_rows || - // note: these are sorted by length excep in unit tests - (1. - (double)msa->seq_lens[msa->seq_no-1] / (double)msa->seq_lens[0] > max_prog_length_diff)) { - abpt->progressive_poa = 0; - } - -#ifdef CACTUS_ABPOA_MSA_DUMP_DIR - // dump the input to file - char abpoa_input_path[1024], abpoa_matrix_path[1024], abpoa_command_path[1024], abpoa_output_path[1024]; - sprintf(abpoa_input_path, "%s/ap_in_%ld.fa", CACTUS_ABPOA_MSA_DUMP_DIR, (int64_t)msa); - sprintf(abpoa_matrix_path, "%s.mat", abpoa_input_path); - sprintf(abpoa_command_path, "%s.cmd", abpoa_input_path); - sprintf(abpoa_output_path, "%s.out", abpoa_input_path); - char* abpoa_command_line = dump_abpoa_input(msa, abpt, bseqs, - abpoa_input_path, abpoa_matrix_path, abpoa_command_path, abpoa_output_path); -#endif - -#ifdef CACTUS_ABPOA_FROM_COMMAND_LINE - // run abpoa from the command line - abpoa_msa_from_command_line(abpoa_command_line, abpoa_output_path, &(msa->msa_seq), &(msa->column_no)); - - int test_cols = 0; - uint8_t** test_msa = NULL; - abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); - // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they - // end up in the ab->abc struct -- we extract them here - test_msa = ab->abc->msa_base; - ab->abc->msa_base = NULL; - test_cols = ab->abc->msa_len; - - // sanity check to make sure we get the same output - assert(msa->column_no == test_cols); - for (int i = 0; i < msa->seq_no; ++i) { - for (int j = 0; j < test_cols; ++j) { - //todo: not sure why this doesn't work anymore !!!! - //assert(test_msa[i][j] == msa->msa_seq[i][j]); - } - free(test_msa[i]); - } - free(test_msa); -#else - // perform abpoa-msa - abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); - // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they - // end up in the ab->abc struct -- we extract them here - msa->msa_seq = ab->abc->msa_base; - ab->abc->msa_base = NULL; - msa->column_no = ab->abc->msa_len; -#endif - -#ifdef CACTUS_ABPOA_MSA_DUMP_DIR - // we got this far without crashing, so delete the dumped file (they can really pile up otherwise) - remove(abpoa_input_path); - remove(abpoa_matrix_path); - remove(abpoa_command_path); - remove(abpoa_output_path); - free(abpoa_command_line); -#endif - - // free abpoa - abpoa_free(ab); - abpoa_free_para(abpt); + run_poa_window(msa, bseqs, poa_parameters, max_prog_rows, max_prog_length_diff); // mask out empty sequences that were phonied in as Ns above for (int64_t i = 0; i < msa->seq_no && emptyCount > 0; ++i) { @@ -762,7 +1150,7 @@ Msa *msa_make_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no Msa **make_consistent_partial_order_alignments(int64_t end_no, int64_t *end_lengths, char ***end_strings, int **end_string_lengths, int64_t **right_end_indexes, int64_t **right_end_row_indexes, int64_t **overlaps, - int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, abpoa_para_t *poa_parameters) { + int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { // Calculate the initial, potentially inconsistent msas and column scores for each msa float *column_scores[end_no]; Msa **msas = st_malloc(sizeof(Msa *) * end_no); @@ -1148,7 +1536,7 @@ int64_t getMaxSequenceLength(End *end) { } stList *make_flower_alignment_poa(Flower *flower, int64_t max_seq_length, int64_t window_size, int64_t mask_filter, - int64_t max_prog_rows, double max_prog_length_diff, abpoa_para_t * poa_parameters) { + int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { End *dominantEnd = getDominantEnd(flower); int64_t seq_no = dominantEnd != NULL ? end_getInstanceNumber(dominantEnd) : -1; if(dominantEnd != NULL && getMaxSequenceLength(dominantEnd) < max_seq_length) { diff --git a/bar/inc/poaBarAligner.h b/bar/inc/poaBarAligner.h index 260d91846..2a481e74a 100644 --- a/bar/inc/poaBarAligner.h +++ b/bar/inc/poaBarAligner.h @@ -20,11 +20,82 @@ */ void bar(stList *flowers, CactusParams *p, CactusDisk *cactusDisk, stList *listOfEndAlignmentFiles); +/* + * Read CACTUS_BAR_DUMP_DIR once, before any alignment runs. When set, every window handed to the + * base aligner is dumped there with a command line that replays it. + */ +void bar_dump_dir_init(void); + /* * Construct a pairwise alignment parameters object parsing the cactus params specified parameters. */ PairwiseAlignmentParameters *pairwiseAlignmentParameters_constructFromCactusParams(CactusParams *params); +/** + * Which engine computes the base-level multiple alignment. + * + * Selected by . Configs written before that attribute + * existed fall back to the older boolean, so nothing needs + * rewriting. + */ +typedef enum { + BASE_ALIGNER_PECAN = 0, + BASE_ALIGNER_ABPOA = 1, + BASE_ALIGNER_MINIPOA = 2 +} BaseAligner; + +/** + * Read the selected engine out of the cactus params. + */ +BaseAligner baseAligner_constructFromCactusParams(CactusParams *params); + +/** + * Name of an engine, for logging and error messages. + */ +const char *baseAligner_toString(BaseAligner engine); + +/** + * Engine-agnostic handle for whatever the chosen aligner needs. + * + * mpt is void* rather than minipoa_para_t* on purpose: this header is included by C code that + * is compiled whether or not minipoa was built (see the minipoa=on|off switch in include.mk), + * and nothing outside poaBarAligner.c ever dereferences it. + */ +typedef struct _PoaParameters { + BaseAligner engine; + abpoa_para_t *abpt; /* engine == BASE_ALIGNER_ABPOA */ + void *mpt; /* engine == BASE_ALIGNER_MINIPOA, a minipoa_para_t * */ + /* + * The same minipoa parameters with progressive ordering forced off, for windows with too many + * rows to afford the guide tree. abPOA applies that cap by deep-copying its params per + * window; minipoa's handle is immutable from here, so keeping a second one is simpler and + * costs a few hundred bytes for the whole run. NULL when progressive is off anyway. + */ + void *mptNoProgressive; + int64_t progressiveMaxRows; + /* + * The same settings again, in plain scalars, purely so a dumped window can carry a command + * line that replays it. minipoa's handle is opaque and abPOA's are spread over abpoa_para_t; + * keeping a copy here is cheaper than accessors on both, and it is written once per run. + */ + int mat[25]; + int gapOpen, gapExt; + int bandConstant; + double bandFraction; + bool seeding; + bool adaptiveBand; + bool progressive; + int minimizerK, minimizerW, anchorWindow; +} PoaParameters; + +/** + * Build the parameters for the given engine from the cactus params. Returns NULL for + * BASE_ALIGNER_PECAN, which has its own parameters. Free with poaParameters_destruct(). + */ +PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, BaseAligner engine); + +void poaParameters_destruct(PoaParameters *poaParameters); + /** * Construct the abpoa parameters object parsing the cactus params specified parameters. * It needs to get freed with abpoa_free_para(abpt); @@ -70,7 +141,7 @@ void msa_print(Msa *msa, FILE *f); * @param window_size Sliding window size which limits length of poa sub-alignments. Memory usage is quardatic in this. * @param max_prog_rows Disable abpoas progressive alignment if there are more than this many rows (avoid quadratic dist mat blowup) * @param max_prog_length_diff Disable abpoa's progresive alignment if the 1 - shortest (last) sequence / longest (first) sequence is more than this - * @param poa_parameters abpoa parameters + * @param poa_parameters base aligner parameters * @return An msa of the strings. */ Msa *msa_make_partial_order_alignment(char **seqs, @@ -79,7 +150,7 @@ Msa *msa_make_partial_order_alignment(char **seqs, int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, - abpoa_para_t *poa_parameters); + PoaParameters *poa_parameters); /** * Takes a set of ends and returns a set of consistent multiple alignments, @@ -102,12 +173,12 @@ Msa *msa_make_partial_order_alignment(char **seqs, * @param window_size Sliding window size which limits length of poa sub-alignments. Memory usage is quardatic in this. * @param max_prog_rows Disable abpoas progressive alignment if there are more than this many rows (avoid quadratic dist mat blowup) * @param max_prog_length_diff Disable abpoa's progresive alignment if the 1 - shortest (last) sequence / longest (first) sequence is more than this - * @param poa_parameters abpoa parameters + * @param poa_parameters base aligner parameters * @return A consistent Msa for each end */ Msa **make_consistent_partial_order_alignments(int64_t end_no, int64_t *end_lengths, char ***end_strings, int **end_string_lengths, int64_t **right_end_indexes, int64_t **right_end_row_indexes, int64_t **overlaps, - int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, abpoa_para_t *poa_parameters); + int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters); /** * Represents a gapless alignment of a set of sequences. @@ -151,7 +222,7 @@ char *get_adjacency_string_and_overlap(Cap *cap, int *length, int64_t *overlap, * @param mask_filter Trim input sequences if encountering this many consecutive soft of hard masked bases (0 = disabled) * @param max_prog_rows Disable abpoa's progressive alignment if there are more than this many rows (avoid quadratic dist mat blowup) * @param max_prog_length_diff Disable abpoa's progresive alignment if the 1 - shortest (last) sequence / longest (first) sequence is more than this - * @param poa_parameters abpoa parameters + * @param poa_parameters base aligner parameters */ stList *make_flower_alignment_poa(Flower *flower, int64_t max_seq_length, @@ -159,7 +230,7 @@ stList *make_flower_alignment_poa(Flower *flower, int64_t mask_filter, int64_t max_prog_rows, double max_prog_length_diff, - abpoa_para_t * poa_parameters); + PoaParameters *poa_parameters); /** * Create a pinch iterator for a list of alignment blocks. diff --git a/bar/tests/poaBarTest.c b/bar/tests/poaBarTest.c index 854d1caa3..64b8124e8 100644 --- a/bar/tests/poaBarTest.c +++ b/bar/tests/poaBarTest.c @@ -1,3 +1,4 @@ +#define _GNU_SOURCE /* mkstemp, under -std=c99 */ /* * Copyright (C) 2009-2011 by Benedict Paten (benedictpaten@gmail.com) * @@ -7,11 +8,51 @@ #include "flowersShared.h" #include "randomSequences.h" #include "poaBarAligner.h" +#ifdef HAVE_MINIPOA +#include "minipoa_c.h" +#endif #include "stCaf.h" #include +#include +#include #include //#define stderr_logging +/* + * Every test below runs against both MSA engines. The parameters mirror what the config would + * produce, modulo the deliberately tiny band the tests use to keep the random cases quick. + */ +static const BaseAligner TEST_ENGINES[] = { + BASE_ALIGNER_ABPOA, +#ifdef HAVE_MINIPOA + BASE_ALIGNER_MINIPOA, +#endif +}; +#define TEST_ENGINE_NO ((int64_t)(sizeof(TEST_ENGINES) / sizeof(TEST_ENGINES[0]))) + +static PoaParameters *test_poa_params(BaseAligner engine) { + PoaParameters *poaParams = st_calloc(1, sizeof(PoaParameters)); + poaParams->engine = engine; + if (engine == BASE_ALIGNER_ABPOA) { + abpoa_para_t *abpt = abpoa_init_para(); + abpt->wb = 10; + abpt->wf = 0.01; + abpoa_post_set_para(abpt); + poaParams->abpt = abpt; + } +#ifdef HAVE_MINIPOA + else { + minipoa_para_t *mpt = minipoa_init_para(); + minipoa_set_band(mpt, 10, 0.01); + minipoa_set_adaptive_band(mpt, 0); + minipoa_set_seeding(mpt, 0, 19, 10, 0); + minipoa_set_progressive(mpt, 0); + poaParams->mpt = mpt; + } +#endif + return poaParams; +} + /** * Validate MSA. Lengths is an array that is populated with the lengths of the * sequences found on the MSA. @@ -34,10 +75,8 @@ void validate_msa(CuTest *testCase, Msa *msa, int64_t *lengths) { * Repeatedly generate random sets of closely related strings and test that returned msa is valid */ void test_make_partial_order_alignment(CuTest *testCase) { - abpoa_para_t *abpt = abpoa_init_para(); - abpt->wb = 10; - abpt->wf = 0.01; - abpoa_post_set_para(abpt); + for (int64_t engine_i = 0; engine_i < TEST_ENGINE_NO; engine_i++) { + PoaParameters *poaParams = test_poa_params(TEST_ENGINES[engine_i]); for(int64_t test=0; test<100; test++) { for (int64_t poa_window_size = 5; poa_window_size < 120; poa_window_size += 15) { #ifdef stderr_logging @@ -65,7 +104,7 @@ void test_make_partial_order_alignment(CuTest *testCase) { } // generate the alignment - Msa *msa = msa_make_partial_order_alignment(seqs, seq_lens, seq_no, poa_window_size, 1000, 0.02, abpt); + Msa *msa = msa_make_partial_order_alignment(seqs, seq_lens, seq_no, poa_window_size, 1000, 0.02, poaParams); // print the msa #ifdef stderr_logging @@ -84,18 +123,17 @@ void test_make_partial_order_alignment(CuTest *testCase) { free(parent_string); } } - abpoa_free_para(abpt); + poaParameters_destruct(poaParams); + } } /** * Repeatedly generate random sets of two ends connected by set of strings, check that the resulting msa is valid */ void test_make_consistent_partial_order_alignments_two_ends(CuTest *testCase) { - abpoa_para_t *abpt = abpoa_init_para(); - abpt->wb = 10; - abpt->wf = 0.01; - abpoa_post_set_para(abpt); - + for (int64_t engine_i = 0; engine_i < TEST_ENGINE_NO; engine_i++) { + PoaParameters *poaParams = test_poa_params(TEST_ENGINES[engine_i]); + for(int64_t test=0; test<100; test++) { #ifdef stderr_logging fprintf(stderr, "Running test_make_consistent_partial_order_alignments_two_ends, test %i\n", (int)test); @@ -146,7 +184,7 @@ void test_make_consistent_partial_order_alignments_two_ends(CuTest *testCase) { // generate the alignments Msa **msas = make_consistent_partial_order_alignments(end_no, end_lengths, end_strings, end_string_lengths, right_end_indexes, right_end_row_indexes, overlaps, - 1000000, 100, 0.02, abpt); + 1000000, 100, 0.02, poaParams); // print the msas #ifdef stderr_logging @@ -175,16 +213,14 @@ void test_make_consistent_partial_order_alignments_two_ends(CuTest *testCase) { free(msas); free(parent_string); } - abpoa_free_para(abpt); + poaParameters_destruct(poaParams); + } } void test_make_flower_alignment_poa(CuTest *testCase) { + for (int64_t engine_i = 0; engine_i < TEST_ENGINE_NO; engine_i++) { setup(testCase); - - abpoa_para_t *abpt = abpoa_init_para(); - abpt->wb = 10; - abpt->wf = 0.01; - abpoa_post_set_para(abpt); + PoaParameters *poaParams = test_poa_params(TEST_ENGINES[engine_i]); #ifdef stderr_logging fprintf(stderr, "There are %i ends in the flower\n", (int)flower_getEndNumber(flower)); #endif @@ -213,7 +249,7 @@ void test_make_flower_alignment_poa(CuTest *testCase) { } flower_destructEndIterator(endIterator); - stList *alignment_blocks = make_flower_alignment_poa(flower, 2, 1000000, 5, 1000, 0.02, abpt); + stList *alignment_blocks = make_flower_alignment_poa(flower, 2, 1000000, 5, 1000, 0.02, poaParams); for(int64_t i=0; iwb = 10; - abpt->wf = 0.01; - abpoa_post_set_para(abpt); - - stList *alignment_blocks = make_flower_alignment_poa(flower, 10000, 1000000, 5, 50, 0.05, abpt); + stList *alignment_blocks = make_flower_alignment_poa(flower, 10000, 1000000, 5, 50, 0.05, poaParams); - abpoa_free_para(abpt); + poaParameters_destruct(poaParams); #ifdef stderr_logging for(int64_t i=0; i is new; is what every config written before it + * says. The fallback between them is the whole backwards-compatibility guarantee of this change, + * so it gets a test of its own -- including the case that actually bites, a config carrying both + * because it was copied from the shipped one and then edited. + */ +static BaseAligner engine_for_config(CuTest *testCase, const char *bar_attrs) { + const char *tmp = getenv("TMPDIR"); + char path[1024]; + snprintf(path, sizeof path, "%s/cactus_baseAlignerTestXXXXXX", (tmp && *tmp) ? tmp : "/tmp"); + int fd = mkstemp(path); + CuAssertTrue(testCase, fd >= 0); + FILE *f = fdopen(fd, "w"); + CuAssertPtrNotNull(testCase, f); + fprintf(f, "\n", bar_attrs); + fclose(f); + CactusParams *params = cactusParams_load(path); + BaseAligner engine = baseAligner_constructFromCactusParams(params); + cactusParams_destruct(params); + remove(path); + return engine; +} + +void test_baseAligner_selection(CuTest *testCase) { + // No baseAligner at all: the old boolean decides, as every pre-existing config expects. + CuAssertIntEquals(testCase, BASE_ALIGNER_ABPOA, engine_for_config(testCase, "partialOrderAlignment=\"1\"")); + CuAssertIntEquals(testCase, BASE_ALIGNER_PECAN, engine_for_config(testCase, "partialOrderAlignment=\"0\"")); + + // baseAligner present: it decides, and all three values resolve. + CuAssertIntEquals(testCase, BASE_ALIGNER_PECAN, + engine_for_config(testCase, "partialOrderAlignment=\"1\" baseAligner=\"pecan\"")); + CuAssertIntEquals(testCase, BASE_ALIGNER_ABPOA, + engine_for_config(testCase, "partialOrderAlignment=\"1\" baseAligner=\"abpoa\"")); + CuAssertIntEquals(testCase, BASE_ALIGNER_MINIPOA, + engine_for_config(testCase, "partialOrderAlignment=\"1\" baseAligner=\"minipoa\"")); + + // Only baseAligner, no legacy boolean. This is what the warning below tells users to write, + // and reading partialOrderAlignment unguarded used to st_errAbort on exactly this config. + CuAssertIntEquals(testCase, BASE_ALIGNER_MINIPOA, engine_for_config(testCase, "baseAligner=\"minipoa\"")); + CuAssertIntEquals(testCase, BASE_ALIGNER_ABPOA, engine_for_config(testCase, "baseAligner=\"abpoa\"")); + CuAssertIntEquals(testCase, BASE_ALIGNER_PECAN, engine_for_config(testCase, "baseAligner=\"pecan\"")); + + // Disagreeing: baseAligner wins (and the C code logs about it). This is the trap -- the + // shipped config now carries baseAligner="abpoa", so someone who copies it and sets + // partialOrderAlignment="0" the documented way would otherwise silently keep abpoa. + CuAssertIntEquals(testCase, BASE_ALIGNER_MINIPOA, + engine_for_config(testCase, "partialOrderAlignment=\"0\" baseAligner=\"minipoa\"")); + CuAssertIntEquals(testCase, BASE_ALIGNER_PECAN, + engine_for_config(testCase, "partialOrderAlignment=\"1\" baseAligner=\"pecan\"")); +} + +/* + * Parse the config cactus actually ships and read every attribute the C side depends on. + * + * Nothing else does this: the build does not look at the XML, and the other tests here write + * their own. A stray "--" inside an XML comment got all the way past a clean build and a green + * suite before an alignment run caught it, which is too late and too indirect. + */ +void test_shipped_config_is_loadable(CuTest *testCase) { + const char *path = "src/cactus/cactus_progressive_config.xml"; + FILE *f = fopen(path, "r"); + if (f == NULL) { + // run from somewhere other than the repo root; nothing to check + return; + } + fclose(f); + CactusParams *params = cactusParams_load((char *)path); + CuAssertPtrNotNull(testCase, params); + + CuAssertIntEquals(testCase, BASE_ALIGNER_ABPOA, baseAligner_constructFromCactusParams(params)); + + // every attribute poaBarAligner.c reads; the getters st_errAbort on a missing one + PoaParameters *abpoa = poaParameters_constructFromCactusParams(params, BASE_ALIGNER_ABPOA); + CuAssertPtrNotNull(testCase, abpoa); + poaParameters_destruct(abpoa); +#ifdef HAVE_MINIPOA + PoaParameters *minipoa = poaParameters_constructFromCactusParams(params, BASE_ALIGNER_MINIPOA); + CuAssertPtrNotNull(testCase, minipoa); + CuAssertTrue(testCase, minipoa->gapOpen > 0); + CuAssertTrue(testCase, minipoa->gapExt > 0); + CuAssertTrue(testCase, minipoa->mat[0] > 0); // A/A must be a match + CuAssertTrue(testCase, minipoa->mat[1] < 0); // A/C must be a mismatch + CuAssertTrue(testCase, minipoa->mat[24] > minipoa->mat[4]); // minipoa requires N/N > N/other + poaParameters_destruct(minipoa); +#endif + cactusParams_destruct(params); } CuSuite* poaBarAlignerTestSuite(void) { CuSuite* suite = CuSuiteNew(); + SUITE_ADD_TEST(suite, test_shipped_config_is_loadable); + SUITE_ADD_TEST(suite, test_baseAligner_selection); SUITE_ADD_TEST(suite, test_make_partial_order_alignment); SUITE_ADD_TEST(suite, test_make_consistent_partial_order_alignments_two_ends); SUITE_ADD_TEST(suite, test_make_flower_alignment_poa); diff --git a/doc/progressive.md b/doc/progressive.md index c3d71edd2..e6b8ff693 100644 --- a/doc/progressive.md +++ b/doc/progressive.md @@ -416,6 +416,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 `--maxMemory/--doubleMem` should be enough to work around this type of issue. * `--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 unless the memory estimate is more than the job can be given (the system memory on a single machine, or `--maxMemory`), in which case the estimate is scaled down and the pages are not kept. 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 much less memory per window, which is why its default window is 100000 rather than 10000; the memory *estimate* for `cactus_consolidated` is still the one fitted to abpoa, so under minipoa it over-requests rather than under-requests. 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`. On a cluster with partitions and/or time limits, make sure to use diff --git a/include.mk b/include.mk index 583bb617c..530212c3d 100644 --- a/include.mk +++ b/include.mk @@ -97,6 +97,25 @@ endif ifeq ($(shell arch || true), arm64) arm=1 endif + +# Control variable for minipoa, the optional third base aligner in bar. +# +# Off means the submodule is not built, -lminipoa is not linked, and selecting +# aborts with a message saying so. abpoa and pecan are unaffected, +# so a tree that cannot build minipoa still builds cactus -- which matters because minipoa is a +# young single-maintainer library and nobody who never selects it should be held up by it. +# +# Off on ARM until minipoa's simde/NEON path has been proven there: upstream has no ARM build and +# no CI that would catch a regression in one. +# +# This has to sit BELOW the arm detection above. make evaluates conditionals as it reads the +# file, so testing `ifdef arm` before those three shell probes assign it left minipoa on for every +# machine that did not pass arm=1 explicitly -- i.e. exactly the ARM machines the detection exists +# to catch. +minipoa = on +ifdef arm + minipoa = off +endif # CPU baseline for the code we generate. A plain "make" is portable. Opt into a # machine-specific build, for binaries that will only ever run on the machine that compiled # them: @@ -265,5 +284,13 @@ endif jemallocSubLibs = $(if ${jemallocLib},-L$(abspath ${LIBDIR}) ${jemallocLib}) # note: the CACTUS_STATIC_LINK_FLAGS below can generally be empty -- it's used by the static builder script only -LDLIBS += ${cactusLibs} ${sonLibLibs} ${LIBS} -L${rootPath}/lib -Wl,-rpath,${rootPath}/lib -labpoa ${jemallocLib} -lz -lbz2 -lpthread -lm -lstdc++ -lm -lxml2 ${CACTUS_STATIC_LINK_FLAGS} +# -lminipoa sits next to -labpoa, i.e. ahead of -lz: ld resolves archives left to right and +# libminipoa.a has undefined references into zlib (its kseq reader and klib's err_* wrappers are +# linked in regardless of whether we call them). Putting it after -lz fails the static release +# link, and fails the ordinary one too wherever --as-needed is the default. +ifeq ($(minipoa),on) + minipoaLib = -lminipoa + CFLAGS += -DHAVE_MINIPOA +endif +LDLIBS += ${cactusLibs} ${sonLibLibs} ${LIBS} -L${rootPath}/lib -Wl,-rpath,${rootPath}/lib -labpoa ${minipoaLib} ${jemallocLib} -lz -lbz2 -lpthread -lm -lstdc++ -lm -lxml2 ${CACTUS_STATIC_LINK_FLAGS} LIBDEPENDS = ${sonLibDir}/sonLib.a ${sonLibDir}/cuTest.a ${jemallocDepends} diff --git a/pipeline/cactus_consolidated.c b/pipeline/cactus_consolidated.c index 82b8d029c..5f9f82069 100644 --- a/pipeline/cactus_consolidated.c +++ b/pipeline/cactus_consolidated.c @@ -624,7 +624,9 @@ int main(int argc, char *argv[]) { stList_length(leafFlowers), time(NULL) - startTime); bar(leafFlowers, params, cactusDisk, NULL); - int64_t usePoa = cactusParams_get_int(params, 2, "bar", "partialOrderAlignment"); + // optional: a config may carry only , with no legacy boolean at all + int64_t usePoa = cactusParams_has(params, 2, "bar", "partialOrderAlignment") + ? cactusParams_get_int(params, 2, "bar", "partialOrderAlignment") : 1; st_logInfo("Ran cactus bar (use poa:%i), %" PRIi64 " seconds have elapsed\n", (int)usePoa, time(NULL) - startTime); stList_destruct(leafFlowers); diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index 3dd11d336..669eba561 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -234,8 +234,16 @@ - + + @@ -243,6 +251,7 @@ + + + + + + + + + + + + + + + + + + diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index d63c155da..e07add42f 100644 --- a/src/cactus/paf/last_scoring.py +++ b/src/cactus/paf/last_scoring.py @@ -97,7 +97,10 @@ def apply_scores_to_config(score_dict, config_xml): # that's applied to the trained gap parameters long_gap_open_factor = int(poa_node.attrib['partialOrderAlignmentTrainedGapOpen2Factor']) long_gap_extend_factor = int(poa_node.attrib['partialOrderAlignmentTrainedGapExtension2Factor']) - apply_long_gap(score_dict, long_gap_extend_factor, long_gap_extend_factor) + # the open factor was being passed the extend factor, so TrainedGapOpen2Factor was read and + # then discarded. Both default to 3, so the shipped config is unaffected; setting them to + # different values silently did the wrong thing. + apply_long_gap(score_dict, long_gap_open_factor, long_gap_extend_factor) poa_node.attrib['partialOrderAlignmentGapOpenPenalty1'] = str(score_dict['GAP-OPEN']) poa_node.attrib['partialOrderAlignmentGapExtensionPenalty1'] = str(score_dict['GAP-EXTEND']) @@ -130,6 +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_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']) + RealtimeLogger.info("Overriding minipoa scores with trained values: GapOpen {}; GapExtend {} (single affine, as trained)".format( + minipoa_node.attrib['minipoaGapOpenPenalty'], + minipoa_node.attrib['minipoaGapExtensionPenalty'])) + RealtimeLogger.info("Overriding abPOA scores with trained values: GapOpen {}; GapExtend {}; GapOpen2 {}; GapExtend2 {}; SubMatrix {}".format( poa_node.attrib['partialOrderAlignmentGapOpenPenalty1'], poa_node.attrib['partialOrderAlignmentGapExtensionPenalty1'], diff --git a/src/cactus/pipeline/cactus_workflow.py b/src/cactus/pipeline/cactus_workflow.py index 3e4827c22..3d1d99c50 100644 --- a/src/cactus/pipeline/cactus_workflow.py +++ b/src/cactus/pipeline/cactus_workflow.py @@ -83,8 +83,22 @@ def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_ma # at once. That ratio is an exponent of 0.43, not 2. Nothing is subtracted for # partialOrderAlignmentMaskFilter even though it matters more, because every alignment in # the fit ran with it disabled -- so enabling it can only make the estimate conservative. - poa_node = findRequiredNode(config_node, 'bar').find('poa') - poa_window = getOptionalAttrib(poa_node, 'partialOrderAlignmentWindow', typeFn=int, default=10000) if poa_node is not None else 10000 + # + # The window comes from whichever base aligner is selected. The exponent above was fitted to + # abPOA runs, so applying it to minipoa's much larger window over-requests -- by roughly 7x at + # a 100kb window. That is deliberate until someone fits an exponent from real minipoa runs and + # records it: over-requesting wastes cluster share, under-requesting gets the largest ancestors + # OOM-killed, and those are the very jobs minipoa is for. + bar_node = findRequiredNode(config_node, 'bar') + base_aligner = getOptionalAttrib(bar_node, 'baseAligner', typeFn=str, default=None) + if base_aligner is None: + base_aligner = 'abpoa' if getOptionalAttrib(bar_node, 'partialOrderAlignment', typeFn=bool, default=True) else 'pecan' + if base_aligner == 'minipoa': + engine_node = bar_node.find('minipoa') + poa_window = getOptionalAttrib(engine_node, 'minipoaWindow', typeFn=int, default=100000) if engine_node is not None else 100000 + else: + poa_node = bar_node.find('poa') + poa_window = getOptionalAttrib(poa_node, 'partialOrderAlignmentWindow', typeFn=int, default=10000) if poa_node is not None else 10000 window_exp = getOptionalAttrib(cons_node, 'memory_poa_window_exponent', typeFn=float, default=0.43) if poa_window > 0 and poa_window != 10000: mem = int(mem * (poa_window / 10000.0) ** window_exp) @@ -111,9 +125,9 @@ def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_ma bytes2human(mem), bytes2human(scaled_mem))) mem = scaled_mem - # abPOA needs a table even for tiny alignments; apply the floor last so neither the window - # nor the core scaling can push a small job below it - if getOptionalAttrib(findRequiredNode(config_node, 'bar'), 'partialOrderAlignment', typeFn=bool, default=True): + # A POA aligner needs a table even for tiny alignments; apply the floor last so neither the + # window nor the core scaling can push a small job below it + if base_aligner != 'pecan': mem = max(mem, int(4e9)) RealtimeLogger.info('Estimating cactus_consolidated({}) memory as {} from {} sequences with total-sequence-size {} and paf-size {} using configuration settings'.format(chrom_name if chrom_name else ancestor_event, bytes2human(mem), len(seq_id_map), bytes2human(total_sequence_size), paf_id.size)) diff --git a/src/cactus/setup/cactus_align.py b/src/cactus/setup/cactus_align.py index d7b5b4849..c65c5efe4 100644 --- a/src/cactus/setup/cactus_align.py +++ b/src/cactus/setup/cactus_align.py @@ -71,7 +71,7 @@ def main(): parser.add_argument("--singleCopySpecies", type=str, help="Filter out all self-alignments in given species") parser.add_argument("--barMaskFilter", type=int, default=None, - help="BAR's POA aligner will ignore softmasked regions greater than this length. (overrides partialOrderAlignmentMaskFilter in config)") + help="BAR's POA aligner will ignore softmasked regions greater than this length. (overrides the mask-filter attribute of whichever base aligner is selected)") parser.add_argument("--pafMaskFilter", type=int, default=None, help="softmasked (query) regions greather than this length will be removed from the input PAF before it is processed") # note: when changing this, make sure to keep version in cactus-pangenome consistent @@ -344,10 +344,17 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): cafNode = findRequiredNode(config_node, "caf") barNode = findRequiredNode(config_node, "bar") poaNode = findRequiredNode(barNode, "poa") + # The engine-specific nodes are named differently, so an override has to know which engine is + # selected -- writing only to would make --barMaskFilter a silent no-op under minipoa, + # and that is the lever that decides whether a repeat-rich run finishes at all. + engineNode, maskFilterAttr, seedingAttr = poaNode, "partialOrderAlignmentMaskFilter", "partialOrderAlignmentDisableSeeding" + if getOptionalAttrib(barNode, "baseAligner", typeFn=str, default="abpoa") == "minipoa": + engineNode = findRequiredNode(barNode, "minipoa") + maskFilterAttr, seedingAttr = "minipoaMaskFilter", "minipoaDisableSeeding" if options.singleCopySpecies: cafNode.attrib["alignmentFilter"] = "singleCopyEvent:{}".format(options.singleCopySpecies) if options.barMaskFilter: - poaNode.attrib["partialOrderAlignmentMaskFilter"] = str(options.barMaskFilter) + engineNode.attrib[maskFilterAttr] = str(options.barMaskFilter) if options.maxLen is None and options.pangenome: # consistent behaviour with cactus-pangenome @@ -368,7 +375,7 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): # turn down minimum block degree to get a fat ancestor barNode.attrib["minimumBlockDegree"] = "1" # turn off POA seeding - poaNode.attrib["partialOrderAlignmentDisableSeeding"] = "1" + engineNode.attrib[seedingAttr] = "1" # import the PAF alignments paf_id = toil.importFile(makeURL(options.pafFile)) diff --git a/submodules/minipoa b/submodules/minipoa new file mode 160000 index 000000000..fd134725e --- /dev/null +++ b/submodules/minipoa @@ -0,0 +1 @@ +Subproject commit fd134725ee16acb5de152f096a0a84e7664d51cc diff --git a/test/evolverTest.py b/test/evolverTest.py index 0d4ee1b39..2922aa627 100644 --- a/test/evolverTest.py +++ b/test/evolverTest.py @@ -1646,6 +1646,47 @@ def testEvolverPOALocal(self): # check the output self._check_maf_accuracy(self._out_hal("local"), delta=(0.0025,0.0075), dataset='primates') + def _write_base_aligner_config(self, base_aligner, star=False): + """ Copy the shipped config, point at the given base aligner, and return the path. + Returns the path so a test can hand it to --configFile. + """ + config_path = 'src/cactus/cactus_progressive_config.xml' + xml_root = ET.parse(config_path).getroot() + bar_elem = xml_root.find("bar") + bar_elem.attrib["baseAligner"] = base_aligner + # keep the legacy boolean consistent, so nothing still reading it disagrees + bar_elem.attrib["partialOrderAlignment"] = "0" if base_aligner == "pecan" else "1" + if star: + # force cactus to accept multifurcation in tree + xml_root.find("multi_cactus").find("decomposition").attrib["allow_multifurcations"] = "1" + out_path = os.path.join(self.tempDir, "config.{}.xml".format(base_aligner)) + with open(out_path, 'w') as out_file: + xmlString = ET.tostring(xml_root, encoding='unicode') + xmlString = minidom.parseString(xmlString).toprettyxml() + out_file.write(xmlString) + return out_path + + def testEvolverMinipoaLocal(self): + """ Same shape as testEvolverPOALocal, but with minipoa as the base aligner instead of + abpoa. Same dataset and same tolerance, so the two are directly comparable: this is the + sharpest accuracy instrument in the suite. + """ + minipoa_config_path = self._write_base_aligner_config("minipoa", star=True) + name = "local" + self._run_evolver_primates_star(name, configFile = minipoa_config_path) + self._check_maf_accuracy(self._out_hal("local"), delta=(0.0025,0.0075), dataset='primates') + + def testEvolverMinipoaMammalsLocal(self): + """ The mammals head-to-head: testEvolverLocal runs the same data through abpoa (the + shipped default), so running it again with baseAligner="minipoa" and the same tolerance + compares the two engines and nothing else. + """ + minipoa_config_path = self._write_base_aligner_config("minipoa") + name = "local" + self._run_evolver(name, configFile = minipoa_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 78bfd111b39ec081af7473b0e4d7c98923c2a0c2 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 16 Sep 2026 08:07:05 -0400 Subject: [PATCH 2/4] bar: group poaBarAligner.c by backend Everything lived in one undifferentiated file, with abPOA's parameter construction at the top, minipoa's five hundred lines further down, and abPOA's actual alignment call inlined in the dispatcher while minipoa's sat in its own function. Adding a third engine meant finding the seams again. No behaviour change. The abPOA call is extracted into run_abpoa_window so the two backends are symmetric, run_poa_window becomes a four-line switch over the engine, and the file is grouped into labelled sections: the shared alphabet and Msa handling, the shared trimming and stitching, the shared dump helper, one section per backend, the selection and dispatch, then the shared windowing, adjacency and block-extraction code. The engine-specific code is now about 270 lines for abPOA and 195 for minipoa against roughly 1250 shared, which is the honest ratio -- most of BAR's base-alignment logic is engine-neutral and always was. Verified as a no-op: no line was added or removed except section banners, and both engines produce byte-for-byte identical MSAs across 229 real BAR windows before and after. Co-Authored-By: Claude Opus 5 (1M context) --- bar/impl/poaBarAligner.c | 1062 ++++++++++++++++++++------------------ 1 file changed, 564 insertions(+), 498 deletions(-) diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index ea68835b6..3864998ae 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -16,6 +16,10 @@ #include #include +/* ============================================================================ */ +/* Shared: the CACTUS_BAR_DUMP_DIR switch */ +/* ============================================================================ */ + /* * Set CACTUS_BAR_DUMP_DIR in the environment and every window handed to the base aligner is * written there as a FASTA, alongside the substitution matrix and a command line that replays it. @@ -47,95 +51,11 @@ void bar_dump_dir_init(void) { #include #endif -abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { - abpoa_para_t *abpt = abpoa_init_para(); - - // output options - abpt->out_msa = 1; // generate Row-Column multiple sequence alignment(RC-MSA), set 0 to disable - abpt->out_cons = 0; // generate consensus sequence, set 0 to disable - - // alignment mode. 0:global alignment, 1:local, 2:extension - // only global works - abpt->align_mode = ABPOA_GLOBAL_MODE; - - // banding parameters - abpt->wb = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentBandConstant"); - abpt->wf = cactusParams_get_float(params, 3, "bar", "poa", "partialOrderAlignmentBandFraction"); - - // gap scoring model - abpt->gap_open1 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapOpenPenalty1"); - abpt->gap_ext1 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapExtensionPenalty1"); - abpt->gap_open2 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapOpenPenalty2"); - abpt->gap_ext2 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapExtensionPenalty2"); - - // seeding paramters - abpt->disable_seeding = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentDisableSeeding"); - assert(abpt->disable_seeding == 0 || abpt->disable_seeding == 1); - abpt->k = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMinimizerK"); - abpt->w = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMinimizerW"); - abpt->min_w = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMinimizerMinW"); - - // progressive toggle - abpt->progressive_poa = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentProgressiveMode"); - - // generate the substitution matrix - abpt->use_score_matrix = 0; - abpoa_post_set_para(abpt); - - // optionally override the substitution matrix - char *submat_string = cactusParams_get_string(params, 3, "bar", "poa", "partialOrderAlignmentSubMatrix"); - if (submat_string && strlen(submat_string) > 0) { - // Note, this will be used to explicitly override abpoa's subsitution matrix just before aligning - abpt->use_score_matrix = 1; - assert(abpt->m == 5); - int count = 0; - for (char* val = strtok(submat_string, " "); val != NULL; val = strtok(NULL, " ")) { - abpt->mat[count++] = atoi(val); - } - assert(count == 25); - int i; abpt->min_mis = 0, abpt->max_mat = 0; - for (i = 0; i < abpt->m * abpt->m; ++i) { - if (abpt->mat[i] > abpt->max_mat) - abpt->max_mat = abpt->mat[i]; - if (-abpt->mat[i] > abpt->min_mis) - abpt->min_mis = -abpt->mat[i]; - } - } - free(submat_string); - return abpt; -} +/* ============================================================================ */ +/* Shared: alphabet, and the Msa the backends fill in */ +/* ============================================================================ */ -// It turns out abpoa can write to these, so we make a quick copy before using -static abpoa_para_t *copy_abpoa_params(abpoa_para_t *abpt) { - abpoa_para_t *abpt_cpy = abpoa_init_para(); - abpt_cpy->out_msa = 1; - abpt_cpy->out_cons = 0; - abpt_cpy->align_mode = abpt->align_mode; - abpt_cpy->wb = abpt->wb; - abpt_cpy->wf = abpt->wf; - abpt_cpy->match = abpt->match; - abpt_cpy->mismatch = abpt->mismatch; - abpt_cpy->gap_mode = abpt->gap_mode; - abpt_cpy->gap_open1 = abpt->gap_open1; - abpt_cpy->gap_ext1 = abpt->gap_ext1; - abpt_cpy->gap_open2 = abpt->gap_open2; - abpt_cpy->gap_ext2 = abpt->gap_ext2; - abpt_cpy->disable_seeding = abpt->disable_seeding; - abpt_cpy->k = abpt->k; - abpt_cpy->w = abpt->w; - abpt_cpy->min_w = abpt->min_w; - abpt_cpy->progressive_poa = abpt->progressive_poa; - abpt_cpy->use_score_matrix = 0; - abpoa_post_set_para(abpt_cpy); - abpt_cpy->use_score_matrix = abpt->use_score_matrix; - if (abpt->use_score_matrix == 1) { - memcpy(abpt_cpy->mat, abpt->mat, abpt->m * abpt->m * sizeof(int)); - } - abpt_cpy->max_mat = abpt->max_mat; - abpt_cpy->min_mis = abpt->min_mis; - return abpt_cpy; -} // char <--> uint8_t conversion copied over from abPOA example // AaCcGgTtNn ==> 0,1,2,3,4 @@ -192,277 +112,110 @@ static inline uint8_t msa_to_rc(uint8_t n) { return rc_table[n]; } -/* - * Write the window as FASTA plus its 5x5 matrix. The matrix file format is the one both abpoa -t - * and minipoa -m read, so either aligner can be pointed straight at it. - */ -static void dump_window_fasta_and_matrix(Msa *msa, uint8_t **bseqs, const int *mat, - const char *input_path, const char *matrix_path) { - FILE *mat_file = fopen(matrix_path, "w"); - if (mat_file != NULL) { - fprintf(mat_file, "\tA\tC\tG\tT\tN\n"); - for (size_t i = 0; i < 5; ++i) { - fprintf(mat_file, "%c", "ACGTN"[i]); - for (size_t j = 0; j < 5; ++j) { - fprintf(mat_file, "\t%d", mat[i * 5 + j]); - } - fprintf(mat_file, "\n"); + + + + +void msa_destruct(Msa *msa) { + for(int64_t i=0; iseq_no; i++) { + if (msa->seqs != NULL) { + free(msa->seqs[i]); } - fclose(mat_file); - } - FILE *fa_file = fopen(input_path, "w"); - if (fa_file == NULL) { - return; + free(msa->msa_seq[i]); } - for (int64_t i = 0; i < msa->seq_no; ++i) { - fprintf(fa_file, ">%" PRIi64 "\n", i); - for (int64_t j = 0; j < msa->seq_lens[i]; ++j) { - fputc(msa_to_base(bseqs[i][j]), fa_file); + free(msa->seqs); + free(msa->msa_seq); + free(msa->seq_lens); + free(msa); +} + +void msa_print(Msa *msa, FILE *f) { + fprintf(f, "MSA. Seq no: %i column no: %i \n", (int)msa->seq_no, (int)msa->column_no); + for(int64_t i=0; iseq_no; i++) { + fprintf(f, "Row:%i [len=%i]\t", (int)i, (int)msa->seq_lens[i]); + for(int64_t j=0; jcolumn_no; j++) { + fprintf(f, "%c", msa_to_base(msa->msa_seq[i][j])); } - fputc('\n', fa_file); + fprintf(f, "\n"); } - fclose(fa_file); + fprintf(f, "\n"); } -/* - * A command line that replays this window through minipoa. Gap penalties are negated: minipoa - * maximises, so penalties are negative there, while cactus and abpoa carry them positive. - */ -static char *dump_minipoa_input(Msa *msa, PoaParameters *pp, uint8_t **bseqs, char *input_path, - char *matrix_path, char *command_path, char *output_path) { - dump_window_fasta_and_matrix(msa, bseqs, pp->mat, input_path, matrix_path); +/* ============================================================================ */ +/* Shared: column scores, trimming and window stitching */ +/* ============================================================================ */ - int f = pp->bandFraction > 0.0 ? (int)(1.0 / pp->bandFraction + 0.5) : 0; - 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->seeding) { - char kw_opts[128]; - sprintf(kw_opts, " -S -k %d -w %d", pp->minimizerK, pp->minimizerW); - strcat(command, kw_opts); - if (pp->anchorWindow > 0) { - sprintf(kw_opts, " -W %d", pp->anchorWindow); - strcat(command, kw_opts); + +/** + * flip msa to its reverse complement (for trimming purposees) + */ +static void flip_msa_seq(Msa* msa) { + if (msa != NULL) { + int64_t middle = msa->column_no / 2; + bool odd = msa->column_no % 2 == 1; + for (int64_t i = 0; i < msa->seq_no; ++i) { + for (int64_t j = 0; j < middle; ++j) { + uint8_t buf = msa->msa_seq[i][j]; + msa->msa_seq[i][j] = msa_to_rc(msa->msa_seq[i][msa->column_no - 1 - j]); + msa->msa_seq[i][msa->column_no - 1 - j] = msa_to_rc(buf); + } + if (odd) { + msa->msa_seq[i][middle] = msa_to_rc(msa->msa_seq[i][middle]); + } } } - // -p and -B are on in the shipped config, and progressive ordering in particular changes the - // alignment, so a replay that omitted them would not reproduce the window it is meant to - // explain. - if (pp->progressive) { - strcat(command, " -p"); - } - if (pp->adaptiveBand) { - strcat(command, " -B"); - } - strcat(command, " > "); - strcat(command, output_path); +} - FILE *cmd_file = fopen(command_path, "w"); - if (cmd_file != NULL) { - fprintf(cmd_file, "%s\n", command); - fclose(cmd_file); +/** + * Returns an array of floats, one for each corresponding column in the MSA. Each float + * is the score of the column in the alignment. + */ +static float *make_column_scores(Msa *msa) { + float *column_scores = st_calloc(msa->column_no, sizeof(float)); + for(int64_t i=0; icolumn_no; i++) { + // Score is simply max(number of aligned bases in the column - 1, 0) + for(int64_t j=0; jseq_no; j++) { + if(msa_to_base(msa->msa_seq[j][i]) != '-') { + column_scores[i]++; + } + } + if(column_scores[i] >= 1.0) { + column_scores[i]--; + } + assert(column_scores[i] >= 0.0); } - return command; + return column_scores; } -// dump the abpoa input to files, and return a command line for running abpoa on them -char* dump_abpoa_input(Msa* msa, abpoa_para_t* abpt, uint8_t **bseqs, char* abpoa_input_path, char* abpoa_matrix_path, - char* abpoa_command_path, char* abpoa_output_path) { - // dump the abpoa input sequences to a FASTA file - FILE* dump_file = fopen(abpoa_input_path, "w"); - for (int64_t i = 0; i < msa->seq_no; ++i) { - int64_t seq_len = msa->seq_lens[i]; - char* buffer = (char*)malloc((seq_len + 1) * sizeof(char)); - for (int64_t j = 0; j < seq_len; ++j) { - buffer[j] = msa_to_base(bseqs[i][j]); +/** + * Fills in cu_column_scores with the cumulative sum of column scores, from left-to-right, of columns + * containing a non-gap character in the given "row". + */ +static void sum_column_scores(int64_t row, Msa *msa, float *column_scores, float *cu_column_scores) { + float cu_score = 0.0; // The cumulative sum of column scores containing bases for the given row + int64_t j=0; // The index in the DNA string for the given row + for(int64_t i=0; icolumn_no; i++) { + if(msa_to_base(msa->msa_seq[row][i]) != '-') { + cu_score += column_scores[i]; + cu_column_scores[j++] = cu_score; } - buffer[msa->seq_lens[i]] = '\0'; - fprintf(dump_file, ">%ld\n%s\n", i, buffer); - free(buffer); } - fclose(dump_file); + assert(msa->seq_lens[row] == j); // We should cover all the bases in the DNA sequence +} - // dump the abpoa input matrix to file - FILE* mat_file = fopen(abpoa_matrix_path, "w"); - fprintf(mat_file, "\tA\tC\tG\tT\tN\n"); - for (size_t i = 0; i < 5; ++i) { - fprintf(mat_file, "%c", "ACGTN"[i]); - for (size_t j = 0; j < 5; ++j) { - fprintf(mat_file, "\t%d", abpt->mat[i * 5 + j]); - } - fprintf(mat_file, "\n"); - } - fclose(mat_file); - - // make a command line - char* abpoa_command = st_malloc(4096 * sizeof(char)); - sprintf(abpoa_command, "abpoa %s -O %d,%d -E %d,%d -b %d -f %lf -t %s -r 1 -m 0", - abpoa_input_path, - abpt->gap_open1, - abpt->gap_open2, - abpt->gap_ext1, - abpt->gap_ext2, - abpt->wb, - abpt->wf, - abpoa_matrix_path); - if (!abpt->disable_seeding) { - strcat(abpoa_command, " -S"); - char kw_opts[128]; - sprintf(kw_opts, " -k %d -w %d -n %d", abpt->k, abpt->w, abpt->min_w); - strcat(abpoa_command, kw_opts); - } - if (abpt->progressive_poa) { - strcat(abpoa_command, " -p"); - } - strcat(abpoa_command, " > "); - strcat(abpoa_command, abpoa_output_path); - - // dump the command line - FILE* cmd_file = fopen(abpoa_command_path, "w"); - fprintf(cmd_file, "%s\n", abpoa_command); - fclose(cmd_file); - - return abpoa_command; -} - -#ifdef CACTUS_ABPOA_FROM_COMMAND_LINE -void abpoa_msa_from_command_line(char* abpoa_command_line, char* abpoa_output_path, uint8_t*** msa_seq, int* col_no) { - // run abpoa - st_system(abpoa_command_line); - - // read the result (ascii alignment) back into memory - size_t n_rows = 0; - size_t n_cols = 0; - FILE* msa_file = fopen(abpoa_output_path, "r"); - int64_t buf_size = 500000; - char* buf = st_malloc(buf_size * sizeof(char)); - - while (benLine(&buf, &buf_size, msa_file) != -1) { - if (strlen(buf) && buf[0] != '>') { - ++n_rows; - } - } - - *msa_seq = st_malloc(n_rows * sizeof(uint8_t*)); - fclose(msa_file); - msa_file = fopen(abpoa_output_path, "r"); - n_rows = 0; - while (benLine(&buf, &buf_size, msa_file) != -1) { - if (strlen(buf) && buf[0] != '>') { - if (n_cols == 0) { - n_cols = strlen(buf); - } else { - assert(n_cols == strlen(buf)); - } - (*msa_seq)[n_rows] = st_malloc(n_cols * sizeof(uint8_t)); - for (size_t i = 0; i < n_cols; ++i) { - (*msa_seq)[n_rows][i] = msa_to_byte(buf[i]); - } - ++n_rows; - } - } - fclose(msa_file); - *col_no = (int)n_cols; - - free(buf); -} -#endif - -void msa_destruct(Msa *msa) { - for(int64_t i=0; iseq_no; i++) { - if (msa->seqs != NULL) { - free(msa->seqs[i]); - } - free(msa->msa_seq[i]); - } - free(msa->seqs); - free(msa->msa_seq); - free(msa->seq_lens); - free(msa); -} - -void msa_print(Msa *msa, FILE *f) { - fprintf(f, "MSA. Seq no: %i column no: %i \n", (int)msa->seq_no, (int)msa->column_no); - for(int64_t i=0; iseq_no; i++) { - fprintf(f, "Row:%i [len=%i]\t", (int)i, (int)msa->seq_lens[i]); - for(int64_t j=0; jcolumn_no; j++) { - fprintf(f, "%c", msa_to_base(msa->msa_seq[i][j])); - } - fprintf(f, "\n"); - } - fprintf(f, "\n"); -} - -/** - * flip msa to its reverse complement (for trimming purposees) - */ -static void flip_msa_seq(Msa* msa) { - if (msa != NULL) { - int64_t middle = msa->column_no / 2; - bool odd = msa->column_no % 2 == 1; - for (int64_t i = 0; i < msa->seq_no; ++i) { - for (int64_t j = 0; j < middle; ++j) { - uint8_t buf = msa->msa_seq[i][j]; - msa->msa_seq[i][j] = msa_to_rc(msa->msa_seq[i][msa->column_no - 1 - j]); - msa->msa_seq[i][msa->column_no - 1 - j] = msa_to_rc(buf); - } - if (odd) { - msa->msa_seq[i][middle] = msa_to_rc(msa->msa_seq[i][middle]); - } - } - } -} - -/** - * Returns an array of floats, one for each corresponding column in the MSA. Each float - * is the score of the column in the alignment. - */ -static float *make_column_scores(Msa *msa) { - float *column_scores = st_calloc(msa->column_no, sizeof(float)); - for(int64_t i=0; icolumn_no; i++) { - // Score is simply max(number of aligned bases in the column - 1, 0) - for(int64_t j=0; jseq_no; j++) { - if(msa_to_base(msa->msa_seq[j][i]) != '-') { - column_scores[i]++; - } - } - if(column_scores[i] >= 1.0) { - column_scores[i]--; - } - assert(column_scores[i] >= 0.0); - } - return column_scores; -} - -/** - * Fills in cu_column_scores with the cumulative sum of column scores, from left-to-right, of columns - * containing a non-gap character in the given "row". - */ -static void sum_column_scores(int64_t row, Msa *msa, float *column_scores, float *cu_column_scores) { - float cu_score = 0.0; // The cumulative sum of column scores containing bases for the given row - int64_t j=0; // The index in the DNA string for the given row - for(int64_t i=0; icolumn_no; i++) { - if(msa_to_base(msa->msa_seq[row][i]) != '-') { - cu_score += column_scores[i]; - cu_column_scores[j++] = cu_score; - } - } - assert(msa->seq_lens[row] == j); // We should cover all the bases in the DNA sequence -} - -/** - * Removes the suffix of the given row from the MSA and updates the column scores. suffix_start is the beginning - * suffix to remove. - */ -static void trim_msa_suffix(Msa *msa, float *column_scores, int64_t row, int64_t suffix_start) { - int64_t seq_index = 0; - for(int64_t i=0; icolumn_no; i++) { - if(msa_to_base(msa->msa_seq[row][i]) != '-') { - if(seq_index++ >= suffix_start) { - msa->msa_seq[row][i] = msa_to_byte('-'); - column_scores[i] = column_scores[i] > 1 ? column_scores[i]-1 : 0; - assert(column_scores[i] >= 0.0); - } +/** + * Removes the suffix of the given row from the MSA and updates the column scores. suffix_start is the beginning + * suffix to remove. + */ +static void trim_msa_suffix(Msa *msa, float *column_scores, int64_t row, int64_t suffix_start) { + int64_t seq_index = 0; + for(int64_t i=0; icolumn_no; i++) { + if(msa_to_base(msa->msa_seq[row][i]) != '-') { + if(seq_index++ >= suffix_start) { + msa->msa_seq[row][i] = msa_to_byte('-'); + column_scores[i] = column_scores[i] > 1 ? column_scores[i]-1 : 0; + assert(column_scores[i] >= 0.0); + } } } } @@ -557,60 +310,324 @@ static void msa_fix_trimmed(Msa* msa) { msa->column_no -= empty_columns; } + +/* ============================================================================ */ +/* Window dumping: shared between the backends */ +/* ============================================================================ */ + /* - * picks the engine. It is read with cactusParams_has so that a config - * predating the attribute -- including any a user has saved -- still selects what it used to via - * the older boolean. + * Write the window as FASTA plus its 5x5 matrix. The matrix file format is the one both abpoa -t + * and minipoa -m read, so either aligner can be pointed straight at it. */ -BaseAligner baseAligner_constructFromCactusParams(CactusParams *params) { - /* - * Both attributes are optional, in both directions. An old config has only - * partialOrderAlignment; a config written from now on may reasonably have only baseAligner -- - * including one produced by following the warning below, which tells the user to delete the - * legacy attribute. Reading either unguarded would st_errAbort on the other's config. - */ - bool hasLegacy = cactusParams_has(params, 2, "bar", "partialOrderAlignment"); - int64_t usePoa = hasLegacy ? cactusParams_get_int(params, 2, "bar", "partialOrderAlignment") : 1; - if (!cactusParams_has(params, 2, "bar", "baseAligner")) { - return usePoa ? BASE_ALIGNER_ABPOA : BASE_ALIGNER_PECAN; +static void dump_window_fasta_and_matrix(Msa *msa, uint8_t **bseqs, const int *mat, + const char *input_path, const char *matrix_path) { + FILE *mat_file = fopen(matrix_path, "w"); + if (mat_file != NULL) { + fprintf(mat_file, "\tA\tC\tG\tT\tN\n"); + for (size_t i = 0; i < 5; ++i) { + fprintf(mat_file, "%c", "ACGTN"[i]); + for (size_t j = 0; j < 5; ++j) { + fprintf(mat_file, "\t%d", mat[i * 5 + j]); + } + fprintf(mat_file, "\n"); + } + fclose(mat_file); } - char *name = cactusParams_get_string(params, 2, "bar", "baseAligner"); - BaseAligner engine; - if (strcmp(name, "pecan") == 0) { - engine = BASE_ALIGNER_PECAN; - } else if (strcmp(name, "abpoa") == 0) { - engine = BASE_ALIGNER_ABPOA; - } else if (strcmp(name, "minipoa") == 0) { - engine = BASE_ALIGNER_MINIPOA; - } else { - st_errAbort("Unknown ; expected pecan, abpoa or minipoa", name); - engine = BASE_ALIGNER_ABPOA; /* not reached */ + FILE *fa_file = fopen(input_path, "w"); + if (fa_file == NULL) { + return; } - free(name); + for (int64_t i = 0; i < msa->seq_no; ++i) { + fprintf(fa_file, ">%" PRIi64 "\n", i); + for (int64_t j = 0; j < msa->seq_lens[i]; ++j) { + fputc(msa_to_base(bseqs[i][j]), fa_file); + } + fputc('\n', fa_file); + } + fclose(fa_file); +} - /* - * Both attributes present and disagreeing is worth saying out loud. partialOrderAlignment="0" - * is how the config has always documented "use pecan", so someone who sets it and gets abpoa - * anyway should not have to discover that from the alignment. - */ - bool poaImplied = engine != BASE_ALIGNER_PECAN; - if (hasLegacy && (usePoa != 0) != poaImplied) { - st_logCritical("Warning: overrides , which asks for the opposite. baseAligner wins; remove the other to silence this.\n", - baseAligner_toString(engine), usePoa); +/* ============================================================================ */ +/* abPOA backend */ +/* ============================================================================ */ + +abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { + abpoa_para_t *abpt = abpoa_init_para(); + + // output options + abpt->out_msa = 1; // generate Row-Column multiple sequence alignment(RC-MSA), set 0 to disable + abpt->out_cons = 0; // generate consensus sequence, set 0 to disable + + // alignment mode. 0:global alignment, 1:local, 2:extension + // only global works + abpt->align_mode = ABPOA_GLOBAL_MODE; + + // banding parameters + abpt->wb = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentBandConstant"); + abpt->wf = cactusParams_get_float(params, 3, "bar", "poa", "partialOrderAlignmentBandFraction"); + + // gap scoring model + abpt->gap_open1 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapOpenPenalty1"); + abpt->gap_ext1 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapExtensionPenalty1"); + abpt->gap_open2 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapOpenPenalty2"); + abpt->gap_ext2 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapExtensionPenalty2"); + + // seeding paramters + abpt->disable_seeding = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentDisableSeeding"); + assert(abpt->disable_seeding == 0 || abpt->disable_seeding == 1); + abpt->k = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMinimizerK"); + abpt->w = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMinimizerW"); + abpt->min_w = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentMinimizerMinW"); + + // progressive toggle + abpt->progressive_poa = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentProgressiveMode"); + + // generate the substitution matrix + abpt->use_score_matrix = 0; + abpoa_post_set_para(abpt); + + // optionally override the substitution matrix + char *submat_string = cactusParams_get_string(params, 3, "bar", "poa", "partialOrderAlignmentSubMatrix"); + if (submat_string && strlen(submat_string) > 0) { + // Note, this will be used to explicitly override abpoa's subsitution matrix just before aligning + abpt->use_score_matrix = 1; + assert(abpt->m == 5); + int count = 0; + for (char* val = strtok(submat_string, " "); val != NULL; val = strtok(NULL, " ")) { + abpt->mat[count++] = atoi(val); + } + assert(count == 25); + int i; abpt->min_mis = 0, abpt->max_mat = 0; + for (i = 0; i < abpt->m * abpt->m; ++i) { + if (abpt->mat[i] > abpt->max_mat) + abpt->max_mat = abpt->mat[i]; + if (-abpt->mat[i] > abpt->min_mis) + abpt->min_mis = -abpt->mat[i]; + } } - return engine; + free(submat_string); + + return abpt; +} + +// It turns out abpoa can write to these, so we make a quick copy before using +static abpoa_para_t *copy_abpoa_params(abpoa_para_t *abpt) { + abpoa_para_t *abpt_cpy = abpoa_init_para(); + abpt_cpy->out_msa = 1; + abpt_cpy->out_cons = 0; + abpt_cpy->align_mode = abpt->align_mode; + abpt_cpy->wb = abpt->wb; + abpt_cpy->wf = abpt->wf; + abpt_cpy->match = abpt->match; + abpt_cpy->mismatch = abpt->mismatch; + abpt_cpy->gap_mode = abpt->gap_mode; + abpt_cpy->gap_open1 = abpt->gap_open1; + abpt_cpy->gap_ext1 = abpt->gap_ext1; + abpt_cpy->gap_open2 = abpt->gap_open2; + abpt_cpy->gap_ext2 = abpt->gap_ext2; + abpt_cpy->disable_seeding = abpt->disable_seeding; + abpt_cpy->k = abpt->k; + abpt_cpy->w = abpt->w; + abpt_cpy->min_w = abpt->min_w; + abpt_cpy->progressive_poa = abpt->progressive_poa; + abpt_cpy->use_score_matrix = 0; + abpoa_post_set_para(abpt_cpy); + abpt_cpy->use_score_matrix = abpt->use_score_matrix; + if (abpt->use_score_matrix == 1) { + memcpy(abpt_cpy->mat, abpt->mat, abpt->m * abpt->m * sizeof(int)); + } + abpt_cpy->max_mat = abpt->max_mat; + abpt_cpy->min_mis = abpt->min_mis; + return abpt_cpy; } -const char *baseAligner_toString(BaseAligner engine) { - switch (engine) { - case BASE_ALIGNER_PECAN: return "pecan"; - case BASE_ALIGNER_ABPOA: return "abpoa"; - case BASE_ALIGNER_MINIPOA: return "minipoa"; +// dump the abpoa input to files, and return a command line for running abpoa on them +char* dump_abpoa_input(Msa* msa, abpoa_para_t* abpt, uint8_t **bseqs, char* abpoa_input_path, char* abpoa_matrix_path, + char* abpoa_command_path, char* abpoa_output_path) { + // dump the abpoa input sequences to a FASTA file + FILE* dump_file = fopen(abpoa_input_path, "w"); + for (int64_t i = 0; i < msa->seq_no; ++i) { + int64_t seq_len = msa->seq_lens[i]; + char* buffer = (char*)malloc((seq_len + 1) * sizeof(char)); + for (int64_t j = 0; j < seq_len; ++j) { + buffer[j] = msa_to_base(bseqs[i][j]); + } + buffer[msa->seq_lens[i]] = '\0'; + fprintf(dump_file, ">%ld\n%s\n", i, buffer); + free(buffer); + } + fclose(dump_file); + + // dump the abpoa input matrix to file + FILE* mat_file = fopen(abpoa_matrix_path, "w"); + fprintf(mat_file, "\tA\tC\tG\tT\tN\n"); + for (size_t i = 0; i < 5; ++i) { + fprintf(mat_file, "%c", "ACGTN"[i]); + for (size_t j = 0; j < 5; ++j) { + fprintf(mat_file, "\t%d", abpt->mat[i * 5 + j]); + } + fprintf(mat_file, "\n"); + } + fclose(mat_file); + + // make a command line + char* abpoa_command = st_malloc(4096 * sizeof(char)); + sprintf(abpoa_command, "abpoa %s -O %d,%d -E %d,%d -b %d -f %lf -t %s -r 1 -m 0", + abpoa_input_path, + abpt->gap_open1, + abpt->gap_open2, + abpt->gap_ext1, + abpt->gap_ext2, + abpt->wb, + abpt->wf, + abpoa_matrix_path); + if (!abpt->disable_seeding) { + strcat(abpoa_command, " -S"); + char kw_opts[128]; + sprintf(kw_opts, " -k %d -w %d -n %d", abpt->k, abpt->w, abpt->min_w); + strcat(abpoa_command, kw_opts); + } + if (abpt->progressive_poa) { + strcat(abpoa_command, " -p"); + } + strcat(abpoa_command, " > "); + strcat(abpoa_command, abpoa_output_path); + + // dump the command line + FILE* cmd_file = fopen(abpoa_command_path, "w"); + fprintf(cmd_file, "%s\n", abpoa_command); + fclose(cmd_file); + + return abpoa_command; +} + +#ifdef CACTUS_ABPOA_FROM_COMMAND_LINE +void abpoa_msa_from_command_line(char* abpoa_command_line, char* abpoa_output_path, uint8_t*** msa_seq, int* col_no) { + // run abpoa + st_system(abpoa_command_line); + + // read the result (ascii alignment) back into memory + size_t n_rows = 0; + size_t n_cols = 0; + FILE* msa_file = fopen(abpoa_output_path, "r"); + int64_t buf_size = 500000; + char* buf = st_malloc(buf_size * sizeof(char)); + + while (benLine(&buf, &buf_size, msa_file) != -1) { + if (strlen(buf) && buf[0] != '>') { + ++n_rows; + } + } + + *msa_seq = st_malloc(n_rows * sizeof(uint8_t*)); + fclose(msa_file); + msa_file = fopen(abpoa_output_path, "r"); + n_rows = 0; + while (benLine(&buf, &buf_size, msa_file) != -1) { + if (strlen(buf) && buf[0] != '>') { + if (n_cols == 0) { + n_cols = strlen(buf); + } else { + assert(n_cols == strlen(buf)); + } + (*msa_seq)[n_rows] = st_malloc(n_cols * sizeof(uint8_t)); + for (size_t i = 0; i < n_cols; ++i) { + (*msa_seq)[n_rows][i] = msa_to_byte(buf[i]); + } + ++n_rows; + } + } + fclose(msa_file); + *col_no = (int)n_cols; + + free(buf); +} +#endif + +/* + * The abPOA backend. Fills msa->msa_seq and msa->column_no for one window. + */ +static void run_abpoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_parameters, + int64_t max_prog_rows, double max_prog_length_diff) { + // init abpoa + abpoa_t *ab = abpoa_init(); + abpoa_para_t *abpt = copy_abpoa_params(poa_parameters->abpt); + if (msa->seq_no > max_prog_rows || + // note: these are sorted by length excep in unit tests + (1. - (double)msa->seq_lens[msa->seq_no-1] / (double)msa->seq_lens[0] > max_prog_length_diff)) { + abpt->progressive_poa = 0; + } + + // dump the input to file, if asked to at run time + char abpoa_input_path[1024], abpoa_matrix_path[1024], abpoa_command_path[1024], abpoa_output_path[1024]; + char *abpoa_command_line = NULL; + if (bar_dump_dir != NULL) { + // The old name keyed off the Msa pointer, so two windows that reused the same freed + // allocation silently overwrote each other. pid + counter is unique for the run. + int64_t dump_id; +#if defined(_OPENMP) +#pragma omp atomic capture +#endif + dump_id = ++bar_dump_counter; + sprintf(abpoa_input_path, "%s/bar_window_%d_%" PRIi64 ".fa", bar_dump_dir, (int)getpid(), dump_id); + sprintf(abpoa_matrix_path, "%s.mat", abpoa_input_path); + sprintf(abpoa_command_path, "%s.cmd", abpoa_input_path); + sprintf(abpoa_output_path, "%s.out", abpoa_input_path); + abpoa_command_line = dump_abpoa_input(msa, abpt, bseqs, + abpoa_input_path, abpoa_matrix_path, abpoa_command_path, abpoa_output_path); + } + +#ifdef CACTUS_ABPOA_FROM_COMMAND_LINE + // run abpoa from the command line + if (abpoa_command_line == NULL) { + st_errAbort("CACTUS_ABPOA_FROM_COMMAND_LINE needs CACTUS_BAR_DUMP_DIR set in the environment"); + } + abpoa_msa_from_command_line(abpoa_command_line, abpoa_output_path, &(msa->msa_seq), &(msa->column_no)); + + int test_cols = 0; + uint8_t** test_msa = NULL; + abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); + // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they + // end up in the ab->abc struct -- we extract them here + test_msa = ab->abc->msa_base; + ab->abc->msa_base = NULL; + test_cols = ab->abc->msa_len; + + // sanity check to make sure we get the same output + assert(msa->column_no == test_cols); + for (int i = 0; i < msa->seq_no; ++i) { + for (int j = 0; j < test_cols; ++j) { + //todo: not sure why this doesn't work anymore !!!! + //assert(test_msa[i][j] == msa->msa_seq[i][j]); + } + free(test_msa[i]); + } + free(test_msa); +#else + // perform abpoa-msa + abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); + // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they + // end up in the ab->abc struct -- we extract them here + msa->msa_seq = ab->abc->msa_base; + ab->abc->msa_base = NULL; + msa->column_no = ab->abc->msa_len; +#endif + + // The dumps are kept. Deleting them on success made the switch useless for its main + // job -- looking at a window that aligned badly rather than one that crashed. + if (abpoa_command_line != NULL) { + free(abpoa_command_line); } - return "unknown"; + + // free abpoa + abpoa_free(ab); + abpoa_free_para(abpt); } +/* ============================================================================ */ +/* minipoa backend */ +/* ============================================================================ */ + /* * minipoa's substitution matrix comes from ; its gap penalties do not. * @@ -707,46 +724,45 @@ static void *minipoaParameters_constructFromCactusParams(CactusParams *params, P } #endif -PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, BaseAligner engine) { - if (engine == BASE_ALIGNER_PECAN) { - return NULL; - } - PoaParameters *poaParameters = st_calloc(1, sizeof(PoaParameters)); - poaParameters->engine = engine; - if (engine == BASE_ALIGNER_ABPOA) { - poaParameters->abpt = abpoaParamaters_constructFromCactusParams(params); - } else { - poaParameters->mpt = minipoaParameters_constructFromCactusParams(params, poaParameters); -#ifdef HAVE_MINIPOA - if (poaParameters->progressive) { - // Same settings with the guide tree off, for windows too wide to afford a dense NxN - // distance matrix. abPOA caps this the same way, via partialOrderAlignmentProgressiveMaxRows. - PoaParameters scratch = *poaParameters; - minipoa_para_t *plain = minipoaParameters_constructFromCactusParams(params, &scratch); - minipoa_set_progressive(plain, 0); - poaParameters->mptNoProgressive = plain; - } -#endif - } - return poaParameters; -} +/* + * A command line that replays this window through minipoa. Gap penalties are negated: minipoa + * maximises, so penalties are negative there, while cactus and abpoa carry them positive. + */ +static char *dump_minipoa_input(Msa *msa, PoaParameters *pp, uint8_t **bseqs, char *input_path, + char *matrix_path, char *command_path, char *output_path) { + dump_window_fasta_and_matrix(msa, bseqs, pp->mat, input_path, matrix_path); -void poaParameters_destruct(PoaParameters *poaParameters) { - if (poaParameters == NULL) { - return; + int f = pp->bandFraction > 0.0 ? (int)(1.0 / pp->bandFraction + 0.5) : 0; + 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->seeding) { + char kw_opts[128]; + sprintf(kw_opts, " -S -k %d -w %d", pp->minimizerK, pp->minimizerW); + strcat(command, kw_opts); + if (pp->anchorWindow > 0) { + sprintf(kw_opts, " -W %d", pp->anchorWindow); + strcat(command, kw_opts); + } } - if (poaParameters->abpt != NULL) { - abpoa_free_para(poaParameters->abpt); + // -p and -B are on in the shipped config, and progressive ordering in particular changes the + // alignment, so a replay that omitted them would not reproduce the window it is meant to + // explain. + if (pp->progressive) { + strcat(command, " -p"); } -#ifdef HAVE_MINIPOA - if (poaParameters->mpt != NULL) { - minipoa_free_para((minipoa_para_t *)poaParameters->mpt); + if (pp->adaptiveBand) { + strcat(command, " -B"); } - if (poaParameters->mptNoProgressive != NULL) { - minipoa_free_para((minipoa_para_t *)poaParameters->mptNoProgressive); + strcat(command, " > "); + strcat(command, output_path); + + FILE *cmd_file = fopen(command_path, "w"); + if (cmd_file != NULL) { + fprintf(cmd_file, "%s\n", command); + fclose(cmd_file); } -#endif - free(poaParameters); + return command; } /* @@ -820,99 +836,139 @@ static void run_minipoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_par #endif } +/* ============================================================================ */ +/* Base aligner selection and dispatch */ +/* ============================================================================ */ + /* - * Align one window with whichever engine was selected, filling msa->msa_seq and msa->column_no. - * - * The contract both engines meet: seq_no rows of column_no bytes, one malloc per row plus one for - * the row array (msa_destruct frees them that way), rows in input order, values 0-4 for ACGTN and - * 5 for a gap. Everything around this -- the sliding window, the empty-sequence hack, trimming, - * stitching, block extraction -- is engine-neutral and shared. - * - * max_prog_rows and max_prog_length_diff are abPOA's progressive-mode guards; minipoa has no - * equivalent knob, so they are inert for it. + * picks the engine. It is read with cactusParams_has so that a config + * predating the attribute -- including any a user has saved -- still selects what it used to via + * the older boolean. */ -static void run_poa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_parameters, - int64_t max_prog_rows, double max_prog_length_diff) { - if (poa_parameters->engine == BASE_ALIGNER_MINIPOA) { - run_minipoa_window(msa, bseqs, poa_parameters); - return; - } - - // init abpoa - abpoa_t *ab = abpoa_init(); - abpoa_para_t *abpt = copy_abpoa_params(poa_parameters->abpt); - if (msa->seq_no > max_prog_rows || - // note: these are sorted by length excep in unit tests - (1. - (double)msa->seq_lens[msa->seq_no-1] / (double)msa->seq_lens[0] > max_prog_length_diff)) { - abpt->progressive_poa = 0; +BaseAligner baseAligner_constructFromCactusParams(CactusParams *params) { + /* + * Both attributes are optional, in both directions. An old config has only + * partialOrderAlignment; a config written from now on may reasonably have only baseAligner -- + * including one produced by following the warning below, which tells the user to delete the + * legacy attribute. Reading either unguarded would st_errAbort on the other's config. + */ + bool hasLegacy = cactusParams_has(params, 2, "bar", "partialOrderAlignment"); + int64_t usePoa = hasLegacy ? cactusParams_get_int(params, 2, "bar", "partialOrderAlignment") : 1; + if (!cactusParams_has(params, 2, "bar", "baseAligner")) { + return usePoa ? BASE_ALIGNER_ABPOA : BASE_ALIGNER_PECAN; } - - // dump the input to file, if asked to at run time - char abpoa_input_path[1024], abpoa_matrix_path[1024], abpoa_command_path[1024], abpoa_output_path[1024]; - char *abpoa_command_line = NULL; - if (bar_dump_dir != NULL) { - // The old name keyed off the Msa pointer, so two windows that reused the same freed - // allocation silently overwrote each other. pid + counter is unique for the run. - int64_t dump_id; -#if defined(_OPENMP) -#pragma omp atomic capture -#endif - dump_id = ++bar_dump_counter; - sprintf(abpoa_input_path, "%s/bar_window_%d_%" PRIi64 ".fa", bar_dump_dir, (int)getpid(), dump_id); - sprintf(abpoa_matrix_path, "%s.mat", abpoa_input_path); - sprintf(abpoa_command_path, "%s.cmd", abpoa_input_path); - sprintf(abpoa_output_path, "%s.out", abpoa_input_path); - abpoa_command_line = dump_abpoa_input(msa, abpt, bseqs, - abpoa_input_path, abpoa_matrix_path, abpoa_command_path, abpoa_output_path); + char *name = cactusParams_get_string(params, 2, "bar", "baseAligner"); + BaseAligner engine; + if (strcmp(name, "pecan") == 0) { + engine = BASE_ALIGNER_PECAN; + } else if (strcmp(name, "abpoa") == 0) { + engine = BASE_ALIGNER_ABPOA; + } else if (strcmp(name, "minipoa") == 0) { + engine = BASE_ALIGNER_MINIPOA; + } else { + st_errAbort("Unknown ; expected pecan, abpoa or minipoa", name); + engine = BASE_ALIGNER_ABPOA; /* not reached */ } + free(name); -#ifdef CACTUS_ABPOA_FROM_COMMAND_LINE - // run abpoa from the command line - if (abpoa_command_line == NULL) { - st_errAbort("CACTUS_ABPOA_FROM_COMMAND_LINE needs CACTUS_BAR_DUMP_DIR set in the environment"); + /* + * Both attributes present and disagreeing is worth saying out loud. partialOrderAlignment="0" + * is how the config has always documented "use pecan", so someone who sets it and gets abpoa + * anyway should not have to discover that from the alignment. + */ + bool poaImplied = engine != BASE_ALIGNER_PECAN; + if (hasLegacy && (usePoa != 0) != poaImplied) { + st_logCritical("Warning: overrides , which asks for the opposite. baseAligner wins; remove the other to silence this.\n", + baseAligner_toString(engine), usePoa); } - abpoa_msa_from_command_line(abpoa_command_line, abpoa_output_path, &(msa->msa_seq), &(msa->column_no)); + return engine; +} - int test_cols = 0; - uint8_t** test_msa = NULL; - abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); - // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they - // end up in the ab->abc struct -- we extract them here - test_msa = ab->abc->msa_base; - ab->abc->msa_base = NULL; - test_cols = ab->abc->msa_len; +const char *baseAligner_toString(BaseAligner engine) { + switch (engine) { + case BASE_ALIGNER_PECAN: return "pecan"; + case BASE_ALIGNER_ABPOA: return "abpoa"; + case BASE_ALIGNER_MINIPOA: return "minipoa"; + } + return "unknown"; +} - // sanity check to make sure we get the same output - assert(msa->column_no == test_cols); - for (int i = 0; i < msa->seq_no; ++i) { - for (int j = 0; j < test_cols; ++j) { - //todo: not sure why this doesn't work anymore !!!! - //assert(test_msa[i][j] == msa->msa_seq[i][j]); - } - free(test_msa[i]); +PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, BaseAligner engine) { + if (engine == BASE_ALIGNER_PECAN) { + return NULL; } - free(test_msa); -#else - // perform abpoa-msa - abpoa_msa(ab, abpt, msa->seq_no, NULL, msa->seq_lens, bseqs, NULL, NULL); - // abpoa's interface has changed a bit -- instead of passing in pointers to the results, they - // end up in the ab->abc struct -- we extract them here - msa->msa_seq = ab->abc->msa_base; - ab->abc->msa_base = NULL; - msa->column_no = ab->abc->msa_len; + PoaParameters *poaParameters = st_calloc(1, sizeof(PoaParameters)); + poaParameters->engine = engine; + if (engine == BASE_ALIGNER_ABPOA) { + poaParameters->abpt = abpoaParamaters_constructFromCactusParams(params); + } else { + poaParameters->mpt = minipoaParameters_constructFromCactusParams(params, poaParameters); +#ifdef HAVE_MINIPOA + if (poaParameters->progressive) { + // Same settings with the guide tree off, for windows too wide to afford a dense NxN + // distance matrix. abPOA caps this the same way, via partialOrderAlignmentProgressiveMaxRows. + PoaParameters scratch = *poaParameters; + minipoa_para_t *plain = minipoaParameters_constructFromCactusParams(params, &scratch); + minipoa_set_progressive(plain, 0); + poaParameters->mptNoProgressive = plain; + } #endif + } + return poaParameters; +} - // The dumps are kept. Deleting them on success made the switch useless for its main - // job -- looking at a window that aligned badly rather than one that crashed. - if (abpoa_command_line != NULL) { - free(abpoa_command_line); +void poaParameters_destruct(PoaParameters *poaParameters) { + if (poaParameters == NULL) { + return; + } + if (poaParameters->abpt != NULL) { + abpoa_free_para(poaParameters->abpt); } +#ifdef HAVE_MINIPOA + if (poaParameters->mpt != NULL) { + minipoa_free_para((minipoa_para_t *)poaParameters->mpt); + } + if (poaParameters->mptNoProgressive != NULL) { + minipoa_free_para((minipoa_para_t *)poaParameters->mptNoProgressive); + } +#endif + free(poaParameters); +} - // free abpoa - abpoa_free(ab); - abpoa_free_para(abpt); +/* + * Align one window with whichever engine was selected. + * + * The contract both backends meet: seq_no rows of column_no bytes, one malloc per row plus one + * for the row array (msa_destruct frees them that way), rows in input order, values 0-4 for ACGTN + * and 5 for a gap. Everything around this -- the sliding window, the empty-sequence hack, + * trimming, stitching, block extraction -- is engine-neutral and shared. + */ +static void run_poa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_parameters, + int64_t max_prog_rows, double max_prog_length_diff) { + switch (poa_parameters->engine) { + case BASE_ALIGNER_ABPOA: + run_abpoa_window(msa, bseqs, poa_parameters, max_prog_rows, max_prog_length_diff); + break; + case BASE_ALIGNER_MINIPOA: + run_minipoa_window(msa, bseqs, poa_parameters); + break; + default: + st_errAbort("run_poa_window called with base aligner %s, which produces no MSA", + baseAligner_toString(poa_parameters->engine)); + } } + + + + + +/* ============================================================================ */ +/* Shared: windowed MSA construction, engine-neutral */ +/* ============================================================================ */ + + Msa *msa_make_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no, int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { @@ -1198,6 +1254,11 @@ void alignmentBlock_destruct(AlignmentBlock *alignmentBlock) { } } +/* ============================================================================ */ +/* Shared: adjacency strings out of the cactus graph */ +/* ============================================================================ */ + + char *get_adjacency_string(Cap *cap, int64_t *length, bool return_string) { assert(!cap_getSide(cap)); Sequence *sequence = cap_getSequence(cap); @@ -1325,6 +1386,11 @@ char *get_adjacency_string_and_overlap(Cap *cap, int *length, int64_t *overlap, return adjacency_string; } +/* ============================================================================ */ +/* Shared: MSA to alignment blocks to pinches */ +/* ============================================================================ */ + + /** * Gets the length and sequences present in the next maximal gapless alignment block. * @param msa The msa to scan From 176d21d70de1a7989ad536626561250b55f108bd Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 16 Sep 2026 09:13:39 -0400 Subject: [PATCH 3/4] bar: minipoa window back to 10000, guide tree off The 100000 window was picked from the plan text before anything was measured, and evolver could never have caught it: BAR adjacencies there have a median length of 29bp, so no window setting between 10kb and 1Mb changes a single alignment. A random sample of 6000 dumped windows from an evolver run found zero that reached even the 10kb cap. Measured on human/chimp chr10 instead, where BAR is 85% of the runtime and adjacencies are real, at 8 cores: window 10kb 15kb 20kb 30kb 50kb+ peak 3.5G 5.9G 9.3G 15.1G OOM-killed above 23G BAR 191s 196s 197s 204s - Memory goes roughly as window x (bandConstant + bandFraction x window), so it is quadratic-ish, not flat, and 100000 is comfortably fatal. BAR time is flat across the whole safe range, so a larger window buys nothing even where it fits. Neither minimizer seeding nor a tighter band rescued 100kb -- both died at the same ceiling, because the first sequences build a full-size graph with no consensus to anchor against. So minipoa gets abpoa's window, and the claim that it "affords a much larger window" is removed from the config, the docs and bar.c. That claim was wrong: minipoa's advantage is per-window cost at the same setting -- 3.5G and 191s against abpoa's 7.9G and 221s on the same data. minipoaProgressiveMode also defaults off now. It was on to match abpoa, on the reasoning that insertion order matters for POA. It does in principle, but this guide tree has not earned its cost on anything measurable: bit-identical output on evolver, and 7% more memory and 2% more time on chr10. Co-Authored-By: Claude Opus 5 (1M context) --- ReleaseNotes.md | 2 +- bar/impl/bar.c | 4 ++-- doc/progressive.md | 2 +- src/cactus/cactus_progressive_config.xml | 22 ++++++++++++++-------- src/cactus/pipeline/cactus_workflow.py | 10 +++++----- 5 files changed, 23 insertions(+), 17 deletions(-) diff --git a/ReleaseNotes.md b/ReleaseNotes.md index 6f6cd4e59..e225e7c88 100644 --- a/ReleaseNotes.md +++ b/ReleaseNotes.md @@ -1,6 +1,6 @@ # Unreleased -- `minipoa` added as an optional third base aligner in BAR, selected with `` in the config. abPOA remains the default. On the evolver tests the two are equivalent: primates 0.9974 against abPOA's 0.9973, mammals 0.8506 against 0.8513 -- at roughly half the peak memory (1.0 GB against 2.3 GB on evolver mammals). +- `minipoa` added as an optional third base aligner in BAR, selected with `` in the config. abPOA remains the default. The two are equivalent on accuracy (evolver primates 0.9974 against abPOA's 0.9973, mammals 0.8506 against 0.8513) and minipoa uses roughly half the memory at the same window size: on human/chimp chr10, 3.5 GB and 191 s of BAR against abPOA's 7.9 GB and 221 s. Note the win is per-window cost, not a larger window -- minipoa's memory is quadratic-ish in the window too, and 50 kb upwards is OOM-killed on that data. - minipoa shares abPOA's substitution matrix (from ``) but has its own gap penalties in ``, because the same nominal values do not produce the same alignment. abPOA's gap model is convex, so its effective extension past ~28bp is 1, not 30; minipoa has a single affine piece and charging it 30 made long gaps ~30x dearer, costing six points of mammals accuracy. `--lastTrain` fits a *single* affine gap model, which is exactly minipoa's model, so the learned gaps are passed to minipoa verbatim; abPOA additionally gets a synthesised second piece it needs for stability, which minipoa does not. - `` replaces the `` boolean as the way to choose between pecan, abpoa and minipoa. Configs with no `baseAligner` still work: the boolean is used as before. - Fix `--lastTrain` ignoring `partialOrderAlignmentTrainedGapOpen2Factor`: `apply_long_gap` was being passed the extension factor in both the open and extension positions, so the open factor was read from the config and discarded. Both default to 3, so results with the shipped config are unchanged; setting them to different values previously did nothing (and setting the extension factor to 1 raised an assertion). diff --git a/bar/impl/bar.c b/bar/impl/bar.c index c0860193e..93bc0c82d 100644 --- a/bar/impl/bar.c +++ b/bar/impl/bar.c @@ -77,8 +77,8 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * bool pruneOutStubAlignments = cactusParams_get_int(params, 3, "bar", "pecan", "pruneOutStubAlignments"); // Poa params. The window and the mask filter come from whichever engine is selected: - // abpoa's memory is quadratic in the window so it is held to 10s of kb, while minipoa is the - // reason to have a second engine at all and can take a much larger one. + // Both engines are held to 10s of kb: memory is quadratic-ish in the window for either. + // minipoa's advantage is per-window cost, not a bigger window. int64_t poaWindow, maskFilter; if (engine == BASE_ALIGNER_MINIPOA) { poaWindow = cactusParams_get_int(params, 3, "bar", "minipoa", "minipoaWindow"); diff --git a/doc/progressive.md b/doc/progressive.md index e6b8ff693..2b35a8518 100644 --- a/doc/progressive.md +++ b/doc/progressive.md @@ -416,7 +416,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 `--maxMemory/--doubleMem` should be enough to work around this type of issue. * `--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 unless the memory estimate is more than the job can be given (the system memory on a single machine, or `--maxMemory`), in which case the estimate is scaled down and the pages are not kept. 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 much less memory per window, which is why its default window is 100000 rather than 10000; the memory *estimate* for `cactus_consolidated` is still the one fitted to abpoa, so under minipoa it over-requests rather than under-requests. 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, `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`. 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 669eba561..f893f8aa8 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -236,9 +236,10 @@ - + + Off by default: insertion order matters for POA in principle, but this tree has not earned its + cost on anything measured - bit-identical output on evolver, and 7% more memory plus 2% more + time on human/chimp chr10. Worth revisiting on more diverged input. --> diff --git a/src/cactus/pipeline/cactus_workflow.py b/src/cactus/pipeline/cactus_workflow.py index 3d1d99c50..d44168a77 100644 --- a/src/cactus/pipeline/cactus_workflow.py +++ b/src/cactus/pipeline/cactus_workflow.py @@ -85,17 +85,17 @@ def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_ma # the fit ran with it disabled -- so enabling it can only make the estimate conservative. # # The window comes from whichever base aligner is selected. The exponent above was fitted to - # abPOA runs, so applying it to minipoa's much larger window over-requests -- by roughly 7x at - # a 100kb window. That is deliberate until someone fits an exponent from real minipoa runs and - # records it: over-requesting wastes cluster share, under-requesting gets the largest ancestors - # OOM-killed, and those are the very jobs minipoa is for. + # abPOA runs. minipoa ships the same 10000 window, so the exponent is applied to the same + # number and nothing changes -- but if anyone raises minipoaWindow, note that minipoa's measured + # curve is far steeper than 0.43 in that range (about 1.35 on human/chimp chr10), so this + # estimate would under-request, which is the dangerous direction. bar_node = findRequiredNode(config_node, 'bar') base_aligner = getOptionalAttrib(bar_node, 'baseAligner', typeFn=str, default=None) if base_aligner is None: base_aligner = 'abpoa' if getOptionalAttrib(bar_node, 'partialOrderAlignment', typeFn=bool, default=True) else 'pecan' if base_aligner == 'minipoa': engine_node = bar_node.find('minipoa') - poa_window = getOptionalAttrib(engine_node, 'minipoaWindow', typeFn=int, default=100000) if engine_node is not None else 100000 + poa_window = getOptionalAttrib(engine_node, 'minipoaWindow', typeFn=int, default=10000) if engine_node is not None else 10000 else: poa_node = bar_node.find('poa') poa_window = getOptionalAttrib(poa_node, 'partialOrderAlignmentWindow', typeFn=int, default=10000) if poa_node is not None else 10000 From 2c110862a706cb1bb1499659fe5e1db7044559a7 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 16 Sep 2026 09:37:51 -0400 Subject: [PATCH 4/4] bar: build the window-dump paths without sprintf CI failed the build: CGL_DEBUG=ultra compiles with -Werror, and six sprintf(path, "%s.mat", other_path) calls between two char[1024] buffers trip -Wformat-overflow, because the compiler cannot prove the directory name leaves room for the suffix. impl/poaBarAligner.c:790:33: error: '.mat' directive writing 4 bytes into a region of size between 1 and 1024 [-Werror=format-overflow=] Both backends had their own copy of the same eight lines, so this fixes the duplication and the warning together: next_dump_paths() builds all four names once, with snprintf for the base name and memcpy for the suffixes. A dump directory too long to name a window now logs and skips the dump rather than writing a truncated path that could collide with another window's files. Local builds never caught this because they do not set CGL_DEBUG. Co-Authored-By: Claude Opus 5 (1M context) --- bar/impl/poaBarAligner.c | 65 +++++++++++++++++++++++++--------------- 1 file changed, 41 insertions(+), 24 deletions(-) diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index 3864998ae..790ec9e26 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -315,6 +315,41 @@ static void msa_fix_trimmed(Msa* msa) { /* Window dumping: shared between the backends */ /* ============================================================================ */ +/* + * Build the four file names for one dumped window. False when dumping is off. + * + * snprintf plus memcpy rather than sprintf: the compiler cannot prove the directory is short + * enough, and CGL_DEBUG=ultra builds with -Werror, so plain sprintf into a fixed buffer fails the + * build. A directory long enough to overflow gets no dump at all rather than a truncated name + * that could collide with another window's files. + */ +#define BAR_DUMP_PATH_MAX 1024 + +static bool next_dump_paths(char *fa, char *mat, char *cmd, char *out) { + if (bar_dump_dir == NULL) { + return false; + } + int64_t dump_id; +#if defined(_OPENMP) +#pragma omp atomic capture +#endif + /* + * The old name keyed off the Msa pointer, so two windows that reused the same freed + * allocation silently overwrote each other. pid + counter is unique for the run. + */ + dump_id = ++bar_dump_counter; + int len = snprintf(fa, BAR_DUMP_PATH_MAX - 5, "%s/bar_window_%d_%" PRIi64 ".fa", + bar_dump_dir, (int)getpid(), dump_id); + if (len < 0 || len >= BAR_DUMP_PATH_MAX - 5) { + st_logCritical("bar: CACTUS_BAR_DUMP_DIR is too long to name window files; not dumping\n"); + return false; + } + memcpy(mat, fa, (size_t)len); memcpy(mat + len, ".mat", 5); + memcpy(cmd, fa, (size_t)len); memcpy(cmd + len, ".cmd", 5); + memcpy(out, fa, (size_t)len); memcpy(out + len, ".out", 5); + return true; +} + /* * Write the window as FASTA plus its 5x5 matrix. The matrix file format is the one both abpoa -t * and minipoa -m read, so either aligner can be pointed straight at it. @@ -559,20 +594,10 @@ static void run_abpoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_param } // dump the input to file, if asked to at run time - char abpoa_input_path[1024], abpoa_matrix_path[1024], abpoa_command_path[1024], abpoa_output_path[1024]; + char abpoa_input_path[BAR_DUMP_PATH_MAX], abpoa_matrix_path[BAR_DUMP_PATH_MAX]; + char abpoa_command_path[BAR_DUMP_PATH_MAX], abpoa_output_path[BAR_DUMP_PATH_MAX]; char *abpoa_command_line = NULL; - if (bar_dump_dir != NULL) { - // The old name keyed off the Msa pointer, so two windows that reused the same freed - // allocation silently overwrote each other. pid + counter is unique for the run. - int64_t dump_id; -#if defined(_OPENMP) -#pragma omp atomic capture -#endif - dump_id = ++bar_dump_counter; - sprintf(abpoa_input_path, "%s/bar_window_%d_%" PRIi64 ".fa", bar_dump_dir, (int)getpid(), dump_id); - sprintf(abpoa_matrix_path, "%s.mat", abpoa_input_path); - sprintf(abpoa_command_path, "%s.cmd", abpoa_input_path); - sprintf(abpoa_output_path, "%s.out", abpoa_input_path); + if (next_dump_paths(abpoa_input_path, abpoa_matrix_path, abpoa_command_path, abpoa_output_path)) { abpoa_command_line = dump_abpoa_input(msa, abpt, bseqs, abpoa_input_path, abpoa_matrix_path, abpoa_command_path, abpoa_output_path); } @@ -778,18 +803,10 @@ static void run_minipoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_par (void)msa; (void)bseqs; (void)poa_parameters; st_errAbort("minipoa was selected but this cactus was built with minipoa=off (see include.mk)"); #else - char input_path[1024], matrix_path[1024], command_path[1024], output_path[1024]; + char input_path[BAR_DUMP_PATH_MAX], matrix_path[BAR_DUMP_PATH_MAX]; + char command_path[BAR_DUMP_PATH_MAX], output_path[BAR_DUMP_PATH_MAX]; char *command_line = NULL; - if (bar_dump_dir != NULL) { - int64_t dump_id; -#if defined(_OPENMP) -#pragma omp atomic capture -#endif - dump_id = ++bar_dump_counter; - sprintf(input_path, "%s/bar_window_%d_%" PRIi64 ".fa", bar_dump_dir, (int)getpid(), dump_id); - sprintf(matrix_path, "%s.mat", input_path); - sprintf(command_path, "%s.cmd", input_path); - sprintf(output_path, "%s.out", input_path); + if (next_dump_paths(input_path, matrix_path, command_path, output_path)) { command_line = dump_minipoa_input(msa, poa_parameters, bseqs, input_path, matrix_path, command_path, output_path); }