From 1bb5353533f789b313a61ab5dc5ed81afcc085e7 Mon Sep 17 00:00:00 2001 From: Sergei Date: Tue, 12 Dec 2023 10:31:23 -0500 Subject: [PATCH] BF updates --- .../SelectionAnalyses/RELAX.bf | 20 ++++++++++++++++--- .../SelectionAnalyses/aBSREL.bf | 8 ++++---- 2 files changed, 21 insertions(+), 7 deletions(-) diff --git a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf index 3a68b4a5d..fcb572573 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf @@ -43,19 +43,21 @@ relax.analysis_description = { Version 3.1 adds LHC + Nedler-Mead initial fit phase and keyword support. Version 3.1.1 adds some bug fixes for better convergence. Version 4.0 adds support for synonymous rate variation. - Version 4.1 adds further support for multiple hit models", + Version 4.1 adds further support for multiple hit models. + Version 4.1.1 adds reporting for convergence diagnostics", terms.io.version : "3.1.1", terms.io.reference : "RELAX: Detecting Relaxed Selection in a Phylogenetic Framework (2015). Mol Biol Evol 32 (3): 820-832", terms.io.authors : "Sergei L Kosakovsky Pond, Ben Murrell, Steven Weaver and Temple iGEM / UCSD viral evolution g", terms.io.contact : "spond@temple.edu", - terms.io.requirements : "in-frame codon alignment and a phylogenetic tree, with at least two groups of branches defined using the {} notation (one group can be defined as all unlabeled branches)" + terms.io.requirements : "in-frame codon alignment and a phylogenetic tree, with at least two groups of branches defined using the {} notation (one group can be defined as all unlabeled branches)", + terms.settings: {} }; relax.json = { terms.json.analysis: relax.analysis_description, terms.json.input: {}, terms.json.fits : {}, terms.json.timers : {}, - terms.json.test_results : {} + terms.json.test_results : {}, }; relax.relaxation_parameter = "relax.K"; @@ -1001,6 +1003,9 @@ function relax.FitMainTestPair (prompt) { if (Abs (relax.best_samples[1] - relax.best_samples[0]) < 5.) { // could be diagnostic of convergence problems + selection.io.json_store_setting (relax.json, "convergence-flat-surface", relax.best_samples[1] - relax.best_samples[0]); + + io.ReportProgressMessageMD("RELAX", "alt-2", "* Potential convergence issues due to flat likelihood surfaces; checking to see whether K > 1 or K < 1 is robustly inferred"); if (relax.fitted.K > 1) { parameters.SetRange (model.generic.GetGlobalParameter (relax.model_object_map ["relax.test"] , terms.relax.k), terms.range01); @@ -1022,6 +1027,8 @@ function relax.FitMainTestPair (prompt) { if (relax.alternative_model.fit.take2 [terms.fit.log_likelihood] > relax.alternative_model.fit[terms.fit.log_likelihood]) { + relax.stash_fitted_K = relax.fitted.K; + io.ReportProgressMessageMD("RELAX", "alt-2", "\n### Potential for highly unreliable K inference due to multiple local maxima in the likelihood function, treat results with caution "); io.ReportProgressMessageMD("RELAX", "alt-2", "> Relaxation parameter reset to opposite mode of evolution from that obtained in the initial optimization."); @@ -1032,6 +1039,11 @@ function relax.FitMainTestPair (prompt) { relax.inferred_distribution = parameters.GetStickBreakingDistribution (models.codon.BS_REL.ExtractMixtureDistribution (relax.model_object_map ["relax.test"])) % 0; selection.io.report_dnds (relax.inferred_distribution); + selection.io.json_store_setting (relax.json, "convergence-unstable", { + {relax.fitted.K, relax.alternative_model.fit.take2 [terms.fit.log_likelihood]} + {relax.stash_fitted_K, relax.alternative_model.fit [terms.fit.log_likelihood]} + + }); io.ReportProgressMessageMD("RELAX", "alt-2", "* The following rate distribution was inferred for **reference** branches"); relax.inferred_distribution_ref = parameters.GetStickBreakingDistribution (models.codon.BS_REL.ExtractMixtureDistribution (relax.model_object_map ["relax.reference"])) % 0; @@ -1182,6 +1194,8 @@ do { if (relax.LRT [terms.LRT] < 0) { io.ReportProgressMessageMD("RELAX", "refit", "* Detected convergence issues (negative LRT). Refitting the alterative/null model pair from a new starting point"); + selection.io.json_store_setting (relax.json, "convergence-negative-lrt", relax.loop_passes); + relax.general_descriptive.fit = relax.null_model.fit; // reset initial conditions for (relax.k = 1; relax.k < relax.numbers_of_tested_groups; relax.k += 1) { // remove constraints on K relax.model_nmsp = relax.model_namespaces[relax.k ]; diff --git a/res/TemplateBatchFiles/SelectionAnalyses/aBSREL.bf b/res/TemplateBatchFiles/SelectionAnalyses/aBSREL.bf index 53270a7db..12f6beef9 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/aBSREL.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/aBSREL.bf @@ -69,11 +69,11 @@ absrel.display_orders = {terms.original_name: -1, absrel.analysis_description = {terms.io.info : "aBSREL (Adaptive branch-site random effects likelihood) uses an adaptive random effects branch-site model framework to test whether each branch has evolved under positive selection, - using a procedure which infers an optimal number of rate categories per branch. v2.3 adds support for SRV. v2.5 adds support for ancestral state reconstruction, identification of sites contributing to selection signal, and some diagnostics. ", + using a procedure which infers an optimal number of rate categories per branch. v2.2 adds support for multiple-hit models. v2.3 adds support for SRV. v2.5 adds support for ancestral state reconstruction, identification of sites contributing to selection signal, and some diagnostics. ", terms.io.version : "2.5", terms.io.reference : " Less Is More: An Adaptive Branch-Site Random Effects Model for Efficient Detection of Episodic Diversifying Selection (2015). - Mol Biol Evol 32 (5): 1342-1353. v2.2 adds support for multiple-hit models. + Mol Biol Evol 32 (5): 1342-1353. ", terms.io.authors : "Sergei L Kosakovsky Pond, Ben Murrell, Steven Weaver and Temple iGEM / UCSD viral evolution group", terms.io.contact : "spond@temple.edu", @@ -703,12 +703,12 @@ for (absrel.branch_id = 0; absrel.branch_id < absrel.branch_count; absrel.branch if (absrel.multi_hit == "Double" || absrel.multi_hit == "Double+Triple") { absrel.report.row [4] = (absrel.final_estimates[utility.getGlobalValue ("terms.parameters.multiple_hit_rate")])[utility.getGlobalValue ("terms.fit.MLE")]; - absrel.branch.delta [absrel.current_branch] = absrel.report.row [3]; + absrel.branch.delta [absrel.current_branch] = absrel.report.row [4]; absrel.offset = 1; } if (absrel.multi_hit == "Double+Triple") { absrel.report.row [5] = (absrel.final_estimates[utility.getGlobalValue ("terms.parameters.triple_hit_rate")])[utility.getGlobalValue ("terms.fit.MLE")]; - absrel.branch.psi [absrel.current_branch] = absrel.report.row [4]; + absrel.branch.psi [absrel.current_branch] = absrel.report.row [5]; absrel.offset += 1; }