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..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" @@ -52,9 +53,201 @@ 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); +} + +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; // 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; +}; + +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->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); + 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)); + 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 (%" PRIi64 " for pairs of genomes), distances on the %s tree\n", + modelNo, keyed, 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]); + } + 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); + free(pairModels); +} + +/* + * cPecan's per-pair callback: the model for the two sequences' genomes, else for the distance between them. + */ +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; + } + 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; +} + 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 +255,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 +266,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 +279,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 feb9d92bb..4995452ae 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,154 @@ const char *baseAligner_toString(BaseAligner engine) { return "unknown"; } +/* + * 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->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); + 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 (%" 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 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) { + 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; + } + // 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]); + if (d < nearest[j]) { + nearest[j] = d, partner[j] = k; + } + if (d < nearest[k]) { + nearest[k] = d, partner[k] = j; + } + 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; + } + } + } + // 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 (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] = eventModels[row_event[i]]; + } + free(events); + free(row_event); + free(nearest); + free(nearestIngroup); + free(partner); + free(ingroupPartner); + free(eventModels); + return models; +} + PoaParameters *poaParameters_constructFromCactusParams(CactusParams *params, BaseAligner engine) { if (engine == BASE_ALIGNER_PECAN) { return NULL; @@ -903,6 +1097,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 +1143,18 @@ 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); + 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); + } #ifdef HAVE_MINIPOA if (poaParameters->mpt != NULL) { minipoa_free_para((minipoa_para_t *)poaParameters->mpt); @@ -968,10 +1175,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 +1201,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 +1308,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 +1520,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 +1560,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 (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); + +/* + * 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 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; + +/* + * 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 462b227b0..5340bf99c 100644 --- a/bar/inc/poaBarAligner.h +++ b/bar/inc/poaBarAligner.h @@ -106,6 +106,27 @@ typedef struct _PoaParameters { /* 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 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; // 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; /** @@ -181,6 +202,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 +249,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. */ 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..af4c63167 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,15 @@ 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 +320,7 @@ baseAligner="abpoa" deepAligner="pecan" deepAlignerDivergence="0" + trainedModels="1" minimumBlockDegree="2" minimumIngroupDegree="1" minimumOutgroupDegree="0" diff --git a/src/cactus/paf/last_scoring.py b/src/cactus/paf/last_scoring.py index 1e810944d..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. @@ -292,7 +445,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 +463,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: diff --git a/src/cactus/paf/local_alignment.py b/src/cactus/paf/local_alignment.py index 8eaabd844..2714bb070 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 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' @@ -658,7 +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} 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 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 +1081,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 +1103,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 +1151,12 @@ 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): + 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] @@ -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), @@ -1152,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 = [] @@ -1171,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 @@ -1185,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, @@ -1198,26 +1227,31 @@ 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))] + # 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 +1262,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..d7d537989 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,28 @@ 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 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_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/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..cf64d7547 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 @@ -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 @@ -418,6 +423,23 @@ 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 + # 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_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 = {} @@ -452,13 +474,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 +553,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) 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 diff --git a/submodules/cPecan b/submodules/cPecan index 55920efb7..2ec01d160 160000 --- a/submodules/cPecan +++ b/submodules/cPecan @@ -1 +1 @@ -Subproject commit 55920efb7ad7bc9e5da9703668e93bdd0585e353 +Subproject commit 2ec01d160bd4a558ef213a168a57748b82c9313d