From c1b24ea6f9bd7bfe2dd6511daeef4e34bb3aba4f Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 29 Sep 2026 22:40:28 -0400 Subject: [PATCH 01/10] last-train: lastdb -c and -C2, LAST's genome-genome recipe Without them each lastal pass spends most of its time aligning tandem repeats every which way. On chr17 (4 threads, 1 Mb sample): CHM13 vs GRCh38 went from 2831 to 207 CPU seconds for the first pass and ~350 to ~60 for each one after; human vs mouse chr11 from ~440 to ~105 a pass. The trained model is the same to within last-train's own iteration jitter. Co-Authored-By: Claude Opus 5.5 --- src/cactus/paf/last_scoring.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index 1e810944d..82f8fe1ba 100644 --- a/src/cactus/paf/last_scoring.py +++ b/src/cactus/paf/last_scoring.py @@ -292,7 +292,17 @@ def last_train(job, config, seq_order, seq_id_map, ref_name=None): # note: there are some specific options for distant genomes that should be # incorporated if/when this ever gets used in progressive cactus - train_cmd = ['last-train', '--revsym', '--matsym', '--gapsym', + # + # -C2 (with lastdb -c below) is LAST's own recipe for genome-genome training: without them each + # lastal pass spends most of its time aligning tandem repeats every which way. lastdb -c stops + # them seeding (its default -R01 marks them with tantan; the input's own soft-masking is not + # used), and -C2 drops gapless alignments nested in two others. Measured with 4 threads and + # last-train's 1 Mb sample against an 83 Mb chr17 database, in lastal CPU seconds per pass: + # CHM13 against GRCh38, centromere included, went from 2831 to 207 for the first pass (run with + # LAST's default HUMSUM scores, far too loose for so close a pair) and from ~350 to ~60 for each + # pass after it; human against mouse chr11 went from ~440 to ~105. The trained model came out + # the same both times, to within the jitter between last-train's own iterations. + train_cmd = ['last-train', '--revsym', '--matsym', '--gapsym', '-C2', '-P', str(job.cores), name1 + '_db', fa2_path] train_file = os.path.join(work_dir, '{}_{}.train'.format(name1, name2)) @@ -300,7 +310,7 @@ def last_train(job, config, seq_order, seq_id_map, ref_name=None): # to converge on too little alignment, or writing something parse_train_file won't accept), # the alignment falls back to the default scores, or to another chromosome's model in batch mode try: - cactus_call(parameters=['lastdb', name1 + '_db', fa1_path, '-P', str(job.cores)]) + cactus_call(parameters=['lastdb', '-c', name1 + '_db', fa1_path, '-P', str(job.cores)]) cactus_call(parameters=train_cmd, outfile=train_file) parse_train_file(train_file) except Exception as e: From bc6777e36b4c10fb706f61e55421c1f2d251a820 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Sun, 4 Oct 2026 13:34:23 -0400 Subject: [PATCH 02/10] abPOA: per-read scoring parameters Pins the fork's main-per-read-params: abPOA as master has it (the fork's main, v1.5.7 with ai-optimization) plus one commit, an optional abpoa_para_t.read_abpt, so a caller can score each read against the graph with its own parameters. Unset, abPOA is unchanged (identical CLI output over its test data). Co-Authored-By: Claude Opus 5.5 --- submodules/abPOA | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/abPOA b/submodules/abPOA index 36110ac80..47361c5d1 160000 --- a/submodules/abPOA +++ b/submodules/abPOA @@ -1 +1 @@ -Subproject commit 36110ac805f6f2d250621f8a08e5ef092e09d48b +Subproject commit 47361c5d1de96ea95bd5bb9a8ab1f03820fbe0a0 From c33b54cf4d15b443a087131508d395b32943a67f Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Wed, 30 Sep 2026 09:13:21 -0400 Subject: [PATCH 03/10] bar: optional per-row abPOA scoring models, chosen by divergence One matrix scores every row of a POA, but the rows are genomes at very different distances: two close ingroups and more distant outgroups. With and children (each a maxDistance plus its own substitution matrix and gap penalties), each row is scored with the model for the distance from its genome to the nearest other genome in the alignment (or the nearest ingroup, with partialOrderAlignmentRowModelDistance="ingroup"), measured on the unscaled reconstructionTree. Off unless configured: with no row models, evolver mammals aligns exactly as before (identical hal2maf output). On evolver mammals it keeps the ingroup-distance model's mouse-rat precision while recovering half the ancestor recall that model loses, but overall it is within 0.1 F1 of a single matrix at the ingroup distance. Evolver primates can't tell them apart. Needs harder data. Co-Authored-By: Claude Opus 5.5 --- bar/impl/poaBarAligner.c | 275 ++++++++++++++++++++++++++++++++++++--- bar/inc/poaBarAligner.h | 44 +++++++ 2 files changed, 299 insertions(+), 20 deletions(-) diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index feb9d92bb..1021c6f15 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -387,7 +387,23 @@ static void dump_window_fasta_and_matrix(Msa *msa, uint8_t **bseqs, const int *m /* abPOA backend */ /* ============================================================================ */ -abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { +// A scoring attribute of , or of its child (a row model) when model is not NULL +static int64_t poa_scoring_int(CactusParams *params, const char *model, const char *name) { + return model == NULL ? cactusParams_get_int(params, 3, "bar", "poa", name) + : cactusParams_get_int(params, 4, "bar", "poa", model, name); +} + +static char *poa_scoring_string(CactusParams *params, const char *model, const char *name) { + return model == NULL ? cactusParams_get_string(params, 3, "bar", "poa", name) + : cactusParams_get_string(params, 4, "bar", "poa", model, name); +} + +/* + * abPOA parameters from , with the gap penalties and substitution matrix taken from its + * child instead when model is not NULL. Everything else (banding, seeding, progressive + * mode) always comes from , so row models differ from the default only in how they score. + */ +static abpoa_para_t *abpoaParameters_construct(CactusParams *params, const char *model) { abpoa_para_t *abpt = abpoa_init_para(); // output options @@ -403,10 +419,10 @@ abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { abpt->wf = cactusParams_get_float(params, 3, "bar", "poa", "partialOrderAlignmentBandFraction"); // gap scoring model - abpt->gap_open1 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapOpenPenalty1"); - abpt->gap_ext1 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapExtensionPenalty1"); - abpt->gap_open2 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapOpenPenalty2"); - abpt->gap_ext2 = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentGapExtensionPenalty2"); + abpt->gap_open1 = poa_scoring_int(params, model, "partialOrderAlignmentGapOpenPenalty1"); + abpt->gap_ext1 = poa_scoring_int(params, model, "partialOrderAlignmentGapExtensionPenalty1"); + abpt->gap_open2 = poa_scoring_int(params, model, "partialOrderAlignmentGapOpenPenalty2"); + abpt->gap_ext2 = poa_scoring_int(params, model, "partialOrderAlignmentGapExtensionPenalty2"); // seeding paramters abpt->disable_seeding = cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentDisableSeeding"); @@ -423,7 +439,7 @@ abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { abpoa_post_set_para(abpt); // optionally override the substitution matrix - char *submat_string = cactusParams_get_string(params, 3, "bar", "poa", "partialOrderAlignmentSubMatrix"); + char *submat_string = poa_scoring_string(params, model, "partialOrderAlignmentSubMatrix"); if (submat_string && strlen(submat_string) > 0) { // Note, this will be used to explicitly override abpoa's subsitution matrix just before aligning abpt->use_score_matrix = 1; @@ -446,6 +462,10 @@ abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { return abpt; } +abpoa_para_t *abpoaParamaters_constructFromCactusParams(CactusParams *params) { + return abpoaParameters_construct(params, NULL); +} + // It turns out abpoa can write to these, so we make a quick copy before using static abpoa_para_t *copy_abpoa_params(abpoa_para_t *abpt) { abpoa_para_t *abpt_cpy = abpoa_init_para(); @@ -584,7 +604,7 @@ void abpoa_msa_from_command_line(char* abpoa_command_line, char* abpoa_output_pa * 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) { + int64_t max_prog_rows, double max_prog_length_diff, int *row_models) { // init abpoa abpoa_t *ab = abpoa_init(); abpoa_para_t *abpt = copy_abpoa_params(poa_parameters->abpt); @@ -593,6 +613,23 @@ static void run_abpoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_param (1. - (double)msa->seq_lens[msa->seq_no-1] / (double)msa->seq_lens[0] > max_prog_length_diff)) { abpt->progressive_poa = 0; } + + // per-row scoring: copies, like abpt's, since abpoa can write to its params and they are shared + // between threads + abpoa_para_t **model_copies = NULL; + if (row_models != NULL && poa_parameters->rowModelNo > 0) { + model_copies = st_calloc(poa_parameters->rowModelNo, sizeof(abpoa_para_t *)); + abpt->read_abpt = st_calloc(msa->seq_no, sizeof(void *)); + for (int64_t i = 0; i < msa->seq_no; ++i) { + int m = row_models[i]; + if (m >= 0) { + if (model_copies[m] == NULL) { + model_copies[m] = copy_abpoa_params(poa_parameters->rowModelAbpt[m]); + } + abpt->read_abpt[i] = model_copies[m]; + } + } + } // 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]; @@ -647,6 +684,15 @@ static void run_abpoa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_param // free abpoa abpoa_free(ab); + if (model_copies != NULL) { + for (int64_t m = 0; m < poa_parameters->rowModelNo; ++m) { + if (model_copies[m] != NULL) { + abpoa_free_para(model_copies[m]); + } + } + free(model_copies); + free(abpt->read_abpt); + } abpoa_free_para(abpt); } @@ -895,6 +941,145 @@ const char *baseAligner_toString(BaseAligner engine) { return "unknown"; } +static void addDistanceNodes(stHash *nodes, stTree *tree) { + if (stTree_getLabel(tree) != NULL) { + stHash_insert(nodes, (void *)stTree_getLabel(tree), tree); + } + for (int64_t i = 0; i < stTree_getChildNumber(tree); i++) { + addDistanceNodes(nodes, stTree_getChild(tree, i)); + } +} + +/* + * and its children .. (see + * PoaParameters). A count rather than probing for the children, because asking CactusParams about + * a child that is not there prints an error. + */ +static void rowModels_load(CactusParams *params, PoaParameters *pp) { + pp->rowModelNo = cactusParams_has(params, 3, "bar", "poa", "partialOrderAlignmentRowModels") ? + cactusParams_get_int(params, 3, "bar", "poa", "partialOrderAlignmentRowModels") : 0; + if (pp->rowModelNo <= 0) { + pp->rowModelNo = 0; + return; + } + if (cactusParams_has(params, 3, "bar", "poa", "partialOrderAlignmentRowModelDistance")) { + char *rule = cactusParams_get_string(params, 3, "bar", "poa", "partialOrderAlignmentRowModelDistance"); + if (strcmp(rule, "ingroup") == 0) { + pp->rowModelToIngroups = true; + } else if (strcmp(rule, "nearest") != 0) { + st_errAbort("partialOrderAlignmentRowModelDistance must be nearest or ingroup, not %s", rule); + } + free(rule); + } + pp->rowModelMaxDistance = st_malloc(sizeof(double) * pp->rowModelNo); + pp->rowModelAbpt = st_malloc(sizeof(abpoa_para_t *) * pp->rowModelNo); + for (int64_t m = 0; m < pp->rowModelNo; ++m) { + char *name = stString_print("rowModel%" PRIi64, m); + pp->rowModelMaxDistance[m] = cactusParams_get_float(params, 4, "bar", "poa", name, "maxDistance"); + if (m > 0 && pp->rowModelMaxDistance[m] < pp->rowModelMaxDistance[m-1]) { + st_errAbort("bar row models must be in ascending maxDistance, but %s has %f after %f", name, + pp->rowModelMaxDistance[m], pp->rowModelMaxDistance[m-1]); + } + pp->rowModelAbpt[m] = abpoaParameters_construct(params, name); + free(name); + } + if (cactusParams_has(params, 2, "reference", "reconstructionTree")) { + char *newick = cactusParams_get_string(params, 2, "reference", "reconstructionTree"); + pp->distanceTree = stTree_parseNewickString(newick); + pp->distanceNodes = stHash_construct3(stHash_stringKey, stHash_stringEqualKey, NULL, NULL); + addDistanceNodes(pp->distanceNodes, pp->distanceTree); + free(newick); + } + st_logInfo("bar: %" PRIi64 " row models, distances to the nearest %s on the %s tree\n", pp->rowModelNo, + pp->rowModelToIngroups ? "ingroup" : "genome", pp->distanceTree != NULL ? "reconstruction" : "event"); +} + +/* + * Path length between two genomes, on the reconstruction tree when it names both, else on the event + * tree (whose branches are the ones scaled for alignment sensitivity). + */ +static double rowModels_distance(PoaParameters *pp, Event *event1, Event *event2) { + stTree *node1 = pp->distanceNodes == NULL ? NULL : stHash_search(pp->distanceNodes, (void *)event_getHeader(event1)); + stTree *node2 = pp->distanceNodes == NULL ? NULL : stHash_search(pp->distanceNodes, (void *)event_getHeader(event2)); + if (node1 != NULL && node2 != NULL) { + double d1 = 0.0; + for (stTree *a = node1; a != NULL; a = stTree_getParent(a)) { + double d2 = 0.0; + for (stTree *b = node2; b != NULL; b = stTree_getParent(b)) { + if (a == b) { + return d1 + d2; + } + d2 += stTree_getBranchLength(b); + } + d1 += stTree_getBranchLength(a); + } + } + Event *ancestor = eventTree_getCommonAncestor(event1, event2); + double d = 0.0; + for (Event *e = event1; e != ancestor; e = event_getParent(e)) { + d += event_getBranchLength(e); + } + for (Event *e = event2; e != ancestor; e = event_getParent(e)) { + d += event_getBranchLength(e); + } + return d; +} + +/* + * The row model for each row of an end's alignment (-1 for the default), from the distance between + * the row's genome and the nearest other genome among the rows (the nearest other ingroup genome with + * rowModelToIngroups). NULL when there are no row models. + */ +static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { + if (pp == NULL || pp->rowModelNo == 0) { + return NULL; + } + // distances between the distinct genomes only: rows are often several copies of a few genomes + Event **events = st_malloc(sizeof(Event *) * row_no); + int64_t *row_event = st_malloc(sizeof(int64_t) * row_no); + int64_t event_no = 0; + for (int64_t i = 0; i < row_no; ++i) { + Event *event = cap_getEvent(caps[i]); + int64_t j = 0; + while (j < event_no && events[j] != event) { + ++j; + } + if (j == event_no) { + events[event_no++] = event; + } + row_event[i] = j; + } + double *nearest = st_malloc(sizeof(double) * event_no); + for (int64_t j = 0; j < event_no; ++j) { + nearest[j] = INFINITY; + } + for (int64_t j = 0; j < event_no; ++j) { + for (int64_t k = j + 1; k < event_no; ++k) { + double d = rowModels_distance(pp, events[j], events[k]); + if (!pp->rowModelToIngroups || !event_isOutgroup(events[k])) { + nearest[j] = d < nearest[j] ? d : nearest[j]; + } + if (!pp->rowModelToIngroups || !event_isOutgroup(events[j])) { + nearest[k] = d < nearest[k] ? d : nearest[k]; + } + } + } + int *models = st_malloc(sizeof(int) * row_no); + for (int64_t i = 0; i < row_no; ++i) { + models[i] = -1; + for (int64_t m = 0; m < pp->rowModelNo; ++m) { + if (nearest[row_event[i]] <= pp->rowModelMaxDistance[m]) { + models[i] = (int)m; + break; + } + } + } + free(events); + free(row_event); + free(nearest); + return models; +} + PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, BaseAligner engine) { if (engine == BASE_ALIGNER_PECAN) { return NULL; @@ -903,6 +1088,7 @@ PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, Bas poaParameters->engine = engine; if (engine == BASE_ALIGNER_ABPOA) { poaParameters->abpt = abpoaParamaters_constructFromCactusParams(params); + rowModels_load(params, poaParameters); } else { poaParameters->mpt = minipoaParameters_constructFromCactusParams(params, poaParameters); #ifdef HAVE_MINIPOA @@ -948,6 +1134,15 @@ void poaParameters_destruct(PoaParameters *poaParameters) { abpoa_free_para(poaParameters->abpt); } free(poaParameters->referenceEvent); + for (int64_t m = 0; m < poaParameters->rowModelNo; ++m) { + abpoa_free_para(poaParameters->rowModelAbpt[m]); + } + free(poaParameters->rowModelAbpt); + free(poaParameters->rowModelMaxDistance); + if (poaParameters->distanceTree != NULL) { + stHash_destruct(poaParameters->distanceNodes); + stTree_destruct(poaParameters->distanceTree); + } #ifdef HAVE_MINIPOA if (poaParameters->mpt != NULL) { minipoa_free_para((minipoa_para_t *)poaParameters->mpt); @@ -968,10 +1163,10 @@ void poaParameters_destruct(PoaParameters *poaParameters) { * trimming, stitching, block extraction -- is engine-neutral and shared. */ static void run_poa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_parameters, - int64_t max_prog_rows, double max_prog_length_diff) { + int64_t max_prog_rows, double max_prog_length_diff, int *row_models) { switch (poa_parameters->engine) { case BASE_ALIGNER_ABPOA: - run_abpoa_window(msa, bseqs, poa_parameters, max_prog_rows, max_prog_length_diff); + run_abpoa_window(msa, bseqs, poa_parameters, max_prog_rows, max_prog_length_diff, row_models); break; case BASE_ALIGNER_MINIPOA: run_minipoa_window(msa, bseqs, poa_parameters); @@ -994,6 +1189,13 @@ static void run_poa_window(Msa *msa, uint8_t **bseqs, PoaParameters *poa_paramet Msa *msa_make_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no, int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { + return msa_make_partial_order_alignment_with_row_models(seqs, seq_lens, seq_no, window_size, max_prog_rows, + max_prog_length_diff, poa_parameters, NULL); +} + +Msa *msa_make_partial_order_alignment_with_row_models(char **seqs, int *seq_lens, int64_t seq_no, int64_t window_size, + int64_t max_prog_rows, double max_prog_length_diff, + PoaParameters *poa_parameters, int *row_models) { assert(seq_no > 0); @@ -1094,7 +1296,7 @@ Msa *msa_make_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no } } - run_poa_window(msa, bseqs, poa_parameters, max_prog_rows, max_prog_length_diff); + run_poa_window(msa, bseqs, poa_parameters, max_prog_rows, max_prog_length_diff, row_models); // mask out empty sequences that were phonied in as Ns above for (int64_t i = 0; i < msa->seq_no && emptyCount > 0; ++i) { @@ -1306,13 +1508,22 @@ static bool poa_should_flip(Flower *flower, Cap **caps, char **strings, int *len Msa *msa_make_oriented_partial_order_alignment(char **seqs, int *seq_lens, int64_t seq_no, bool flip, int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters) { + return msa_make_oriented_partial_order_alignment_with_row_models(seqs, seq_lens, seq_no, flip, window_size, + max_prog_rows, max_prog_length_diff, + poa_parameters, NULL); +} + +Msa *msa_make_oriented_partial_order_alignment_with_row_models(char **seqs, int *seq_lens, int64_t seq_no, bool flip, + int64_t window_size, int64_t max_prog_rows, + double max_prog_length_diff, + PoaParameters *poa_parameters, int *row_models) { if (flip) { for (int64_t i = 0; i < seq_no; i++) { reverse_complement_in_place(seqs[i], seq_lens[i]); } } - Msa *msa = msa_make_partial_order_alignment(seqs, seq_lens, seq_no, window_size, max_prog_rows, - max_prog_length_diff, poa_parameters); + Msa *msa = msa_make_partial_order_alignment_with_row_models(seqs, seq_lens, seq_no, window_size, max_prog_rows, + max_prog_length_diff, poa_parameters, row_models); if (flip) { flip_msa_seq(msa); for (int64_t i = 0; i < seq_no; i++) { @@ -1337,14 +1548,26 @@ Msa **make_consistent_partial_order_alignments(int64_t end_no, int64_t *end_leng 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, PoaParameters *poa_parameters, bool *flips) { + return make_consistent_partial_order_alignments_with_row_models(end_no, end_lengths, end_strings, end_string_lengths, + right_end_indexes, right_end_row_indexes, overlaps, + window_size, max_prog_rows, max_prog_length_diff, + poa_parameters, flips, NULL); +} + +Msa **make_consistent_partial_order_alignments_with_row_models(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, PoaParameters *poa_parameters, + bool *flips, int **row_models) { // 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); for(int64_t i=0; i and , the event header REFERENCE goes by (NULL if none) */ PoaOrientation orientation; char *referenceEvent; + /* + * Per-row scoring, abPOA only. switches it on, and its + * children .. each carry a maxDistance and their own + * partialOrderAlignmentSubMatrix and gap penalties, in ascending maxDistance. A row is scored + * with the first whose maxDistance is at least the distance from its genome to the nearest other + * genome in the same alignment; rows beyond the last, or with no other genome, keep abpt. + * Distances are on (the unscaled tree), falling back to the event + * tree for a genome it does not name. rowModelNo is 0 when the config has no row models. + * With rowModelToIngroups (partialOrderAlignmentRowModelDistance="ingroup") the distance is to the + * nearest other ingroup genome instead, so an outgroup is scored for how far it is from what is + * being reconstructed rather than from another outgroup. + */ + int64_t rowModelNo; + bool rowModelToIngroups; + double *rowModelMaxDistance; + abpoa_para_t **rowModelAbpt; + stTree *distanceTree; + stHash *distanceNodes; // node label -> stTree node of distanceTree } PoaParameters; /** @@ -181,6 +199,23 @@ Msa *msa_make_oriented_partial_order_alignment(char **seqs, int *seq_lens, int64 int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters); +/** + * As msa_make_partial_order_alignment, with row_models[i] naming the poa_parameters row model + * that scores row i (-1 for the default). row_models may be NULL. + */ +Msa *msa_make_partial_order_alignment_with_row_models(char **seqs, int *seq_lens, int64_t seq_no, int64_t window_size, + int64_t max_prog_rows, double max_prog_length_diff, + PoaParameters *poa_parameters, int *row_models); + +/** + * Both: msa_make_oriented_partial_order_alignment with row models. Reverse complementing leaves each sequence + * in its row, so row_models applies either way round. + */ +Msa *msa_make_oriented_partial_order_alignment_with_row_models(char **seqs, int *seq_lens, int64_t seq_no, bool flip, + int64_t window_size, int64_t max_prog_rows, + double max_prog_length_diff, + PoaParameters *poa_parameters, int *row_models); + /** * Takes a set of ends and returns a set of consistent multiple alignments, * one for each of them. @@ -211,6 +246,15 @@ Msa **make_consistent_partial_order_alignments(int64_t end_no, int64_t *end_leng int64_t window_size, int64_t max_prog_rows, double max_prog_length_diff, PoaParameters *poa_parameters, bool *flips); +/** + * As make_consistent_partial_order_alignments, with row_models[i] the row models for end i (see + * msa_make_partial_order_alignment_with_row_models). row_models, or any entry of it, may be NULL. + */ +Msa **make_consistent_partial_order_alignments_with_row_models(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, PoaParameters *poa_parameters, + bool *flips, int **row_models); + /** * Represents a gapless alignment of a set of sequences. */ From e287ea4b345337c0e640856ba8fcc375e0250c67 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Thu, 1 Oct 2026 12:41:09 -0400 Subject: [PATCH 04/10] bar row models: under the ingroup rule, fall back to the nearest genome With partialOrderAlignmentRowModelDistance="ingroup", a row with no other ingroup genome among the rows of its alignment (only one child present, or only outgroups) was given the default scoring, which mixes the shipped matrix into an otherwise trained alignment. It now takes the model for its nearest genome of any kind instead. Co-Authored-By: Claude Opus 5.5 --- bar/impl/poaBarAligner.c | 24 +++++++++++++++++++----- bar/inc/poaBarAligner.h | 3 ++- 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index 1021c6f15..905e112c5 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -1049,18 +1049,31 @@ static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { } row_event[i] = j; } + // nearest other genome, and nearest other ingroup genome double *nearest = st_malloc(sizeof(double) * event_no); + double *nearestIngroup = st_malloc(sizeof(double) * event_no); for (int64_t j = 0; j < event_no; ++j) { - nearest[j] = INFINITY; + nearest[j] = nearestIngroup[j] = INFINITY; } for (int64_t j = 0; j < event_no; ++j) { for (int64_t k = j + 1; k < event_no; ++k) { double d = rowModels_distance(pp, events[j], events[k]); - if (!pp->rowModelToIngroups || !event_isOutgroup(events[k])) { - nearest[j] = d < nearest[j] ? d : nearest[j]; + nearest[j] = d < nearest[j] ? d : nearest[j]; + nearest[k] = d < nearest[k] ? d : nearest[k]; + if (!event_isOutgroup(events[k])) { + nearestIngroup[j] = d < nearestIngroup[j] ? d : nearestIngroup[j]; } - if (!pp->rowModelToIngroups || !event_isOutgroup(events[j])) { - nearest[k] = d < nearest[k] ? d : nearest[k]; + if (!event_isOutgroup(events[j])) { + nearestIngroup[k] = d < nearestIngroup[k] ? d : nearestIngroup[k]; + } + } + } + // with rowModelToIngroups, rows that have another ingroup among them are measured to the nearest one, + // and the rest (only one ingroup genome, or none, in this alignment) to the nearest genome of any kind + if (pp->rowModelToIngroups) { + for (int64_t j = 0; j < event_no; ++j) { + if (nearestIngroup[j] != INFINITY) { + nearest[j] = nearestIngroup[j]; } } } @@ -1077,6 +1090,7 @@ static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { free(events); free(row_event); free(nearest); + free(nearestIngroup); return models; } diff --git a/bar/inc/poaBarAligner.h b/bar/inc/poaBarAligner.h index 61a7d5fc5..089ab5d46 100644 --- a/bar/inc/poaBarAligner.h +++ b/bar/inc/poaBarAligner.h @@ -116,7 +116,8 @@ typedef struct _PoaParameters { * tree for a genome it does not name. rowModelNo is 0 when the config has no row models. * With rowModelToIngroups (partialOrderAlignmentRowModelDistance="ingroup") the distance is to the * nearest other ingroup genome instead, so an outgroup is scored for how far it is from what is - * being reconstructed rather than from another outgroup. + * being reconstructed rather than from another outgroup; a row with no other ingroup genome in the + * alignment falls back to the nearest genome of any kind. */ int64_t rowModelNo; bool rowModelToIngroups; From 66e46eb05b6c51bae23fef00e18bee60da2b289f Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Thu, 1 Oct 2026 15:15:59 -0400 Subject: [PATCH 05/10] bar: optional per-pair pecan models, chosen by divergence pecan aligns the sequences of an end in pairs, so each pair can be scored with the model for the distance between its two genomes. With and children (each a maxDistance and an hmm: a five-state pair HMM in the JSON hmm_jsonParse reads), a pair is aligned with the first model whose maxDistance covers its distance on the reconstruction tree. Pairs within one genome, and pairs beyond the last model, keep pecan's built-in model. Off unless configured. The genome distances the abPOA row models used are now shared (BarDistances, in endAligner). cPecan moves to pecan-params, which adds the per-pair callback (makeAlignmentWithPairStateMachines) to the multiple aligner. Re-aligning two evolver mammals nodes from fixed lastz output, against the truth: the ingroup-distance model alone gives the best mouse-rat precision (+0.29 points over pecan's own) but loses 0.46 points of ancestor recall; per pair keeps the recall (0.9478 vs 0.9481) and +0.17 of the precision. Co-Authored-By: Claude Opus 5.5 --- bar/impl/bar.c | 5 +- bar/impl/endAligner.c | 156 ++++++++++++++++++++++++++++++++++++++- bar/impl/flowerAligner.c | 10 +-- bar/impl/poaBarAligner.c | 58 ++------------- bar/inc/endAligner.h | 41 ++++++++++ bar/inc/flowerAligner.h | 5 +- bar/inc/poaBarAligner.h | 3 +- submodules/cPecan | 2 +- 8 files changed, 216 insertions(+), 64 deletions(-) diff --git a/bar/impl/bar.c b/bar/impl/bar.c index 59a54e268..ecdb24c6e 100644 --- a/bar/impl/bar.c +++ b/bar/impl/bar.c @@ -108,6 +108,8 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * float matchGamma = cactusParams_get_float(params, 3, "bar", "pecan", "matchGamma"); PairwiseAlignmentParameters *pairwiseAlignmentParameters = pairwiseAlignmentParameters_constructFromCactusParams(params); StateMachine *sM = stateMachine5_construct(fiveState); + // per-pair pecan models, if the config has any (see endAligner.h) + PecanPairModels *pecanPairModels = usePoa ? NULL : pecanPairModels_constructFromCactusParams(params); bool pruneOutStubAlignments = cactusParams_get_int(params, 3, "bar", "pecan", "pruneOutStubAlignments"); // Poa params. The window and the mask filter come from whichever engine is selected: @@ -233,7 +235,7 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * poaMaxProgRows, poaMaxLenDiff, poaParameters); st_logDebug("Created the poa alignments: %" PRIi64 " poa alignment blocks for flower\n", stList_length(alignments)); } else { - alignments = makeFlowerAlignment3(sM, flower, listOfEndAlignmentFiles, spanningTrees, maximumLength, + alignments = makeFlowerAlignment3(sM, pecanPairModels, flower, listOfEndAlignmentFiles, spanningTrees, maximumLength, useProgressiveMerging, matchGamma, pairwiseAlignmentParameters, pruneOutStubAlignments); st_logDebug("Created the alignment: %" PRIi64 " pairs for flower\n", stSortedSet_size(alignments)); @@ -356,6 +358,7 @@ void bar(stList *flowers, CactusParams *params, CactusDisk *cactusDisk, stList * pairwiseAlignmentBandingParameters_destruct(pairwiseAlignmentParameters); stateMachine_destruct(sM); + pecanPairModels_destruct(pecanPairModels); poaParameters_destruct(poaParameters); } diff --git a/bar/impl/endAligner.c b/bar/impl/endAligner.c index 646e23c30..0abf14b68 100644 --- a/bar/impl/endAligner.c +++ b/bar/impl/endAligner.c @@ -52,9 +52,158 @@ int alignedPair_cmpFn(const AlignedPair *alignedPair1, const AlignedPair *aligne return i; } +struct _barDistances { + stTree *tree; // the reconstruction tree, or NULL + stHash *nodes; // its nodes by label +}; + +static void barDistances_addNodes(stHash *nodes, stTree *tree) { + if (stTree_getLabel(tree) != NULL) { + stHash_insert(nodes, (void *)stTree_getLabel(tree), tree); + } + for (int64_t i = 0; i < stTree_getChildNumber(tree); i++) { + barDistances_addNodes(nodes, stTree_getChild(tree, i)); + } +} + +BarDistances *barDistances_constructFromCactusParams(CactusParams *params) { + BarDistances *distances = st_calloc(1, sizeof(BarDistances)); + if (cactusParams_has(params, 2, "reference", "reconstructionTree")) { + char *newick = cactusParams_get_string(params, 2, "reference", "reconstructionTree"); + distances->tree = stTree_parseNewickString(newick); + distances->nodes = stHash_construct3(stHash_stringKey, stHash_stringEqualKey, NULL, NULL); + barDistances_addNodes(distances->nodes, distances->tree); + free(newick); + } + return distances; +} + +bool barDistances_hasReconstructionTree(BarDistances *distances) { + return distances->tree != NULL; +} + +double barDistances_get(BarDistances *distances, Event *event1, Event *event2) { + stTree *node1 = distances->nodes == NULL ? NULL : stHash_search(distances->nodes, (void *)event_getHeader(event1)); + stTree *node2 = distances->nodes == NULL ? NULL : stHash_search(distances->nodes, (void *)event_getHeader(event2)); + if (node1 != NULL && node2 != NULL) { + double d1 = 0.0; + for (stTree *a = node1; a != NULL; a = stTree_getParent(a)) { + double d2 = 0.0; + for (stTree *b = node2; b != NULL; b = stTree_getParent(b)) { + if (a == b) { + return d1 + d2; + } + d2 += stTree_getBranchLength(b); + } + d1 += stTree_getBranchLength(a); + } + } + Event *ancestor = eventTree_getCommonAncestor(event1, event2); + double d = 0.0; + for (Event *e = event1; e != ancestor; e = event_getParent(e)) { + d += event_getBranchLength(e); + } + for (Event *e = event2; e != ancestor; e = event_getParent(e)) { + d += event_getBranchLength(e); + } + return d; +} + +void barDistances_destruct(BarDistances *distances) { + if (distances->tree != NULL) { + stHash_destruct(distances->nodes); + stTree_destruct(distances->tree); + } + free(distances); +} + +struct _pecanPairModels { + int64_t modelNo; + double *maxDistance; + StateMachine **stateMachines; + BarDistances *distances; +}; + +PecanPairModels *pecanPairModels_constructFromCactusParams(CactusParams *params) { + int64_t modelNo = cactusParams_has(params, 3, "bar", "pecan", "pecanPairModels") ? + cactusParams_get_int(params, 3, "bar", "pecan", "pecanPairModels") : 0; + if (modelNo <= 0) { + return NULL; + } + PecanPairModels *pairModels = st_calloc(1, sizeof(PecanPairModels)); + pairModels->modelNo = modelNo; + pairModels->maxDistance = st_malloc(sizeof(double) * modelNo); + pairModels->stateMachines = st_malloc(sizeof(StateMachine *) * modelNo); + for (int64_t m = 0; m < modelNo; m++) { + char *name = stString_print("pairModel%" PRIi64, m); + pairModels->maxDistance[m] = cactusParams_get_float(params, 4, "bar", "pecan", name, "maxDistance"); + if (m > 0 && pairModels->maxDistance[m] < pairModels->maxDistance[m - 1]) { + st_errAbort("pecan pair models must be in ascending maxDistance, but %s has %f after %f", name, + pairModels->maxDistance[m], pairModels->maxDistance[m - 1]); + } + char *json = cactusParams_get_string(params, 4, "bar", "pecan", name, "hmm"); + Hmm *hmm = hmm_jsonParse(json, strlen(json)); + if (hmm->type != fiveState) { + st_errAbort("pecan pair model %s is a type %i HMM; bar's pecan uses the five-state model (type %i)", + name, (int)hmm->type, (int)fiveState); + } + pairModels->stateMachines[m] = hmm_getStateMachine(hmm); + hmm_destruct(hmm); + free(json); + free(name); + } + pairModels->distances = barDistances_constructFromCactusParams(params); + st_logInfo("bar: %" PRIi64 " pecan pair models, distances on the %s tree\n", modelNo, + barDistances_hasReconstructionTree(pairModels->distances) ? "reconstruction" : "event"); + return pairModels; +} + +void pecanPairModels_destruct(PecanPairModels *pairModels) { + if (pairModels == NULL) { + return; + } + for (int64_t m = 0; m < pairModels->modelNo; m++) { + stateMachine_destruct(pairModels->stateMachines[m]); + } + free(pairModels->stateMachines); + free(pairModels->maxDistance); + barDistances_destruct(pairModels->distances); + free(pairModels); +} + +/* + * cPecan's per-pair callback: the model for the distance between the two sequences' genomes. + */ +typedef struct { + PecanPairModels *pairModels; + stList *events; // the genome of each sequence, by its index in the alignment +} PairModelChoice; + +static StateMachine *choosePairModel(int64_t seqX, int64_t seqY, void *extraArgs) { + PairModelChoice *choice = extraArgs; + Event *eventX = stList_get(choice->events, seqX), *eventY = stList_get(choice->events, seqY); + if (eventX == eventY) { + return NULL; + } + double distance = barDistances_get(choice->pairModels->distances, eventX, eventY); + for (int64_t m = 0; m < choice->pairModels->modelNo; m++) { + if (distance <= choice->pairModels->maxDistance[m]) { + return choice->pairModels->stateMachines[m]; + } + } + return NULL; +} + stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees, int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { + return makeEndAlignmentWithPairModels(sM, NULL, end, spanningTrees, maxSequenceLength, useProgressiveMerging, gapGamma, + pairwiseAlignmentBandingParameters); +} + +stSortedSet *makeEndAlignmentWithPairModels(StateMachine *sM, PecanPairModels *pairModels, End *end, int64_t spanningTrees, + int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma, + PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { //Make an alignment of the sequences in the ends //Get the adjacency sequences to be aligned. @@ -62,6 +211,7 @@ stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees, End_InstanceIterator *it = end_getInstanceIterator(end); stList *sequences = stList_construct3(0, (void (*)(void *))adjacencySequence_destruct); stList *seqFrags = stList_construct3(0, (void (*)(void *))seqFrag_destruct); + stList *events = stList_construct(); stHash *endInstanceNumbers = stHash_construct2(NULL, free); while((cap = end_getNext(it)) != NULL) { if(cap_getSide(cap)) { @@ -72,6 +222,7 @@ stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees, assert(cap_getAdjacency(cap) != NULL); End *otherEnd = end_getPositiveOrientation(cap_getEnd(cap_getAdjacency(cap))); stList_append(seqFrags, seqFrag_construct(adjacencySequence->string, 0, end_getName(otherEnd))); + stList_append(events, cap_getEvent(cap)); //Increase count of seqfrags with a given end. int64_t *c = stHash_search(endInstanceNumbers, otherEnd); if(c == NULL) { @@ -84,7 +235,10 @@ stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees, end_destructInstanceIterator(it); //Get the alignment. - MultipleAlignment *mA = makeAlignment(sM, seqFrags, spanningTrees, 100000000, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters); + PairModelChoice choice = { pairModels, events }; + MultipleAlignment *mA = makeAlignmentWithPairStateMachines(sM, pairModels == NULL ? NULL : choosePairModel, &choice, + seqFrags, spanningTrees, 100000000, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters); + stList_destruct(events); //Build an array of weights to reweight pairs in the alignment. int64_t *pairwiseAlignmentsPerSequenceNonCommonEnds = st_calloc(stList_length(seqFrags), sizeof(int64_t)); diff --git a/bar/impl/flowerAligner.c b/bar/impl/flowerAligner.c index 390d9bc3e..b93a1b849 100644 --- a/bar/impl/flowerAligner.c +++ b/bar/impl/flowerAligner.c @@ -528,7 +528,7 @@ static stSortedSet *getEndsToAlign(Flower *flower, int64_t maxSequenceLength) { * then call the makeFlowerAlignment2 consistency generating function. */ -static void computeMissingEndAlignments(StateMachine *sM, Flower *flower, stHash *endAlignments, int64_t spanningTrees, +static void computeMissingEndAlignments(StateMachine *sM, PecanPairModels *pairModels, Flower *flower, stHash *endAlignments, int64_t spanningTrees, int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters) { /* @@ -546,7 +546,7 @@ static void computeMissingEndAlignments(StateMachine *sM, Flower *flower, stHash stHash_insert( endAlignments, end, - makeEndAlignment(sM, end, spanningTrees, maxSequenceLength, + makeEndAlignmentWithPairModels(sM, pairModels, end, spanningTrees, maxSequenceLength, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters)); } else { @@ -562,7 +562,7 @@ stSortedSet *makeFlowerAlignment(StateMachine *sM, Flower *flower, int64_t spann bool useProgressiveMerging, float gapGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments) { stHash *endAlignments = stHash_construct2(NULL, (void(*)(void *)) stSortedSet_destruct); - computeMissingEndAlignments(sM, flower, endAlignments, spanningTrees, maxSequenceLength, + computeMissingEndAlignments(sM, NULL, flower, endAlignments, spanningTrees, maxSequenceLength, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters); return makeFlowerAlignment2(flower, endAlignments, pruneOutStubAlignments); } @@ -583,14 +583,14 @@ static void loadEndAlignments(Flower *flower, stHash *endAlignments, stList *lis } } -stSortedSet *makeFlowerAlignment3(StateMachine *sM, Flower *flower, stList *listOfEndAlignmentFiles, int64_t spanningTrees, +stSortedSet *makeFlowerAlignment3(StateMachine *sM, PecanPairModels *pairModels, Flower *flower, stList *listOfEndAlignmentFiles, int64_t spanningTrees, int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments) { stHash *endAlignments = stHash_construct2(NULL, (void(*)(void *)) stSortedSet_destruct); if(listOfEndAlignmentFiles != NULL) { loadEndAlignments(flower, endAlignments, listOfEndAlignmentFiles); } - computeMissingEndAlignments(sM, flower, endAlignments, spanningTrees, maxSequenceLength, + computeMissingEndAlignments(sM, pairModels, flower, endAlignments, spanningTrees, maxSequenceLength, useProgressiveMerging, gapGamma, pairwiseAlignmentBandingParameters); return makeFlowerAlignment2(flower, endAlignments, pruneOutStubAlignments); } diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index 905e112c5..55a75ea21 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -941,15 +941,6 @@ const char *baseAligner_toString(BaseAligner engine) { return "unknown"; } -static void addDistanceNodes(stHash *nodes, stTree *tree) { - if (stTree_getLabel(tree) != NULL) { - stHash_insert(nodes, (void *)stTree_getLabel(tree), tree); - } - for (int64_t i = 0; i < stTree_getChildNumber(tree); i++) { - addDistanceNodes(nodes, stTree_getChild(tree, i)); - } -} - /* * and its children .. (see * PoaParameters). A count rather than probing for the children, because asking CactusParams about @@ -983,46 +974,10 @@ static void rowModels_load(CactusParams *params, PoaParameters *pp) { pp->rowModelAbpt[m] = abpoaParameters_construct(params, name); free(name); } - if (cactusParams_has(params, 2, "reference", "reconstructionTree")) { - char *newick = cactusParams_get_string(params, 2, "reference", "reconstructionTree"); - pp->distanceTree = stTree_parseNewickString(newick); - pp->distanceNodes = stHash_construct3(stHash_stringKey, stHash_stringEqualKey, NULL, NULL); - addDistanceNodes(pp->distanceNodes, pp->distanceTree); - free(newick); - } + pp->distances = barDistances_constructFromCactusParams(params); st_logInfo("bar: %" PRIi64 " row models, distances to the nearest %s on the %s tree\n", pp->rowModelNo, - pp->rowModelToIngroups ? "ingroup" : "genome", pp->distanceTree != NULL ? "reconstruction" : "event"); -} - -/* - * Path length between two genomes, on the reconstruction tree when it names both, else on the event - * tree (whose branches are the ones scaled for alignment sensitivity). - */ -static double rowModels_distance(PoaParameters *pp, Event *event1, Event *event2) { - stTree *node1 = pp->distanceNodes == NULL ? NULL : stHash_search(pp->distanceNodes, (void *)event_getHeader(event1)); - stTree *node2 = pp->distanceNodes == NULL ? NULL : stHash_search(pp->distanceNodes, (void *)event_getHeader(event2)); - if (node1 != NULL && node2 != NULL) { - double d1 = 0.0; - for (stTree *a = node1; a != NULL; a = stTree_getParent(a)) { - double d2 = 0.0; - for (stTree *b = node2; b != NULL; b = stTree_getParent(b)) { - if (a == b) { - return d1 + d2; - } - d2 += stTree_getBranchLength(b); - } - d1 += stTree_getBranchLength(a); - } - } - Event *ancestor = eventTree_getCommonAncestor(event1, event2); - double d = 0.0; - for (Event *e = event1; e != ancestor; e = event_getParent(e)) { - d += event_getBranchLength(e); - } - for (Event *e = event2; e != ancestor; e = event_getParent(e)) { - d += event_getBranchLength(e); - } - return d; + pp->rowModelToIngroups ? "ingroup" : "genome", + barDistances_hasReconstructionTree(pp->distances) ? "reconstruction" : "event"); } /* @@ -1057,7 +1012,7 @@ static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { } for (int64_t j = 0; j < event_no; ++j) { for (int64_t k = j + 1; k < event_no; ++k) { - double d = rowModels_distance(pp, events[j], events[k]); + double d = barDistances_get(pp->distances, events[j], events[k]); nearest[j] = d < nearest[j] ? d : nearest[j]; nearest[k] = d < nearest[k] ? d : nearest[k]; if (!event_isOutgroup(events[k])) { @@ -1153,9 +1108,8 @@ void poaParameters_destruct(PoaParameters *poaParameters) { } free(poaParameters->rowModelAbpt); free(poaParameters->rowModelMaxDistance); - if (poaParameters->distanceTree != NULL) { - stHash_destruct(poaParameters->distanceNodes); - stTree_destruct(poaParameters->distanceTree); + if (poaParameters->distances != NULL) { + barDistances_destruct(poaParameters->distances); } #ifdef HAVE_MINIPOA if (poaParameters->mpt != NULL) { diff --git a/bar/inc/endAligner.h b/bar/inc/endAligner.h index 850b2989a..4d57e2174 100644 --- a/bar/inc/endAligner.h +++ b/bar/inc/endAligner.h @@ -51,6 +51,47 @@ stSortedSet *makeEndAlignment(StateMachine *sM, End *end, int64_t spanningTrees, bool useProgressiveMerging, float gapGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters); +/* + * Path lengths between genomes, for choosing scoring models by how diverged two genomes are: on + * (the unscaled tree, as ancestral reconstruction uses it) when it names + * both genomes, else on the event tree, whose branches are the ones scaled for alignment sensitivity. + */ +typedef struct _barDistances BarDistances; + +BarDistances *barDistances_constructFromCactusParams(CactusParams *params); + +double barDistances_get(BarDistances *distances, Event *event1, Event *event2); + +/* + * Whether the distances are on the reconstruction tree (else the event tree). + */ +bool barDistances_hasReconstructionTree(BarDistances *distances); + +void barDistances_destruct(BarDistances *distances); + +/* + * Per-pair pecan models. switches them on, and its children .. + * each carry a maxDistance and an hmm (a five-state pair HMM in the JSON hmm_jsonParse reads), + * in ascending maxDistance. Each pair of sequences is aligned with the first model whose maxDistance is at + * least the distance between their genomes (see BarDistances). Pairs from one genome (paralogs, whose + * divergence the tree does not give) and pairs beyond the last model keep the default machine. + */ +typedef struct _pecanPairModels PecanPairModels; + +/* + * NULL when the config has no pair models. + */ +PecanPairModels *pecanPairModels_constructFromCactusParams(CactusParams *params); + +void pecanPairModels_destruct(PecanPairModels *pairModels); + +/* + * As makeEndAlignment, aligning each pair of sequences with its model from pairModels (may be NULL). + */ +stSortedSet *makeEndAlignmentWithPairModels(StateMachine *sM, PecanPairModels *pairModels, End *end, int64_t spanningTrees, + int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma, + PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters); + /* * Writes an end alignment to the given file. */ diff --git a/bar/inc/flowerAligner.h b/bar/inc/flowerAligner.h index 44d663c94..e83ab85e8 100644 --- a/bar/inc/flowerAligner.h +++ b/bar/inc/flowerAligner.h @@ -15,6 +15,7 @@ #define FLOWER_ALIGNER_H_ #include "pairwiseAligner.h" +#include "endAligner.h" /* * Constructs an alignment for the flower by constructing an alignment for each end @@ -28,9 +29,9 @@ stSortedSet *makeFlowerAlignment(StateMachine *sM, Flower *flower, int64_t spann PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments); /* - * As above, but including alignments from disk. + * As above, but including alignments from disk, and with per-pair pecan models (pairModels may be NULL). */ -stSortedSet *makeFlowerAlignment3(StateMachine *sM, Flower *flower, stList *listOfEndAlignmentFiles, int64_t spanningTrees, +stSortedSet *makeFlowerAlignment3(StateMachine *sM, PecanPairModels *pairModels, Flower *flower, stList *listOfEndAlignmentFiles, int64_t spanningTrees, int64_t maxSequenceLength, bool useProgressiveMerging, float gapGamma, PairwiseAlignmentParameters *pairwiseAlignmentBandingParameters, bool pruneOutStubAlignments); diff --git a/bar/inc/poaBarAligner.h b/bar/inc/poaBarAligner.h index 089ab5d46..c77f1d44d 100644 --- a/bar/inc/poaBarAligner.h +++ b/bar/inc/poaBarAligner.h @@ -123,8 +123,7 @@ typedef struct _PoaParameters { bool rowModelToIngroups; double *rowModelMaxDistance; abpoa_para_t **rowModelAbpt; - stTree *distanceTree; - stHash *distanceNodes; // node label -> stTree node of distanceTree + BarDistances *distances; // NULL when there are no row models } PoaParameters; /** diff --git a/submodules/cPecan b/submodules/cPecan index 55920efb7..ccdbc32ca 160000 --- a/submodules/cPecan +++ b/submodules/cPecan @@ -1 +1 @@ -Subproject commit 55920efb7ad7bc9e5da9703668e93bdd0585e353 +Subproject commit ccdbc32caa6f1768867a23020dd8034015a2bef3 From 4a57b45ac2077774cc050961f39fd01b0d3bdb4b Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Fri, 2 Oct 2026 08:56:19 -0400 Subject: [PATCH 06/10] bar: pin cPecan at its deterministic multiple aligner cPecan 2ec01d1 makes the multiple aligner a function of its inputs alone: the tie-breaking jitter and the spanning-tree tie-break are hashes rather than draws from sonLib's generator, and nothing is ordered by address. Cactus already reseeds the generator per flower, so its runs were reproducible before. What changes is that an end's alignment no longer depends on how many draws the ends aligned before it in the flower used -- which a change of scoring model alters -- so comparing two configurations is not confounded by it. On evolver mammals with stock pecan, two runs are identical on every genome; against the previous pin, 17 of 4.5 million leaf-pair alignments change, and accuracy is the same to four decimals. Co-Authored-By: Claude Opus 5.5 --- submodules/cPecan | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/submodules/cPecan b/submodules/cPecan index ccdbc32ca..2ec01d160 160000 --- a/submodules/cPecan +++ b/submodules/cPecan @@ -1 +1 @@ -Subproject commit ccdbc32caa6f1768867a23020dd8034015a2bef3 +Subproject commit 2ec01d160bd4a558ef213a168a57748b82c9313d From 6b8eacfb0511ddd2007f76b897cfca2af4c32b67 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Mon, 5 Oct 2026 14:10:11 -0400 Subject: [PATCH 07/10] bar: row and pair models for named pairs of genomes A row model or pair model can now name the pair of genomes it is for (genomes="A B"), as one trained on them would, in place of or as well as a maxDistance. A pecan pair of sequences uses the model for their two genomes if there is one, else the first whose maxDistance covers their distance, as before. An abPOA row's partner is now kept along with its distance -- the nearest other genome, or ingroup -- and the row uses the model for it and its partner, else the distance bins. The models with a maxDistance must still come in ascending order of it. Configs with only maxDistance models behave as they did. Co-Authored-By: Claude Opus 5.5 --- bar/impl/endAligner.c | 68 ++++++++++++++++++++++++----- bar/impl/poaBarAligner.c | 94 +++++++++++++++++++++++++++++----------- bar/inc/endAligner.h | 21 +++++++-- bar/inc/poaBarAligner.h | 25 ++++++----- 4 files changed, 156 insertions(+), 52 deletions(-) diff --git a/bar/impl/endAligner.c b/bar/impl/endAligner.c index 0abf14b68..7365e7840 100644 --- a/bar/impl/endAligner.c +++ b/bar/impl/endAligner.c @@ -5,6 +5,7 @@ */ #include "endAligner.h" +#include #include "multipleAligner.h" #include "adjacencySequences.h" #include "pairwiseAligner.h" @@ -117,9 +118,31 @@ void barDistances_destruct(BarDistances *distances) { free(distances); } +void barModel_readGenomes(CactusParams *params, const char *section, const char *name, char **genomes) { + genomes[0] = genomes[1] = NULL; + if (!cactusParams_has(params, 4, "bar", section, name, "genomes")) { + return; + } + char *value = cactusParams_get_string(params, 4, "bar", section, name, "genomes"); + stList *tokens = stString_split(value); + if (stList_length(tokens) != 2) { + st_errAbort("bar model %s/%s has genomes=\"%s\", which is not two genome names", section, name, value); + } + genomes[0] = stString_copy(stList_get(tokens, 0)); + genomes[1] = stString_copy(stList_get(tokens, 1)); + stList_destruct(tokens); + free(value); +} + +bool barModel_isForGenomes(char **genomes, const char *a, const char *b) { + return genomes[0] != NULL && ((strcmp(genomes[0], a) == 0 && strcmp(genomes[1], b) == 0) || + (strcmp(genomes[0], b) == 0 && strcmp(genomes[1], a) == 0)); +} + struct _pecanPairModels { int64_t modelNo; - double *maxDistance; + double *maxDistance; // NAN for a model with none + char **genomes; // model m's two genomes at 2m and 2m+1, NULL for a model with none StateMachine **stateMachines; BarDistances *distances; }; @@ -133,13 +156,24 @@ PecanPairModels *pecanPairModels_constructFromCactusParams(CactusParams *params) PecanPairModels *pairModels = st_calloc(1, sizeof(PecanPairModels)); pairModels->modelNo = modelNo; pairModels->maxDistance = st_malloc(sizeof(double) * modelNo); + pairModels->genomes = st_malloc(sizeof(char *) * 2 * modelNo); pairModels->stateMachines = st_malloc(sizeof(StateMachine *) * modelNo); + int64_t keyed = 0; + double lastMaxDistance = -INFINITY; for (int64_t m = 0; m < modelNo; m++) { char *name = stString_print("pairModel%" PRIi64, m); - pairModels->maxDistance[m] = cactusParams_get_float(params, 4, "bar", "pecan", name, "maxDistance"); - if (m > 0 && pairModels->maxDistance[m] < pairModels->maxDistance[m - 1]) { - st_errAbort("pecan pair models must be in ascending maxDistance, but %s has %f after %f", name, - pairModels->maxDistance[m], pairModels->maxDistance[m - 1]); + barModel_readGenomes(params, "pecan", name, pairModels->genomes + 2 * m); + keyed += pairModels->genomes[2 * m] != NULL; + pairModels->maxDistance[m] = NAN; + if (cactusParams_has(params, 4, "bar", "pecan", name, "maxDistance")) { + pairModels->maxDistance[m] = cactusParams_get_float(params, 4, "bar", "pecan", name, "maxDistance"); + if (pairModels->maxDistance[m] < lastMaxDistance) { + st_errAbort("pecan pair models must be in ascending maxDistance, but %s has %f after %f", name, + pairModels->maxDistance[m], lastMaxDistance); + } + lastMaxDistance = pairModels->maxDistance[m]; + } else if (pairModels->genomes[2 * m] == NULL) { + st_errAbort("pecan pair model %s has neither genomes nor a maxDistance", name); } char *json = cactusParams_get_string(params, 4, "bar", "pecan", name, "hmm"); Hmm *hmm = hmm_jsonParse(json, strlen(json)); @@ -153,8 +187,8 @@ PecanPairModels *pecanPairModels_constructFromCactusParams(CactusParams *params) free(name); } pairModels->distances = barDistances_constructFromCactusParams(params); - st_logInfo("bar: %" PRIi64 " pecan pair models, distances on the %s tree\n", modelNo, - barDistances_hasReconstructionTree(pairModels->distances) ? "reconstruction" : "event"); + st_logInfo("bar: %" PRIi64 " pecan pair models (%" PRIi64 " for pairs of genomes), distances on the %s tree\n", + modelNo, keyed, barDistances_hasReconstructionTree(pairModels->distances) ? "reconstruction" : "event"); return pairModels; } @@ -165,6 +199,10 @@ void pecanPairModels_destruct(PecanPairModels *pairModels) { for (int64_t m = 0; m < pairModels->modelNo; m++) { stateMachine_destruct(pairModels->stateMachines[m]); } + for (int64_t i = 0; i < 2 * pairModels->modelNo; i++) { + free(pairModels->genomes[i]); + } + free(pairModels->genomes); free(pairModels->stateMachines); free(pairModels->maxDistance); barDistances_destruct(pairModels->distances); @@ -172,7 +210,7 @@ void pecanPairModels_destruct(PecanPairModels *pairModels) { } /* - * cPecan's per-pair callback: the model for the distance between the two sequences' genomes. + * cPecan's per-pair callback: the model for the two sequences' genomes, else for the distance between them. */ typedef struct { PecanPairModels *pairModels; @@ -185,10 +223,16 @@ static StateMachine *choosePairModel(int64_t seqX, int64_t seqY, void *extraArgs if (eventX == eventY) { return NULL; } - double distance = barDistances_get(choice->pairModels->distances, eventX, eventY); - for (int64_t m = 0; m < choice->pairModels->modelNo; m++) { - if (distance <= choice->pairModels->maxDistance[m]) { - return choice->pairModels->stateMachines[m]; + PecanPairModels *pairModels = choice->pairModels; + for (int64_t m = 0; m < pairModels->modelNo; m++) { + if (barModel_isForGenomes(pairModels->genomes + 2 * m, event_getHeader(eventX), event_getHeader(eventY))) { + return pairModels->stateMachines[m]; + } + } + double distance = barDistances_get(pairModels->distances, eventX, eventY); + for (int64_t m = 0; m < pairModels->modelNo; m++) { + if (distance <= pairModels->maxDistance[m]) { // false for a NAN maxDistance + return pairModels->stateMachines[m]; } } return NULL; diff --git a/bar/impl/poaBarAligner.c b/bar/impl/poaBarAligner.c index 55a75ea21..4995452ae 100644 --- a/bar/impl/poaBarAligner.c +++ b/bar/impl/poaBarAligner.c @@ -963,27 +963,58 @@ static void rowModels_load(CactusParams *params, PoaParameters *pp) { free(rule); } pp->rowModelMaxDistance = st_malloc(sizeof(double) * pp->rowModelNo); + pp->rowModelGenomes = st_malloc(sizeof(char *) * 2 * pp->rowModelNo); pp->rowModelAbpt = st_malloc(sizeof(abpoa_para_t *) * pp->rowModelNo); + int64_t keyed = 0; + double lastMaxDistance = -INFINITY; for (int64_t m = 0; m < pp->rowModelNo; ++m) { char *name = stString_print("rowModel%" PRIi64, m); - pp->rowModelMaxDistance[m] = cactusParams_get_float(params, 4, "bar", "poa", name, "maxDistance"); - if (m > 0 && pp->rowModelMaxDistance[m] < pp->rowModelMaxDistance[m-1]) { - st_errAbort("bar row models must be in ascending maxDistance, but %s has %f after %f", name, - pp->rowModelMaxDistance[m], pp->rowModelMaxDistance[m-1]); + barModel_readGenomes(params, "poa", name, pp->rowModelGenomes + 2 * m); + keyed += pp->rowModelGenomes[2 * m] != NULL; + pp->rowModelMaxDistance[m] = NAN; + if (cactusParams_has(params, 4, "bar", "poa", name, "maxDistance")) { + pp->rowModelMaxDistance[m] = cactusParams_get_float(params, 4, "bar", "poa", name, "maxDistance"); + if (pp->rowModelMaxDistance[m] < lastMaxDistance) { + st_errAbort("bar row models must be in ascending maxDistance, but %s has %f after %f", name, + pp->rowModelMaxDistance[m], lastMaxDistance); + } + lastMaxDistance = pp->rowModelMaxDistance[m]; + } else if (pp->rowModelGenomes[2 * m] == NULL) { + st_errAbort("bar row model %s has neither genomes nor a maxDistance", name); } pp->rowModelAbpt[m] = abpoaParameters_construct(params, name); free(name); } pp->distances = barDistances_constructFromCactusParams(params); - st_logInfo("bar: %" PRIi64 " row models, distances to the nearest %s on the %s tree\n", pp->rowModelNo, - pp->rowModelToIngroups ? "ingroup" : "genome", + st_logInfo("bar: %" PRIi64 " row models (%" PRIi64 " for pairs of genomes), distances to the nearest %s on the %s tree\n", + pp->rowModelNo, keyed, pp->rowModelToIngroups ? "ingroup" : "genome", barDistances_hasReconstructionTree(pp->distances) ? "reconstruction" : "event"); } /* - * The row model for each row of an end's alignment (-1 for the default), from the distance between - * the row's genome and the nearest other genome among the rows (the nearest other ingroup genome with - * rowModelToIngroups). NULL when there are no row models. + * The row model for a genome and its partner at the given distance (no partner: NULL, and INFINITY): the + * model for the pair, else the first whose maxDistance covers the distance, else -1, the default. + */ +static int rowModels_choose(PoaParameters *pp, Event *event, Event *partner, double distance) { + if (partner != NULL) { + for (int64_t m = 0; m < pp->rowModelNo; ++m) { + if (barModel_isForGenomes(pp->rowModelGenomes + 2 * m, event_getHeader(event), event_getHeader(partner))) { + return (int)m; + } + } + } + for (int64_t m = 0; m < pp->rowModelNo; ++m) { + if (distance <= pp->rowModelMaxDistance[m]) { // false for a NAN maxDistance + return (int)m; + } + } + return -1; +} + +/* + * The row model for each row of an end's alignment (-1 for the default), for the row's genome and its + * partner, the nearest other genome among the rows (the nearest other ingroup genome with + * rowModelToIngroups): see rowModels_choose. NULL when there are no row models. */ static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { if (pp == NULL || pp->rowModelNo == 0) { @@ -1004,22 +1035,30 @@ static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { } row_event[i] = j; } - // nearest other genome, and nearest other ingroup genome + // nearest other genome, and nearest other ingroup genome, and which they are (-1 for none; ties go to + // the first met) double *nearest = st_malloc(sizeof(double) * event_no); double *nearestIngroup = st_malloc(sizeof(double) * event_no); + int64_t *partner = st_malloc(sizeof(int64_t) * event_no); + int64_t *ingroupPartner = st_malloc(sizeof(int64_t) * event_no); for (int64_t j = 0; j < event_no; ++j) { nearest[j] = nearestIngroup[j] = INFINITY; + partner[j] = ingroupPartner[j] = -1; } for (int64_t j = 0; j < event_no; ++j) { for (int64_t k = j + 1; k < event_no; ++k) { double d = barDistances_get(pp->distances, events[j], events[k]); - nearest[j] = d < nearest[j] ? d : nearest[j]; - nearest[k] = d < nearest[k] ? d : nearest[k]; - if (!event_isOutgroup(events[k])) { - nearestIngroup[j] = d < nearestIngroup[j] ? d : nearestIngroup[j]; + if (d < nearest[j]) { + nearest[j] = d, partner[j] = k; + } + if (d < nearest[k]) { + nearest[k] = d, partner[k] = j; } - if (!event_isOutgroup(events[j])) { - nearestIngroup[k] = d < nearestIngroup[k] ? d : nearestIngroup[k]; + if (!event_isOutgroup(events[k]) && d < nearestIngroup[j]) { + nearestIngroup[j] = d, ingroupPartner[j] = k; + } + if (!event_isOutgroup(events[j]) && d < nearestIngroup[k]) { + nearestIngroup[k] = d, ingroupPartner[k] = j; } } } @@ -1027,25 +1066,26 @@ static int *rowModels_assign(PoaParameters *pp, Cap **caps, int64_t row_no) { // and the rest (only one ingroup genome, or none, in this alignment) to the nearest genome of any kind if (pp->rowModelToIngroups) { for (int64_t j = 0; j < event_no; ++j) { - if (nearestIngroup[j] != INFINITY) { - nearest[j] = nearestIngroup[j]; + if (ingroupPartner[j] >= 0) { + nearest[j] = nearestIngroup[j], partner[j] = ingroupPartner[j]; } } } + int *eventModels = st_malloc(sizeof(int) * event_no); + for (int64_t j = 0; j < event_no; ++j) { + eventModels[j] = rowModels_choose(pp, events[j], partner[j] >= 0 ? events[partner[j]] : NULL, nearest[j]); + } int *models = st_malloc(sizeof(int) * row_no); for (int64_t i = 0; i < row_no; ++i) { - models[i] = -1; - for (int64_t m = 0; m < pp->rowModelNo; ++m) { - if (nearest[row_event[i]] <= pp->rowModelMaxDistance[m]) { - models[i] = (int)m; - break; - } - } + models[i] = eventModels[row_event[i]]; } free(events); free(row_event); free(nearest); free(nearestIngroup); + free(partner); + free(ingroupPartner); + free(eventModels); return models; } @@ -1108,6 +1148,10 @@ void poaParameters_destruct(PoaParameters *poaParameters) { } free(poaParameters->rowModelAbpt); free(poaParameters->rowModelMaxDistance); + for (int64_t i = 0; i < 2 * poaParameters->rowModelNo; ++i) { + free(poaParameters->rowModelGenomes[i]); + } + free(poaParameters->rowModelGenomes); if (poaParameters->distances != NULL) { barDistances_destruct(poaParameters->distances); } diff --git a/bar/inc/endAligner.h b/bar/inc/endAligner.h index 4d57e2174..6d2fe1ba3 100644 --- a/bar/inc/endAligner.h +++ b/bar/inc/endAligner.h @@ -69,12 +69,25 @@ bool barDistances_hasReconstructionTree(BarDistances *distances); void barDistances_destruct(BarDistances *distances); +/* + * A scoring model's optional genomes attribute (
): its two genome names, copied + * into genomes[0] and genomes[1], or NULLs when it has none. + */ +void barModel_readGenomes(CactusParams *params, const char *section, const char *name, char **genomes); + +/* + * Whether a model's genomes (as barModel_readGenomes gives them) are the given pair, either way round. + */ +bool barModel_isForGenomes(char **genomes, const char *a, const char *b); + /* * Per-pair pecan models. switches them on, and its children .. - * each carry a maxDistance and an hmm (a five-state pair HMM in the JSON hmm_jsonParse reads), - * in ascending maxDistance. Each pair of sequences is aligned with the first model whose maxDistance is at - * least the distance between their genomes (see BarDistances). Pairs from one genome (paralogs, whose - * divergence the tree does not give) and pairs beyond the last model keep the default machine. + * each carry an hmm (a five-state pair HMM in the JSON hmm_jsonParse reads) and a genomes (two + * genome names: a model for that pair of genomes, as one trained on them) or a maxDistance, or both. Each pair + * of sequences is aligned with the model for its two genomes if there is one, else with the first, in order, + * whose maxDistance is at least the distance between them (see BarDistances; the ones that have a maxDistance + * must come in ascending order of it). Pairs from one genome (paralogs, whose divergence the tree does not + * give) and pairs with no model keep the default machine. */ typedef struct _pecanPairModels PecanPairModels; diff --git a/bar/inc/poaBarAligner.h b/bar/inc/poaBarAligner.h index c77f1d44d..5340bf99c 100644 --- a/bar/inc/poaBarAligner.h +++ b/bar/inc/poaBarAligner.h @@ -108,20 +108,23 @@ typedef struct _PoaParameters { char *referenceEvent; /* * Per-row scoring, abPOA only. switches it on, and its - * children .. each carry a maxDistance and their own - * partialOrderAlignmentSubMatrix and gap penalties, in ascending maxDistance. A row is scored - * with the first whose maxDistance is at least the distance from its genome to the nearest other - * genome in the same alignment; rows beyond the last, or with no other genome, keep abpt. - * Distances are on (the unscaled tree), falling back to the event - * tree for a genome it does not name. rowModelNo is 0 when the config has no row models. - * With rowModelToIngroups (partialOrderAlignmentRowModelDistance="ingroup") the distance is to the - * nearest other ingroup genome instead, so an outgroup is scored for how far it is from what is - * being reconstructed rather than from another outgroup; a row with no other ingroup genome in the - * alignment falls back to the nearest genome of any kind. + * children .. each carry their own partialOrderAlignmentSubMatrix and gap + * penalties, and a genomes (two genome names: a model for that pair of genomes, as one trained on them) + * or a maxDistance (a model for a distance), or both. A row's partner is the nearest other genome in + * the same alignment, and the row is scored with the model for it and its partner if there is one, else + * with the first, in order, whose maxDistance is at least the distance between them (the ones that have + * a maxDistance must come in ascending order of it). Rows with no model, or no other genome, keep abpt. + * Distances are on (the unscaled tree), falling back to the event tree + * for a genome it does not name. rowModelNo is 0 when the config has no row models. With + * rowModelToIngroups (partialOrderAlignmentRowModelDistance="ingroup") the partner is the nearest other + * ingroup genome instead, so an outgroup is scored for how far it is from what is being reconstructed + * rather than from another outgroup; a row with no other ingroup genome in the alignment falls back to + * the nearest genome of any kind. */ int64_t rowModelNo; bool rowModelToIngroups; - double *rowModelMaxDistance; + double *rowModelMaxDistance; // NAN for a model with none + char **rowModelGenomes; // model m's two genomes at 2m and 2m+1, NULL for a model with none abpoa_para_t **rowModelAbpt; BarDistances *distances; // NULL when there are no row models } PoaParameters; From 3f089e45d0253ab649eb53cb63fda64e61ab6d43 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Mon, 5 Oct 2026 14:10:20 -0400 Subject: [PATCH 08/10] Score bar with the models the blast step trains per pair of genomes With , the blast step's last-train run on each pair of genomes it aligns (see , which it now runs whether or not lastz uses the models) also makes the pair's model for bar: abPOA's scoring, on HOXD70's scale like the shipped matrix, with the trained gap as the first piece and a long-gap piece made as the pangenome makes one; and pecan's five-state pair HMM, from the trained probabilities. The consolidated job writes them into bar's config, in place of any row or pair models there: abPOA scores each row with the model for its genome and the nearest ingroup (always a pair blast trained), and pecan each pair of sequences with the model for its two genomes. Pairs blast did not train, as pecan's outgroup pairs, take the model trained nearest their distance in log terms on the reconstruction tree. In one cactus run the models pass straight from blast to consolidated; cactus-blast leaves them next to its PAF, as .models.json, for cactus-align. On evolver mammals (pooled leaf-pair precision/recall), against the shipped config's 0.816/0.739: bar abPOA trained, lastz HOXD70 0.824/0.736 (hand-built per-row table: 0.827/0.734) bar pecan trained, lastz HOXD70 0.841/0.754 (hand-built per-pair table: 0.841/0.755) lastz and abPOA trained 0.836/0.742 lastz and pecan trained 0.855/0.760 (hand-built pecan table: 0.855/0.760) cactus-prepare's steps and a single cactus run both applied every ancestor's models; they agreed to within 0.0014. Co-Authored-By: Claude Opus 5.5 --- src/cactus/blast/cactus_blast.py | 15 +- src/cactus/cactus_progressive_config.xml | 14 +- src/cactus/paf/last_scoring.py | 165 ++++++++++++++++++- src/cactus/paf/local_alignment.py | 83 +++++++--- src/cactus/pipeline/cactus_workflow.py | 28 +++- src/cactus/progressive/cactus_progressive.py | 26 ++- src/cactus/setup/cactus_align.py | 18 +- 7 files changed, 294 insertions(+), 55 deletions(-) diff --git a/src/cactus/blast/cactus_blast.py b/src/cactus/blast/cactus_blast.py index a1c809f46..94bfd452d 100644 --- a/src/cactus/blast/cactus_blast.py +++ b/src/cactus/blast/cactus_blast.py @@ -23,6 +23,7 @@ from cactus.progressive.cactus_prepare import human2bytesN from cactus.paf.local_alignment import sanitize_then_make_paf_alignments +from cactus.paf.last_scoring import bar_train_enabled, trained_models_path from toil.job import Job from toil.common import Toil @@ -157,10 +158,18 @@ def runCactusBlastOnly(options): input_seq_id_map[genome] = toil.importFile(seq) paf_id = toil.start(Job.wrapJobFn(sanitize_then_make_paf_alignments, NXNewick().writeString(spanning_tree), - input_seq_id_map, options.root, config_node, options.outputFile, walltime=cactus_walltime())) - - # export the alignments + input_seq_id_map, options.root, config_node, options.outputFile, + return_trained_models=bar_train_enabled(config_node), + walltime=cactus_walltime())) + + # export the alignments, and with the scoring models trained for them, where + # cactus-align looks for them (see trained_models_path) + models_id = None + if isinstance(paf_id, dict): + paf_id, models_id = paf_id['paf'], paf_id['models'] toil.exportFile(paf_id, makeURL(options.outputFile)) + if models_id is not None: + toil.exportFile(models_id, makeURL(trained_models_path(options.outputFile))) if __name__ == '__main__': diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index 8f4a4389b..f952d2edd 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -92,9 +92,10 @@ + lastzScoreModel table, else HOXD70. With the same models score bar too. It does not apply + to KegAlign, minimap2 or FastGA. On evolver mammals it raised pooled leaf-pair precision from 0.816 to 0.830 + and recall from 0.740 to 0.744, for about 0.3% of the 577-way's lastz CPU; set it to 0 for lastz's HOXD70 + everywhere --> @@ -298,6 +299,12 @@ So an ancestor of two ancestors can reach 0.5 however deep it is. 0 turns it off. Over the 576 ancestors of the VGP 577-way tree with branchScale 2 the median is 0.17 and 15% are 0.5 or more: MammalsAnc0 0.58, MammalsAnc1 0.50, RayFinnedFishesAnc00 0.81, RayFinnedFishesAnc01 0.53, PrimatesAnc0 0.46, HominidaeAnc4 0.001. --> + @@ -310,6 +317,7 @@ baseAligner="abpoa" deepAligner="pecan" deepAlignerDivergence="0" + trainedModels="0" minimumBlockDegree="2" minimumIngroupDegree="1" minimumOutgroupDegree="0" diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index 82f8fe1ba..dd511ebd4 100644 --- a/src/cactus/paf/last_scoring.py +++ b/src/cactus/paf/last_scoring.py @@ -16,6 +16,9 @@ import os import re import math +import json +import xml.etree.ElementTree as ET +from sonLib.nxnewick import NXNewick from cactus.shared.common import cactus_call, getOptionalAttrib from cactus.shared.common import cactus_clamp_memory from cactus.shared.common import findRequiredNode @@ -217,7 +220,8 @@ def excess(lam): def parse_train_model(train_file_path): """ last-train's final model at the scale it trains at (parse_train_file reads the integer one it ends with, whose scores are too coarse to rescale): matrix rows and columns A, C, G, T, the gap existence and - extension costs, the query's base frequencies and the identity of the last pass's alignments. The + extension costs, the substitution and gap probabilities they come from (as the last pass counted them), + the query's base frequencies, and the identity of the last pass's alignments. The training scale differs from one run to the next (last-train raises it 10% at a time when rounding gets too coarse), so the numbers are only comparable between runs once rescaled. """ with open(train_file_path) as train_file: @@ -235,19 +239,27 @@ def values(key): if not found: raise RuntimeError('no {} in {}'.format(key, train_file_path)) return found + probs = {} + for line in txt.split('# probability matrix (query letters = columns, reference letters = rows):')[-1].splitlines()[2:6]: + toks = line.split() + probs[toks[1]] = [float(x) for x in toks[2:6]] # the last identity printed is the one the final integer scores imply, whose rounding can move it by # several points; the one before is the last pass's, from the alignments it counted - return {'matrix': [rows[b] for b in 'ACGT'], + return {'matrix': [rows[b] for b in 'ACGT'], 'probs': [probs[b] for b in 'ACGT'], 'open': int(values('delExistCost')[-1]), 'extend': int(values('delExtendCost')[-1]), + 'open_prob': float(values('delOpenProb')[-1]), 'extend_prob': float(values('delExtendProb')[-1]), 'freqs': [float(x) / 100 for x in re.findall(r'# qry letter %: (.*)', txt)[-1].split()], 'identity': float(values('substitution percent identity')[-2])} -def lastz_scores_from_train(train_file_path, trained_gaps=True): - """ the model last-train fitted, as lastz scores on HOXD70's scale: a dict of a 's +def hoxd70_scale(model): + """ the factor that puts a parse_train_model model on HOXD70's scale (lambda) """ + return implied_lambda(model['matrix'], model['freqs']) / implied_lambda(HOXD70, HOXD70_FREQS) + +def lastz_scores_from_model(model, trained_gaps=True): + """ a parse_train_model model as lastz scores on HOXD70's scale: a dict of a 's matrix, gapOpen and gapExtend (see local_alignment.write_lastz_scores), plus the trained identity. Without trained_gaps the gap penalties are HOXD70's, 400 and 30, as lastz has them by default. """ - model = parse_train_model(train_file_path) - scale = implied_lambda(model['matrix'], model['freqs']) / implied_lambda(HOXD70, HOXD70_FREQS) + scale = hoxd70_scale(model) matrix = [int(round(s * scale)) for row in model['matrix'] for s in row] if any(matrix[5 * i] <= 0 for i in range(4)): raise RuntimeError('trained matrix {} has a non-positive match score'.format(matrix)) @@ -256,6 +268,147 @@ def lastz_scores_from_train(train_file_path, trained_gaps=True): 'gapExtend': str(max(1, int(round(model['extend'] * scale)))) if trained_gaps else '30', 'identity': model['identity']} +def lastz_scores_from_train(train_file_path, trained_gaps=True): + """ lastz_scores_from_model for a last-train output file """ + return lastz_scores_from_model(parse_train_model(train_file_path), trained_gaps) + +def poa_scores_from_model(model, open2_factor, extend2_factor): + """ a parse_train_model model as abPOA's scoring attributes, on HOXD70's scale, which is also that of the + shipped matrix: the substitution matrix (with N scored as apply_scores_to_config scores it), the + trained gap as the first piece of abPOA's convex gap penalty, and a second, long-gap piece made the way + apply_long_gap makes one for the pangenome, open2_factor times dearer to open and extend2_factor times + cheaper to extend (, ...Extension2Factor). """ + scale = hoxd70_scale(model) + m = [[int(round(s * scale)) for s in row] for row in model['matrix']] + gap_open, gap_extend = int(round(model['open'] * scale)), max(1, int(round(model['extend'] * scale))) + mismatch = max(m[i][j] for i in range(4) for j in range(4) if i != j) + match = min(m[i][i] for i in range(4)) + sub = [] + for i in range(4): + sub += m[i] + [mismatch] + sub += [mismatch] * 4 + [match] + return {'partialOrderAlignmentSubMatrix': ' '.join(str(x) for x in sub), + 'partialOrderAlignmentGapOpenPenalty1': str(gap_open), + 'partialOrderAlignmentGapExtensionPenalty1': str(gap_extend), + 'partialOrderAlignmentGapOpenPenalty2': str(gap_open * open2_factor), + 'partialOrderAlignmentGapExtensionPenalty2': str(max(1, int(round(gap_extend / extend2_factor))))} + +# cPecan's built-in five-state model (stateMachine5_construct), whose short:long gap ratios pecan_hmm_from_model keeps +PECAN_DEFAULT = {'short_open': 0.0129868352330243, 'short_extend': 0.7126062401851738, 'short_switch': 0.0073673675173412815, + 'long_open': 0.001821479941473, 'long_extend': 0.99656342579062, 'long_switch': 0.001821479941473} + +def pecan_hmm_from_model(model): + """ a parse_train_model model as cPecan's symmetric five-state pair HMM (match, short gap x/y, long gap x/y), + in the JSON hmm_jsonParse reads. Match emissions are last-train's substitution probabilities, symmetrised, + and gap emissions the letter frequencies. last-train fits one affine gap where pecan has a short and a long + one: its gap-open probability is split between them in the ratio of pecan's own (trained) defaults, so the + total is last-train's; short gaps extend as last-train's do, long ones as pecan's; and switching between an + insertion and a deletion scales with the open probability. """ + d = PECAN_DEFAULT + f = model['open_prob'] / (d['short_open'] + d['long_open']) + short_open, long_open = d['short_open'] * f, d['long_open'] * f + short_extend, long_extend = model['extend_prob'], d['long_extend'] + short_switch = min(d['short_switch'] * f, (1 - short_extend) / 2) + long_switch = min(d['long_switch'] * f, (1 - long_extend) / 2) + M, SX, SY, LX, LY = range(5) + t = [[0.0] * 5 for _ in range(5)] + t[M][M] = 1 - 2 * (short_open + long_open) + t[M][SX] = t[M][SY] = short_open + t[M][LX] = t[M][LY] = long_open + for g, o, ext, sw in ((SX, SY, short_extend, short_switch), (SY, SX, short_extend, short_switch), + (LX, LY, long_extend, long_switch), (LY, LX, long_extend, long_switch)): + t[g][g], t[g][o], t[g][M] = ext, sw, 1 - ext - sw + if any(abs(sum(row) - 1) > 1e-9 or min(row) < 0 for row in t): + raise RuntimeError('trained gap probabilities open {} extend {} make no pecan model'.format( + model['open_prob'], model['extend_prob'])) + total = sum(sum(row) for row in model['probs']) + probs = [[(model['probs'][x][y] + model['probs'][y][x]) / (2 * total) for y in range(4)] for x in range(4)] + q = model['freqs'] + emissions = [] + for state in range(5): + for x in range(4): + for y in range(4): + emissions.append(probs[x][y] if state == M else (q[x] if state in (SX, LX) else q[y]) / 4) + return json.dumps({'type': 0, 'transitions': [v for row in t for v in row], 'emissions': emissions, 'likelihood': 0.0}) + +def scoring_models_from_train(train_file_path, config_node, trained_gaps=True): + """ everything one last-train run on a pair of genomes scores: the pair's model for lastz + (lastz_scores_from_model, with or without its trained gaps), for abPOA's rows (poa_scores_from_model, with + the long-gap factors) and for pecan (pecan_hmm_from_model), and the identity it was trained to """ + model = parse_train_model(train_file_path) + poa_node = findRequiredNode(findRequiredNode(config_node, 'bar'), 'poa') + return {'lastz': lastz_scores_from_model(model, trained_gaps), + 'poa': poa_scores_from_model(model, int(poa_node.attrib['partialOrderAlignmentTrainedGapOpen2Factor']), + int(poa_node.attrib['partialOrderAlignmentTrainedGapExtension2Factor'])), + 'pecan': pecan_hmm_from_model(model), + 'identity': model['identity']} + +def trained_models_path(paf_path): + """ where cactus-blast leaves the scoring models it trained for a PAF, for cactus-align to give bar """ + return paf_path + '.models.json' + +def bar_train_enabled(config_node): + """ whether bar scores with the models trained per pair of genomes in the blast step: """ + return getOptionalAttrib(findRequiredNode(config_node, 'bar'), 'trainedModels', typeFn=bool, default=False) + +def tree_distances(newick, names): + """ path lengths between the named nodes of a newick tree, by unordered pair (a frozenset), for those it has """ + tree = NXNewick().parseString(newick) + nodes = {tree.getName(n): n for n in tree.postOrderTraversal() if tree.getName(n) in names} + def ancestors(n): + path, depth = [n], [0.0] + while tree.hasParent(path[-1]): + parent = tree.getParent(path[-1]) + depth.append(depth[-1] + tree.getWeight(parent, path[-1])) + path.append(parent) + return dict(zip(path, depth)) + up = {name: ancestors(n) for name, n in nodes.items()} + distances = {} + for a in nodes: + for b in nodes: + if a < b: + common = min((d + up[b][n] for n, d in up[a].items() if n in up[b]), default=None) + if common is not None: + distances[frozenset((a, b))] = common + return distances + +def apply_trained_models_to_config(config_node, trained_models, reconstruction_tree=None): + """ give bar the models trained per pair of genomes in the blast step, in place of any it has: each as + a row model and a pair model for its pair of genomes, abPOA's measured to the nearest ingroup + (partialOrderAlignmentRowModelDistance="ingroup", so every row's pair is one blast trained). Pairs blast + did not train (pecan aligns outgroups to each other) take the model trained nearest their distance in log + terms: each model also covers the distances up to the geometric mean of its pair's and the next one's, on + reconstruction_tree, the tree bar measures distances on (the models of pairs it does not name only match + their own pair). trained_models is the list the blast step writes: dicts of genomes (two names), poa + and pecan. """ + bar_node = findRequiredNode(config_node, 'bar') + poa_node, pecan_node = findRequiredNode(bar_node, 'poa'), findRequiredNode(bar_node, 'pecan') + for node, prefix in (poa_node, 'rowModel'), (pecan_node, 'pairModel'): + for child in [c for c in node if c.tag.startswith(prefix)]: + node.remove(child) + names = {g for m in trained_models for g in m['genomes']} + distances = tree_distances(reconstruction_tree, names) if reconstruction_tree else {} + # in order of distance, the pairs the tree does not give last; a bin each for the distinct distances + def distance(m): + return distances.get(frozenset(m['genomes'])) + ordered = sorted(trained_models, key=lambda m: (distance(m) is None, distance(m) or 0.0, sorted(m['genomes']))) + known = sorted({distance(m) for m in ordered if distance(m) is not None}) + bound = {d: math.sqrt(d * known[i + 1]) if i + 1 < len(known) else 1e9 for i, d in enumerate(known)} + binned = set() + for i, m in enumerate(ordered): + d = distance(m) + attrib = {'genomes': ' '.join(m['genomes'])} + if d is not None and d not in binned: # the first model at a distance covers it for unnamed pairs + attrib['maxDistance'] = repr(bound[d]) + binned.add(d) + poa_node.append(ET.Element('rowModel{}'.format(i), dict(attrib, **m['poa']))) + pecan_node.append(ET.Element('pairModel{}'.format(i), dict(attrib, hmm=m['pecan']))) + poa_node.set('partialOrderAlignmentRowModels', str(len(ordered))) + poa_node.set('partialOrderAlignmentRowModelDistance', 'ingroup') + pecan_node.set('pecanPairModels', str(len(ordered))) + RealtimeLogger.info('bar scores with {} models trained per pair of genomes, {} also by distance'.format( + len(ordered), len(binned))) + def last_train(job, config, seq_order, seq_id_map, ref_name=None): """ run last_train on a pair of fasta files, using the first as the database. diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py index 8eaabd844..400af6bb9 100755 --- a/src/cactus/paf/local_alignment.py +++ b/src/cactus/paf/local_alignment.py @@ -14,13 +14,14 @@ from toil.realtimeLogger import RealtimeLogger from sonLib.bioio import newickTreeParser import glob +import json import os import shutil import math import copy from Bio import SeqIO from cactus.paf.paf import get_event_pairs, get_leaves, get_node, get_distances -from cactus.paf.last_scoring import lastz_scores_from_train +from cactus.paf.last_scoring import scoring_models_from_train, bar_train_enabled from cactus.shared.common import cactus_call, getOptionalAttrib, zip_gz, cactus_walltime from cactus.preprocessor.checkUniqueHeaders import sanitize_fasta_headers from cactus.preprocessor.unmasking import unmask_contigs_all @@ -105,10 +106,11 @@ def last_train_has_tolerance(): return False -def train_lastz_score_model(job, event_a, genome_a, event_b, genome_b, params): - """ A scoring model for lastz to align a pair of genomes with, trained on the pair with last-train, as - lastz_scores_from_train makes it -- or None, leaving the pair to the lastzScoreModel table or HOXD70, - when the smaller genome has under 500kb to train on or training fails. """ +def train_pair_scoring_models(job, event_a, genome_a, event_b, genome_b, params): + """ Scoring models for a pair of genomes, trained on the pair with last-train: for lastz to align them + with (), and for bar's abPOA rows and pecan pairs (), as + scoring_models_from_train makes them -- or None when the smaller genome has under 500kb to train on or + training fails, which leaves lastz the lastzScoreModel table or HOXD70, and bar its own settings. """ # lastdb is most of the work, so the smaller genome is the database (an ancestor's fasta is gzipped) genomes = [] for event, genome in (event_a, genome_a), (event_b, genome_b): @@ -140,18 +142,30 @@ def train_lastz_score_model(job, event_a, genome_a, event_b, genome_b, params): with open(os.path.join(work_dir, 'db.prj')) as prj: db_letters = int(next(line for line in prj if line.startswith('numofletters=')).split('=')[1]) if db_letters < 500000: - RealtimeLogger.warning("Not training lastz scores for {} and {}: {} has {} bases, too few to train on".format( + RealtimeLogger.warning("Not training scoring for {} and {}: {} has {} bases, too few to train on".format( event_a, event_b, db_event, db_letters)) return None cactus_call(parameters=train_cmd, outfile=train_file, work_dir=work_dir) - model = lastz_scores_from_train(train_file, trained_gaps) + models = scoring_models_from_train(train_file, params, trained_gaps) except Exception as e: - RealtimeLogger.warning("Training lastz scores for {} and {} failed, so the pair falls back to the " - "lastzScoreModel table or HOXD70: {}".format(event_a, event_b, e)) + RealtimeLogger.warning("Training scoring for {} and {} failed, so lastz falls back to the lastzScoreModel " + "table or HOXD70, and bar to its own settings: {}".format(event_a, event_b, e)) return None - RealtimeLogger.info("Trained lastz scores for {} and {} ({} against {}): identity {:.1f}%, matrix {}, gap {}+{}k".format( - event_a, event_b, query_event, db_event, model['identity'], model['matrix'], model['gapOpen'], model['gapExtend'])) - return model + RealtimeLogger.info("Trained scoring for {} and {} ({} against {}): identity {:.1f}%, lastz matrix {}, gap {}+{}k".format( + event_a, event_b, query_event, db_event, models['identity'], models['lastz']['matrix'], + models['lastz']['gapOpen'], models['lastz']['gapExtend'])) + return models + + +def write_trained_models(job, trained): + """ The models train_pair_scoring_models made for the pairs of genomes it could train, by lastz_pair_key, + as the JSON list cactus-align, and progressive cactus's consolidated, give bar (see + last_scoring.apply_trained_models_to_config): its file ID """ + models = [dict(model, genomes=key.split(' ')) for key, model in sorted(trained.items()) if model is not None] + path = job.fileStore.getLocalTempFile() + with open(path, 'w') as out_file: + json.dump(models, out_file, indent=1) + return job.fileStore.writeGlobalFile(path) def get_lastz_walltime(distance, params, chunk_a, chunk_b): @@ -183,7 +197,7 @@ def get_lastz_walltime(distance, params, chunk_a, chunk_b): def run_lastz(job, name_A, genome_A, name_B, genome_B, distance, params, trained_score_model=None): - """ trained_score_model is the pair's model from train_lastz_score_model, if made one """ + """ trained_score_model is the pair's lastz model from train_pair_scoring_models, with """ # Create a local temporary file to put the alignments in. work_dir = job.fileStore.getLocalTempDir() alignment_file = os.path.join(work_dir, '{}_{}.paf'.format(name_A, name_B)) @@ -625,7 +639,8 @@ def chunked_alignment_walltime(genome_a, genome_b): def make_chunked_alignments(job, event_a, genome_a, event_b, genome_b, distance, params, score_model=None): - """ score_model is the pair's trained lastz model, if made one """ + """ score_model is the pair's trained models, if it has them (see train_pair_scoring_models): lastz scores + with them with """ lastz_params_node = params.find("blast") gpu = getOptionalAttrib(lastz_params_node, 'gpu', typeFn=int, default=0) fastga = getOptionalAttrib(lastz_params_node, 'mapper', typeFn=str) == 'fastga' @@ -658,7 +673,8 @@ def make_chunks(genome): for j, chunk_b in enumerate(chunks_b): mappers = { "lastz":run_lastz, "minimap2":run_minimap2, "fastga":run_fastga} mappingFn = mappers[params.find("blast").attrib["mapper"]] - mapper_args = {'trained_score_model': score_model} if score_model is not None and mappingFn is run_lastz else {} + mapper_args = {'trained_score_model': score_model['lastz']} if score_model is not None and \ + mappingFn is run_lastz and lastz_train_enabled(params) else {} memory = lastz_memory if lastz_memory else max(200000000, 15*(chunk_a.size+chunk_b.size)) chunked_alignment_files.append(job.addChildJobFn(mappingFn, '{}_{}'.format(event_a, i), chunk_a, '{}_{}'.format(event_b, j), chunk_b, distance, params, @@ -1066,10 +1082,17 @@ def merge_processed_alignments(job, processed_file_ids): def sanitize_then_make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancestor_event_string, params, - output_path=None): + output_path=None, return_trained_models=False): + """ As make_paf_alignments, with the output gzipped if output_path ends in .gz """ sanitize_job = job.addChildJobFn(sanitize_fasta_headers, event_names_to_sequences, walltime=cactus_walltime()) paf_job = sanitize_job.addFollowOnJobFn(make_paf_alignments, event_tree_string, sanitize_job.rv(), - ancestor_event_string, params, walltime=cactus_walltime()) + ancestor_event_string, params, return_trained_models=return_trained_models, + walltime=cactus_walltime()) + if return_trained_models: + paf, models = paf_job.rv('paf'), paf_job.rv('models') + if output_path and output_path.endswith('.gz'): + paf = paf_job.addFollowOnJobFn(zip_gz, output_path, paf, walltime=cactus_walltime(900)).rv() + return {'paf': paf, 'models': models} # gzip the output if requested if output_path and output_path.endswith('.gz'): # gzipping the final paf: one per run in the 577-way blast logs, p50 165s, p99 575s, max @@ -1081,7 +1104,11 @@ def sanitize_then_make_paf_alignments(job, event_tree_string, event_names_to_seq return paf_job.rv() -def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancestor_event_string, params): +def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancestor_event_string, params, + return_trained_models=False): + """ The PAF of the alignments for an ancestor, or with return_trained_models a dict of it ('paf') and of the + file of the scoring models trained for its pairs of genomes ('models', see write_trained_models; None when + none were trained) """ # a job should never set its own follow-on, so we hang everything off the root_job here to encapsulate root_job = Job(walltime=cactus_walltime()) job.addChild(root_job) @@ -1125,11 +1152,11 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest for i,ingroup in enumerate(ingroups): event_names_to_sequences[ingroup] = unmask_job.rv(i) - # train lastz's scores on each pair of genomes it will align (see ), on the whole genomes - # (the trimmed ingroups aligned to later outgroups are what earlier ones left) and alongside the unmasking, - # which only changes case: last-train masks repeats its own way + # train scoring on each pair of genomes to be aligned (see and ), on the + # whole genomes (the trimmed ingroups aligned to later outgroups are what earlier ones left) and alongside the + # unmasking, which only changes case: last-train masks repeats its own way score_models = {} - if lastz_train_enabled(params): + if lastz_train_enabled(params) or bar_train_enabled(params): train_cores = getOptionalAttrib(lastz_params_node, 'lastzTrainCores', typeFn=int, default=4) pairs = [(ingroup, ingroup2) for ingroup, ingroup2, _ in get_event_pairs(ancestor_event, ingroup_events)] + \ [(ingroup, outgroup) for ingroup in ingroup_events for outgroup in outgroup_events] @@ -1142,7 +1169,7 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest # fish genome and 34 on a mouse one (split into volumes), and training took 2 to 17 minutes. db_letters = min(letters[event_a.iD], letters[event_b.iD]) score_models[lastz_pair_key(event_a.iD, event_b.iD)] = root_job.addChildJobFn( - train_lastz_score_model, event_a.iD, genome_a, event_b.iD, genome_b, params, + train_pair_scoring_models, event_a.iD, genome_a, event_b.iD, genome_b, params, cores=cactus_clamp_cores(train_cores), memory=cactus_clamp_memory(int(1e9 + 5.5 * db_letters)), disk=int(2e9 + 5 * db_letters + genome_a.size + genome_b.size), @@ -1203,21 +1230,26 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest # for better logs outgroup_alignment_names = ['{}-og_{}'.format(ancestor_event_string, i) for i in range(len(outgroup_alignments))] + # the trained models, for bar + models = root_job.addChildJobFn(write_trained_models, score_models, walltime=cactus_walltime()).rv() \ + if return_trained_models and score_models else None + # Now do the chaining, either getting ingroup primary alignments separately to primary ingroup to outgroup # alignments if the following option is true, otherwise get every sequence to pick its primary alignment without # regard to whether the other sequence is an ingroup or an outgroup if int(params.find("blast").attrib["pickIngroupPrimaryAlignmentsSeparatelyToOutgroups"]): - return root_job.addFollowOnJobFn(chain_alignments_splitting_ingroups_and_outgroups, + paf = root_job.addFollowOnJobFn(chain_alignments_splitting_ingroups_and_outgroups, ingroup_alignments, ingroup_alignment_names, outgroup_alignments, outgroup_alignment_names, ancestor_event_string, params, total_sequence_size=total_sequence_size, walltime=cactus_walltime()).rv() + return {'paf': paf, 'models': models} if return_trained_models else paf # Delete the unmasked fastas (todo: should we do the unmasking somewhere further upstream?) for ingroup in ingroup_events: root_job.addFollowOnJobFn(clean_if_different, event_names_to_sequences[ingroup.iD], input_sequence_map[ingroup.iD], walltime=cactus_walltime()) - return root_job.addFollowOnJobFn(chain_alignments, ingroup_alignments + outgroup_alignments, + paf = root_job.addFollowOnJobFn(chain_alignments, ingroup_alignments + outgroup_alignments, ingroup_alignment_names + outgroup_alignment_names, ancestor_event_string, params, total_sequence_size=total_sequence_size, @@ -1228,6 +1260,7 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest # and split_file at the p99 paf-to-sequence ratio of about 7 walltime=cactus_walltime(600 * total_sequence_size / 1e9, io_bytes=8 * total_sequence_size)).rv() + return {'paf': paf, 'models': models} if return_trained_models else paf def trim_unaligned_sequences(job, sequences, alignments, params, has_resources=False): diff --git a/src/cactus/pipeline/cactus_workflow.py b/src/cactus/pipeline/cactus_workflow.py index 8fa71ee53..4bcf0de11 100644 --- a/src/cactus/pipeline/cactus_workflow.py +++ b/src/cactus/pipeline/cactus_workflow.py @@ -10,6 +10,7 @@ import os import copy import sys +import json from toil.lib.bioio import system from toil.lib.bioio import getLogLevelString from toil.realtimeLogger import RealtimeLogger @@ -21,6 +22,7 @@ from cactus.shared.common import findRequiredNode, getOptionalAttrib from cactus.shared.common import cactus_clamp_memory from cactus.shared.common import cactus_walltime +from cactus.paf.last_scoring import bar_train_enabled, apply_trained_models_to_config ############################################################ ############################################################ @@ -73,8 +75,9 @@ def bar_base_aligner(tree, ancestor_event, config_node, name=None): def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_map, og_map, paf_id, cons_cores = None, cons_memory = None, intermediate_results_url = None, chrom_name = None, - cons_retain_pages = None): - ''' run cactus_consolidated as a child job, requesting resources based on input sizes ''' + cons_retain_pages = None, trained_models = None): + ''' run cactus_consolidated as a child job, requesting resources based on input sizes. trained_models is the + file of scoring models the blast step trained for the ancestor's pairs of genomes, if any (see cactus_cons) ''' cons_node = findRequiredNode(config_node, 'consolidated') name = chrom_name if chrom_name else ancestor_event @@ -248,18 +251,31 @@ def cactus_cons_with_resources(job, tree, ancestor_event, config_node, seq_id_ma cons_job = job.addChildJobFn(cactus_cons, tree, ancestor_event, config_node, seq_id_map, og_map, paf_id, intermediate_results_url=intermediate_results_url, chrom_name=chrom_name, cores = cons_cores, memory=cactus_clamp_memory(mem), disk=disk, retain_pages=retain_pages, - poa_window=poa_window, base_aligner=base_aligner, walltime=cactus_walltime(walltime_secs)) + poa_window=poa_window, base_aligner=base_aligner, trained_models=trained_models, + walltime=cactus_walltime(walltime_secs)) return cons_job.rv() def cactus_cons(job, tree, ancestor_event, config_node, seq_id_map, og_map, paf_id, intermediate_results_url = None, chrom_name = None, retain_pages = None, - poa_window = None, base_aligner = None): - ''' run cactus_consolidated ''' + poa_window = None, base_aligner = None, trained_models = None): + ''' run cactus_consolidated. With , bar scores with trained_models, the file of + scoring models the blast step trained for the ancestor's pairs of genomes (see + last_scoring.apply_trained_models_to_config); without it, or without the file, bar keeps its settings ''' # cactus_consolidated reads its settings from the config, so the resolved page retention # goes into the copy it is given (this job's copy of the node, so nothing else sees it) - if retain_pages is not None or poa_window is not None or base_aligner is not None: + use_trained_models = bar_train_enabled(config_node) + if use_trained_models and trained_models is None: + RealtimeLogger.warning('cactus_consolidated({}): is set, but no models were trained for ' + 'this ancestor, so bar keeps its settings'.format(ancestor_event)) + if retain_pages is not None or poa_window is not None or base_aligner is not None or \ + (use_trained_models and trained_models is not None): config_node = copy.deepcopy(config_node) + if use_trained_models and trained_models is not None: + with open(job.fileStore.readGlobalFile(trained_models)) as models_file: + models = json.load(models_file) + apply_trained_models_to_config(config_node, models, findRequiredNode(config_node, 'reference').attrib.get( + 'reconstructionTree')) if retain_pages is not None: findRequiredNode(config_node, 'consolidated').set('retain_pages', str(retain_pages)) # only the workflow applies partialOrderAlignmentWindowBigGenome, so the window it resolved diff --git a/src/cactus/progressive/cactus_progressive.py b/src/cactus/progressive/cactus_progressive.py index 9d497570f..5168d39f8 100755 --- a/src/cactus/progressive/cactus_progressive.py +++ b/src/cactus/progressive/cactus_progressive.py @@ -45,6 +45,7 @@ from cactus.preprocessor.checkPreprocessedSequence import preprocessed_fasta_id from cactus.hal.cactus_validate import fail_hal_validation, validate_hal_export from cactus.paf.local_alignment import make_paf_alignments, trim_unaligned_sequences +from cactus.paf.last_scoring import bar_train_enabled from cactus.shared.configWrapper import ConfigWrapper from cactus.progressive.multiCactusTree import MultiCactusTree from cactus.shared.common import setupBinaries, importSingularityImage @@ -189,30 +190,36 @@ def progressive_step(job, options, config_node, seq_id_map, tree, og_map, event) # get the spanning tree (which is what consolidated wants) spanning_tree = get_spanning_subtree(tree, event, ConfigWrapper(config_node), og_map) - # do the blast + # do the blast, keeping the scoring models it trains for bar if it is to use them + trained = bar_train_enabled(config_node) paf_job = job.addChildJobFn(make_paf_alignments, NXNewick().writeString(spanning_tree), - subtree_eventmap, event, config_node, walltime=cactus_walltime()).encapsulate() + subtree_eventmap, event, config_node, return_trained_models=trained, + walltime=cactus_walltime()).encapsulate() + paf = paf_job.rv('paf') if trained else paf_job.rv() + trained_models = paf_job.rv('models') if trained else None outgroups = og_map[event] if event in og_map else [] # trim the outgroups if outgroups and int(config_node.find("blast").attrib["trimOutgroups"]): # Trim the outgroup sequences trim_sequences = paf_job.addChildJobFn(trim_unaligned_sequences, - [subtree_eventmap[i] for i in outgroups], paf_job.rv(), config_node, + [subtree_eventmap[i] for i in outgroups], paf, config_node, walltime=cactus_walltime()) cons_job = paf_job.addFollowOnJobFn(progressive_step_2, trim_sequences.rv(), options, config_node, subtree_eventmap, - spanning_tree, og_map, event, walltime=cactus_walltime()) + spanning_tree, og_map, event, trained_models=trained_models, + walltime=cactus_walltime()) else: # Without outgroup trimming (or if there are no outgroups to trim) cons_job = paf_job.addChildJobFn(cactus_cons_with_resources, spanning_tree, event, config_node, subtree_eventmap, - og_map, paf_job.rv(), cons_cores=options.consCores, cons_memory=options.consMemory, cons_retain_pages=getattr(options, 'consRetainPages', None), - intermediate_results_url=options.intermediateResultsUrl, walltime=cactus_walltime()) + og_map, paf, cons_cores=options.consCores, cons_memory=options.consMemory, cons_retain_pages=getattr(options, 'consRetainPages', None), + intermediate_results_url=options.intermediateResultsUrl, + trained_models=trained_models, walltime=cactus_walltime()) # erase the paf since its now longer needed - cons_job.addFollowOnJobFn(clean_jobstore_files, file_ids=[paf_job.rv()], walltime=cactus_walltime()) + cons_job.addFollowOnJobFn(clean_jobstore_files, file_ids=[paf], walltime=cactus_walltime()) return cons_job.rv() def progressive_step_2(job, trimmed_outgroups_and_alignments, options, config_node, subtree_eventmap, - spanning_tree, og_map, event): + spanning_tree, og_map, event, trained_models=None): trimmed_outgroup_seqs, pafs = trimmed_outgroups_and_alignments # unpack the pafs and outgroup sequences # set the outgroup seqs for outgroup, sequence in zip(og_map[event] if event in og_map else [], trimmed_outgroup_seqs): @@ -221,7 +228,8 @@ def progressive_step_2(job, trimmed_outgroups_and_alignments, options, config_no # now do consolidated return job.addChildJobFn(cactus_cons_with_resources, spanning_tree, event, config_node, subtree_eventmap, og_map, pafs, cons_cores=options.consCores, cons_memory=options.consMemory, cons_retain_pages=getattr(options, 'consRetainPages', None), - intermediate_results_url=options.intermediateResultsUrl, walltime=cactus_walltime()).rv() + intermediate_results_url=options.intermediateResultsUrl, + trained_models=trained_models, walltime=cactus_walltime()).rv() def export_hal(job, mc_tree, config_node, seq_id_map, og_map, results, event=None, cacheBytes=None, diff --git a/src/cactus/setup/cactus_align.py b/src/cactus/setup/cactus_align.py index bce3e5f3b..615a601e1 100644 --- a/src/cactus/setup/cactus_align.py +++ b/src/cactus/setup/cactus_align.py @@ -37,7 +37,7 @@ from cactus.refmap.cactus_graphmap import filter_paf from cactus.refmap.cactus_minigraph import check_sample_names from cactus.preprocessor.checkUniqueHeaders import sanitize_fasta_headers -from cactus.paf.last_scoring import parse_train_file, apply_scores_to_config +from cactus.paf.last_scoring import parse_train_file, apply_scores_to_config, bar_train_enabled, trained_models_path from toil.job import Job from toil.common import Toil @@ -418,6 +418,16 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): # import the PAF alignments paf_id = toil.importFile(makeURL(options.pafFile)) + + # and with , the scoring models cactus-blast trained for them, which it leaves next to them + trained_models_id = None + if bar_train_enabled(config_wrapper.xmlRoot): + models_url = makeURL(trained_models_path(options.pafFile)) + try: + trained_models_id = toil.importFile(models_url) + except Exception as e: + logger.warning(' is set, but there are no trained models at {} ({}), so bar keeps its ' + 'settings: were the alignments made by cactus-blast with it set?'.format(models_url, e)) #import the sequences input_seq_id_map = {} @@ -452,13 +462,14 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): do_filter_paf=options.pangenome, chrom_name=chrom_name, scores_id=scores_id, + trained_models_id=trained_models_id, branch_scale=options.branchScale, validate=getattr(options, 'validate', False), walltime=cactus_walltime()) return align_job def cactus_align(job, config_wrapper, mc_tree, input_seq_map, input_seq_id_map, paf_id, paf_path, root_name, og_map, checkpointInfo, doVG, doGFA, delay=0, referenceEvents=None, pafMaskFilter=None, paf2Stable=False, cons_cores = None, cons_memory = None, cons_retain_pages = None, do_filter_paf=False, chrom_name=None, scores_id=None, branch_scale=1.0, - validate=False): + validate=False, trained_models_id=None): head_job = Job(walltime=cactus_walltime()) job.addChild(head_job) @@ -530,7 +541,8 @@ def cactus_align(job, config_wrapper, mc_tree, input_seq_map, input_seq_id_map, # run consolidated cons_job = head_job.addFollowOnJobFn(cactus_cons_with_resources, spanning_tree, root_name, config_wrapper.xmlRoot, new_seq_id_map, og_map, paf_id, cons_cores = cons_cores, cons_memory=cons_memory, chrom_name=chrom_name, - cons_retain_pages=cons_retain_pages, walltime=cactus_walltime()) + cons_retain_pages=cons_retain_pages, trained_models=trained_models_id, + walltime=cactus_walltime()) results = {root_name : (cons_job.rv(1), cons_job.rv(2))} # get the immediate subtree (which is all export_hal can use) From 3651e76d60f9b70e10330b552dd25bded235cfe4 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Mon, 5 Oct 2026 16:59:06 -0400 Subject: [PATCH 09/10] bar: score with the trained models by default is now on, so with (on in the lastz scoring branch) progressive cactus trains a model per pair of genomes and scores both lastz and bar with it. cactus-align no longer looks for trained models next to pangenome alignments, which come from minigraph, and the consolidated job no longer warns for an ancestor without them: cactus-align says when a models file it expected is missing. Whether lastz uses the models is decided once per ancestor, so a KegAlign run warns once rather than for every chunk. On evolver mammals the shipped config now gives the same alignment, byte for byte, as the earlier run with both switched on: pooled leaf-pair precision/recall 0.836/0.742, against 0.830/0.744 with lastz alone trained, and 0.816/0.739 with neither. abPOA's trained rows cost up to 0.9 points of ancestor recall there, against 0.830/0.744's; with pecan, training both gave 0.855/0.760. Co-Authored-By: Claude Opus 5.5 --- src/cactus/cactus_progressive_config.xml | 7 +++++-- src/cactus/paf/local_alignment.py | 18 ++++++++++-------- src/cactus/pipeline/cactus_workflow.py | 3 --- src/cactus/setup/cactus_align.py | 3 ++- 4 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/cactus/cactus_progressive_config.xml b/src/cactus/cactus_progressive_config.xml index f952d2edd..af4c63167 100644 --- a/src/cactus/cactus_progressive_config.xml +++ b/src/cactus/cactus_progressive_config.xml @@ -304,7 +304,10 @@ with the model for its genome and the nearest ingroup (, measured "ingroup"), and pecan each pair of sequences with the model for their two genomes (); a pair blast did not train, as pecan's outgroup pairs, takes the model trained nearest its distance. They replace any row or pair - models in the config. cactus-blast leaves them next to its PAF, as .models.json, for cactus-align. --> + models in the config. cactus-blast leaves them next to its PAF, as .models.json, for cactus-align; pangenome + alignments have none. On by default: on evolver mammals, with lastz trained too, abPOA went from pooled leaf-pair + precision/recall 0.830/0.744 (lastz alone) to 0.836/0.742, losing up to 0.9 points of ancestor recall, and pecan + from 0.838/0.752 (neither trained) to 0.855/0.760. 0 keeps bar's own settings --> @@ -317,7 +320,7 @@ baseAligner="abpoa" deepAligner="pecan" deepAlignerDivergence="0" - trainedModels="0" + trainedModels="1" minimumBlockDegree="2" minimumIngroupDegree="1" minimumOutgroupDegree="0" diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py index 400af6bb9..2714bb070 100755 --- a/src/cactus/paf/local_alignment.py +++ b/src/cactus/paf/local_alignment.py @@ -639,8 +639,8 @@ def chunked_alignment_walltime(genome_a, genome_b): def make_chunked_alignments(job, event_a, genome_a, event_b, genome_b, distance, params, score_model=None): - """ score_model is the pair's trained models, if it has them (see train_pair_scoring_models): lastz scores - with them with """ + """ score_model is the pair's trained models, if lastz is to score with them (see train_pair_scoring_models, + and ) """ lastz_params_node = params.find("blast") gpu = getOptionalAttrib(lastz_params_node, 'gpu', typeFn=int, default=0) fastga = getOptionalAttrib(lastz_params_node, 'mapper', typeFn=str) == 'fastga' @@ -673,8 +673,7 @@ def make_chunks(genome): for j, chunk_b in enumerate(chunks_b): mappers = { "lastz":run_lastz, "minimap2":run_minimap2, "fastga":run_fastga} mappingFn = mappers[params.find("blast").attrib["mapper"]] - mapper_args = {'trained_score_model': score_model['lastz']} if score_model is not None and \ - mappingFn is run_lastz and lastz_train_enabled(params) else {} + mapper_args = {'trained_score_model': score_model['lastz']} if score_model is not None and mappingFn is run_lastz else {} memory = lastz_memory if lastz_memory else max(200000000, 15*(chunk_a.size+chunk_b.size)) chunked_alignment_files.append(job.addChildJobFn(mappingFn, '{}_{}'.format(event_a, i), chunk_a, '{}_{}'.format(event_b, j), chunk_b, distance, params, @@ -1156,7 +1155,8 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest # whole genomes (the trimmed ingroups aligned to later outgroups are what earlier ones left) and alongside the # unmasking, which only changes case: last-train masks repeats its own way score_models = {} - if lastz_train_enabled(params) or bar_train_enabled(params): + lastz_trained = lastz_train_enabled(params) + if lastz_trained or bar_train_enabled(params): train_cores = getOptionalAttrib(lastz_params_node, 'lastzTrainCores', typeFn=int, default=4) pairs = [(ingroup, ingroup2) for ingroup, ingroup2, _ in get_event_pairs(ancestor_event, ingroup_events)] + \ [(ingroup, outgroup) for ingroup in ingroup_events for outgroup in outgroup_events] @@ -1179,6 +1179,8 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest new_root_job = Job(walltime=cactus_walltime()) root_job.addFollowOn(new_root_job) root_job = new_root_job + # the models lastz aligns with, if it is to use them (bar's go by write_trained_models below) + lastz_models = score_models if lastz_trained else {} # for each pair of ingroups make alignments ingroup_alignments = [] @@ -1198,7 +1200,7 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest disk=2*total_sequence_size, walltime=chunked_alignment_walltime(input_sequence_map[ingroup.iD], input_sequence_map[ingroup2.iD]), - score_model=score_models.get(lastz_pair_key(ingroup.iD, ingroup2.iD))).rv()) + score_model=lastz_models.get(lastz_pair_key(ingroup.iD, ingroup2.iD))).rv()) ingroup_alignment_names.append('{}-{}_vs_{}'.format(ancestor_event_string, ingroup.iD, ingroup2.iD)) distances = get_distances(event_tree) # Distances between all pairs of nodes @@ -1212,7 +1214,7 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest if int(params.find("blast").attrib["trimIngroups"]): # Trim the ingroup sequences outgroup_alignments = [root_job.addChildJobFn(make_ingroup_to_outgroup_alignments_0, ingroup, outgroup_events, dict(event_names_to_sequences), distances, params, walltime=cactus_walltime(), - score_models=score_models).rv() + score_models=lastz_models).rv() for ingroup in ingroup_events] if len(outgroup_events) > 0 else [] else: outgroup_alignments = [root_job.addChildJobFn(make_chunked_alignments, @@ -1225,7 +1227,7 @@ def make_paf_alignments(job, event_tree_string, event_names_to_sequences, ancest disk=2*total_sequence_size, walltime=chunked_alignment_walltime(input_sequence_map[ingroup.iD], input_sequence_map[outgroup.iD]), - score_model=score_models.get(lastz_pair_key(ingroup.iD, outgroup.iD))).rv() + score_model=lastz_models.get(lastz_pair_key(ingroup.iD, outgroup.iD))).rv() for ingroup in ingroup_events for outgroup in outgroup_events] # for better logs outgroup_alignment_names = ['{}-og_{}'.format(ancestor_event_string, i) for i in range(len(outgroup_alignments))] diff --git a/src/cactus/pipeline/cactus_workflow.py b/src/cactus/pipeline/cactus_workflow.py index 4bcf0de11..d7d537989 100644 --- a/src/cactus/pipeline/cactus_workflow.py +++ b/src/cactus/pipeline/cactus_workflow.py @@ -265,9 +265,6 @@ def cactus_cons(job, tree, ancestor_event, config_node, seq_id_map, og_map, paf_ # cactus_consolidated reads its settings from the config, so the resolved page retention # goes into the copy it is given (this job's copy of the node, so nothing else sees it) use_trained_models = bar_train_enabled(config_node) - if use_trained_models and trained_models is None: - RealtimeLogger.warning('cactus_consolidated({}): is set, but no models were trained for ' - 'this ancestor, so bar keeps its settings'.format(ancestor_event)) if retain_pages is not None or poa_window is not None or base_aligner is not None or \ (use_trained_models and trained_models is not None): config_node = copy.deepcopy(config_node) diff --git a/src/cactus/setup/cactus_align.py b/src/cactus/setup/cactus_align.py index 615a601e1..7cbbc1339 100644 --- a/src/cactus/setup/cactus_align.py +++ b/src/cactus/setup/cactus_align.py @@ -420,8 +420,9 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): paf_id = toil.importFile(makeURL(options.pafFile)) # and with , the scoring models cactus-blast trained for them, which it leaves next to them + # (pangenome alignments come from minigraph, with nothing of the kind) trained_models_id = None - if bar_train_enabled(config_wrapper.xmlRoot): + if bar_train_enabled(config_wrapper.xmlRoot) and not options.pangenome: models_url = makeURL(trained_models_path(options.pafFile)) try: trained_models_id = toil.importFile(models_url) From a3044b8a7d415da65b7c23fd57f4231745c593b0 Mon Sep 17 00:00:00 2001 From: Glenn Hickey Date: Tue, 6 Oct 2026 07:53:01 -0400 Subject: [PATCH 10/10] cactus-prepare: give cactus-align blast's trained models by option cactus-prepare's WDL and Toil workflows gave cactus-align everything cactus-blast wrote whose name starts with the PAF's, as alignments. With on by default, that took in the trained models file, .models.json, and cactus-align stopped with a usage error. cactus-align now takes the models file with --trainedModels (still the file next to the PAF by default), and both workflows give it that way. The WDL's blast task names its two outputs, out_paf and out_models: cromwell puts the two in different directories, so cactus-align could not have found the models next to the PAF by itself. Co-Authored-By: Claude Opus 5.5 --- src/cactus/progressive/cactus_prepare.py | 33 ++++++++++++++++++------ src/cactus/setup/cactus_align.py | 15 +++++++++-- 2 files changed, 38 insertions(+), 10 deletions(-) diff --git a/src/cactus/progressive/cactus_prepare.py b/src/cactus/progressive/cactus_prepare.py index 3d3530745..43ece964b 100644 --- a/src/cactus/progressive/cactus_prepare.py +++ b/src/cactus/progressive/cactus_prepare.py @@ -36,6 +36,7 @@ from cactus.pipeline.cactus_workflow import cons_core_scale from cactus.shared.common import cactus_clamp_memory from cactus.shared.common import cactus_walltime +from cactus.paf.last_scoring import trained_models_path from toil.job import Job from toil.common import Toil @@ -1039,7 +1040,11 @@ def wdl_task_blast(options): s += ' docker: \"{}\"\n'.format(options.dockerImage) s += ' zones: \"{}\"\n'.format(options.zone) s += ' }\n' - s += ' output {\n Array[File] out_files=glob(\"${out_name}*\")\n }\n' + # the alignments, and the scoring models trained for bar (see ) if any were + s += ' output {\n' + s += ' File out_paf=\"${out_name}\"\n' + s += ' Array[File] out_models=glob(\"{}\")\n'.format(trained_models_path('${out_name}')) + s += ' }\n' s += '}\n' return s @@ -1132,7 +1137,8 @@ def wdl_task_align(options): s += ' Array[String] in_fa_names\n' s += ' Array[File] in_fa_files\n' s += ' Array[String]? in_fa_urls\n' - s += ' Array[File] in_blast_files\n' + s += ' File in_blast_file\n' + s += ' Array[File] in_blast_models\n' s += ' String in_root\n' s += ' File? in_config_file\n' s += ' File? in_chrom_info_file\n' @@ -1144,7 +1150,10 @@ def wdl_task_align(options): # two commands run here, so without set -e a failed cactus-align would be # masked by whatever cactus-hal2fasta returns afterwards s += 'set -eo pipefail\n ' - s += 'cactus-align {} ${{in_seq_file}} ${{sep=\" \" in_blast_files}} ${{out_hal_name}} --root ${{in_root}}'.format(get_jobstore(options, 'align')) + s += 'cactus-align {} ${{in_seq_file}} ${{in_blast_file}} ${{out_hal_name}} --root ${{in_root}}'.format(get_jobstore(options, 'align')) + # and --trainedModels, if blast trained any. Not with prefix(): cromwell 49 can't put the empty array it makes + # into a command, and a String declared from a File is its path before localization + s += ' ${true=\"--trainedModels \" false=\"\" length(in_blast_models) > 0}${sep=\" \" in_blast_models}' s += ' --pathOverrides ${{sep=\" \" in_fa_files}} ${{sep=\" \" in_fa_urls}} --pathOverrideNames ${{sep=\" \" in_fa_names}} {}'.format(options.cactusOptions) s += ' ${\"--chromInfo \" + in_chrom_info_file}' s += ' {} {}{} ${{\"--configFile \" + in_config_file}} ${{in_options}}'.format(get_toil_resource_opts(options, 'align'), @@ -1201,7 +1210,8 @@ def wdl_call_align(options, in_seq_file, mc_tree, og_map, event, cigar_name, hal s += ' in_fa_files=[{}],'.format(', '.join(input_fas)) if input_urls: s += ' in_fa_urls=[{}],'.format(', '.join(input_urls)) - s += ' in_blast_files={}.out_files,'.format(blast_call_name(event)) + s += ' in_blast_file={}.out_paf,'.format(blast_call_name(event)) + s += ' in_blast_models={}.out_models,'.format(blast_call_name(event)) s += ' in_root=\"{}\",'.format(event) s += ' in_config_file=config_file,' s += ' in_chrom_info_file=chrom_info_file,' @@ -1222,11 +1232,17 @@ def toil_call_align(job, options, seq_file, mc_tree, og_map, event, cigar_name, with open(seq_file_path, 'w') as sf: sf.write(str(seq_file)) - # download the blast output from the file store + # download the blast output from the file store: the alignments, and any scoring models trained for bar + # (see ), which cactus-align takes by option blast_files = [] + models_files = [] for blast_file_name, blast_file_id in blast_output: - blast_files.append(os.path.join(work_dir, blast_file_name)) - job.fileStore.readGlobalFile(blast_file_id, blast_files[-1]) + local_path = os.path.join(work_dir, blast_file_name) + job.fileStore.readGlobalFile(blast_file_id, local_path) + if blast_file_name == os.path.basename(trained_models_path(cigar_name)): + models_files.append(local_path) + else: + blast_files.append(local_path) # read the fasta files assert len(dep_names) == len(dep_fa_ids) @@ -1239,7 +1255,8 @@ def toil_call_align(job, options, seq_file, mc_tree, og_map, event, cigar_name, cactus_call(parameters=['cactus-align', os.path.join(work_dir, 'js'), seq_file_path] + blast_files + [out_hal_path, '--root', event, '--pathOverrides'] + fa_paths + ['--pathOverrideNames'] + dep_names + - ['--workDir', work_dir, '--maxCores', str(int(job.cores)), '--maxDisk', bytes2humanN(job.disk), '--maxMemory', bytes2humanN(job.memory)] + options.cactusOptions.strip().split(' ') + (['--chromInfo', options.chromInfo] if options.chromInfo else []) + (['--branchScale', str(options.branchScale)] if options.branchScale else []) + ([] if options.noValidate else ['--validate'])) + ['--workDir', work_dir, '--maxCores', str(int(job.cores)), '--maxDisk', bytes2humanN(job.disk), '--maxMemory', bytes2humanN(job.memory)] + options.cactusOptions.strip().split(' ') + (['--chromInfo', options.chromInfo] if options.chromInfo else []) + (['--branchScale', str(options.branchScale)] if options.branchScale else []) + ([] if options.noValidate else ['--validate']) + + (['--trainedModels', models_files[0]] if models_files else [])) out_hal_id = job.fileStore.writeGlobalFile(out_hal_path) diff --git a/src/cactus/setup/cactus_align.py b/src/cactus/setup/cactus_align.py index 7cbbc1339..cf64d7547 100644 --- a/src/cactus/setup/cactus_align.py +++ b/src/cactus/setup/cactus_align.py @@ -98,6 +98,9 @@ def main(): help = "File containing scoring parameters (output of last-train, as made by cactus-minigraph by default)") parser.add_argument("--scoresFromChromfile", action="store_true", default=False, help = "Load scoring parameters from the 4th column of chromfile (as made by cactus-minigraph --batch)") + parser.add_argument("--trainedModels", type=str, + help = "Scoring models cactus-blast trained for the alignments, for bar (see in " + "the config) [default=.models.json, where cactus-blast leaves them]") parser.add_argument("--singleCopySpecies", type=str, help="Filter out all self-alignments in given species") @@ -192,6 +195,8 @@ def main(): raise RuntimeError('--scoresFile is only currently supported with --pangenome') if options.scoresFromChromfile and (not options.pangenome or not options.batch or options.scoresFile): raise RuntimeError('--scoresFromChromfile can only be used with --batch --pangenome and without --scoresFile') + if options.trainedModels and options.pangenome: + raise RuntimeError('--trainedModels is not supported with --pangenome') options.buildHal = True options.buildFasta = True @@ -420,15 +425,21 @@ def make_align_job(options, toil, config_wrapper=None, chrom_name=None): paf_id = toil.importFile(makeURL(options.pafFile)) # and with , the scoring models cactus-blast trained for them, which it leaves next to them - # (pangenome alignments come from minigraph, with nothing of the kind) + # unless --trainedModels says where (pangenome alignments come from minigraph, with nothing of the kind). + # getattr as cactus-pangenome comes through here with its own options trained_models_id = None + trained_models_file = getattr(options, 'trainedModels', None) if bar_train_enabled(config_wrapper.xmlRoot) and not options.pangenome: - models_url = makeURL(trained_models_path(options.pafFile)) + models_url = makeURL(trained_models_file if trained_models_file else trained_models_path(options.pafFile)) try: trained_models_id = toil.importFile(models_url) except Exception as e: + if trained_models_file: + raise logger.warning(' is set, but there are no trained models at {} ({}), so bar keeps its ' 'settings: were the alignments made by cactus-blast with it set?'.format(models_url, e)) + elif trained_models_file: + logger.warning('Ignoring --trainedModels, as is not set') #import the sequences input_seq_id_map = {}