[med-svn] [Git][med-team/hyphy][upstream] New upstream version 2.5.68+dfsg
Étienne Mollier (@emollier)
gitlab at salsa.debian.org
Sat Feb 15 18:26:45 GMT 2025
Étienne Mollier pushed to branch upstream at Debian Med / hyphy
Commits:
62b8190b by Étienne Mollier at 2025-02-15T19:18:01+01:00
New upstream version 2.5.68+dfsg
- - - - -
26 changed files:
- CMakeLists.txt
- + Dockerfile
- README.md
- res/TemplateBatchFiles/F_ST.bf
- res/TemplateBatchFiles/SelectionAnalyses/FEL.bf
- res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
- res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf
- res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf
- res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf
- + res/TemplateBatchFiles/libv3/models/codon/MG_GTR.bf
- res/TemplateBatchFiles/libv3/models/codon/MG_REV_PROPERTIES.bf
- res/TemplateBatchFiles/libv3/models/codon/MSS.bf
- res/TemplateBatchFiles/libv3/models/frequencies.bf
- res/TemplateBatchFiles/libv3/models/parameters.bf
- res/TemplateBatchFiles/libv3/tasks/estimators.bf
- + res/TemplateBatchFiles/libv3/tasks/mapping.bf
- res/TemplateBatchFiles/libv3/tasks/trees.bf
- src/core/fisher_exact.cpp
- src/core/global_things.cpp
- src/core/likefunc.cpp
- src/core/likefunc2.cpp
- src/core/matrix.cpp
- src/core/tree_evaluator.cpp
- tests/hbltests/libv3/FEL.wbf
- tests/hbltests/libv3/MEME.wbf
- tests/hbltests/libv3/algae-mtDNA.wbf
Changes:
=====================================
CMakeLists.txt
=====================================
@@ -7,6 +7,7 @@ enable_testing()
set(CMAKE_CONFIGURATION_TYPES Release)
+
#-------------------------------------------------------------------------------
# OPTIONS
#-------------------------------------------------------------------------------
@@ -170,12 +171,14 @@ set_source_files_properties(${SRC_CORE} ${SRC_NEW} {SRC_UTILS} PROPERTIES COMPIL
set(DEFAULT_WARNING_FLAGS " -w -Weffc++ -Wextra -Wall ")
set(DEFAULT_DEBUG_WARNING_FLAGS "-Wall -Wno-int-to-pointer-cast -Wno-conversion-null -Wno-sign-compare")
-set(ADDITIONAL_FLAGS "")
+set(ADDITIONAL_FLAGS "-g")
if (DEBUGFLAGS)
set(ADDITIONAL_FLAGS "-g")
endif()
+#include_directories("/Library/Developer/CommandLineTools/SDKs/MacOSX15.2.sdk/usr/include/c++/v1")
+
if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
execute_process(
@@ -371,6 +374,7 @@ endif(NOT NOZLIB)
+
#-------------------------------------------------------------------------------
# uninstall target
#-------------------------------------------------------------------------------
@@ -678,10 +682,12 @@ makeLink(HYPHYMP hyphy hyphy)
# TESTS
#-------------------------------------------------------------------------------
+set(HYPHY_ENV_ARG "ENV=TOLERATE_NUMERICAL_ERRORS=1;")
+
add_test (NAME UNIT-TESTS COMMAND bash run_unit_tests.sh)
add_test (CODON HYPHYMP tests/hbltests/SimpleOptimizations/SmallCodon.bf)
-add_test (PROTEIN HYPHYMP CPU=1 tests/hbltests/SimpleOptimizations/IntermediateProtein.bf)
+add_test (PROTEIN HYPHYMP tests/hbltests/SimpleOptimizations/IntermediateProtein.bf)
add_test (MTCODON HYPHYMP tests/hbltests/libv3/mtDNA-code.wbf)
add_test (ALGAE HYPHYMP tests/hbltests/libv3/algae-mtDNA.wbf)
add_test (CILIATES HYPHYMP tests/hbltests/libv3/ciliate-code.wbf)
@@ -697,8 +703,12 @@ add_test (RELAX HYPHYMP tests/hbltests/libv3/RELAX.wbf)
add_test (FUBAR HYPHYMP tests/hbltests/libv3/FUBAR.wbf)
add_test (BGM HYPHYMP tests/hbltests/libv3/BGM.wbf)
add_test (CONTRAST-FEL HYPHYMP tests/hbltests/libv3/CFEL.wbf)
-add_test (GARD HYPHYMP tests/hbltests/libv3/GARD.wbf)
add_test (FADE HYPHYMP tests/hbltests/libv3/FADE.wbf)
+add_test (NAME GARD
+ COMMAND HYPHYMP
+ tests/hbltests/libv3/GARD.wbf
+ ${HYPHY_ENV_ARG}
+ )
add_test (ABSREL HYPHYMP tests/hbltests/libv3/ABSREL.wbf)
#add_test (FMM HYPHYMP tests/hbltests/libv3/FMM.wbf)
#add_test (NAME ABSREL-MH COMMAND HYPHYMP tests/hbltests/libv3/ABSREL-MH.wbf)
=====================================
Dockerfile
=====================================
@@ -0,0 +1,26 @@
+FROM ubuntu:24.04
+
+# Install build dependencies
+RUN apt-get update && \
+ apt-get install --no-install-recommends -y build-essential cmake
+
+# Add a non-root user
+RUN groupadd -r hyphyuser && useradd -r -g hyphyuser hyphyuser
+
+# Create a directory for the project
+WORKDIR /hyphy
+
+# Copy project files
+COPY ./cmake /hyphy/cmake
+COPY ./src /hyphy/src
+COPY ./contrib /hyphy/contrib
+COPY ./res /hyphy/res
+COPY CMakeLists.txt /hyphy
+COPY ./tests /hyphy/tests
+
+RUN chown -R hyphyuser:hyphyuser .
+
+# Install project
+RUN cmake . && make -j install
+
+USER hyphyuser
=====================================
README.md
=====================================
@@ -9,6 +9,20 @@ HyPhy is an open-source software package for the analysis of genetic sequences u
#### Install
`conda install -c bioconda hyphy`
+#### Running with Docker
+You can also run HyPhy without having to install it on your system using the provided Dockerfile. Following the below instructions starts an interactive Docker container where HyPhy is already available.
+
+If you don't have Docker installed, first look at [Docker Desktop](https://docs.docker.com/desktop/).
+
+**Please note you must change this code snippet to point to the appropriate location for your input data. This will be made available inside the container at `/hyphy/data`.**
+
+```
+git clone https://github.com/veg/hyphy.git
+cd hyphy
+docker build -t hyphy:latest .
+docker run --rm -v [path-to-your-input-data]:/hyphy/data -it hyphy:latest
+```
+
#### Run with Command Line Arguments
`hyphy <method_name> --alignment <path_to_alignment_file> <additional_method_specific_arguments>`
+ _`<method_name>` is the name of the analysis you wish to run (can be: absrel, bgm, busted, fade, fel, fubar, gard, meme, relax or slac)_
@@ -86,6 +100,27 @@ And then run make install to install the software
+ libhyphy_mp.(so/dylib/dll) will be installed at `/location/of/choice/lib`
+ HyPhy's standard library of batchfiles will go into `/location/of/choice/lib/hyphy`
+#### Building HyPhy for Web Assembly
+
+Use [Emscripten](https://emscripten.org/docs/compiling/Building-Projects.html) to produce web assembly files, which can be run entirely within a modern browser. See https://observablehq.com/@spond/hyphy-biowasm for an example.
+
+```
+emcmake cmake -DCMAKE_EXE_LINKER_FLAGS="-sTOTAL_STACK=2097152 -02 -sASSERTIONS=1 -sMODULARIZE=1 -sALLOW_MEMORY_GROWTH -sFORCE_FILESYSTEM=1 -sEXIT_RUNTIME=0 -s EXPORTED_RUNTIME_METHODS=["callMain","FS","PROXYFS","WORKERFS","UTF8ToString","getValue","AsciiToString"] -lworkerfs.js -lproxyfs.js -s INVOKE_RUN=0 -s ENVIRONMENT="web,worker" ${EM_FLAGS//-s /-s} -fwasm-exceptions --preload-file res@/hyphy --preload-file tests/hbltests@/tests"
+```
+
+```
+emmake make -j hyphy
+```
+
+This creates
+
+```
+hyphy.js
+hyphy.wasm
+hyphy.data
+```
+
+Which should be served from the same directory.
#### Testing
=====================================
res/TemplateBatchFiles/F_ST.bf
=====================================
@@ -242,42 +242,32 @@ if (distanceChoice == "Full likelihood") {
{
for (j = 0; j<=i; j = j+1)
{
- if (dataType)
- {
+ if (dataType == "Codon") {
DataSetFilter twoSpecFilter = CreateFilter (filteredData,3,"",(speciesIndex==i+1)||(speciesIndex==j),GeneticCodeExclusions);
}
- else
- {
+ else {
DataSetFilter twoSpecFilter = CreateFilter (filteredData,1,"",(speciesIndex==i+1)||(speciesIndex==j));
}
- if (FREQUENCY_SENSITIVE)
- {
- if (USE_POSITION_SPECIFIC_FREQS)
- {
+ if (FREQUENCY_SENSITIVE) {
+ if (USE_POSITION_SPECIFIC_FREQS) {
HarvestFrequencies (vectorOfFrequencies,filteredData,3,1,1);
}
- else
- {
+ else {
HarvestFrequencies (vectorOfFrequencies,twoSpecFilter,1,1,0);
}
}
- if (FREQUENCY_SENSITIVE)
- {
+ if (FREQUENCY_SENSITIVE) {
MULTIPLY_BY_FREQS = PopulateModelMatrix ("modelMatrix",vectorOfFrequencies);
- if (dataType)
- {
+ if (dataType == "Codon") {
codonFrequencies = BuildCodonFrequencies (vectorOfFrequencies);
Model pairModel = (modelMatrix, codonFrequencies, MULTIPLY_BY_FREQS);
}
- else
- {
+ else {
Model pairModel = (modelMatrix, vectorOfFrequencies, MULTIPLY_BY_FREQS);
}
}
- else
- {
- if (i+j==0)
- {
+ else {
+ if (i+j==0) {
MULTIPLY_BY_FREQS = PopulateModelMatrix ("modelMatrix",equalFreqs);
Model pairModel = (modelMatrix, equalFreqs, MULTIPLY_BY_FREQS);
}
=====================================
res/TemplateBatchFiles/SelectionAnalyses/FEL.bf
=====================================
@@ -610,7 +610,6 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
GetString (lfInfo, ^lf,-1);
ExecuteCommands (filter_data);
__make_filter ((lfInfo["Datafilters"])[0]);
-
utility.SetEnvVariable ("USE_LAST_RESULTS", TRUE);
@@ -724,9 +723,6 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
}
}
- //Export (lfe, ^lf);
- //console.log (lfe);
- //assert (0);
Optimize (results, ^lf
, {
@@ -737,6 +733,7 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
);
+
if (^"fel.ci") {
if (!sim_mode) {
if (^"fel.srv") {
@@ -776,10 +773,13 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
}
}
- ^"fel.alpha_scaler" = (^"fel.alpha_scaler" + 3*^"fel.beta_scaler_test")/4;
+ ^"fel.alpha_scaler" = (Min (^"fel.alpha_scaler", 100) + 3*Min(^"fel.beta_scaler_test",100))/4;
parameters.SetConstraint ("fel.beta_scaler_test","fel.alpha_scaler", "");
- Optimize (results, ^lf);
+
+ Optimize (results, ^lf, {
+ "OPTIMIZATION_METHOD" : "coordinate-wise"
+ });
if (sim_mode) {
=====================================
res/TemplateBatchFiles/SelectionAnalyses/MEME.bf
=====================================
@@ -284,7 +284,7 @@ if (meme.run_full_mg94) {
terms.run_options.retain_lf_object: TRUE,
terms.run_options.apply_user_constraints: meme.define_constraints,
terms.run_options.optimization_settings: {
- "OPTIMIZATION_METHOD" : "coordinate-wise"
+ //"OPTIMIZATION_METHOD" : "coordinate-wise"
}
}, meme.partitioned_mg_results);
@@ -946,6 +946,7 @@ lfunction meme.apply_initial_estimates (spec) {
for (t,d;in;spec) {
if (None != d["init"]) {
if (Type (d["scaler"]) == "String") {
+ //console.log (d["scaler"] + " = " + d["init"]);
^(d["scaler"]) = d["init"];
}
}
@@ -983,6 +984,7 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
rate.fel.alpha = (((scaler_mapping['BG']))[^"terms.parameters.synonymous_rate"])["scaler"];
rate.fel.beta_fg = (((scaler_mapping['FEL-FG']))[^"terms.parameters.nonsynonymous_rate"])["scaler"];
rate.fel.beta_bg = (((scaler_mapping['BG']))[^"terms.parameters.nonsynonymous_rate"])["scaler"];
+
meme.apply_initial_estimates (scaler_mapping["BG"]);
meme.apply_initial_estimates (scaler_mapping["FEL-FG"]);
@@ -1002,20 +1004,22 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
bsrel_tree_id = (lfInfo["Trees"])[0];
+ init_alpha = ^rate.fel.alpha;
+
character_map = None;
if (rate_count == 2) {
// separate clauses for initial parameter values
if ( ^rate.fel.alpha < ^rate.fel.beta_fg) { // FEL estimated beta > alpha
- meme.set_ith_omega (scaler_mapping, 0, 0.25, 0.25);
+ meme.set_ith_omega (scaler_mapping, 0, 0.0, 0.25);
} else { // FEL estimated beta <= alpha
- if (^rate.fel.alpha > 0) {
+ if (^rate.fel.alpha > 1e-5) {
omega_rate = ^rate.fel.beta_fg / ^rate.fel.alpha;
} else {
omega_rate = 1;
}
meme.set_ith_omega (scaler_mapping, 0, omega_rate, 0.75);
- meme.set_ith_omega (scaler_mapping, 1, 1.5*^rate.fel.alpha , None);
+ meme.set_ith_omega (scaler_mapping, 1, Max (0.1, 1.5*^rate.fel.alpha) , None);
}
rate.meme.mixture_weight1 = ((scaler_mapping["OMEGA"])[^"terms.parameters.weights"])[0];
@@ -1092,14 +1096,14 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
meme.set_ith_omega (scaler_mapping, 1, 0.5, 0.5);
meme.set_ith_omega (scaler_mapping, 2, 0.5, None);
} else { // FEL estimated beta <= alpha
- if (^rate.fel.alpha > 0) {
+ if (^rate.fel.alpha > 1e-5) {
omega_rate = ^rate.fel.beta_fg / ^rate.fel.alpha;
} else {
omega_rate = 1;
}
meme.set_ith_omega (scaler_mapping, 0, omega_rate, 0.5);
meme.set_ith_omega (scaler_mapping, 1, Min (1, Max (0.5,omega_rate*2)), 0.5);
- meme.set_ith_omega (scaler_mapping, 2, 1.5*^rate.fel.alpha , None);
+ meme.set_ith_omega (scaler_mapping, 2, Max (0.1, 1.5*^rate.fel.alpha) , None);
}
rate.meme.mixture_weight1 = ((scaler_mapping["OMEGA"])[^"terms.parameters.weights"])[0];
@@ -1191,7 +1195,7 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
meme.set_ith_omega (scaler_mapping, 2, 0.5, 0.5);
meme.set_ith_omega (scaler_mapping, 3, ^rate.fel.beta_fg, None);
} else { // FEL estimated beta <= alpha
- if (^rate.fel.alpha > 0) {
+ if (^rate.fel.alpha > 1e-5) {
omega_rate = ^rate.fel.beta_fg / ^rate.fel.alpha;
} else {
omega_rate = 1;
@@ -1199,7 +1203,7 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
meme.set_ith_omega (scaler_mapping, 0, omega_rate, 0.5);
meme.set_ith_omega (scaler_mapping, 1, Min (1, omega_rate*0.25), 0.5);
meme.set_ith_omega (scaler_mapping, 2, Min (1, 0.5*omega_rate) , 0.5);
- meme.set_ith_omega (scaler_mapping, 3, 1.5*^rate.fel.alpha , None);
+ meme.set_ith_omega (scaler_mapping, 3, Max (0.1, 1.5*^rate.fel.alpha) , None);
}
rate.meme.mixture_weight1 = ((scaler_mapping["OMEGA"])[^"terms.parameters.weights"])[0];
@@ -1309,7 +1313,8 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
- ^rate.fel.alpha = ^rate.fel.alpha;
+ ^rate.fel.alpha = Min (100,^rate.fel.alpha);
+
// SLKP 20201028 : without this, there are occasional initialization issues with
// the likelihood function below
@@ -1317,9 +1322,10 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
Optimize (results, ^lf_bsrel, {
"OPTIMIZATION_METHOD" : "nedler-mead",
"OPTIMIZATION_START_GRID" : initial_guess_grid
- //"OPTIMIZATION_PRECISION" : 1e-5,
+ //,"OPTIMIZATION_PRECISION" : 1e-5,
});
+
/*after_opt = {
"alpha" : Eval("meme.site_alpha"),
@@ -1337,7 +1343,14 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
alternative = estimators.ExtractMLEsOptions (lf_bsrel, model_mapping, {^"terms.globals_only" : TRUE});
alternative [utility.getGlobalValue("terms.fit.log_likelihood")] = results[1][0];
-
+ /*
+ if (^rate.fel.alpha / init_alpha > 5) {
+ console.log ("\n" + ^rate.fel.alpha / init_alpha);
+ console.log (fel["global"]);
+ console.log (alternative["global"]);
+ assert (0);
+ }
+ */
/*
Export (lfe, ^lf_bsrel);
console.log (lfe);
@@ -1422,7 +1435,7 @@ lfunction meme.handle_a_site (lf_fel, lf_bsrel, filter_data, partition_index, pa
}
^rate.fel.beta_fg := ^rate.fel.alpha;
- Optimize (results, ^lf_bsrel, {"OPTIMIZATION_METHOD" : "nedler-mead"});
+ Optimize (results, ^lf_bsrel, {"OPTIMIZATION_METHOD" : "coordinate-wise"});
if (sim_mode) {
return lrt - 2*results[1][0];
=====================================
res/TemplateBatchFiles/SelectionAnalyses/RELAX.bf
=====================================
@@ -276,7 +276,7 @@ relax.do_srv = io.SelectAnOption (
if (relax.do_srv == "Branch-site") {
relax.do_bs_srv = TRUE;
relax.do_srv = TRUE;
- (relax.json[terms.json.analysis])[terms.settings] = "Branch-site synonymous rate variation";
+ selection.io.json_store_setting (relax.json, "SRV type", "Branch-site synonymous rate variation");
} else {
if (relax.do_srv == "Yes") {
relax.do_bs_srv = FALSE;
=====================================
res/TemplateBatchFiles/SelectionAnalyses/error-filter.bf
=====================================
@@ -274,10 +274,14 @@ for (s;in;TipName (efilter.T,-1)) {
fprintf (efilter.path, ">", s, "\n", efilter.sequences[s], "\n");
}
+utility.ToggleEnvVariable ("INCLUDE_MODEL_SPECS", TRUE);
+
if (efilter.input[terms.json.partition_count] == 1) {
fprintf (efilter.path, "\n", Format (efilter.T,1,0));
}
+utility.ToggleEnvVariable ("INCLUDE_MODEL_SPECS", None);
+
io.ReportProgressMessageMD ("efilter","results","\nMasked a total of **`efilter.total_filtered`** or `Format (efilter.total_filtered/efilter.input[terms.json.sequences]/efilter.input[terms.json.sites]*100,8,3)`% sites");
=====================================
res/TemplateBatchFiles/SelectionAnalyses/modules/shared-load-file.bf
=====================================
@@ -266,11 +266,13 @@ function load_file (prefix) {
},
...
*/
-
- for (_key_, _value_; in; partitions_and_trees) {
- (partitions_and_trees[_key_])[utility.getGlobalValue("terms.data.filter_string")] =
- selection.io.adjust_partition_string (_value_[utility.getGlobalValue("terms.data.filter_string")], 3*codon_data_info[utility.getGlobalValue("terms.data.sites")]);
- }
+
+
+ for (_key_, _value_; in; partitions_and_trees) {
+ (partitions_and_trees[_key_])[utility.getGlobalValue("terms.data.filter_string")] =
+ selection.io.adjust_partition_string (_value_[utility.getGlobalValue("terms.data.filter_string")], 3*codon_data_info[utility.getGlobalValue("terms.data.sites")]);
+ }
+
}
@@ -453,7 +455,8 @@ function doGTR (prefix) {
terms.nucleotideRate ("A","C") : { utility.getGlobalValue ("terms.fit.MLE") : 0.25} ,
terms.nucleotideRate ("A","T") : { utility.getGlobalValue ("terms.fit.MLE") : 0.25} ,
terms.nucleotideRate ("C","G") : { utility.getGlobalValue ("terms.fit.MLE") : 0.25} ,
- terms.nucleotideRate ("G","T") : { utility.getGlobalValue ("terms.fit.MLE") : 0.25}
+ terms.nucleotideRate ("G","T") : { utility.getGlobalValue ("terms.fit.MLE") : 0.25} ,
+ terms.nucleotideRate ("C","T") : { utility.getGlobalValue ("terms.fit.MLE") : 1.00}
}
});
@@ -525,7 +528,12 @@ function doGTR (prefix) {
for (i = 0; i < partition_count; i+=1) {
(partitions_and_trees[i])[^"terms.data.tree"] = trees[i];
}
+
store_tree_information ();
+
+ for (index, tree; in; trees) {
+ trees.AdjustZerosToNearZeros (tree, (gtr_results[^"terms.branch_length"])[index]);
+ }
} else {
if (kill0 == "Constrain") {
zero_branch_length_constrain = ^"namespace_tag"+".constrain_zero_branches";
=====================================
res/TemplateBatchFiles/libv3/models/codon/MG_GTR.bf
=====================================
@@ -0,0 +1,114 @@
+RequireVersion ("2.5.60");
+
+LoadFunctionLibrary("../codon.bf");
+LoadFunctionLibrary("../DNA.bf");
+LoadFunctionLibrary("../parameters.bf");
+LoadFunctionLibrary("../frequencies.bf");
+LoadFunctionLibrary("../../UtilityFunctions.bf");
+LoadFunctionLibrary("MG_REV.bf");
+LoadFunctionLibrary("../protein.bf");
+
+/** @module models.codon.GTR */
+//----------------------------------------------------------------------------------------------------------------
+
+
+lfunction model.codon.MG_GTR.prompt_and_define (type, code) {
+ console.log (type);
+ return models.codon.MG_GTR.ModelDescription(type, code);
+}
+
+lfunction models.codon.MG_GTR.ModelDescription(type, code) {
+
+
+ io.CheckAssertion ("`&type` == '`^'terms.global'`'", "MG_GTR only supports `^'terms.global'` type");
+
+ // piggyback on the standard MG_REV model for most of the code
+
+ mg_base = models.codon.MG_REV.ModelDescription (type, code);
+ mg_base[utility.getGlobalValue("terms.description")] = "The Muse-Gaut 94 codon-substitution model coupled with the general time reversible (GTR) model of nucleotide substitution, which defines a separate dN/dS ratio for each pair of amino-acids";
+ mg_base[utility.getGlobalValue("terms.model.q_ij")] = "models.codon.MG_GTR._GenerateRate";
+ mg_base[utility.getGlobalValue("terms.model.post_definition")] = "models.codon.MG_GTR.post_definition";
+ return mg_base;
+}
+
+lfunction models.codon.MG_GTR._GenerateRate(fromChar, toChar, namespace, model_type, model) {
+
+ return models.codon.MG_GTR._GenerateRate_generic (fromChar, toChar, namespace, model_type,
+ model[utility.getGlobalValue("terms.translation_table")],
+ "alpha", utility.getGlobalValue("terms.parameters.synonymous_rate"),
+ "omega", utility.getGlobalValue("terms.parameters.omega_ratio"),
+ );
+}
+
+/**
+ * @name models.codon.MG_GTR._GenerateRate
+ * @param {Number} fromChar
+ * @param {Number} toChar
+ * @param {String} namespace
+ * @param {String} model_type
+ * @param {Matrix} _tt - translation table
+ */
+
+
+lfunction models.codon.MG_GTR._GenerateRate_generic (fromChar, toChar, namespace, model_type, _tt, alpha, alpha_term, omega, omega_term) {
+
+ _GenerateRate.p = {};
+ _GenerateRate.diff = models.codon.diff.complete(fromChar, toChar);
+ diff_count = utility.Array1D (_GenerateRate.diff);
+
+ if (diff_count == 1) {
+
+ _GenerateRate.p[model_type] = {};
+ _GenerateRate.p[utility.getGlobalValue("terms.global")] = {};
+
+ nuc_rate = "";
+
+ for (i = 0; i < diff_count; i += 1) {
+ if ((_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")] > (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")]) {
+ nuc_p = "theta_" + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")] + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")];
+ } else {
+ nuc_p = "theta_" + (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")] +(_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")];
+ }
+ nuc_p = parameters.ApplyNameSpace(nuc_p, namespace);
+ (_GenerateRate.p[utility.getGlobalValue("terms.global")])[terms.nucleotideRateReversible((_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.from")], (_GenerateRate.diff[i])[utility.getGlobalValue("terms.diff.to")])] = nuc_p;
+
+ nuc_rate = parameters.AppendMultiplicativeTerm (nuc_rate, nuc_p);
+ }
+
+ rate_entry = nuc_rate;
+
+ if (_tt[fromChar] != _tt[toChar]) {
+
+ if (_tt[fromChar] < _tt[toChar]) {
+ aa_pair = _tt[fromChar] + " and " + _tt[toChar];
+ aa_tag = "_" + _tt[fromChar] + _tt[toChar];
+ } else {
+ aa_pair = _tt[toChar] + " and " + _tt[fromChar];
+ aa_tag = "_" + _tt[toChar] + _tt[fromChar];
+ }
+
+ omega_p = parameters.ApplyNameSpace (omega + aa_tag , namespace);
+ (_GenerateRate.p[model_type]) [^"terms.parameters.omega_ratio" + " for " + aa_pair] = omega_p;
+
+ rate_entry = parameters.AppendMultiplicativeTerm (rate_entry, omega_p);
+ } else {
+ _GenerateRate.p[utility.getGlobalValue("terms.model.rate_entry")] = nuc_rate;
+ }
+
+ _GenerateRate.p[utility.getGlobalValue("terms.model.rate_entry")] = rate_entry;
+ }
+ return _GenerateRate.p;
+}
+
+lfunction models.codon.MG_GTR.post_definition (model) {
+
+ for (id; in ; model.GetParameters_RegExp (model, ^"terms.parameters.omega_ratio")) {
+ parameters.SetValue(id, 0.25);
+ }
+
+ models.generic.post.definition (model);
+}
+
+lfunction models.codon.MG_GTR.set_branch_length(model, value, parameter) {
+ return models.codon.MG_REV.set_branch_length(model,value,parameter);
+}
=====================================
res/TemplateBatchFiles/libv3/models/codon/MG_REV_PROPERTIES.bf
=====================================
@@ -252,7 +252,211 @@ terms.model.MG_REV_PROPERTIES.predefined = {
{0.2306464200526206}
}
},
- "Random" : {
+ "Random-2" : {
+ "Random Factor 1":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ },
+ "Random Factor 2":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ }
+ },
+ "Random-3" : {
+ "Random Factor 1":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ },
+ "Random Factor 2":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ },
+ "Random Factor 3":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ }
+ },
+ "Random-4" : {
+ "Random Factor 1":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ },
+ "Random Factor 2":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ },
+ "Random Factor 3":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ },
+ "Random Factor 4":{
+ "A":Random(-2,2),
+ "C":Random(-2,2),
+ "D":Random(-2,2),
+ "E":Random(-2,2),
+ "F":Random(-2,2),
+ "G":Random(-2,2),
+ "H":Random(-2,2),
+ "I":Random(-2,2),
+ "K":Random(-2,2),
+ "L":Random(-2,2),
+ "M":Random(-2,2),
+ "N":Random(-2,2),
+ "P":Random(-2,2),
+ "Q":Random(-2,2),
+ "R":Random(-2,2),
+ "S":Random(-2,2),
+ "T":Random(-2,2),
+ "V":Random(-2,2),
+ "W":Random(-2,2),
+ "Y":Random(-2,2)
+ }
+ },
+ "Random-5" : {
"Random Factor 1":{
"A":Random(-2,2),
"C":Random(-2,2),
@@ -376,7 +580,10 @@ lfunction model.codon.MG_REV_PROPERTIES.prompt_and_define (type, code) {
{
"Atchley":"Use the five properties derived from a factor analysis of 500 amino-acid properties [Table 2 in PNAS (2005) 102(18) 6395-6400 doi: 10.1073/pnas.0408677102]",
"LCAP":"Use the five properties defined in the Conant and Stadler LCAP model [Mol Biol Evol (2009) 26 (5): 1155-1161. doi: 10.1093/molbev/msp031]",
- "Random" : "Random properties (for null hypothesis testing)",
+ "Random-2" : "Two random properties (for null hypothesis testing)",
+ "Random-3" : "Three random properties (for null hypothesis testing)",
+ "Random-4" : "Four random properties (for null hypothesis testing)",
+ "Random-5" : "Five random properties (for null hypothesis testing)",
"Custom":"Load the set of properties from a file"
},
"The set of properties to use in the model");
@@ -408,6 +615,7 @@ lfunction models.codon.MG_REV_PROPERTIES.ModelDescription(type, code, properties
mg_base[utility.getGlobalValue("terms.model.residue_name_map")] = parameters.ValidateIDs (utility.Keys (mg_base[utility.getGlobalValue("terms.model.residue_properties")]));
mg_base[utility.getGlobalValue("terms.model.post_definition")] = "models.codon.MG_REV_PROPERTIES.post_definition";
mg_base[utility.getGlobalValue("terms.model.set_branch_length")] = "models.codon.MG_REV_PROPERTIES.set_branch_length";
+
return mg_base;
}
=====================================
res/TemplateBatchFiles/libv3/models/codon/MSS.bf
=====================================
@@ -14,8 +14,11 @@ terms.model.MSS.syn_rate_within = " within codon class ";
terms.model.MSS.syn_rate_between = " between codon classes ";
terms.model.MSS.between = "synonymous rate between codon classes";
terms.model.MSS.neutral = "neutral reference";
+terms.model.MSS.empirical = "empirical";
terms.model.MSS.codon_classes = "codon classes";
+terms.model.MSS.codon_pairs = "codon pairs";
terms.model.MSS.normalize = "normalize rates";
+terms.model.MSS.non_syn_reference = "non-synonymous reference";
//----------------------------------------------------------------------------------------------------------------
@@ -32,7 +35,9 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) {
{"SynREV2g", "Each pair of synonymous codons mapping to the same amino-acid class and separated by a transition have a separate substitution rate (Valine == neutral). All between-class synonymous substitutions share a rate."}
{"SynREVCodon", "Each codon pair that is exchangeable gets its own substitution rate (fully estimated, mean = 1)"}
{"Random", "Random partition (specify how many classes; largest class = neutral)"}
+ {"Empirical", "Load a TSV file with an empirical rate estimate for each codon pair"}
{"File", "Load a TSV partition from file (prompted for neutral class)"}
+ {"Codon-file", "Load a TSV partition for pairs of codons from a file (prompted for neutral class)"}
},
"Synonymous Codon Class Definitions");
@@ -69,6 +74,17 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) {
}
if (partitioning_option == "SynREV" || partitioning_option == "SynREVFull" || partitioning_option == "SynREVCodon" ) {
+
+ KeywordArgument ("mss-reference-rate", "Normalize relative to these rates", "synonymous");
+
+ reference_option = io.SelectAnOption (
+ {
+ {"synonymous", "Synonymous rates have mean 1."},
+ {"non-synonymous", "Non-synonymous rates are set to 1."}
+ },
+ "Rate normalization"
+ );
+
bins = {};
mapping = {};
mapping_codon = {};
@@ -82,9 +98,9 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) {
return models.codon.MSS.ModelDescription(type, code,
{
^"terms.model.MSS.codon_classes" : mapping,
- ^"terms.model.MSS.normalize" : TRUE
- //^"terms.model.MSS.neutral" : "V"
- }
+ ^"terms.model.MSS.normalize" : reference_option == "synonymous",
+ ^"terms.model.MSS.non_syn_reference" : reference_option != "synonymous"
+ }
);
}
if (partitioning_option == "SynREVFull") {
@@ -92,8 +108,14 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) {
{^"terms.model.MSS.codon_classes" : mapping}
);
}
+
+
return models.codon.MSS.ModelDescription(type, code,
- {^"terms.model.MSS.codon_classes" : mapping_codon, ^"terms.model.MSS.normalize" : TRUE}
+ {
+ ^"terms.model.MSS.codon_classes" : mapping_codon,
+ ^"terms.model.MSS.normalize" : reference_option == "synonymous",
+ ^"terms.model.MSS.non_syn_reference" : reference_option != "synonymous"
+ }
);
}
@@ -182,6 +204,18 @@ lfunction model.codon.MSS.prompt_and_define_freq (type, code, freq) {
return models.codon.MSS.ModelDescription(type, code, models.codon.MSS.LoadClasses (null));
}
+ if (partitioning_option == "Codon-file") {
+ KeywordArgument ("mss-file", "File defining the model partition");
+ KeywordArgument ("mss-neutral", "Designation for the neutral substitution rate");
+ return models.codon.MSS.ModelDescription(type, code, models.codon.MSS.LoadClassesCodon (null));
+ }
+
+ if (partitioning_option == "Empirical") {
+ KeywordArgument ("mss-file", "File defining empirical rates for each pair of codons");
+ return models.codon.MSS.ModelDescription(type, code, models.codon.MSS.LoadEmpiricalRates (null));
+ }
+
+
return {};
}
@@ -209,6 +243,10 @@ lfunction models.codon.MSS.ModelDescription(type, code, codon_classes) {
m[utility.getGlobalValue("terms.model.q_ij")] = "models.codon.MSS._GenerateRate";
m[utility.getGlobalValue("terms.model.MSS.codon_classes")] = codon_classes [^"terms.model.MSS.codon_classes"];
m[utility.getGlobalValue("terms.model.MSS.neutral")] = codon_classes [^"terms.model.MSS.neutral"];
+ m[utility.getGlobalValue("terms.model.MSS.empirical")] = codon_classes [^"terms.model.MSS.empirical"];
+ m[utility.getGlobalValue("terms.model.MSS.codon_pairs")] = codon_classes [^"terms.model.MSS.codon_pairs"];
+ m[utility.getGlobalValue("terms.model.MSS.non_syn_reference")] = codon_classes [^"terms.model.MSS.non_syn_reference"];
+
if (codon_classes/utility.getGlobalValue("terms.model.MSS.between")) {
m[^"terms.model.MSS.between"] = codon_classes [^"terms.model.MSS.between"];
}
@@ -225,6 +263,9 @@ lfunction models.codon.MSS.ModelDescription(type, code, codon_classes) {
lfunction models.codon.MSS.post_definition (model) {
// TBD
+ if (model[^"terms.model.MSS.non_syn_reference"]) {
+ return;
+ }
rates = model.GetParameters_RegExp (model,"^" + utility.getGlobalValue ("terms.parameters.synonymous_rate"));
D = utility.Array1D (rates);
w = 1 / D;
@@ -257,6 +298,9 @@ lfunction models.codon.MSS._GenerateRate (fromChar, toChar, namespace, model_typ
alpha_term = utility.getGlobalValue ("terms.parameters.synonymous_rate");
beta_term = utility.getGlobalValue ("terms.parameters.nonsynonymous_rate");
nr = model[utility.getGlobalValue("terms.model.MSS.neutral")];
+ empirical = model[utility.getGlobalValue("terms.model.MSS.empirical")];
+ codon_pairs = model[utility.getGlobalValue("terms.model.MSS.codon_pairs")];
+ non_syn_ref = model[utility.getGlobalValue("terms.model.MSS.non_syn_reference")];
omega = "omega";
alpha = "alpha";
beta = "beta";
@@ -288,73 +332,98 @@ lfunction models.codon.MSS._GenerateRate (fromChar, toChar, namespace, model_typ
if (_tt[fromChar] != _tt[toChar]) {
- if (model_type == utility.getGlobalValue("terms.global")) {
- aa_rate = parameters.ApplyNameSpace(omega, namespace);
- (_GenerateRate.p[model_type])[omega_term] = aa_rate;
- } else {
- aa_rate = beta;
- (_GenerateRate.p[model_type])[beta_term] = aa_rate;
- }
- rate_entry += "*" + aa_rate;
+ if (!non_syn_ref) {
+ if (model_type == utility.getGlobalValue("terms.global")) {
+ aa_rate = parameters.ApplyNameSpace(omega, namespace);
+ (_GenerateRate.p[model_type])[omega_term] = aa_rate;
+ } else {
+ aa_rate = beta;
+ (_GenerateRate.p[model_type])[beta_term] = aa_rate;
+ }
+ rate_entry += "*" + aa_rate;
+ }
} else {
- class_from = (model[^"terms.model.MSS.codon_classes"])[fromChar];
- class_to = (model[^"terms.model.MSS.codon_classes"])[toChar];
-
-
- if ((Abs (class_from) && Abs (class_to)) == FALSE) {
- class_from = (model[^"terms.model.MSS.codon_classes"])[fromChar+toChar];
- class_to = class_from;
- }
-
- assert (Abs (class_from) && Abs (class_to), "The neutral class for `fromChar` to `toChar` is not specified in the model definition");
-
- if (class_from == class_to) {
- if (class_from == nr) {
- if (model_type == utility.getGlobalValue("terms.local")) {
- codon_rate = alpha + "_" + class_from;
- (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_within" + class_from] = codon_rate;
- rate_entry += "*" + codon_rate;
- } else {
- rate_entry = nuc_rate;
- }
+ if (empirical) {
+ if (fromChar < toChar) {
+ key = fromChar + "|" + toChar;
} else {
- if (model_type == utility.getGlobalValue("terms.local")) {
- codon_rate = alpha + "_" + class_from;
+ key = toChar + "|" + fromChar;
+ }
+
+ assert (model[^"terms.model.MSS.codon_classes"] / key, "Rate for codon pair `key` was missing from the empirical file definition");
+ rate_entry += "*" + (model[^"terms.model.MSS.codon_classes"])[key];
+
+ } else {
+ if (codon_pairs) {
+ if (fromChar < toChar) {
+ key = fromChar + "|" + toChar;
} else {
- codon_rate = parameters.ApplyNameSpace(alpha + "_" + class_from, namespace);
+ key = toChar + "|" + fromChar;
}
- (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_within" + class_from] = codon_rate;
- rate_entry += "*" + codon_rate;
+ assert (model[^"terms.model.MSS.codon_classes"] / key, "Rate for codon pair `key` was missing from the empirical file definition");
+ class_from = (model[^"terms.model.MSS.codon_classes"])[key];
+ class_to = class_from;
+ } else {
+ class_from = (model[^"terms.model.MSS.codon_classes"])[fromChar];
+ class_to = (model[^"terms.model.MSS.codon_classes"])[toChar];
}
- } else {
- if (class_from > class_to) {
- codon_rate = class_to;
+
+
+ if ((Abs (class_from) && Abs (class_to)) == FALSE) {
+ class_from = (model[^"terms.model.MSS.codon_classes"])[fromChar+toChar];
class_to = class_from;
- class_from = codon_rate;
}
- if (Abs (between_rate)) {
- if (model_type == utility.getGlobalValue("terms.local")) {
- codon_rate = between_rate;
+
+ assert (Abs (class_from) && Abs (class_to), "The neutral class for `fromChar` to `toChar` is not specified in the model definition");
+
+ if (class_from == class_to) {
+ if (class_from == nr) {
+ if (model_type == utility.getGlobalValue("terms.local")) {
+ codon_rate = alpha + "_" + class_from;
+ (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_within" + class_from] = codon_rate;
+ rate_entry += "*" + codon_rate;
+ } else {
+ rate_entry = nuc_rate;
+ }
} else {
- codon_rate = parameters.ApplyNameSpace(between_rate, namespace);
+ if (model_type == utility.getGlobalValue("terms.local")) {
+ codon_rate = alpha + "_" + class_from;
+ } else {
+ codon_rate = parameters.ApplyNameSpace(alpha + "_" + class_from, namespace);
+ }
+ (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_within" + class_from] = codon_rate;
+ rate_entry += "*" + codon_rate;
}
- (_GenerateRate.p[model_type])[^"terms.model.MSS.between"] = codon_rate;
-
} else {
- if (class_from + class_to == nr) {
- //console.log ("NEUTRAL");
- codon_rate = 1;
- } else {
+ if (class_from > class_to) {
+ codon_rate = class_to;
+ class_to = class_from;
+ class_from = codon_rate;
+ }
+ if (Abs (between_rate)) {
if (model_type == utility.getGlobalValue("terms.local")) {
- codon_rate = alpha + "_" + class_from + "_" + class_to;
+ codon_rate = between_rate;
} else {
- codon_rate = parameters.ApplyNameSpace(alpha + "_" + class_from + "_" + class_to, namespace);
+ codon_rate = parameters.ApplyNameSpace(between_rate, namespace);
+ }
+ (_GenerateRate.p[model_type])[^"terms.model.MSS.between"] = codon_rate;
+
+ } else {
+ if (class_from + class_to == nr) {
+ //console.log ("NEUTRAL");
+ codon_rate = 1;
+ } else {
+ if (model_type == utility.getGlobalValue("terms.local")) {
+ codon_rate = alpha + "_" + class_from + "_" + class_to;
+ } else {
+ codon_rate = parameters.ApplyNameSpace(alpha + "_" + class_from + "_" + class_to, namespace);
+ }
+ (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_between" + class_from + " and " + class_to] = codon_rate;
}
- (_GenerateRate.p[model_type])[alpha_term + ^"terms.model.MSS.syn_rate_between" + class_from + " and " + class_to] = codon_rate;
}
+ rate_entry += "*" + codon_rate;
}
- rate_entry += "*" + codon_rate;
}
}
@@ -363,6 +432,27 @@ lfunction models.codon.MSS._GenerateRate (fromChar, toChar, namespace, model_typ
return _GenerateRate.p;
}
+
+//----------------------------------------------------------------------------------------------------------------
+
+lfunction models.codon.MSS.LoadEmpiricalRates (file) {
+
+ SetDialogPrompt ("A TSV file with three columns (Codon1, Codon2, Empirical Rate) which is used to define relative synonymous substitution rates");
+ classes = io.ReadDelimitedFile (file, "\t", TRUE);
+ headers = utility.Array1D(classes[^'terms.io.header']);
+ io.CheckAssertion("`&headers`>=3", "Expected a TSV file with at least 3 columns; 2nd column is the codon, 3rd is the class for this codon");
+ codon_pairs = {};
+ for (_record_; in; classes [^"terms.io.rows"]) {
+ if (_record_[0] < _record_[1]) {
+ key = _record_[0] + "|" + _record_[1];
+ } else {
+ key = _record_[1] + "|" + _record_[0];
+ }
+ codon_pairs [key] = Eval (_record_[2]);
+ }
+
+ return {^"terms.model.MSS.codon_classes" : codon_pairs, ^"terms.model.MSS.empirical" : TRUE};
+}
//----------------------------------------------------------------------------------------------------------------
@@ -391,3 +481,37 @@ lfunction models.codon.MSS.LoadClasses (file) {
return {^"terms.model.MSS.codon_classes" : codons_by_class, ^"terms.model.MSS.neutral" : nr};
}
+
+//----------------------------------------------------------------------------------------------------------------
+
+lfunction models.codon.MSS.LoadClassesCodon (file) {
+
+ SetDialogPrompt ("A TSV file with three columns (Codon1, Codon2, Class) which is used to partition synonymous substitutions into groups");
+ classes = io.ReadDelimitedFile (file, "\t", TRUE);
+ headers = utility.Array1D(classes[^'terms.io.header']);
+ io.CheckAssertion("`&headers`>=3", "Expected a TSV file with at least 3 columns; 2nd column is the codon, 3rd is the class for this codon");
+ codon_pairs = {};
+ for (_record_; in; classes [^"terms.io.rows"]) {
+ if (_record_[0] < _record_[1]) {
+ key = _record_[0] + "|" + _record_[1];
+ } else {
+ key = _record_[1] + "|" + _record_[0];
+ }
+ codon_pairs [key] = _record_[2];
+ }
+
+ classes = utility.UniqueValues(codon_pairs);
+ class_count = utility.Array1D(classes);
+ io.CheckAssertion("`&class_count`>=2", "Expected at least 2 codon classes");
+
+ choices = {class_count,2};
+ for (i = 0; i < class_count; i += 1) {
+ choices[i][0] = classes[i];
+ choices[i][1] = "Codon class " + classes[i];
+ }
+
+ nr= io.SelectAnOption (choices, "Select the codon class which will serve as the neutral rate reference (relative rate = 1)");
+
+
+ return {^"terms.model.MSS.codon_classes" : codon_pairs, ^"terms.model.MSS.neutral" : nr, ^"terms.model.MSS.codon_pairs" : TRUE};
+}
=====================================
res/TemplateBatchFiles/libv3/models/frequencies.bf
=====================================
@@ -230,9 +230,34 @@ lfunction frequencies._aux.validate (model) {
}
+lfunction frequencies._aux.pad_zeros (freqs) {
+ nr = Rows (freqs);
+ nc = Columns (freqs);
+ F = freqs;
+
+ for (c = 0; c < nc; c += 1) {
+ s = 0;
+ for (r = 0; r < nr; r += 1) {
+ if (freqs [r][c] > 0) {
+ s += freqs [r][c];
+ F[r][c] = freqs [r][c];
+ } else {
+ s += 1e-8;
+ F[r][c] = 1e-8;
+ }
+ }
+ if (s > 1) {
+ for (r = 0; r < nr; r += 1) {
+ F[r][c] = F[r][c]/s;
+ }
+ }
+ }
+
+ return F;
+}
+
lfunction frequencies.empirical.codon_from_nuc (model, nuc_dict) {
-
result = {model.Dimension (model),1};
corrector = 1;
utility.ForEach (model[^"terms.stop_codons"], "codon",
@@ -267,6 +292,7 @@ lfunction frequencies.empirical.F3x4(model, namespace, datafilter) {
__f = model [^"terms.efv_estimate"];
}
+ __f = frequencies._aux.pad_zeros (__f );
__alphabet = model[^"terms.bases"];
nuc_dict = {};
@@ -303,6 +329,8 @@ lfunction frequencies.empirical.F1x4(model, namespace, datafilter) {
} else {
__f = model [^"terms.efv_estimate"];
}
+
+ __f = frequencies._aux.pad_zeros (__f );
__alphabet = model[^"terms.bases"];
nuc_dict = {};
@@ -338,6 +366,8 @@ lfunction frequencies.empirical.corrected.CF3x4(model, namespace, datafilter) {
__f = model [^"terms.efv_estimate"];
}
+ __f = frequencies._aux.pad_zeros (__f );
+
//TODO
__alphabet = model[^"terms.alphabet"];
__estimates = frequencies._aux.CF3x4(__f, model[^"terms.bases"], __alphabet, model[^"terms.stop_codons"]);
@@ -462,6 +492,7 @@ function frequencies._aux.empirical.singlechar(model, namespace, datafilter) {
} else {
__f = model [^"terms.efv_estimate"];
}
+ __f = frequencies._aux.pad_zeros (__f );
model[terms.efv_estimate] = __f;
return model;
=====================================
res/TemplateBatchFiles/libv3/models/parameters.bf
=====================================
@@ -853,18 +853,23 @@ lfunction parameters.helper.tree_lengths_to_initial_values(dict, type) {
for (i = 0; i < components; i += 1) {
this_component = {};
+ for (_branch_name_ , _branch_length_; in; (dict[i])[ utility.getGlobalValue("terms.branch_length")]) {
+ this_component [_branch_name_] = {
+ utility.getGlobalValue('terms.fit.MLE') :
+ factor*_branch_length_
+ };
+ }
-
- utility.ForEachPair((dict[i])[ utility.getGlobalValue("terms.branch_length")], "_branch_name_", "_branch_length_",
- "
- `&this_component`[_branch_name_] = {utility.getGlobalValue('terms.fit.MLE') : `&factor`*_branch_length_}
- ");
- result[i] = this_component;
+ if (Abs (this_component)) {
+ result[i] = this_component;
+ }
}
- return { utility.getGlobalValue("terms.branch_length"): result
+ return {
+ utility.getGlobalValue("terms.branch_length"): result
};
+
}
/**
=====================================
res/TemplateBatchFiles/libv3/tasks/estimators.bf
=====================================
@@ -48,10 +48,10 @@ lfunction estimators.RestoreLFStateFromSnapshot(lf_id, snapshot) {
lfunction estimators.ConstrainAndRunLRT (lf_id, constraint) {
savedMLES = estimators.TakeLFStateSnapshot (lf_id);
currentLL = estimators.ComputeLF (lf_id);
-
+
df = Call (constraint, TRUE);
Optimize (res, ^lf_id);
-
+
if (df >= 0) {
lrt = math.DoLRT (res[1][0],currentLL,df);
} else {
=====================================
res/TemplateBatchFiles/libv3/tasks/mapping.bf
=====================================
@@ -0,0 +1,519 @@
+LoadFunctionLibrary("../IOFunctions.bf");
+LoadFunctionLibrary("../all-terms.bf");
+
+
+/**
+ * Define scoring parameters for alignment tasks
+ * @name mapping.DefineScoringMatrices
+ * @param code - genetic code definition
+ * @returns {Dictionary} a structure with various scoring options (opaque)
+ */
+
+namespace mapping {
+
+ lfunction DefineScoringMatrices (code) {
+
+ score_matrix = {
+ { 6, -3, -4, -4, -2, -2, -2, -1, -3, -3, -3, -2, -2, -4, -2, 0, -1, -5, -3, -1, -4, -2, -2, -7}
+ { -3, 8, -2, -4, -6, 0, -2, -5, -2, -6, -4, 1, -3, -5, -4, -2, -3, -5, -3, -5, -3, -1, -2, -7}
+ { -4, -2, 8, 0, -5, -1, -2, -2, 0, -6, -6, -1, -4, -5, -4, 0, -1, -7, -4, -5, 6, -2, -2, -7}
+ { -4, -4, 0, 8, -6, -2, 0, -3, -3, -5, -6, -2, -6, -6, -3, -1, -3, -7, -6, -6, 6, 0, -3, -7}
+ { -2, -6, -5, -6, 10, -5, -7, -5, -5, -3, -3, -6, -3, -5, -5, -2, -2, -4, -4, -2, -5, -6, -4, -7}
+ { -2, 0, -1, -2, -5, 8, 1, -4, 0, -6, -4, 0, -1, -6, -3, -1, -2, -3, -3, -4, -1, 6, -2, -7}
+ { -2, -2, -2, 0, -7, 1, 7, -4, -1, -6, -5, 0, -4, -6, -3, -1, -2, -5, -4, -4, 0, 6, -2, -7}
+ { -1, -5, -2, -3, -5, -4, -4, 7, -4, -7, -6, -3, -5, -5, -4, -2, -4, -4, -5, -6, -2, -4, -4, -7}
+ { -3, -2, 0, -3, -5, 0, -1, -4, 10, -6, -5, -2, -3, -3, -4, -2, -4, -5, 0, -6, -1, -1, -3, -7}
+ { -3, -6, -6, -5, -3, -6, -6, -7, -6, 6, 0, -5, 0, -1, -5, -5, -2, -5, -3, 2, -5, -6, -2, -7}
+ { -3, -4, -6, -6, -3, -4, -5, -6, -5, 0, 6, -5, 1, -1, -5, -5, -3, -3, -3, 0, -6, -5, -2, -7}
+ { -2, 1, -1, -2, -6, 0, 0, -3, -2, -5, -5, 7, -3, -6, -2, -1, -2, -5, -3, -4, -2, 0, -2, -7}
+ { -2, -3, -4, -6, -3, -1, -4, -5, -3, 0, 1, -3, 9, -1, -5, -3, -2, -3, -3, 0, -5, -2, -1, -7}
+ { -4, -5, -5, -6, -5, -6, -6, -5, -3, -1, -1, -6, -1, 8, -6, -4, -4, 0, 1, -3, -6, -6, -3, -7}
+ { -2, -4, -4, -3, -5, -3, -3, -4, -4, -5, -5, -2, -5, -6, 9, -2, -3, -6, -5, -4, -4, -3, -4, -7}
+ { 0, -2, 0, -1, -2, -1, -1, -2, -2, -5, -5, -1, -3, -4, -2, 7, 0, -5, -3, -4, -1, -1, -2, -7}
+ { -1, -3, -1, -3, -2, -2, -2, -4, -4, -2, -3, -2, -2, -4, -3, 0, 7, -4, -3, -1, -2, -2, -2, -7}
+ { -5, -5, -7, -7, -4, -3, -5, -4, -5, -5, -3, -5, -3, 0, -6, -5, -4, 12, 0, -6, -7, -4, -4, -7}
+ { -3, -3, -4, -6, -4, -3, -4, -5, 0, -3, -3, -3, -3, 1, -5, -3, -3, 0, 9, -3, -5, -3, -2, -7}
+ { -1, -5, -5, -6, -2, -4, -4, -6, -6, 2, 0, -4, 0, -3, -4, -4, -1, -6, -3, 6, -6, -4, -2, -7}
+ { -4, -3, 6, 6, -5, -1, 0, -2, -1, -5, -6, -2, -5, -6, -4, -1, -2, -7, -5, -6, 7, -1, -3, -7}
+ { -2, -1, -2, 0, -6, 6, 6, -4, -1, -6, -5, 0, -2, -6, -3, -1, -2, -4, -3, -4, -1, 7, -2, -7}
+ { -2, -2, -2, -3, -4, -2, -2, -4, -3, -2, -2, -2, -1, -3, -4, -2, -2, -4, -2, -2, -3, -2, -2, -7}
+ { -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, -7, 1}
+ };
+
+ max_score = Max (score_matrix,0);
+ penalty = Max (max_score, -Min (score_matrix,0));
+
+
+ options = {
+ "SEQ_ALIGN_CHARACTER_MAP" : "ACGTN",
+ "SEQ_ALIGN_GAP_OPEN" : 1.5*penalty,
+ "SEQ_ALIGN_GAP_OPEN2" : 1.5*penalty,
+ "SEQ_ALIGN_GAP_EXTEND" : 0.1*penalty,
+ "SEQ_ALIGN_GAP_EXTEND2" : 0.1*penalty,
+ "SEQ_ALIGN_FRAMESHIFT" : 2*penalty,
+ "SEQ_ALIGN_NO_TP" : 1,
+ "SEQ_ALIGN_AFFINE" : 1,
+ "SEQ_ALIGN_CODON_ALIGN" : 1
+ };
+
+ base_frequencies = {
+ { 0.069352301}
+ { 0.021000333}
+ { 0.049862283}
+ { 0.026563029}
+ { 0.033026253}
+ { 0.1054545}
+ { 0.008970554}
+ { 0.036648952}
+ { 0.036329907}
+ { 0.065164842}
+ { 0.021805663}
+ { 0.032992805}
+ { 0.0348093}
+ { 0.036818766}
+ { 0.054168098}
+ { 0.12149678}
+ { 0.082464001}
+ { 0.053564744}
+ { 0.038383113}
+ { 0.07112377}
+ };
+
+
+
+ options["SEQ_ALIGN_SCORE_MATRIX"] = pSM2cSM(score_matrix, "ACDEFGHIKLMNPQRSTVWY", code["code"], code["ordering"]);
+ shift_penalty = computeExpectedPerBaseScore (.5,score_matrix,base_frequencies);
+
+ _cdnaln_partialScoreMatrices = cSM2partialSMs(options["SEQ_ALIGN_SCORE_MATRIX"],
+ {{shift_penalty__*1.5,shift_penalty__,shift_penalty__,shift_penalty*1.5}});
+
+
+ options ["SEQ_ALIGN_PARTIAL_3x1_SCORES"] = _cdnaln_partialScoreMatrices["3x1"];
+ options ["SEQ_ALIGN_PARTIAL_3x2_SCORES"] = _cdnaln_partialScoreMatrices["3x2"];
+ options ["SEQ_ALIGN_PARTIAL_3x4_SCORES"] = _cdnaln_partialScoreMatrices["3x4"];
+ options ["SEQ_ALIGN_PARTIAL_3x5_SCORES"] = _cdnaln_partialScoreMatrices["3x5"];
+
+ options ["MATCH"] = computeExpectedPerBaseScore (1,score_matrix,base_frequencies);
+ options ["E"] = Max (0.1, computeExpectedPerBaseScore (.4,score_matrix,base_frequencies));
+
+ return options;
+ }
+
+
+ lfunction cSM2partialSMs(_cdnScoreMatrix, penalties) {
+
+ m3x5 = { 126, 1250 };
+ m3x4 = { 126, 500 };
+ m3x2 = { 126, 75 };
+ m3x1 = { 126, 15 };
+
+
+ if (utility.Array1D (penalties) == 4) {
+ p3x1 = penalties [0];
+ p3x2 = penalties [1];
+ p3x4 = penalties [2];
+ p3x5 = penalties [3];
+ } else {
+ p3x5 = 0;
+ p3x4 = 0;
+ p3x2 = 0;
+ p3x1 = 0;
+ }
+
+ for ( thisCodon = 0; thisCodon < 64; thisCodon += 1 ) {
+ for ( d1 = 0; d1 < 5; d1 += 1 ) {
+ max100 = -1e100;
+ max010 = -1e100;
+ max001 = -1e100;
+
+ for ( d2 = 0; d2 < 5; d2 += 1 ) {
+ partialCodon = 5 * d1 + d2;
+ max110 = -1e100;
+ max101 = -1e100;
+ max011 = -1e100;
+
+ for ( d3 = 0; d3 < 5; d3 += 1 ) {
+ thisCodon2 = 5 * partialCodon + d3;
+ thisScore = _cdnScoreMatrix[ thisCodon ][ thisCodon2 ];
+
+ // this is the trivial and stupid way of doing it, but it should work
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 0 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 1 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 2 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 3 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 4 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 5 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 6 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 7 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 8 ] = thisScore - p3x5;
+ m3x5[ thisCodon ][ 10 * thisCodon2 + 9 ] = thisScore - p3x5;
+
+ m3x4[ thisCodon ][ 4 * thisCodon2 + 0 ] = thisScore - p3x4;
+ m3x4[ thisCodon ][ 4 * thisCodon2 + 1 ] = thisScore - p3x4;
+ m3x4[ thisCodon ][ 4 * thisCodon2 + 2 ] = thisScore - p3x4;
+ m3x4[ thisCodon ][ 4 * thisCodon2 + 3 ] = thisScore - p3x4;
+
+ // d1 is 1
+ max100 = Max( max100, _cdnScoreMatrix[ thisCodon ][ 25 * d1 + 5 * d2 + d3 ] );
+ max010 = Max( max010, _cdnScoreMatrix[ thisCodon ][ 25 * d2 + 5 * d1 + d3 ] );
+ max001 = Max( max001, _cdnScoreMatrix[ thisCodon ][ 25 * d2 + 5 * d3 + d1 ] );
+
+ // d1 and d2 are 1
+ max110 = Max( max110, _cdnScoreMatrix[ thisCodon ][ 25 * d1 + 5 * d2 + d3 ] );
+ max101 = Max( max101, _cdnScoreMatrix[ thisCodon ][ 25 * d1 + 5 * d3 + d2 ] );
+ max011 = Max( max011, _cdnScoreMatrix[ thisCodon ][ 25 * d3 + 5 * d1 + d2 ] );
+ }
+
+ m3x2[ thisCodon ][ 3 * partialCodon + 0 ] = max110 - p3x2;
+ m3x2[ thisCodon ][ 3 * partialCodon + 1 ] = max101 - p3x2;
+ m3x2[ thisCodon ][ 3 * partialCodon + 2 ] = max011 - p3x2;
+ }
+
+ m3x1[ thisCodon ][ 3 * d1 + 0 ] = max100 - p3x1;
+ m3x1[ thisCodon ][ 3 * d1 + 1 ] = max010 - p3x1;
+ m3x1[ thisCodon ][ 3 * d1 + 2 ] = max001 - p3x1;
+ }
+ }
+
+
+ return { "3x1": m3x1, "3x2": m3x2, "3x4": m3x4, "3x5": m3x5 };
+ }
+
+
+
+ lfunction _digits (index) {
+ return {{index__$25, index__ % 25 $ 5, index__ % 5}};
+ };
+
+ lfunction _hazN (digits) {
+ return digits [0] == 4 || digits [1] == 4 || digits [2] == 4;
+ };
+
+ lfunction _map_to_nuc (digits) {
+ return digits [0] * 16 + digits [1] * 4 + digits[2];
+ };
+
+ lfunction _generate_resolutions (digits) {
+ resolutions = {};
+ for (k = 0; k < 4; k += 1) {
+ try = digits;
+ if (digits[0] == 4) {
+ try [0] = k;
+ }
+ for (k2 = 0; k2 < 4; k2 += 1) {
+ if (digits[1] == 4) {
+ try [1] = k;
+ }
+ for (k3 = 0; k3 < 4; k3 += 1) {
+ if (digits[2] == 4) {
+ try [2] = k;
+ }
+ resolutions[_map_to_nuc (try)] = 1;
+ }
+ }
+ }
+ return Rows(resolutions);
+ };
+
+ // -------------------------------------------------------------------------- //
+
+ lfunction pSM2cSM (_scorematrix, _letters, code, ordering) {
+
+ _cdnScoreMatrix = { 126 , 126 };
+ _mapping = utility.MapStrings ( ordering, _letters );
+
+
+ for ( _k = 0; _k < 125; _k += 1 ) {
+
+ letters1 = _digits (_k);
+
+ if (_hazN (letters1)) {
+ letters1 = _generate_resolutions (letters1);
+ for ( _k2 = _k; _k2 < 125 ; _k2 += 1 ) {
+ letters2 = _digits (_k2);
+ if (_hazN (letters2) == 0) {
+ codon2 = _map_to_nuc(letters2);
+ _mappedK2 = _mapping[ code[ codon2 ] ];
+ _aScore = -1e4;
+ if (_mappedK2 >= 0) {
+
+ res_count = utility.Array1D (letters1);
+ for (r = 0; r < res_count; r += 1) {
+ resolution_codon = 0 + letters1[r];
+ resolution_aa = _mapping[ code[ resolution_codon ] ];
+ if (resolution_aa >= 0) {
+ try_score = _scorematrix[ resolution_aa ][ _mappedK2 ] - 1;
+ if (resolution_aa == _mappedK2 && codon2 != resolution_codon) {
+ try_score = try_score - 1;
+ }
+ _aScore = Max (_aScore, try_score);
+ }
+ }
+ }
+
+
+ _cdnScoreMatrix[ _k ][ _k2 ] = _aScore;
+ _cdnScoreMatrix[ _k2 ][ _k ] = _aScore;
+ }
+ }
+ } else {
+ codon1 = _map_to_nuc(letters1);
+
+ _mappedK = _mapping[ code[ codon1 ] ];
+
+ if ( _mappedK >= 0) {
+ for ( _k2 = _k; _k2 < 125 ; _k2 += 1 ) {
+ letters2 = _digits (_k2);
+
+ if (_hazN (letters2)) {
+ _aScore = -1e4;
+ letters2 = _generate_resolutions (letters2);
+ res_count = utility.Array1D (letters2);
+ for (r = 0; r < res_count; r += 1) {
+ resolution_codon = 0 + letters2[r];
+ resolution_aa = _mapping[ code[ resolution_codon ] ];
+ if (resolution_aa >= 0) {
+ try_score = _scorematrix[ _mappedK ][ resolution_aa ] - 1;
+ if (resolution_aa == _mappedK && codon1 != resolution_codon) {
+ try_score = try_score - 1;
+ }
+ _aScore = Max (_aScore, try_score);
+ }
+ }
+ _cdnScoreMatrix[ _k ][ _k2 ] = _aScore;
+ _cdnScoreMatrix[ _k2 ][ _k ] = _aScore;
+ continue;
+ }
+
+ codon2 = _map_to_nuc(letters2);
+
+ _mappedK2 = _mapping[ code[ codon2 ] ];
+ if ( _mappedK2 >= 0 ) {
+ _aScore = _scorematrix[ _mappedK ][ _mappedK2 ];
+ if ( _mappedK == _mappedK2 && _k2 > _k ) {
+ _aScore = _aScore - 1; // synonymous match
+ }
+ } else {
+ // stop codons don't match anything
+ _aScore = -1e4;
+ }
+ _cdnScoreMatrix[ _k ][ _k2 ] = _aScore;
+ _cdnScoreMatrix[ _k2 ][ _k ] = _aScore;
+ }
+ } else { // stop codons here
+ for ( _k2 = _k; _k2 < 125; _k2 += 1 ) {
+
+ letters2 = _digits (_k2);
+
+ if (_hazN (letters2)) {
+ continue;
+ }
+
+ codon2 = _map_to_nuc(letters2);
+
+ _mappedK2 = _mapping[ code[ codon2 ] ];
+
+ if ( _mappedK2 < 0 ) {
+ // don't penalize stop codons matching themselves
+ _cdnScoreMatrix[ _k ][ _k2 ] = 0;
+ _cdnScoreMatrix[ _k2 ][ _k ] = 0;
+ } else {
+ _cdnScoreMatrix[ _k ][ _k2 ] = -1e4;
+ _cdnScoreMatrix[ _k2 ][ _k ] = -1e4;
+ }
+ }
+ }
+ }
+ }
+
+ return _cdnScoreMatrix;
+ }
+
+ // -------------------------------------------------------------------------- //
+
+ lfunction computeExpectedPerBaseScore( _expectedIdentity, _cdnaln_scorematrix, _cdnaln_base_freqs ) {
+ meanScore = 0;
+
+ for (_aa1 = 0; _aa1 < 20; _aa1 += 1) {
+ for (_aa2 = 0; _aa2 < 20; _aa2 += 1) {
+ if ( _aa1 != _aa2 ) {
+ meanScore += ( 1 - _expectedIdentity ) * _cdnaln_scorematrix[_aa1][_aa2] * _cdnaln_base_freqs[_aa1] * _cdnaln_base_freqs[_aa2];
+ } else {
+ meanScore += _expectedIdentity * _cdnaln_scorematrix[_aa1][_aa1] * _cdnaln_base_freqs[_aa1] * _cdnaln_base_freqs[_aa1];
+ }
+ }
+ }
+
+ return meanScore;
+ }
+
+ // -------------------------------------------------------------------------- //
+
+ lfunction computeCorrection (str) {
+ result = {1,2};
+
+ result[0] = (str$"^\\-+")[1]+1;
+ result[1] = (str$"\\-+$")[0];
+
+ if (result[1] >= 0) {
+ result[1] = Abs(str)-result[1];
+ }
+ else {
+ result[1] = 0;
+ }
+ return result;
+ }
+
+
+ // -------------------------------------------------------------------------- //
+
+ lfunction igg_alignment_cleanup (reference, query, offset_nuc, code) {
+
+ too_short = 0;
+ too_long = 0;
+ span = 0; // how many nucleotides in the reference were covered by non-gaps
+ _seqL = Abs (reference);
+
+
+ ref_cleaned = ""; ref_cleaned * 128;
+ qry_cleaned = ""; qry_cleaned * 128;
+
+ _codon_in_reference = 0;
+
+ for ( _rcidx = 0; _rcidx < _seqL; _rcidx += 1 ) {
+ _del1 = reference [_rcidx] != (reference [_rcidx]&&1);
+ if (_del1) {
+ too_short += 1;
+ _codon_in_reference += 1;
+ ref_cleaned * (reference [_rcidx]&&1);
+ qry_cleaned * (query [_rcidx]&&1);
+ } else {
+ _del1 = query [_rcidx] != (query [_rcidx]&&1);
+ if (_del1) {
+ if (_seqL-_rcidx < 3 && _codon_in_reference % 3 == 0) {
+ break;
+ }
+ too_long += 1;
+ } else {
+ ref_cleaned * (reference [_rcidx]&&1);
+ qry_cleaned * (query [_rcidx]&&1);
+ span += 1;
+ _codon_in_reference +=1;
+ }
+ }
+ }
+ ref_cleaned * 0; qry_cleaned * 0;
+
+ return {"REF": ref_cleaned, "QRY": qry_cleaned, "TOO_SHORT" : too_short, "TOO_LONG": too_long, "SPAN": span, "OFFSET_AA" : offset_nuc$3 + (offset_nuc % 3 > 0),"OFFSET" : offset_nuc, "AA" : alignments.TranslateCodonsToAminoAcids (qry_cleaned, (3-offset_nuc%3)%3, code), "AA_REF" : alignments.TranslateCodonsToAminoAcids (ref_cleaned, (3-offset_nuc%3)%3, code)};
+ }
+
+ // -------------------------------------------------------------------------- //
+
+ lfunction correctReadUsingCodonAlignedData (ref, qry, code) {
+ reference_shifts = computeCorrection(ref);
+
+ /*reference_shifts is the starting,ending nucleotide on the reference relative to the read. if reference is longer than the read, then both are 0*/
+
+ offsetFrom = (qry$"^\\-+")[1]+1;
+ offsetTo = (qry$"\\-+$")[0]-1;
+
+ /* the $ looks for the regular expression in bestAl[2] and returns a 2x1 array with the starting and ending 0-based positions of the regular expression. in this case multiple indels, -. returns -1 for both if the regular expression is not found.
+ i.e. 0-based index leading indels start at (bestAl[2]$"^\\-+")[0] and end at (bestAl[2]$"^\\-+")[1]; trailing indels start at (bestAl[2]$"\\-+$")[0] and end at (bestAl[2]$"\\-+$")[0];
+
+ so offSetFrom to offSetTo will return the reference sequence co-ordinates overlapping with the read.
+ */
+
+
+ if (offsetTo < 0) {
+ offsetTo = Abs(qry)-1; /*if no trailing indels then to end of read*/
+ }
+
+ // check to see if the prefix in REF has some out-of-frame indels so that we can start offsetFrom in the correct frame */
+
+ if (offsetFrom > 0) {
+ frame_skips = 0;
+ for (i = 0; i < offsetFrom; i += 1) {
+ if ((ref[i] && 1) != ref[i]) {
+ frame_skips += 1;
+ }
+ }
+ offsetFrom += -frame_skips;
+ }
+
+ seqOffset = offsetFrom; /*set the offset of the read relative to the reference. ie the number of indels needed on the read to align to the reference */
+ offsetFrom += reference_shifts[0]; /*if the read starts before the reference then shift to start of reference ie. by reference_shifts[0] */
+ offsetTo = offsetTo - reference_shifts[1]; /*if the read extends beyond the reference then shift to end of reference ie. by reference_shifts[1] */
+
+ theSeq = qry;
+ theSeq = qry[reference_shifts[0]][Abs(qry)-reference_shifts[1]-1]; /*the nucleotide sequence of the read that overlaps with the reference sequence */
+
+ nucSeq = qry[offsetFrom][offsetTo]; /*read sequence pruned to exactly overlapping region*/
+ nucSeqRef = ref[offsetFrom][offsetTo]; /*reference sequence pruned to exactly overlapping region*/
+
+ extra_keys = {};
+
+ if (reference_shifts[0] > 0) { // 'qry' has a prefix that's not aligned to the reference
+ extra_keys['PREFIX'] = qry[0][reference_shifts[0]-1];
+ }
+ if (reference_shifts[1] > 0) {
+ l = Abs (qry);
+ extra_keys['SUFFIX'] = qry[l-reference_shifts[1]][l-1];
+ }
+
+
+ return utility.Extend (igg_alignment_cleanup (nucSeqRef, nucSeq,seqOffset, code), extra_keys);
+ }
+
+
+ // -------------------------------------------------------------------------- //
+
+
+ lfunction align_sequence_to_reference_set (seq, references, alignment_settings, seq_name) {
+
+
+ input = {1,2};
+ input [1] = seq;
+
+ overall = {'SCORE' : -1e100, 'REFNAME' : ''};
+
+ for (id, seq; in; references) {
+ input [0] = seq;
+ AlignSequences (result, input, alignment_settings);
+
+ result = result[0];
+
+ if (result [0] >= overall['SCORE']) {
+
+ if (alignment_settings ["REPORT_VARIANTS"] && result [0] == overall['SCORE'] && Abs (overall['REFNAME'])) {
+ overall['REFNAME'] += "|" + id;
+ } else {
+ overall['REFNAME'] = id;
+ }
+ overall['SCORE'] = result[0];
+ overall['RAW-REF'] = result[1];
+ overall['RAW-QRY'] = result[2];
+ }
+ }
+
+
+ if (Type (overall['RAW-REF']) == "String") {
+
+ computed_score = (overall["SCORE"] - 30 * alignment_settings["MATCH"] * Exp (-Abs(seq)/3) ) / Abs (seq) * 3 ;
+
+ if (alignment_settings["E"] <= computed_score) {
+ utility.Extend (overall, correctReadUsingCodonAlignedData (overall['RAW-REF'], overall['RAW-QRY'], alignment_settings["code"]));
+
+
+ if (alignment_settings["SEQ_ALIGN_CODON_ALIGN"] == TRUE && overall["SPAN"] <= 3) {
+ //assert (0, "Internal error in align_sequence_to_reference_set" + overall + "\nComputed score `computed_score`; expected score " + alignment_settings["E"] + "; match score " + alignment_settings["MATCH"] + "\nInput sequence: `seq`");
+ return None;
+ }
+ return overall;
+ }
+ }
+
+
+ return None;
+ }
+}
+
=====================================
res/TemplateBatchFiles/libv3/tasks/trees.bf
=====================================
@@ -436,6 +436,28 @@ lfunction trees.extract_paml_annotation (tree_string) {
return result;
}
+/**
+ * @name trees.AdjustZerosToNearZeros
+ * Given a tree dict and a list of matching parameter estimates
+ * this traverses the estimates and converts terminal branch lengths that are exactly 0 to something small (1e-6)
+ * internal branches removed and modifies in in place
+ * @param tree - dict of the tree object
+ * @param estimates - dict with branch length estimates
+ * @return nothing; estimates modified in place
+ */
+
+lfunction trees.AdjustZerosToNearZeros (tree, estimates) {
+
+ for (branch, value; in; tree [^"terms.trees.partitioned"]) {
+ if (value != ^"terms.tree_attributes.internal") {
+ if (estimates / branch) {
+ if ((estimates[branch])[^"terms.fit.MLE"] < 1e-10) {
+ (estimates[branch])[^"terms.fit.MLE"] = 1e-6;
+ }
+ }
+ }
+ }
+}
/**
@@ -722,9 +744,10 @@ lfunction trees.ParsimonyLabel(tree_id, given_labels) {
scores = {}; // node name -> score of optimal labeling staring at this now given parent state
optimal_labeling = {}; // node name -> current node labeling which achieves this score
resulting_labels = {}; // internal nodes -> label
-
+
// pass 1 to fill in the score matrix
for (k = 0; k < Abs (tree_avl) ; k += 1) {
+ //console.log (node_name);
node_name = (tree_avl[k])["Name"];
node_children = (tree_avl[k])["Children"];
c_count = Abs (node_children);
@@ -788,7 +811,7 @@ lfunction trees.ParsimonyLabel(tree_id, given_labels) {
}
- // pass 2 to choose the best state for subtree parents
+ // pass 2 to choose the best state for subtree parents
total_score = 0;
@@ -826,6 +849,8 @@ lfunction trees.ParsimonyLabel(tree_id, given_labels) {
}
}
+ console.log (labels);
+
tree_avl = (^tree_id) ^ 1;
for (k = 2; k < Abs (tree_avl); k += 1) {
node_name = (tree_avl[k])["Name"];
=====================================
src/core/fisher_exact.cpp
=====================================
@@ -18,6 +18,7 @@ Function prototype is declared in matrix.h
*/
+#include <stdlib.h>
#include <math.h>
#include "global_things.h"
=====================================
src/core/global_things.cpp
=====================================
@@ -93,6 +93,7 @@ namespace hy_global {
has_terminal_stdout = true,
has_terminal_stderr = true,
ignore_kw_defaults = false,
+ exiting_with_error = false,
force_verbosity_from_cli = false;
hyFile *hy_error_log_file = NULL,
@@ -122,7 +123,7 @@ namespace hy_global {
kErrorStringMatrixExportError ("Export matrix called with a non-polynomial matrix argument"),
kErrorStringNullOperand ("Attempting to operate on an undefined value; this is probably the result of an earlier 'soft' error condition"),
kErrorNumerical ("To treat numerical errors as warnings, please specify \"ENV=TOLERATE_NUMERICAL_ERRORS=1;\" as the command line argument. This often resolves the issue, which is indicative of numerical instability."),
- kHyPhyVersion = _String ("2.5.63"),
+ kHyPhyVersion = _String ("2.5.68"),
kNoneToken = "None",
kNullToken = "null",
@@ -745,6 +746,7 @@ namespace hy_global {
//____________________________________________________________________________________
void HandleApplicationError (const _String & message, bool force_exit, bool dump_core) {
+
if (!force_exit && currentExecutionList && currentExecutionList->errorHandlingMode == HY_BL_ERROR_HANDLING_SOFT) {
currentExecutionList->ReportAnExecutionError(message, true);
return;
@@ -785,9 +787,14 @@ namespace hy_global {
#else
errMsg = message;
#endif
+ if (exiting_with_error) {
+ errMsg = errMsg & "\nAn additional error occurred while handling the first error. Error messages may be incomplete\n";
+ } else {
+ exiting_with_error = true;
errMsg = ConstructAnErrorMessage (errMsg);
- StringToConsole(errMsg);
- #endif
+ }
+ StringToConsole(errMsg);
+ #endif
#if defined __UNIX__
if (hy_drop_into_debug_mode)
=====================================
src/core/likefunc.cpp
=====================================
@@ -3105,7 +3105,7 @@ inline hyFloat sqr (hyFloat x)
_Matrix * eq_freqs = GetIthFrequencies(partition_index);
unsigned long tip_count = filter->NumberSpecies();
- if (filter->GetData()->GetTT()->IsStandardNucleotide() && filter->IsNormalFilter() && tip_count<150 && eq_freqs->IsIndependent()) {
+ if (filter->GetData()->GetTT()->IsStandardNucleotide() && filter->IsNormalFilter() && tip_count < 300 && eq_freqs->IsIndependent()) {
// if not - use distance estimates
if (tree->IsDegenerate()) {
@@ -4341,15 +4341,6 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
_variables_changed_during_last_compute = new _SimpleList ();
variables_changed_during_last_compute = new _AVLList (_variables_changed_during_last_compute);
-#ifdef __HYPHYMPI__
- if (hy_mpi_node_rank == 0) {
-#endif
- BenchmarkThreads (this);
-#ifdef __HYPHYMPI__
- }
-
-
-#endif
//ObjectToConsole(variables_changed_during_last_compute);
//NLToConsole();
@@ -4489,7 +4480,14 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
}
-
+#ifdef __HYPHYMPI__
+ if (hy_mpi_node_rank == 0) {
+#endif
+ BenchmarkThreads (this);
+#ifdef __HYPHYMPI__
+ }
+#endif
+
_OptimiztionProgress progress_tracker;
//checkParameter (optimizationMethod,optMethodP,4.0);
@@ -4604,11 +4602,12 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
hyFloat grad_precision;
if (maxSoFar > -INFINITY) {
- grad_precision = MIN (1000., -maxSoFar / 500.);
+ grad_precision = MIN (1000., -maxSoFar * 0.002);
} else {
- grad_precision = -0.005;
+ grad_precision = -0.0025;
}
+
if (gradientBlocks.nonempty()) {
for (long b = 0; b < gradientBlocks.lLength; b++) {
_SimpleList * this_block = (_SimpleList*)gradientBlocks(b);
@@ -4930,7 +4929,7 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
_Matrix bestMSoFar;
GetAllIndependent (bestMSoFar);
hyFloat prec = Minimum (diffs[0], diffs[1]),
- grad_precision = Maximum(diffs[0],diffs[1]);
+ grad_precision = Maximum (precision, Maximum(diffs[0],diffs[1]));
prec = Minimum (Maximum (prec, precision), 1.);
@@ -5119,9 +5118,9 @@ _Matrix* _LikelihoodFunction::Optimize (_AssociativeList const * options)
//verbosity_level = 101;
if (use_adaptive_step) {
if (convergenceMode < 2) {
- LocateTheBump (current_index,precisionStep, maxSoFar, bestVal, go2Bound, precisionStep);
+ LocateTheBump (current_index,precisionStep, maxSoFar, bestVal, go2Bound, precisionStep * 0.5);
} else {
- LocateTheBump (current_index,precisionStep, maxSoFar, bestVal, go2Bound, precisionStep);//convergenceMode == 2? precisionStep*0.25: precisionStep*0.0625);
+ LocateTheBump (current_index,precisionStep, maxSoFar, bestVal, go2Bound, convergenceMode == 2? precisionStep*0.25: precisionStep*0.0625);
}
} else {
LocateTheBump (current_index,brackStep, maxSoFar, bestVal, go2Bound);
@@ -6516,7 +6515,7 @@ void _LikelihoodFunction::ComputeGradient (_Matrix& gradient, hyFloat& gradi
hyFloat currentValue = GetIthIndependent(index),
ub = GetIthIndependentBound(index,false)-currentValue,
lb = currentValue-GetIthIndependentBound(index,true),
- testStep = MAX(currentValue * gradientStep,gradientStep);
+ testStep = MAX(fabs(currentValue) * gradientStep,gradientStep);
//printf ("%ld %s %20.18g\n", index, GetIthIndependentName (index)->get_str(), currentValue);
@@ -6691,6 +6690,10 @@ hyFloat _LikelihoodFunction::ConjugateGradientDescent (hyFloat step_precision
//printf ("\n\n_LikelihoodFunction::ConjugateGradientDescent ==> %d (%lg)\n", usedCachedResults, maxSoFar);
+ if (maxSoFar == -INFINITY) {
+ return maxSoFar;
+ }
+
if (min_improvement_to_contuinue < 0.) {
min_improvement_to_contuinue = initial_value * min_improvement_to_contuinue;
}
=====================================
src/core/likefunc2.cpp
=====================================
@@ -1254,8 +1254,11 @@ hyFloat mapParameterToInverval (hyFloat in, char type, bool inverse) {
break;
case _hyphyIntervalMapSqueeze:
if (inverse) {
- return in/(1.-in);
+ in = in/(1.-in);
+ //return in;
+ return in * in;
} else {
+ in = sqrt (in);
return in/(1.+in);
}
break;
=====================================
src/core/matrix.cpp
=====================================
@@ -72,8 +72,8 @@ _String MATRIX_AGREEMENT = "CONVERT_TO_POLYNOMIALS",
int _Matrix::precisionArg = 0;
-int _Matrix::storageIncrement = 16;
-int _Matrix::switchThreshold = 35;
+int _Matrix::storageIncrement = 20;
+int _Matrix::switchThreshold = 30;
hyFloat _Matrix::truncPrecision = 1e-16;
#define MatrixMemAllocate(X) MemAllocate(X, false, 64)
=====================================
src/core/tree_evaluator.cpp
=====================================
@@ -3519,8 +3519,12 @@ hyFloat _TheTree::ComputeTreeBlockByBranch ( _SimpleList
result = -INFINITY;
#pragma omp critical
{
- //printf ("BAILING WITH INFINITY %ld / %ld (%ld)\n", siteID, siteOrdering.list_data[siteID], siteOrdering.lLength);
+ //printf ("BAILING WITH INFINITY %ld / %ld (%ld); (%ld);\n", siteID, siteOrdering.list_data[siteID], siteOrdering.lLength, theFilter->theMap .get(siteOrdering.list_data[siteID]*3));
hy_global::ReportWarning (_String("Site ") & (1L+siteOrdering.list_data[siteID]) & " evaluated to a 0 probability in ComputeTreeBlockByBranch");
+ /*long lfID = -1;
+ IsLinkedToALF (lfID);
+ ((_LikelihoodFunction*)likeFuncList (lfID))
+ ->_TerminateAndDump (_String("Site ") & (1L+siteOrdering.list_data[siteID]) & " evaluated to a 0 probability in ComputeTreeBlockByBranch", true);*/
}
break;
}
=====================================
tests/hbltests/libv3/FEL.wbf
=====================================
@@ -3,7 +3,7 @@ GetString (version, HYPHY_VERSION, 0);
if (+version >= 2.4) {
LoadFunctionLibrary ("SelectionAnalyses/FEL.bf",
{"--code" : "Universal", "--alignment" : PATH_TO_CURRENT_BF + "data/CD2.nex"/*,
- "--branches" : "All",
+ "--pvalue" : "0.2",
"--limit-to-sites" : "144", "--save-lf-for-sites" : "144"*/
});
} else {
=====================================
tests/hbltests/libv3/MEME.wbf
=====================================
@@ -1,8 +1,12 @@
GetString (version, HYPHY_VERSION, 0);
+p_cut = 0.05;
if (+version >= 2.4) {
- LoadFunctionLibrary ("SelectionAnalyses/MEME.bf", {"--code" : "Universal", "--alignment" : PATH_TO_CURRENT_BF + "data/CD2.nex", "--branches" : "All"});
+ LoadFunctionLibrary ("SelectionAnalyses/MEME.bf", {"--code" : "Universal",
+ "--alignment" : PATH_TO_CURRENT_BF + "data/CD2.nex",
+ "--branches" : "All",
+ "--pvalue" : "" + p_cut});
} else {
LoadFunctionLibrary ("SelectionAnalyses/MEME.bf", {"0" : "Universal", "1" : PATH_TO_CURRENT_BF + "data/CD2.nex", "2" : "All", "3" : "0.1"});
@@ -44,16 +48,27 @@ branches = (((meme.json["MLE"])["content"])["0"])[-1][7];
| 149 | 1 | 0.000 | 1.968 | 1.000 | 5.229 | Yes, p = 0.0335 | 0 |
| 174 | 1 | 0.737 | 25.356 | 0.176 | 3.851 | Yes, p = 0.0684 | 1 |
+
+| Codon | Partition | alpha |non-syn rate (beta) distribution, rates : weights| LRT |Episodic selection detected?| # branches | List of most common codon substitutions at this site |
+|:----------:|:----------:|:----------:|:-----------------------------------------------:|:----------:|:--------------------------:|:----------:|:---------------------------------------------------------------------:|
+| 43 | 1 | 0.000 | 0.00/15.54 : 0.83/0.17 | 4.878 | Yes, p = 0.0402 | 1 | [1]CTg>AAg |
+| 67 | 1 | 0.000 | 0.00/99.81 : 0.79/0.21 | 5.675 | Yes, p = 0.0266 | 1 | [1]TcC>AcT,TCc>GAc |
+| 76 | 1 | 0.000 | 0.00/48.40 : 0.70/0.30 | 10.326 | Yes, p = 0.0025 | 2 | [1]aCc>aTc,aGA>aCC,AGa>CTa |
+| 98 | 1 | 0.002 | 0.00/26.81 : 0.55/0.45 | 4.573 | Yes, p = 0.0470 | 0 | [1]aCa>aAa,aCa>aTa,ACA>GAC,Aca>Gca |
+| 113 | 1 | 0.002 | 0.00/8.41 : 0.34/0.66 | 5.871 | Yes, p = 0.0241 | 2 | [1]Tac>Cac,tAc>tCc,TCc>AAc,TCc>GAc |
+| 117 | 1 | 0.000 | 0.00/4.81 : 0.00/1.00 | 5.246 | Yes, p = 0.0332 | 0 | [2]aCc>aTc|[1]aCc>aGc,Acc>Ccc,Acc>Gcc |
+| 149 | 1 | 0.000 | 0.00/4.71 : 0.00/1.00 | 5.229 | Yes, p = 0.0335 | 0 | [1]aCc>aAc,ACc>GTc,Gtc>Ctc,Gtc>Ttc |
+
*/
test.lrt_sum = 0;
test.branch_sum = 0;
-test.expected_positives = utility.MatrixToDict({{9,34,43,55,64,67,76,81,98,113,117,122,141,149,174}});
+test.expected_positives = utility.MatrixToDict({{43,67,76,98,113,117,149}});
function confirm_site (site, p, dict) {
- if (p <= 0.1) {
+ if (p <= p_cut) {
if (dict/(site+1)) {
test.lrt_sum += lrts[site];
test.branch_sum += branches[site];
@@ -75,7 +90,7 @@ utility.ForEachPair (p_values,"_index_", "_p_",
");
assert (check_value (
- test.lrt_sum, 72.175, 0.05), "More than 5% difference in cumulative LRT for positively selected sites");
+ test.lrt_sum, 41.798, 0.05), "More than 5% difference in cumulative LRT for positively selected sites");
assert (check_value (
- test.branch_sum, 11, 0.0001), "Incorrect total # of branches with high EBF");
+ test.branch_sum, 6, 0.0001), "Incorrect total # of branches with high EBF");
=====================================
tests/hbltests/libv3/algae-mtDNA.wbf
=====================================
@@ -14,7 +14,7 @@ assert (check_value (
(fitter.results)["LogL"], -2345.80, 0.05), "Incorrect log-likelihood for the codon model");
assert (check_value (
- (((fitter.results)["global"])["non-synonymous/synonymous rate ratio"])["MLE"], 0.0040, 0.01), "Incorrect dN/dS ratio");
+ (((fitter.results)["global"])["non-synonymous/synonymous rate ratio"])["MLE"], 0.0040, 0.05), "Incorrect dN/dS ratio");
bl_sum = 0;
for (bl; in; (fitter.results["branch length"])["0"]){
View it on GitLab: https://salsa.debian.org/med-team/hyphy/-/commit/62b8190b71685492a2d15f7cbd212d3a3b80bf45
--
View it on GitLab: https://salsa.debian.org/med-team/hyphy/-/commit/62b8190b71685492a2d15f7cbd212d3a3b80bf45
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20250215/90689acf/attachment-0001.htm>
More information about the debian-med-commit
mailing list