Integrate SPRTA of CMAPLE into IQ-TREE

main/phyloanalysis.cpp

@@ -4609,6 +4609,7 @@ bool runCMaple(Params &params)
     {
         // record the start time
         auto start = getRealTime();
+        auto start_cpu = getCPUTime();
 
         // Dummy variables
         std::string aln_path(params.aln_file);
@@ -4686,6 +4687,15 @@ bool runCMaple(Params &params)
             const std::string input_treefile(params.user_file ? params.user_file : "");
             cmaple::Tree tree(&aln, &model, input_treefile, (params.fixed_branch_length == BRLEN_FIX), cmaple::ParamsBuilder().build());
 
+            // transfer SPRTA options if any
+            if (params.compute_SPRTA)
+            {
+                tree.params->compute_SPRTA = params.compute_SPRTA;
+                tree.params->compute_SPRTA_zero_length_branches = params.SPRTA_zero_branches;
+                tree.params->print_SPRTA_less_info_seqs = params.SPRTA_zero_branches;
+                tree.params->output_alternative_spr = params.out_alter_spr;
+            }
+
             // Infer a phylogenetic tree
             const cmaple::Tree::TreeSearchType tree_search_type = cmaple::Tree::parseTreeSearchType(params.tree_search_type_str);
             std::ostream null_stream(0);
@@ -4721,7 +4731,23 @@ bool runCMaple(Params &params)
             ofstream out = ofstream(output_treefile);
             out << tree.exportNewick(tree_format);
             out.close();
-
+
+            // Write tree file in NEXUS format (if computing SPRTA)
+            if (params.compute_SPRTA)
+            {
+                ofstream out = ofstream(output_treefile + ".nexus");
+                out << tree.exportNexus(tree_format);
+                out.close();
+            }
+
+            // export a TSV file if SPRTA is computed and we output alternative SPRs
+            if (params.compute_SPRTA && params.out_alter_spr)
+            {
+                ofstream out = ofstream(output_treefile + ".tsv");
+                out << tree.exportTSV();
+                out.close();
+            }
+
             // Show model parameters
             if (cmaple::verbose_mode > cmaple::VB_QUIET)
             {
@@ -4736,14 +4762,20 @@ bool runCMaple(Params &params)
             // Show information about output files
             std::cout << "Analysis results written to:" << std::endl;
             std::cout << "Maximum-likelihood tree:       " << output_treefile << std::endl;
-            /*if (params.aLRT_replicates)
-                std::cout << "Tree with aLRT-SH values:      " << prefix + ".aLRT_SH.treefile" << std::endl;*/
+            if (params.compute_SPRTA)
+                std::cout << "Tree in NEXUS format:          " << output_treefile + ".nexus" << std::endl;
+            if (params.compute_SPRTA && params.out_alter_spr)
+                std::cout << "Meta data in TSV format:       " << output_treefile + ".tsv" << std::endl;
             std::cout << "Screen log file:               " << prefix + ".log" << std::endl << std::endl;
 
             // show runtime
             auto end = getRealTime();
-            if (cmaple::verbose_mode > cmaple::VB_QUIET)
-                cout << "CMAPLE Runtime: " << end - start << "s" << endl;
+            auto end_cpu = getCPUTime();
+            cout << "Total CPU time used: " << fixed << end_cpu - start_cpu << " sec (" <<
+                convert_time(end_cpu-start_cpu) << ")" << endl;
+            cout << "Total wall-clock time used: " << fixed << end - start << " sec (" <<
+                convert_time(end-start) << ")" << endl;
+            cout << endl;
         }
     }
     catch (std::invalid_argument e)

utils/tools.cpp

@@ -1526,6 +1526,11 @@ void parseArg(int argc, char *argv[], Params &params) {
     params.mutation_file = "";
     params.site_starting_index = 0;
 
+    // ----------- SPRTA ----------
+    params.compute_SPRTA = false;
+    params.SPRTA_zero_branches = false;
+    params.out_alter_spr = false;
+
     // store original params
     for (cnt = 1; cnt < argc; cnt++) {
         params.original_params = params.original_params + argv[cnt] + " ";
@@ -3123,6 +3128,24 @@ void parseArg(int argc, char *argv[], Params &params) {
                 params.inference_alg = ALG_CMAPLE;
                 continue;
             }
+            if (strcmp(argv[cnt], "--sprta") == 0 ||
+                strcmp(argv[cnt], "-sprta") == 0) {
+              params.compute_SPRTA = true;
+
+              continue;
+            }
+            if (strcmp(argv[cnt], "--zero-branch-supp") == 0 ||
+                strcmp(argv[cnt], "-zero-branch-supp") == 0) {
+              params.SPRTA_zero_branches = true;
+
+              continue;
+            }
+            if (strcmp(argv[cnt], "--out-alter-spr") == 0 ||
+                strcmp(argv[cnt], "-out-alter-spr") == 0) {
+              params.out_alter_spr = true;
+
+              continue;
+            }
             if (strcmp(argv[cnt], "--out-csv") == 0) {
                 params.output_format = FORMAT_CSV;
                 continue;
@@ -5917,6 +5940,12 @@ void usage(char* argv[]) {
     cout << "                       divergence is low, otherwise, apply IQ-TREE algorithm." << endl;
     cout << "  --pathogen-force     Apply CMAPLE tree search algorithm regardless" << endl;
     cout << "                       of sequence divergence." << endl;
+    cout << "  --alrt <num_rep>     Specify number of replicates to compute SH-aLRT." << endl;
+    cout << "  --sprta              Compute SPRTA (DeMaio et al., 2024) branch supports." << endl;
+    cout << "  --zero-branch-supp   Compute SPRTA supports for zero-length branches." << endl;
+    cout << "  --out-alter-spr      Output alternative SPRs and their SPRTA supports." << endl;
+    cout << "  -T <num_thread>      Specify number of threads used for computing" << endl;
+    cout << "                       branch supports (SH-aLRT or SPRTA)." << endl;
     cout << endl;
     //cout << "HIDDEN OPTIONS: see the source code file pda.cpp::parseArg()" << endl;
 
@@ -6235,6 +6264,12 @@ void usage_iqtree(char* argv[], bool full_command) {
     << "                       divergence is low, otherwise, apply IQ-TREE algorithm." << endl
     << "  --pathogen-force     Apply CMAPLE tree search algorithm regardless" << endl
     << "                       of sequence divergence." << endl
+    << "  --alrt <num_rep>     Specify number of replicates to compute SH-aLRT." << endl
+    << "  --sprta              Compute SPRTA (DeMaio et al., 2024) branch supports." << endl
+    << "  --zero-branch-supp   Compute SPRTA supports for zero-length branches." << endl
+    << "  --out-alter-spr      Output alternative SPRs and their SPRTA supports." << endl
+    << "  -T <num_thread>      Specify number of threads used for computing" << endl
+    << "                       branch supports (SH-aLRT or SPRTA)." << endl
 
 
 

utils/tools.h

@@ -2776,6 +2776,21 @@ class Params {
      * TRUE to make the processes of outputting->re-inputting a tree consistent
      */
     bool make_consistent;
+
+    /**
+     * TRUE to compute the SPRTA branch supports
+     */
+    bool compute_SPRTA;
+
+    /**
+     * TRUE to compute the SPRTA for zero-length branches
+     */
+    bool SPRTA_zero_branches;
+
+    /**
+     * TRUE to output the alternative SPRs with their supports in the tree
+     */
+    bool out_alter_spr;
 
     /**
     *  Mutation file that specifies pre-defined mutations occurs at nodes