Support Exporting Protein, Binary and Morphological Data to MrBayes

main/phyloanalysis.cpp

@@ -2602,6 +2602,82 @@ void printTrees(vector<string> trees, Params &params, string suffix) {
     treesOut.close();
 }
 
+void printMrBayesBlockFile(const char* filename, IQTree* &iqtree, bool inclParams) {
+    ofstream out;
+    try {
+        out.exceptions(ios::failbit | ios::badbit);
+        out.open(filename);
+
+        // Write Warning
+        out << "#nexus" << endl << endl
+            <<"[This MrBayes Block Declaration provides the retrieved values from the IQTree Run.]" << endl
+            << "[Note that MrBayes does not support a large collection of models, so defaults of 'nst=6' for DNA and 'wag' for Protein will be used if a model that does not exist in MrBayes is selected.]" << endl
+            << "[Furthermore, the Model Parameter '+R' will be replaced by '+G'.]" << endl
+            << "[This should be used as a Template Only.]" << endl << endl;
+
+        // Begin File, Print Charsets
+        out << "begin mrbayes;" << endl;
+    } catch (ios::failure &) {
+        outError(ERR_WRITE_OUTPUT, filename);
+    }
+
+    if (!iqtree->isSuperTree()) {
+        // Set Outgroup (if available)
+        if (!iqtree->rooted) out << "  outgroup " << iqtree->root->name << ";" << endl << endl;
+
+        out << "  [Using Model '" << iqtree->getModelName() << "']" << endl;
+        iqtree->getModel()->printMrBayesModelText(iqtree->getRate(), out, "all", "", false, inclParams);
+
+        out << endl << "end;" << endl;
+        out.close();
+        return;
+    }
+
+    auto superTree = (PhyloSuperTree*) iqtree;
+    auto saln = (SuperAlignment*) superTree->aln;
+    auto size = superTree->size();
+
+    for (int part = 0; part < size; part++) {
+        string name = saln->partitions[part]->name;
+        replace(name.begin(), name.end(), '+', '_');
+        out << "  charset " << name << " = ";
+
+        string pos = saln->partitions[part]->position_spec;
+        replace(pos.begin(), pos.end(), ',' , ' ');
+        out << pos << ";" << endl;
+    }
+
+    // Create Partition
+    out << endl << "  partition iqtree = " << size << ": ";
+    for (int part = 0; part < size; part++) {
+        if (part != 0) out << ", ";
+
+        string name = saln->partitions[part]->name;
+        replace(name.begin(), name.end(), '+', '_');
+        out << name;
+    }
+    out << ";" << endl;
+
+    // Set Partition for Use
+    out << "  set partition = iqtree;" << endl;
+
+    // Set Outgroup (if available)
+    if (!superTree->rooted) out << "  outgroup " << superTree->root->name << ";" << endl << endl;
+
+    // Partition-Specific Model Information
+    for (int part = 0; part < size; part++) {
+        PhyloTree* currentTree = superTree->at(part);
+
+        // MrBayes Partitions are 1-indexed
+        out << "  [Partition No. " << convertIntToString(part + 1) << ", Using Model '" << currentTree->getModelName() << "']" << endl;
+        currentTree->getModel()->printMrBayesModelText(currentTree->getRate(), out,
+                convertIntToString(part + 1), saln->partitions[part]->name, true, inclParams);
+        out << endl;
+    }
+    out << "end;" << endl;
+    out.close();
+}
+
 /************************************************************
  *  MAIN TREE RECONSTRUCTION
  ***********************************************************/
@@ -3111,8 +3187,10 @@ void runTreeReconstruction(Params &params, IQTree* &iqtree) {
         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 << endl << "Writing MrBayes Block Files..." << endl;
+        printMrBayesBlockFile((string(params.out_prefix) + ".mr_bayes_scheme.nex").c_str(), iqtree, false);
+        printMrBayesBlockFile((string(params.out_prefix) + ".mr_bayes_model.nex").c_str(), iqtree, true);
+        cout << endl;
     }
 
     cout << "BEST SCORE FOUND : " << iqtree->getCurScore() << endl;

model/modelbin.cpp

@@ -76,3 +76,56 @@ string ModelBIN::getNameParams(bool show_fixed_params) {
     }
     return retname.str();
 }
+
+void ModelBIN::printMrBayesModelText(RateHeterogeneity* rate, ofstream& out, string partition, string charset, bool isSuperTree, bool inclParams) {
+    // MrBayes does not support Invariable Modifier for Binary data
+    if (rate->isFreeRate() || rate->getPInvar() > 0.0) {
+        warnLogStream("MrBayes does not support Invariable Sites with Binary Data! +I has been ignored!", out);
+    }
+
+    // Lset Parameters
+    out << "  lset applyto=(" << partition << ") rates=";
+
+    // Free Rate should be substituted by +G (+I not supported)
+    bool hasGamma = rate->getGammaShape() != 0.0 || rate->isFreeRate();
+    if (hasGamma) {
+        // Rate Categories + Gamma
+        out << "gamma ngammacat=" << rate->getNRate();
+    } else
+        out << "equal";
+
+    out << ";" << endl;
+
+    if (!inclParams) {
+        if (freq_type == FREQ_EQUAL) out << "  prset applyto=(" << partition << ") statefreqpr=fixed(equal);" << endl;
+        return;
+    }
+
+    // Prset Parameters
+    out << "  prset applyto=(" << partition << ")";
+
+    // Freerate (+R)
+    // Get replacement Gamma Shape
+    if (rate->isFreeRate()) {
+        printMrBayesFreeRateReplacement(rate, charset, out, false);
+    }
+
+    // Gamma Distribution (+G/+R)
+    // Dirichlet is not available here, use fixed
+    if (rate->getGammaShape() > 0.0)
+        out << " shapepr=fixed(" << minValueCheckMrBayes(rate->getGammaShape()) << ")";
+
+    // State Frequencies
+    if (freq_type == FREQ_EQUAL)
+        out << " statefreqpr=fixed(equal)";
+    else {
+        out << " statefreqpr=dirichlet(";
+        for (int i = 0; i < num_states; ++i) {
+            if (i != 0) out << ", ";
+            out << minValueCheckMrBayes(state_freq[i]);
+        }
+        out << ")";
+    }
+
+    out << ";" << endl;
+}

model/modelbin.h

@@ -55,6 +55,18 @@ class ModelBIN : public ModelMarkov
     */
     virtual void startCheckpoint();
 
+    /**
+     * Print the model information in a format that can be accepted by MrBayes, using lset and prset.<br>
+     * By default, it simply prints a warning to the log and to the stream, stating that this model is not supported by MrBayes.
+     * @param rate the rate information
+     * @param out the ofstream to print the result to
+     * @param partition the partition to apply lset and prset to
+     * @param charset the current partition's charset. Useful for getting information from the checkpoint file
+     * @param isSuperTree whether the tree is a super tree. Useful for retrieving information from the checkpoint file, which has different locations for PhyloTree and PhyloSuperTree
+     * @param inclParams whether to include IQTree optimized parameters for the model
+     */
+    virtual void printMrBayesModelText(RateHeterogeneity* rate, ofstream& out, string partition, string charset, bool isSuperTree, bool inclParams);
+
 };
 
 #endif

model/modeldna.cpp

@@ -565,6 +565,100 @@ void ModelDNA::setVariables(double *variables) {
 //                      }
 }
 
-string ModelDNA::getModelDNACode() {
-    return param_spec;
+void ModelDNA::printMrBayesModelText(RateHeterogeneity* rate, ofstream& out, string partition, string charset, bool isSuperTree, bool inclParams) {
+    bool equalFreq = freq_type == FREQ_EQUAL;
+    short nst = 6;
+
+    // Find NST value (1 for JC/JC69/F81, 2 for K80/K2P/HKY/HKY85, 6 for SYM/GTR)
+    // NST is set to 6 by default (above), so check for JC/JC69/F81 and K80/K2P/HKY/HKY85
+    // If it returns a valid dna code, we can extract information from there (needed for user-defined models)
+    if (!param_spec.empty()) {
+        // if all equal
+        if (param_spec.find_first_not_of(param_spec[0]) == string::npos) {
+            nst = 1;
+            // if A-G=C-T and everything else is equal to first part of dna code
+        } else if (param_spec[1] == param_spec[4] &&
+                param_spec[0] == param_spec[2] && param_spec[0] == param_spec[3] && param_spec[0] == param_spec[5]) {
+            nst = 2;
+        }
+    } else if (!name.empty()) {
+        // Check the name of the model
+        if (strcmp(name.c_str(), "JC") == 0 || strcmp(name.c_str(), "JC69") == 0 || strcmp(name.c_str(), "F81") == 0) {
+            nst = 1;
+        } else if (strcmp(name.c_str(), "K80") == 0 || strcmp(name.c_str(), "K2P") == 0 || strcmp(name.c_str(), "HKY") == 0 || strcmp(name.c_str(), "HKY85") == 0) {
+            nst = 2;
+        }
+    }
+
+    // NucModel = 4By4: DNA Nucleotides (4 Options, A, C, G, T)
+    out << "  lset applyto=(" << partition << ") nucmodel=4by4 nst=" << nst << " rates=";
+
+    // RHAS Specification
+    // Free Rate should be substituted by +G+I
+    bool hasGamma = rate->getGammaShape() != 0.0 || rate->isFreeRate();
+    bool hasInvariable = rate->getPInvar() != 0.0 || rate->isFreeRate();
+    if (hasGamma) {
+        if (hasInvariable)
+            out << "invgamma";
+        else
+            out << "gamma";
+    } else if (hasInvariable)
+        out << "propinv";
+    else
+        out << "equal";
+
+    // Rate Categories
+    if (hasGamma)
+        out << " ngammacat=" << rate->getNRate();
+
+    out << ";" << endl;
+
+    // Priors (apart from Fixed Freq)
+    if (!inclParams) {
+        // If not include other params, simply set fixed frequency and return
+        if (!equalFreq) return;
+
+        out << "  prset applyto=(" << partition << ") statefreqpr=fixed(equal);" << endl;
+        return;
+    }
+
+    out << "  prset applyto=(" << partition << ")";
+
+    // Freerate (+R)
+    // Get replacement Gamma Shape + Invariable Sites
+    if (rate->isFreeRate()) {
+        printMrBayesFreeRateReplacement(isSuperTree, charset, out);
+    }
+
+    // Gamma Distribution (+G/+R)
+    // Dirichlet is not available here, use fixed
+    if (rate->getGammaShape() > 0.0)
+        out << " shapepr=fixed(" << minValueCheckMrBayes(rate->getGammaShape()) << ")";
+
+    // Invariable Sites (+I)
+    // Dirichlet is not available here, use fixed
+    if (rate->getPInvar() > 0.0)
+        out << " pinvarpr=fixed(" << minValueCheckMrBayes(rate->getPInvar()) << ")";
+
+    // Frequency of Nucleotides (+F)
+    if (equalFreq)
+        out << " statefreqpr=fixed(equal)";
+    else {
+        out << " statefreqpr=dirichlet(";
+        for (int i = 0; i < num_states; ++i) {
+            if (i != 0) out << ", ";
+            out << minValueCheckMrBayes(state_freq[i]);
+        }
+        out << ")";
+    }
+
+    // Reversible Rate Matrix
+    out << " revmatpr=dirichlet(";
+    for (int i = 0; i < getNumRateEntries(); ++i) {
+        if (i != 0) out << ", ";
+        out << minValueCheckMrBayes(rates[i]);
+    }
+
+    // Close revmatpr + prset
+    out << ");" << endl;
 }

model/modeldna.h

@@ -115,10 +115,16 @@ 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)
+     * Print the model information in a format that can be accepted by MrBayes, using lset and prset.<br>
+     * By default, it simply prints a warning to the log and to the stream, stating that this model is not supported by MrBayes.
+     * @param rate the rate information
+     * @param out the ofstream to print the result to
+     * @param partition the partition to apply lset and prset to
+     * @param charset the current partition's charset. Useful for getting information from the checkpoint file
+     * @param isSuperTree whether the tree is a super tree. Useful for retrieving information from the checkpoint file, which has different locations for PhyloTree and PhyloSuperTree
+     * @param inclParams whether to include IQTree optimized parameters for the model
      */
-    virtual string getModelDNACode();
+    virtual void printMrBayesModelText(RateHeterogeneity* rate, ofstream& out, string partition, string charset, bool isSuperTree, bool inclParams);
 
 protected:
 

model/modelmorphology.cpp

@@ -150,3 +150,54 @@ void ModelMorphology::writeInfo(ostream &out) {
         out << endl;
     }
 }
+
+void ModelMorphology::printMrBayesModelText(RateHeterogeneity* rate, ofstream& out, string partition, string charset, bool isSuperTree, bool inclParams) {
+    warnLogStream("MrBayes only supports Morphological Data with states from {0-9}!", out);
+    warnLogStream("Morphological Data with states {A-Z} may cause errors!", out);
+    warnLogStream("Use Morphological Models in MrBayes with Caution!", out);
+
+    // MrBayes does not support Invariable Modifier for Morph data
+    if (rate->isFreeRate() || rate->getPInvar() > 0.0) {
+        warnLogStream("MrBayes does not support Invariable Sites with Morphological Data! +I has been ignored!", out);
+    }
+    // MrBayes does not support State Frequency for Morph data
+    if (freq_type != FREQ_EQUAL) {
+        warnLogStream("MrBayes does not support non-equal frequencies for Morphological Data! Frequencies are left as the default! (All equal)", out);
+    }
+
+    // Lset Parameters
+    out << "  lset applyto=(" << partition << ") rates=";
+
+    // Free Rate should be substituted by +G (+I not supported)
+    bool hasGamma = rate->getGammaShape() != 0.0 || rate->isFreeRate();
+    if (hasGamma) {
+        // Rate Categories + Gamma
+        out << "gamma ngammacat=" << rate->getNRate();
+    } else
+        out << "equal";
+
+    out << ";" << endl;
+
+    // ctype (ordered or not)
+    if (strcmp(name.c_str(), "ORDERED") == 0)
+        out << "  ctype ordered: " << charset << ";" << endl;
+
+    // Since morphological doesn't have state frequencies, if there are no parameters to prset, return
+    if (!inclParams || !rate->isFreeRate() && rate->getGammaShape() <= 0.0) return;
+
+    // Prset Parameters
+    out << "  prset applyto=(" << partition << ")";
+
+    // Freerate (+R)
+    // Get replacement Gamma Shape
+    if (rate->isFreeRate()) {
+        printMrBayesFreeRateReplacement(rate, charset, out, false);
+    }
+
+    // Gamma Distribution (+G/+R)
+    // Dirichlet is not available here, use fixed
+    if (rate->getGammaShape() > 0.0)
+        out << " shapepr=fixed(" << minValueCheckMrBayes(rate->getGammaShape()) << ")";
+
+    out << ";" << endl;
+}

model/modelmorphology.h

@@ -78,6 +78,18 @@ class ModelMorphology: public ModelMarkov {
     virtual void readRates(istream &in) noexcept(false);
 
     virtual ~ModelMorphology();
+
+    /**
+     * Print the model information in a format that can be accepted by MrBayes, using lset and prset.<br>
+     * By default, it simply prints a warning to the log and to the stream, stating that this model is not supported by MrBayes.
+     * @param rate the rate information
+     * @param out the ofstream to print the result to
+     * @param partition the partition to apply lset and prset to
+     * @param charset the current partition's charset. Useful for getting information from the checkpoint file
+     * @param isSuperTree whether the tree is a super tree. Useful for retrieving information from the checkpoint file, which has different locations for PhyloTree and PhyloSuperTree
+     * @param inclParams whether to include IQTree optimized parameters for the model
+     */
+    virtual void printMrBayesModelText(RateHeterogeneity* rate, ofstream& out, string partition, string charset, bool isSuperTree, bool inclParams);
 };
 
 #endif /* MODELMORPHOLOGY_H_ */