Support Exporting DNA/RNA Analysis to MrBayes Block Files (iqtree#260)

* Support Exporting DNA/RNA Analysis to MrBayes Block Files

* Fix Formatting Issues

* Move Functions to Supplementary, Fix +R Remapping

Author: IntegerLimit

.gitignore

@@ -14,6 +14,7 @@ tags
 test_scripts/iqtree_test_cmds.txt
 pllrepo/
 build/
+cmake-build-debug/
 /Default-clang
 .idea/
 .idea/codeStyleSettings.xml

main/phyloanalysis.cpp

@@ -1156,6 +1156,12 @@ void printOutfilesInfo(Params &params, IQTree &tree) {
     if (params.optimize_linked_gtr) {
         cout << "  GTRPMIX nex file:              " << params.out_prefix << ".GTRPMIX.nex" << endl;
     }
+
+    if (params.mr_bayes_output) {
+        cout << endl << "MrBayes block written to:" << endl;
+        cout << "  Base template file:            " << params.out_prefix << ".mr_bayes_scheme.nex" << endl;
+        cout << "  Template file with parameters: " << params.out_prefix << ".mr_bayes_model.nex" << endl;
+    }
     cout << endl;
 
 }
@@ -3100,8 +3106,13 @@ void runTreeReconstruction(Params &params, IQTree* &iqtree) {
 
     }
     if (iqtree->isSuperTree()) {
-        ((PhyloSuperTree*) iqtree)->computeBranchLengths();
-        ((PhyloSuperTree*) iqtree)->printBestPartitionParams((string(params.out_prefix) + ".best_model.nex").c_str());
+        auto superTree = (PhyloSuperTree*) iqtree;
+        superTree->computeBranchLengths();
+        superTree->printBestPartitionParams((string(params.out_prefix) + ".best_model.nex").c_str());
+    }
+    if (params.mr_bayes_output) {
+        iqtree->printMrBayesBlock((string(params.out_prefix) + ".mr_bayes_scheme.nex").c_str(), false);
+        iqtree->printMrBayesBlock((string(params.out_prefix) + ".mr_bayes_model.nex").c_str(), true);
     }
 
     cout << "BEST SCORE FOUND : " << iqtree->getCurScore() << endl;

main/phylotesting.h

@@ -390,5 +390,12 @@ string convertSeqTypeToSeqTypeName(SeqType seq_type);
 
 string detectSeqTypeName(string model_name);
 
+/**
+ * get string name from a SeqType object
+ * @param seq_type input sequence type
+ * @return name
+ */
+string getSeqTypeName(SeqType seq_type);
+
 
 #endif /* PHYLOTESTING_H_ */

model/modeldna.cpp

@@ -564,3 +564,7 @@ void ModelDNA::setVariables(double *variables) {
 //                              j++;
 //                      }
 }
+
+string ModelDNA::getModelDNACode() {
+    return param_spec;
+}

model/modeldna.h

@@ -114,6 +114,12 @@ class ModelDNA : public ModelMarkov
      */
     virtual void computeTipLikelihood(PML::StateType state, double *state_lk);
 
+    /**
+     * Get the Model DNA 'code', in form 'abcdef', used with ModelDNA model
+     * Returns empty string by default (this is not a dna specific model)
+     */
+    virtual string getModelDNACode();
+
 protected:
 
 	/**

model/modelsubst.h

@@ -392,6 +392,12 @@ class ModelSubst: public Optimization, public CheckpointFactory
     */
     virtual ModelSubst *getMutationModel() { return this; }
 
+    /**
+     * Get the Model DNA 'code', in form 'abcdef', used with ModelDNA model
+     * Returns empty string by default (this is not a dna specific model)
+     */
+     virtual string getModelDNACode() { return ""; }
+
 	/*****************************************************
 		Checkpointing facility
 	*****************************************************/

tree/phylosupertree.cpp

@@ -23,6 +23,7 @@
 #include "main/phylotesting.h"
 #include "model/partitionmodel.h"
 #include "utils/MPIHelper.h"
+#include "utils/tools.h"
 
 PhyloSuperTree::PhyloSuperTree()
  : IQTree()
@@ -1519,29 +1520,37 @@ void PhyloSuperTree::writeBranches(ostream &out) {
     }
 }
 
-void PhyloSuperTree::printBestPartitionParams(const char *filename) {
+void printCharsets(SuperAlignment* saln, size_t size, ofstream &out)
+{
+    for (int part = 0; part < size; part++) {
+        string name = saln->partitions[part]->name;
+        replace(name.begin(), name.end(), '+', '_');
+        out << "  charset " << name << " = ";
+        if (!saln->partitions[part]->aln_file.empty()) out << saln->partitions[part]->aln_file << ": ";
+        /*if (saln->partitions[part]->seq_type == SEQ_CODON)
+            out << "CODON, ";*/
+        if (!saln->partitions[part]->sequence_type.empty())
+            out << saln->partitions[part]->sequence_type << ", ";
+        string pos = saln->partitions[part]->position_spec;
+        replace(pos.begin(), pos.end(), ',' , ' ');
+        out << pos << ";" << endl;
+    }
+}
+
+void PhyloSuperTree::printBestPartitionParams(const char *filename)
+{
     try {
         ofstream out;
         out.exceptions(ios::failbit | ios::badbit);
         out.open(filename);
         out << "#nexus" << endl
         << "begin sets;" << endl;
-        int part;
-        SuperAlignment *saln = (SuperAlignment*)aln;
-        for (part = 0; part < size(); part++) {
-            string name = saln->partitions[part]->name;
-            replace(name.begin(), name.end(), '+', '_');
-            out << "  charset " << name << " = ";
-            if (!saln->partitions[part]->aln_file.empty()) out << saln->partitions[part]->aln_file << ": ";
-            /*if (saln->partitions[part]->seq_type == SEQ_CODON)
-                out << "CODON, ";*/
-            out << saln->partitions[part]->sequence_type << ", ";
-            string pos = saln->partitions[part]->position_spec;
-            replace(pos.begin(), pos.end(), ',' , ' ');
-            out << pos << ";" << endl;
-        }
+        auto *saln = (SuperAlignment*)aln;
+
+        printCharsets(saln, size(), out);
+
         out << "  charpartition mymodels =" << endl;
-        for (part = 0; part < size(); part++) {
+        for (int part = 0; part < size(); part++) {
             string name = saln->partitions[part]->name;
             replace(name.begin(), name.end(), '+', '_');
             if (part > 0) out << "," << endl;
@@ -1555,3 +1564,45 @@ void PhyloSuperTree::printBestPartitionParams(const char *filename) {
         outError(ERR_WRITE_OUTPUT, filename);
     }
 }
+
+void PhyloSuperTree::printMrBayesFreeRateReplacement(RateHeterogeneity *rate, string &charset, ofstream &out) {
+    // Call PhyloTree's function inside the correct checkpoint struct
+    checkpoint->startStruct("PartitionModelPlen");
+    checkpoint->startStruct(charset);
+    PhyloTree::printMrBayesFreeRateReplacement(rate, charset, out);
+    checkpoint->endStruct();
+    checkpoint->endStruct();
+}
+
+void PhyloSuperTree::printMrBayesBlock(const char *filename, bool inclParams)
+{
+    ofstream out = startMrBayesBlock(filename);
+    auto *saln = (SuperAlignment*)aln;
+    printCharsets(saln, size(), out);
+    out << endl;
+
+    // Create Partition
+    size_type listSize = size();
+    out << "  partition iqtree = " << listSize << ": ";
+    for (int part = 0; part < listSize; part++) {
+        string name = saln->partitions[part]->name;
+        replace(name.begin(), name.end(), '+', '_');
+        out << name;
+        if (part != listSize - 1) out << ", ";
+        else out << ";" << endl;
+    }
+
+    // Set Partition for Use
+    out << "  set partition = iqtree;" << endl;
+
+    // Set Outgroup (if available)
+    if (!rooted) out << "  outgroup " << root->name << ";" << endl << endl;
+
+    // Partition-Specific Model Information
+    for (int part = 0; part < listSize; part++) {
+        // MrBayes Partitions are 1-indexed
+        printMrBayesModelText(this, at(part), out, convertIntToString(part + 1), saln->partitions[part]->name, inclParams);
+    }
+    out << "end;" << endl;
+    out.close();
+}

tree/phylosupertree.h

@@ -475,12 +475,25 @@ class PhyloSuperTree : public IQTree, public vector<PhyloTree* >
      */
     void printBestPartitionParams(const char *filename);
 
-    
+    /**
+        print mr bayes block with partition information and best model parameters
+        @param filename output file name
+        @param inclParams whether to include iqtree params
+      */
+    virtual void printMrBayesBlock(const char *filename, bool inclParams);
+
+    /**
+     * Prints the replacement prior settings for +R or +R+I.
+     * @param rate the heterogeneity rate
+     * @param charset the (original) charset of the current partition. An empty string if not a partitioned tree
+     * @param out the ofstream to print to
+     */
+    virtual void printMrBayesFreeRateReplacement(RateHeterogeneity* rate, string &charset, ofstream &out);
+
     /** True when mixed codon with other data type */
     bool rescale_codon_brlen;
-    
-    int totalNNIs, evalNNIs;
 
+    int totalNNIs, evalNNIs;
 };
 
 #endif