Merge pull request #701 from veg/beta

2.3.7
veg · Nov 17, 2017 · 1586b66 · 1586b66
2 parents c8917f7 + 4782b04
commit 1586b66
Show file tree

Hide file tree

Showing 23 changed files with 3,031 additions and 1,386 deletions.
diff --git a/res/TemplateBatchFiles/GARD.bf b/res/TemplateBatchFiles/GARD.bf
diff --git a/res/TemplateBatchFiles/LEISR.bf b/res/TemplateBatchFiles/LEISR.bf
diff --git a/res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf b/res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf
@@ -39,6 +39,7 @@ namespace terms.fubar {
 };
 
 fubar.json = {
+                    terms.json.analysis:    {},
                     terms.json.input:       {},
                     terms.json.fits :       {},
                     terms.json.timers :     {},
@@ -97,6 +98,8 @@ fubar.analysis_description = {terms.io.info :
 
 io.DisplayAnalysisBanner ( fubar.analysis_description );
 
+fubar.json[terms.json.analysis] = fubar.analysis_description;
+
 namespace fubar {
     LoadFunctionLibrary ("modules/shared-load-file.bf");
     load_file ({utility.getGlobalValue("terms.prefix"): "fubar", utility.getGlobalValue("terms.settings") : {utility.getGlobalValue("terms.settings.branch_selector") : "selection.io.SelectAllBranches"}});

diff --git a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf
@@ -75,12 +75,12 @@ slac.json = {
     /**
         The dictionary of results to be written to JSON at the end of the run
     */
-    
+
 slac.display_orders =   {terms.original_name: -1,
                         terms.json.nucleotide_gtr: 0,
                         terms.json.global_mg94xrev: 1
                        };
-                       
+
 
 selection.io.startTimer (slac.json [terms.json.timers], "Total time", 0);
 
@@ -264,13 +264,13 @@ for (slac.i = 0; slac.i < Abs (slac.filter_specification); slac.i += 1) {
     }
 
     slac.branch_attributes = selection.substitution_mapper (slac.ancestors["MATRIX"], slac.ancestors["TREE_AVL"], slac.ancestors["AMBIGS"], slac.counts, slac.ancestors ["MAPPING"], slac.codon_data_info[terms.code]);
-    
+
     /*
     selection.io.json_store_branch_attribute(slac.json, terms.original_name, terms.json.node_label, 0,
                                              slac.i,
                                              slac.name_mapping);
     */
-    
+
     selection.io.json_store_branch_attribute(slac.json, terms.codon, terms.json.node_label, 0,
                                              slac.i,
                                              slac.branch_attributes[terms.codon]);
@@ -329,7 +329,11 @@ if (slac.samples > 0) {
         slac.table_output_options_samplers[terms.table_options.header] = TRUE;
         slac.sample.results = {};
 
-        slac.queue = mpi.CreateQueue ({"LikelihoodFunctions": {{slac.partitioned_mg_results[terms.likelihood_function]}}});
+        slac.queue = mpi.CreateQueue ({"LikelihoodFunctions": {{slac.partitioned_mg_results[terms.likelihood_function]}},
+                                        "Headers" : {{"libv3/all-terms.bf"}},
+                                        "Variables" : {{"slac.by_site","slac.AVERAGED","slac.RESOLVED","slac.by_branch"}}
+                                        });
+
 
         for (slac.s = 0; slac.s < slac.samples; slac.s+=1) {
 

diff --git a/res/TemplateBatchFiles/libv3/all-terms.bf b/res/TemplateBatchFiles/libv3/all-terms.bf
@@ -84,6 +84,9 @@ namespace terms{
 
 
     /* Term functions */
+    function characterFrequency (character) {
+        return "Equilibrium frequency for " + character;
+    }
     function nucleotideRate(fromC, toC) {
         return "Substitution rate from nucleotide " + fromC + " to nucleotide " + toC;
     }
@@ -172,6 +175,7 @@ namespace terms{
         equal      = "Equal frequencies";
         CF3x4      = "Corrected 3x4 frequency estimator";
         _20x1      = "Protein 20x1 estimator";
+        MLE        = "Maximum likelihood frequency estimator";
         predefined = "Based on a training set";
         binary     = "Binary character frequency estimator"; 
     }

diff --git a/res/TemplateBatchFiles/libv3/models/frequencies.bf b/res/TemplateBatchFiles/libv3/models/frequencies.bf
@@ -70,6 +70,44 @@ function frequencies.empirical.protein (model, namespace, datafilter) {
     return model;
 }
 
+/**
+ * Sets model's equilibrium frequency estimator to ML for protein data
+ * @name frequencies.empirical.protein
+ * @param {Dictionary} model
+ * @param {String} namespace
+ * @param {DataSetFilter} datafilter
+ * @returns {Dictionary} updated model
+ */
+function frequencies.ML.protein (model, namespace, datafilter) {
+    model = frequencies._aux.empirical.singlechar(model, namespace, datafilter);
+    // define 20 frequency parameters
+    // initialize to empirical freqs
+    // add to model parameter manifest
+    frequencies.ML.protein.emp = model[terms.efv_estimate];
+    model[terms.efv_estimate] = {20,1};
+    frequencies.ML.protein.variables = {20,1};
+    frequencies.ML.protein.scaler = namespace + ".frequency_scaler";
+
+    utility.ForEachPair (model[terms.alphabet], "_index", "_letter",
+                         '  
+                            _idx = _index[1];
+                            _freq_parameter = namespace + ".equilibrium_frequency_of." + _letter;
+                            frequencies.ML.protein.variables [_idx] = _freq_parameter;
+                            (model[terms.efv_estimate]) [_idx] = _freq_parameter + "/" + frequencies.ML.protein.scaler;
+                            parameters.DeclareGlobalWithRanges (_freq_parameter, frequencies.ML.protein.emp[_idx], 0, 1);
+                            model.generic.AddGlobal (model, _freq_parameter, terms.characterFrequency (_letter));
+                         '
+                         );
+
+    // constrain pi_A + ... + pi_A = 1, 
+
+    parameters.SetConstraint ( frequencies.ML.protein.scaler, Join ("+", frequencies.ML.protein.variables), "global");
+
+    model[terms.model.efv_estimate_name] = terms.frequencies.ml;
+    (model[terms.parameters])[terms.model.empirical] = -1; // correct for the restricted sum
+    return model;
+}
+
 /**
  * @name frequencies.empirical.corrected.CF3x4
  * @param {Dictionary} model

diff --git a/res/TemplateBatchFiles/libv3/models/model_functions.bf b/res/TemplateBatchFiles/libv3/models/model_functions.bf
@@ -228,7 +228,7 @@ function model.generic.DefineModel (model_spec, id, arguments, data_filter, esti
 
 
     // Set EFV field
-	Call (model.generic.DefineModel.model [terms.model.frequency_estimator], model.generic.DefineModel.model,
+	model.generic.DefineModel.model = Call (model.generic.DefineModel.model [terms.model.frequency_estimator], model.generic.DefineModel.model,
 													    id,
 													    data_filter);
 
@@ -240,10 +240,14 @@ function model.generic.DefineModel (model_spec, id, arguments, data_filter, esti
 	model.generic.DefineModel.model [terms.id] = id;
 
 
-
 	parameters.StringMatrixToFormulas (model.generic.DefineModel.model [terms.model.matrix_id],model.generic.DefineModel.model[terms.model.rate_matrix]);
-	utility.SetEnvVariable (model.generic.DefineModel.model [terms.model.efv_id], model.generic.DefineModel.model[terms.efv_estimate]);
-
+
+	if (Type ((model.generic.DefineModel.model[terms.efv_estimate])[0]) == "String") {
+	    parameters.StringMatrixToFormulas (model.generic.DefineModel.model [terms.model.efv_id],model.generic.DefineModel.model[terms.efv_estimate]);
+	} else {
+	    utility.SetEnvVariable (model.generic.DefineModel.model [terms.model.efv_id], model.generic.DefineModel.model[terms.efv_estimate]);
+	}
+
 
 	model.define_from_components (id, 	model.generic.DefineModel.model [terms.model.matrix_id], model.generic.DefineModel.model [terms.model.efv_id], model.generic.DefineModel.model [terms.model.canonical]);
 

diff --git a/res/TemplateBatchFiles/libv3/models/parameters.bf b/res/TemplateBatchFiles/libv3/models/parameters.bf
@@ -209,7 +209,16 @@ function parameters.SetValues(set) {
  * @returns nothing
  */
 lfunction parameters.ConstrainMeanOfSet (set, mean, namespace) {
-    unscaled = utility.Map (utility.Values (set), "_name_", "_name_ + '_scaler_variable'");
+    if (Type (set) == "AssociativeList") {
+        unscaled = utility.Map (utility.Values (set), "_name_", "_name_ + '_scaler_variable'");
+    } else {
+        if (Type (set) == "Matrix") {
+         unscaled = utility.Map (set, "_name_", "_name_ + '_scaler_variable'");
+        }
+        else {
+            return;
+        }
+    }
     global_scaler = namespace + ".scaler_variable";
     parameters.SetConstraint (global_scaler, Join ("+", unscaled), "global");
     utility.ForEach (set, "_name_", '
@@ -287,18 +296,18 @@ lfunction parameters.AddMultiplicativeTerm(matrix, term, do_empties) {
  */
 function parameters.StringMatrixToFormulas(id, matrix) {
     __N = Rows(matrix);
+    __M = Columns(matrix);
 
-    ExecuteCommands("`id` = {__N,__N}");
+    ExecuteCommands("`id` = {__N,__M}");
 
     for (__r = 0; __r < __N; __r += 1) {
-        for (__c = 0; __c < __N; __c += 1) {
+        for (__c = 0; __c < __M; __c += 1) {
 
-            if (__r != __c && Abs(matrix[__r][__c])) {
+            if (Abs(matrix[__r][__c])) {
                 ExecuteCommands("`id`[__r][__c] := " + matrix[__r][__c]);
             }
         }
     }
-
 }
 
 /**
@@ -399,6 +408,7 @@ lfunction parameters.GetConstraint(parameter) {
 function parameters.SetConstraint(id, value, global_tag) {
     if (Type(id) == "String") {
         if (Abs(id)) {
+            //console.log ("`global_tag` `id` := " + value);
             ExecuteCommands("`global_tag` `id` := " + value);
         }
     } else {

diff --git a/res/TemplateBatchFiles/libv3/models/protein.bf b/res/TemplateBatchFiles/libv3/models/protein.bf
@@ -14,7 +14,12 @@ models.protein.empirical_models = {{"LG",  "Empirical model of protein evolution
                                    {"WAG", "Empirical model of protein evolution from Whelan and Goldman (2001). Ref: https://doi.org/10.1093/oxfordjournals.molbev.a003851"},
                                    {"JTT", "Empirical model of protein evolution from Jones, Taylor, and Thornton (1996). Ref: https://doi.org/10.1093/bioinformatics/8.3.275"},
                                    {"JC69", "Empirical model of protein evolution with equal exchangeability rates among all amino acids, also known as JC69."},
-                                   {"mtMAM", "Empirical model of protein evolution for mammalian mitochondrial genomes from Yang N, Nielsen R, and Hasegawa M. (1998). Ref: http://dx.doi.org/10.1093/oxfordjournals.molbev.a025888"}};
+                                   {"mtMAM", "Empirical model of protein evolution for mammalian mitochondrial genomes from Yang N, Nielsen R, and Hasegawa M. (1998). Ref: http://dx.doi.org/10.1093/oxfordjournals.molbev.a025888"},
+                                   {"cpREV", "Empirical model of protein evolution for chloroplast genomes from  from Adachi et al. (2000). Ref: https://www.ncbi.nlm.nih.gov/pubmed/10795826"},
+                                   {"HIVBm", "Empirical model of protein evolution for between-host HIV sequences from Nickle et al. (2007). Ref: https://doi.org/10.1371/journal.pone.0000503"},
+                                   {"HIVWm", "Empirical model of protein evolution for within-host HIV sequences from Nickle et al. (2007). Ref: https://doi.org/10.1371/journal.pone.0000503"},
+                                   {"AB", "Empirical model of protein evolution for antibody sequences from Mirsky et al. (2015). Ref: https://doi.org/10.1093/molbev/msu340"}
+                                  };
 
 
 models.protein.dimensions = 20;