Merge pull request #1695 from veg/develop

2.5.60 RC

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.12)
 project(HyPhy)
-cmake_policy(VERSION 3.0.0)
+#cmake_policy(VERSION 3.0.0)
 
 enable_testing()
 #enable_language(Fortran)
@@ -316,6 +316,11 @@ endif()
 MESSAGE ("Set compiler flags to ${DEFAULT_COMPILE_FLAGS}")
 MESSAGE ("Set linker flags to ${DEFAULT_LINK_FLAGS}")
 
+set(BLA_VENDOR Apple)
+find_package(BLAS)
+
+
+
 #-------------------------------------------------------------------------------
 # OpenMP support
 #-------------------------------------------------------------------------------
@@ -422,6 +427,10 @@ include_directories(
     contrib/SQLite-3.8.2 # SQLite
 )
 
+if (${BLAS_FOUND}) 
+	add_definitions (-D_SLKP_USE_APPLE_BLAS)
+	set(DEFAULT_LIBRARIES "${DEFAULT_LIBRARIES}")
+endif()
 
 #-------------------------------------------------------------------------------
 # shared hyphy hbl stdlib target
@@ -442,6 +451,10 @@ else(${OPENMP_FOUND})
     target_link_libraries(hyphy PRIVATE ${DEFAULT_LIBRARIES})
 endif(${OPENMP_FOUND})
 
+if(${BLAS_FOUND})
+    target_link_libraries(hyphy PRIVATE ${BLAS_LIBRARIES})
+endif ()
+
 add_custom_target(MP DEPENDS hyphy)
 
 install(
@@ -477,6 +490,10 @@ if(${MPI_FOUND})
         target_link_libraries(HYPHYMPI ${DEFAULT_LIBRARIES} ${MPI_LIBRARIES})
     endif(${OPENMP_FOUND})
 
+	if(${BLAS_FOUND})
+    	target_link_libraries(HYPHYMPI ${BLAS_LIBRARIES})
+	endif ()
+
     install(
         TARGETS HYPHYMPI
         RUNTIME DESTINATION bin

res/TemplateBatchFiles/SelectionAnalyses/BUSTED.bf

@@ -1,4 +1,4 @@
-RequireVersion ("2.5.38");
+RequireVersion ("2.5.51");
 
 
 LoadFunctionLibrary("libv3/all-terms.bf"); // must be loaded before CF3x4
@@ -22,17 +22,22 @@ utility.SetEnvVariable ("NORMALIZE_SEQUENCE_NAMES", TRUE);
 
 
 busted.analysis_description = {
-                               terms.io.info : 
-"BUSTED (branch-site unrestricted statistical test of episodic diversification) uses a random effects branch-site model fitted 
-jointly to all or a subset of tree branches in order to test for alignment-wide evidence of episodic diversifying selection. 
-Assuming there is evidence of positive selection (i.e. there is an omega > 1),  BUSTED will also perform a quick evidence-ratio 
-style analysis to explore which individual sites may have been subject to selection. v2.0 adds support for synonymous rate variation, 
-and relaxes the test statistic to 0.5 (chi^2_0 + chi^2_2). Version 2.1 adds a grid search for the initial starting point.
-Version 2.2 changes the grid search to LHC, and adds an initial search phase to use adaptive Nedler-Mead. Version 3.0 implements the option
-for branch-site variation in synonymous substitution rates. Version 3.1 adds HMM auto-correlation option for SRV, and binds SRV distributions for multiple branch sets.
+
+                               terms.io.info : "BUSTED (branch-site unrestricted statistical test of episodic diversification) uses a random effects branch-site model fitted  jointly to all or a subset of tree branches in order to test for alignment-wide evidence of episodic diversifying selection. Assuming there is evidence of positive selection (i.e. there is an omega > 1),  BUSTED will also perform a quick evidence-ratio style analysis to explore which individual sites may have been subject to selection. v2.0 adds support for synonymous rate variation, and relaxes the test statistic to 0.5 (chi^2_0 + chi^2_2).
+
+Version 2.1 adds a grid search for the initial starting point.
+
+Version 2.2 changes the grid search to LHC, and adds an initial search phase to use adaptive Nedler-Mead. 
+
+Version 3.0 implements the option for branch-site variation in synonymous substitution rates. 
+
+Version 3.1 adds HMM auto-correlation option for SRV, and binds SRV distributions for multiple branch sets.
+
 Version 4.0 adds support for multiple hits (MH), ancestral state reconstruction saved to JSON, and profiling of branch-site level support for selection / multiple hits.
+
 Version 4.2 adds calculation of MH-attributable fractions of substitutions.
-Version 4.5 adds an "error absorption" component [experimental]
+
+Version 4.5 adds an 'error absorption' component [experimental]
 ",
                                terms.io.version : "4.5",
                                terms.io.reference : "*Gene-wide identification of episodic selection*, Mol Biol Evol. 32(5):1365-71, *Synonymous Site-to-Site Substitution Rate Variation Dramatically Inflates False Positive Rates of Selection Analyses: Ignore at Your Own Peril*, Mol Biol Evol. 37(8):2430-2439",
@@ -42,6 +47,8 @@ Version 4.5 adds an "error absorption" component [experimental]
                                terms.io.requirements : "in-frame codon alignment and a phylogenetic tree (optionally annotated with {})"
                               };
 
+
+console.log (busted.analysis_description);
 io.DisplayAnalysisBanner (busted.analysis_description);
 
 busted.FG = "Test";
@@ -548,8 +555,8 @@ for (busted.partition_index = 0; busted.partition_index < busted.partition_count
     busted.global_scaler_list [busted.partition_index] = "busted.bl.scaler_" + busted.partition_index;
     parameters.DeclareGlobalWithRanges (busted.global_scaler_list [busted.partition_index], 3, 0, 1000);
     busted.initial_ranges      [busted.global_scaler_list [busted.partition_index]] = {
-                terms.lower_bound : 2.5,
-                terms.upper_bound : 3.5
+                terms.lower_bound : 2.0,
+                terms.upper_bound : 4.0
             };
 }
 
@@ -562,7 +569,7 @@ busted.initial_grid         = estimators.LHC (busted.initial_ranges,busted.initi
 //estimators.CreateInitialGrid (busted.initial_grid, busted.initial_grid.N, busted.initial_grid_presets);
 
 busted.initial_grid = utility.Map (busted.initial_grid, "_v_", 
-    'busted._renormalize (_v_, "busted.distribution", busted.initial.test_mean, busted.error_sink)'
+    'busted._renormalize_with_weights (_v_, "busted.distribution", busted.initial.test_mean, busted.error_sink)'
 );
 
 if (busted.has_background) { //GDD rate category
@@ -604,24 +611,29 @@ debug.checkpoint = utility.GetEnvVariable ("DEBUG_CHECKPOINT");
 if (Type (debug.checkpoint) != "String") {
 
     // constrain nucleotide rate parameters
-    busted.tmp_fixed = models.FixParameterSetRegExp (terms.nucleotideRatePrefix,busted.test.bsrel_model);
+    // copy global nucleotide parameter estimates to busted.test.bsrel_model)
+
+
+    busted.tmp_fixed = models.FixParameterSetRegExpFromReference (terms.nucleotideRatePrefix,busted.test.bsrel_model, busted.final_partitioned_mg_results[terms.global]);
+
+    busted.grid_search.results =  estimators.FitLF (busted.filter_names, busted.trees, busted.model_map, busted.final_partitioned_mg_results, busted.model_object_map,                    
+        {
+            terms.run_options.retain_lf_object: TRUE,
+            terms.run_options.proportional_branch_length_scaler : 
+                                                    busted.global_scaler_list,
+
+            terms.run_options.optimization_settings : 
+                {
+                    "OPTIMIZATION_METHOD" : "nedler-mead",
+                    "MAXIMUM_OPTIMIZATION_ITERATIONS" : 500,
+                    "OPTIMIZATION_PRECISION" : busted.nm.precision
+                } ,
+
+            terms.search_grid     : busted.initial_grid,
+            terms.search_restarts : busted.N.initial_guesses
+        }
+    );
 
-
-    busted.grid_search.results =  estimators.FitLF (busted.filter_names, busted.trees, busted.model_map, busted.final_partitioned_mg_results, busted.model_object_map, {
-        "retain-lf-object": TRUE,
-        terms.run_options.proportional_branch_length_scaler : 
-                                                busted.global_scaler_list,
-
-        terms.run_options.optimization_settings : 
-            {
-                "OPTIMIZATION_METHOD" : "nedler-mead",
-                "MAXIMUM_OPTIMIZATION_ITERATIONS" : 500,
-                "OPTIMIZATION_PRECISION" : busted.nm.precision
-            } ,
-
-        terms.search_grid     : busted.initial_grid,
-        terms.search_restarts : busted.N.initial_guesses
-    });
 
     parameters.RemoveConstraint (busted.tmp_fixed );
 
@@ -1247,6 +1259,64 @@ lfunction busted._renormalize (v, distro, mean, skip_first) {
 
 //------------------------------------------------------------------------------
 
+lfunction busted._renormalize_with_weights (v, distro, mean, skip_first) {
+
+    parameters.SetValues (v);
+    m = parameters.GetStickBreakingDistribution (^distro);
+    d = Rows (m);
+
+    if (skip_first) {
+        m0 = m[0][0]*m[0][1];
+    } else {
+        m0 = 0;
+    }
+
+    over_one = m[d-1][0] * m[d-1][1];
+
+    if (over_one >= mean*0.95) {
+       //console.log ("OVERAGE");
+       new_weight = mean * Random (0.9, 0.95) / m[d-1][0];
+       diff = (m[d-1][1] - new_weight)/(d-1);
+       for (k = (skip_first != 0); k < d-1; k += 1) {
+            m[k][1] += diff;
+       }
+       m[d-1][1] = new_weight;
+    }
+
+    over_one = m[d-1][0] * m[d-1][1];
+    under_one = (+(m[-1][0] $ m[-1][1])) / (1-m[d-1][1]); // current mean
+
+    for (i = (skip_first != 0); i < d-1; i+=1) {
+        m[i][0] = m[i][0] * mean / under_one;
+    }
+
+    under_one = +(m[-1][0] $ m[-1][1]);
+
+    for (i = (skip_first != 0); i < d; i+=1) {
+        m[i][0] = m[i][0] * mean / under_one;
+    }
+
+    m = m%0;
+    wts = parameters.SetStickBreakingDistributionWeigths (m[-1][1]);
+
+    //console.log (v);
+
+    for (i = (skip_first != 0); i < d; i+=1) {
+        (v[((^distro)["rates"])[i]])[^"terms.fit.MLE"] = m[i][0];
+    }
+    for (i = (skip_first != 0); i < d-1; i+=1) {
+        (v[((^distro)["weights"])[i]])[^"terms.fit.MLE"] = wts[i];
+    }
+
+    //console.log (v);
+
+    //assert (0);
+    return v;
+
+}
+
+//------------------------------------------------------------------------------
+
 lfunction busted.get_multi_hit (model_fit) {
     params = selection.io.extract_global_MLE_re (model_fit, '^' + ^'terms.parameters.multiple_hit_rate');
     for (k,v; in; selection.io.extract_global_MLE_re (model_fit, '^' + ^'terms.parameters.triple_hit_rate')) {

res/TemplateBatchFiles/SelectionAnalyses/FEL.bf

@@ -736,7 +736,6 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
            }
     );
 
-
 
     if (^"fel.ci") {
         if (!sim_mode) {
@@ -761,6 +760,8 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
         ci = None;
     }
 
+
+
     if (sim_mode) {
         lrt = 2*results[1][0];    
     } else {
@@ -777,9 +778,10 @@ 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;
     parameters.SetConstraint ("fel.beta_scaler_test","fel.alpha_scaler", "");
-
+ 
     Optimize (results, ^lf);
 
+
     if (sim_mode) {
         return lrt - 2*results[1][0];
     } else {
@@ -792,11 +794,11 @@ lfunction fel.handle_a_site (lf, filter_data, partition_index, pattern_info, mod
         tree_name = (lfInfo["Trees"])[0];
         sum = (BranchLength (^tree_name, -1)*^"fel.selected_branches_index")[0];
         if (^"fel.resample") {
+
             N = ^"fel.resample";
             sims = {};
             GetDataInfo (fi, ^((lfInfo["Datafilters"])[0]), "PARAMETERS");
             //Export (lfe, ^lf);
-            //console.log (lfe);
             for (i = 0; i < N; i+=1) {        
                 DataSet null_sim = SimulateDataSet (^lf);
                 DataSetFilter null_filter = CreateFilter (null_sim,3,,,fi["EXCLUSIONS"]);

res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf

@@ -650,12 +650,13 @@ return fubar.json;
 //----------------------------------------------------------------------------
 
 lfunction fubar.ComputeENFP_CI (p_i,sig_level) {
+
     N = Abs (p_i);
-
     PDF_old = {{1-p_i[0],p_i[0]}};
     PDF = PDF_old;
 
     for (i = 1; i < N; i+=1) {
+
         PDF = {1,i+2};
         PDF[0] = PDF_old[0] * (1-p_i[i]);
         for (X = 1; X < i+1; X+=1) {
@@ -678,6 +679,7 @@ lfunction fubar.ComputeENFP_CI (p_i,sig_level) {
         sum += PDF[_idx];
     }
 
+
     return {{lb__, _idx__}}
 }