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 8993e78eb..92fc5af8f 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} @@ -209,6 +212,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 @@ -246,7 +257,7 @@ evolver_test_primates_pangenome_resume_local: all ${CWD}/test/primates-truth.maf evolver_test_primates_pangenome_extend_local: all ${CWD}/test/primates-truth.maf PYTHONPATH="${CWD}/submodules/" CACTUS_BINARIES_MODE=local CACTUS_DOCKER_MODE=0 ${PYTHON} -m pytest ${pytestOpts} -s test/evolverTest.py::TestCase::testEvolverPrimatesPangenomeExtendLocal -evolver_test_all_local: evolver_test_local evolver_test_prepare_toil evolver_test_decomposed_local evolver_test_prepare_no_outgroup_local evolver_test_poa_local evolver_test_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 @@ -331,6 +342,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..e225e7c88 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. 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). +- `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..93bc0c82d 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: + // 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"); + 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..790ec9e26 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -6,15 +6,44 @@ #include #include "abpoa.h" +#ifdef HAVE_MINIPOA +#include "minipoa_c.h" +#endif #include "poaBarAligner.h" #include "flowerAligner.h" #include #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. + * 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; -// 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) +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: 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 @@ -22,6 +51,341 @@ #include #endif + +/* ============================================================================ */ +/* Shared: alphabet, and the Msa the backends fill in */ +/* ============================================================================ */ + + +// char <--> uint8_t conversion copied over from abPOA example +// AaCcGgTtNn ==> 0,1,2,3,4 +static unsigned char nst_nt4_table[256] = { + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5 /*'-'*/, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 +}; + +// 65,97=>A, 67,99=>C, 71,103=>G, 84,85,116,117=>T, else=>N +static const char nst_nt256_table[256] = { + 'A', 'C', 'G', 'T', 'N', '-', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', '-', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'A', 'N', 'C', 'N', 'N', 'N', 'G', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'T', 'T', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'A', 'N', 'C', 'N', 'N', 'N', 'G', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'T', 'T', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', + 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N' +}; + +char msa_to_base(uint8_t n) { + return (char)nst_nt256_table[n]; +} + +uint8_t msa_to_byte(char c) { + return nst_nt4_table[(int)c]; +} + +static uint8_t rc_table[6] = { 3, 2, 1, 0, 4, 5 }; + +static inline uint8_t msa_to_rc(uint8_t n) { + return rc_table[n]; +} + + + + + +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"); +} + +/* ============================================================================ */ +/* Shared: column scores, trimming and window stitching */ +/* ============================================================================ */ + + +/** + * 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); + } + } + } +} + +/** + * Used to make two MSAs consistent with each other for a shared sequence + */ +static void trim(int64_t row1, Msa *msa1, float *column_scores1, + int64_t row2, Msa *msa2, float *column_scores2, int64_t overlap) { + if(overlap == 0) { // There is no overlap, so no need to trim either MSA + return; + } + assert(overlap > 0); // Otherwise the overlap must be positive + + int64_t seq_len1 = msa1->seq_lens[row1]; // The prefix length of the forward complement sequence in the first MSA + int64_t seq_len2 = msa2->seq_lens[row2]; // The prefix length of the reverse complement sequence in the second MSA + // They can be different if either MSA does not include the whole sequence + assert(overlap <= seq_len1); // The overlap must be less than the length of the prefixes + assert(overlap <= seq_len2); + + // Get the cumulative cut scores for the columns containing the shared sequence + float *cu_column_scores1 = st_malloc(msa1->column_no * sizeof(float)); + float *cu_column_scores2 = st_malloc(msa2->column_no * sizeof(float)); + sum_column_scores(row1, msa1, column_scores1, cu_column_scores1); + sum_column_scores(row2, msa2, column_scores2, cu_column_scores2); + + // The score if we cut all of the overlap in msa1 and keep all of the overlap in msa2 + assert(seq_len2 <= msa2->column_no); + float max_cut_score = cu_column_scores2[seq_len2-1]; + if(overlap < seq_len1) { // The overlap is less than the length of the first sequence + assert(seq_len1-overlap-1 >= 0); + max_cut_score += cu_column_scores1[seq_len1-overlap-1]; // We will keep everything before the overlap + } + int64_t max_overlap_cut_point = 0; // the length of the prefix of the overlap of msa1 to keep + + // Now walk through each possible cut point within the overlap + for(int64_t i=0; i= 0); // Sanity check + float cut_score = cu_column_scores1[seq_len1-overlap+i] + cu_column_scores2[seq_len2-i-2]; // The score if we keep prefix up to + // and including column i of MSA1's overlap, and the prefix of msa2 up to and including column seq_len-i-2 + if(cut_score > max_cut_score) { + max_overlap_cut_point = i + 1; + max_cut_score = cut_score; + } + } + + // The score if we cut all of msa2's overlap and keep all of msa1's + float f = cu_column_scores1[seq_len1-1]; + if(overlap < seq_len2) { + assert(seq_len2-overlap-1 >= 0); + f += cu_column_scores2[seq_len2-overlap-1]; + } + + if(f > max_cut_score) { + max_cut_score = f; + max_overlap_cut_point = overlap; + } + + // Now trim back the two MSAs + assert(max_overlap_cut_point <= overlap); + trim_msa_suffix(msa1, column_scores1, row1, seq_len1 - overlap + max_overlap_cut_point); + trim_msa_suffix(msa2, column_scores2, row2, seq_len2 - max_overlap_cut_point); + + free(cu_column_scores1); + free(cu_column_scores2); +} + +/** + * recompute the seq_lens of a trimmed msa and clip off empty suffix columns + * (todo: can this be built into trimming code?) + */ +static void msa_fix_trimmed(Msa* msa) { + for (int64_t i = 0; i < msa->seq_no; ++i) { + // recompute the seq_len + msa->seq_lens[i] = 0; + for (int64_t j = 0; j < msa->column_no; ++j) { + if (msa_to_base(msa->msa_seq[i][j]) != '-') { + ++msa->seq_lens[i]; + } + } + } + // trim empty columns + int64_t empty_columns = 0; + for (bool still_empty = true; empty_columns < msa->column_no; ++empty_columns) { + for (int64_t i = 0; i < msa->seq_no && still_empty; ++i) { + still_empty = msa_to_base(msa->msa_seq[i][msa->column_no - 1 - empty_columns]) == '-'; + } + if (!still_empty) { + break; + } + } + msa->column_no -= empty_columns; +} + + +/* ============================================================================ */ +/* 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. + */ +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); +} + +/* ============================================================================ */ +/* abPOA backend */ +/* ============================================================================ */ + abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { abpoa_para_t *abpt = abpoa_init_para(); @@ -86,88 +450,32 @@ 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 -static unsigned char nst_nt4_table[256] = { - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5 /*'-'*/, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 -}; - -// 65,97=>A, 67,99=>C, 71,103=>G, 84,85,116,117=>T, else=>N -static const char nst_nt256_table[256] = { - 'A', 'C', 'G', 'T', 'N', '-', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', '-', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'A', 'N', 'C', 'N', 'N', 'N', 'G', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'T', 'T', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'A', 'N', 'C', 'N', 'N', 'N', 'G', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'T', 'T', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', - 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N', 'N' -}; - -char msa_to_base(uint8_t n) { - return (char)nst_nt256_table[n]; -} - -uint8_t msa_to_byte(char c) { - return nst_nt4_table[(int)c]; -} - -static uint8_t rc_table[6] = { 3, 2, 1, 0, 4, 5 }; - -static inline uint8_t msa_to_rc(uint8_t n) { - return rc_table[n]; + 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; } -#ifdef CACTUS_ABPOA_MSA_DUMP_DIR // 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 +535,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) { @@ -272,197 +579,415 @@ void abpoa_msa_from_command_line(char* abpoa_command_line, char* abpoa_output_pa } #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]); +/* + * 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[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 (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); } - 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"); +#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"); } - fprintf(f, "\n"); -} + 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 -/** - * 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]); - } - } + // 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); } -/** - * Returns an array of floats, one for each corresponding column in the MSA. Each float - * is the score of the column in the alignment. +/* ============================================================================ */ +/* minipoa backend */ +/* ============================================================================ */ + +/* + * 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. */ -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]++; - } +#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(column_scores[i] >= 1.0) { - column_scores[i]--; + if (count != 25) { + st_errAbort(" needs 25 values, got %d", count); } - assert(column_scores[i] >= 0.0); + minipoa_set_score_matrix(mpt, mat); + memcpy(out->mat, mat, sizeof(mat)); } - return column_scores; + 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; } - -/** - * 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". +#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 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 +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 -/** - * Removes the suffix of the given row from the MSA and updates the column scores. suffix_start is the beginning - * suffix to remove. +/* + * 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 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); - } +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; } -/** - * Used to make two MSAs consistent with each other for a shared sequence +/* + * 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 trim(int64_t row1, Msa *msa1, float *column_scores1, - int64_t row2, Msa *msa2, float *column_scores2, int64_t overlap) { - if(overlap == 0) { // There is no overlap, so no need to trim either MSA - return; +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[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 (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); } - assert(overlap > 0); // Otherwise the overlap must be positive - - int64_t seq_len1 = msa1->seq_lens[row1]; // The prefix length of the forward complement sequence in the first MSA - int64_t seq_len2 = msa2->seq_lens[row2]; // The prefix length of the reverse complement sequence in the second MSA - // They can be different if either MSA does not include the whole sequence - assert(overlap <= seq_len1); // The overlap must be less than the length of the prefixes - assert(overlap <= seq_len2); - - // Get the cumulative cut scores for the columns containing the shared sequence - float *cu_column_scores1 = st_malloc(msa1->column_no * sizeof(float)); - float *cu_column_scores2 = st_malloc(msa2->column_no * sizeof(float)); - sum_column_scores(row1, msa1, column_scores1, cu_column_scores1); - sum_column_scores(row2, msa2, column_scores2, cu_column_scores2); - // The score if we cut all of the overlap in msa1 and keep all of the overlap in msa2 - assert(seq_len2 <= msa2->column_no); - float max_cut_score = cu_column_scores2[seq_len2-1]; - if(overlap < seq_len1) { // The overlap is less than the length of the first sequence - assert(seq_len1-overlap-1 >= 0); - max_cut_score += cu_column_scores1[seq_len1-overlap-1]; // We will keep everything before the overlap + 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; } - int64_t max_overlap_cut_point = 0; // the length of the prefix of the overlap of msa1 to keep - // Now walk through each possible cut point within the overlap - for(int64_t i=0; i= 0); // Sanity check - float cut_score = cu_column_scores1[seq_len1-overlap+i] + cu_column_scores2[seq_len2-i-2]; // The score if we keep prefix up to - // and including column i of MSA1's overlap, and the prefix of msa2 up to and including column seq_len-i-2 - if(cut_score > max_cut_score) { - max_overlap_cut_point = i + 1; - max_cut_score = cut_score; + 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 +} - // The score if we cut all of msa2's overlap and keep all of msa1's - float f = cu_column_scores1[seq_len1-1]; - if(overlap < seq_len2) { - assert(seq_len2-overlap-1 >= 0); - f += cu_column_scores2[seq_len2-overlap-1]; +/* ============================================================================ */ +/* Base aligner selection and dispatch */ +/* ============================================================================ */ + +/* + * 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; +} - if(f > max_cut_score) { - max_cut_score = f; - max_overlap_cut_point = overlap; +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"; +} - // Now trim back the two MSAs - assert(max_overlap_cut_point <= overlap); - trim_msa_suffix(msa1, column_scores1, row1, seq_len1 - overlap + max_overlap_cut_point); - trim_msa_suffix(msa2, column_scores2, row2, seq_len2 - max_overlap_cut_point); +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; +} - free(cu_column_scores1); - free(cu_column_scores2); +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); } -/** - * recompute the seq_lens of a trimmed msa and clip off empty suffix columns - * (todo: can this be built into trimming code?) +/* + * 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 msa_fix_trimmed(Msa* msa) { - for (int64_t i = 0; i < msa->seq_no; ++i) { - // recompute the seq_len - msa->seq_lens[i] = 0; - for (int64_t j = 0; j < msa->column_no; ++j) { - if (msa_to_base(msa->msa_seq[i][j]) != '-') { - ++msa->seq_lens[i]; - } - } - } - // trim empty columns - int64_t empty_columns = 0; - for (bool still_empty = true; empty_columns < msa->column_no; ++empty_columns) { - for (int64_t i = 0; i < msa->seq_no && still_empty; ++i) { - still_empty = msa_to_base(msa->msa_seq[i][msa->column_no - 1 - empty_columns]) == '-'; - } - if (!still_empty) { +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)); } - msa->column_no -= empty_columns; } + + + + + +/* ============================================================================ */ +/* 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, abpoa_para_t *poa_parameters) { + int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { assert(seq_no > 0); @@ -562,71 +1087,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 +1223,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); @@ -810,6 +1271,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); @@ -937,6 +1403,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 @@ -1148,7 +1619,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..2b35a8518 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 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/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..f893f8aa8 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -234,8 +234,17 @@ - + + @@ -243,6 +252,7 @@ + + + + + + + + + + + + + + + + + + diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index b770d3f87..534db52fa 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..d44168a77 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. 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=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 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 3eb4c5560..9d3090a88 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 98ec8f82c..ea70b72c6 100644 --- a/test/evolverTest.py +++ b/test/evolverTest.py @@ -1820,6 +1820,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 """