[WIP: Implementation of a New Algorithm] Extension of MixtureFinder to morphological data

main/phylotesting.cpp

@@ -6207,6 +6207,16 @@ void addModel(string model_str, string& new_model_str, string new_subst) {
     }
 }
 
+// initialise model frequency set in MixtureFinder for different sequence types
+char* initFreqSet(SeqType seq_type) {
+    switch (seq_type) {
+        case SEQ_CODON:   return ",F1X4,F3X4";
+        case SEQ_MORPH:   return "FQ";
+        case SEQ_PROTEIN: return ",FO";
+        default:          return "FO";
+    }
+}
+
 /**
  @brief Find the best component of a Q-mixture model while fixing other components
  @note This function is similar to runModelFinder, but it is adapted for optimisation of Q-Mixture model
@@ -6338,6 +6348,8 @@ CandidateModel findMixtureComponent(Params &params, IQTree &iqtree, ModelCheckpo
 
         skip_all_when_drop = false;
     } else if (mixture_action == MA_NUMBER_CLASS) {
+      // 2025/05/14 @HS6986
+      // I don't know when MA_NUMBER_CLASS is used and what should be expected in it for morphology!
         params.ratehet_set = iqtree.getModelFactory()->site_rate->name;
         // generate candidate models for the possible mixture models
         multi_class_str = "";
@@ -6355,8 +6367,10 @@ CandidateModel findMixtureComponent(Params &params, IQTree &iqtree, ModelCheckpo
             }
         }
         skip_all_when_drop = true;
-    } else if (mixture_action == MA_FIND_CLASS) {
-        char init_state_freq_set[] = "FO";
+      } else if (mixture_action == MA_FIND_CLASS) {
+        // 2025/05/14 @HS6986
+        // I don't know when MA_FIND_CLASS is used and what should be expected in it for morphology!
+        char* init_state_freq_set = initFreqSet(iqtree.aln->seq_type);
         if (!params.state_freq_set) {
             params.state_freq_set = init_state_freq_set;
         }
@@ -6404,6 +6418,7 @@ CandidateModel findMixtureComponent(Params &params, IQTree &iqtree, ModelCheckpo
 
         getModelSubst(iqtree.aln->seq_type, iqtree.aln->isStandardGeneticCode(), params.model_name,
                       params.model_set, params.model_subset, model_names);
+        StrVector model_names_orig = model_names;
 
         if (model_names.empty())
             return best_model; // empty
@@ -6432,11 +6447,54 @@ CandidateModel findMixtureComponent(Params &params, IQTree &iqtree, ModelCheckpo
                 }
             }
         }
-
-        for (i=0; i<model_names.size(); i++) {
-            string new_model_str;
-            addModel(model_str, new_model_str, model_names[i]);
-            candidate_models.push_back(CandidateModel(new_model_str, best_rate_name, iqtree.aln, 0));
+
+        if(getClassNum(model_str) == 1 && params.morph_mix_finder) {
+            auto isGTRXIncluded = [&]() -> bool {
+                return std::any_of(model_names_orig.begin(), model_names_orig.end(),
+                    [](const std::string& s) { return s == "GTRX" || s == "GTR"; });
+            };
+            auto isFOIncluded = [&]() -> bool {
+                return std::any_of(freq_names.begin(), freq_names.end(),
+                    [](const std::string& s) { return s == "+FO"; });
+            };
+            size_t new_model_strs_two_classes_len = (isGTRXIncluded() && isFOIncluded()) ? 3 : 1;
+            const char *new_model_strs_two_classes_GTRX_FO[] = {"MIX{MK+FO,MK+FO}", "MIX{GTRX+FQ,GTRX+FQ}", "MIX{GTRX+FO,GTRX+FO}"};
+            const char *new_model_strs_two_classes_GTRX[] = {"MIX{GTRX+FQ,GTRX+FQ}"};
+            const char *new_model_strs_two_classes_FO[] = {"MIX{MK+FO,MK+FO}"};
+            const char **new_model_strs_two_classes =
+                (isGTRXIncluded() && isFOIncluded())
+                    ? new_model_strs_two_classes_GTRX_FO
+                    : (isGTRXIncluded())
+                        ? new_model_strs_two_classes_GTRX
+                        : new_model_strs_two_classes_FO;
+
+            for (i=0; i<new_model_strs_two_classes_len; i++) {
+                string new_model_str = new_model_strs_two_classes[i];
+                candidate_models.push_back(CandidateModel(new_model_str, best_rate_name, iqtree.aln, 0));
+            }
+        } else {
+            for (i=0; i<model_names.size(); i++) {
+                string new_model_str;
+                addModel(model_str, new_model_str, model_names[i]);
+                if (params.morph_mix_finder) {
+                    int newModelCountGTRX = 0;
+                    size_t newModelPosGTRX = 0;
+                    while (((newModelCountGTRX = new_model_str.find("GTR", newModelPosGTRX)) != string::npos)) {
+                        ++newModelCountGTRX;
+                        newModelPosGTRX += 3;
+                    }
+                    int newModelCountFO = 0;
+                    size_t newModelPosFO = 0;
+                    while (((newModelPosFO = new_model_str.find("+FO", newModelPosFO)) != string::npos)) {
+                        ++newModelCountFO;
+                        newModelPosFO += 3;
+                    }
+                    if (newModelCountGTRX == 1 || newModelCountFO == 1) {
+                        continue;
+                    }
+                }
+                candidate_models.push_back(CandidateModel(new_model_str, best_rate_name, iqtree.aln, 0));
+            }
         }
 
         skip_all_when_drop = false;
@@ -6526,8 +6584,8 @@ void runMixtureFinderMain(Params &params, IQTree* &iqtree, ModelCheckpoint &mode
     bool better_model;
     double LR, df_diff, pvalue;
     string criteria_str;
-
-    char init_state_freq_set[] = "FO";
+    
+    char* init_state_freq_set = initFreqSet(iqtree->aln->seq_type);
     if (!params.state_freq_set) {
         params.state_freq_set = init_state_freq_set;
     }
@@ -6538,18 +6596,71 @@ void runMixtureFinderMain(Params &params, IQTree* &iqtree, ModelCheckpoint &mode
 
     // Step 0: (reorder candidate DNA models when -mset is used) build the nest-relationship network
     map<string, vector<string> > nest_network;
-    if (iqtree->aln->seq_type == SEQ_DNA) {
-        StrVector model_names, freq_names;
-        getModelSubst(iqtree->aln->seq_type, iqtree->aln->isStandardGeneticCode(), params.model_name,
-                      params.model_set, params.model_subset, model_names);
-        getStateFreqs(iqtree->aln->seq_type, params.state_freq_set, freq_names);
+    StrVector model_names, freq_names;
+    getModelSubst(iqtree->aln->seq_type, iqtree->aln->isStandardGeneticCode(), params.model_name,
+                  params.model_set, params.model_subset, model_names);
+    getStateFreqs(iqtree->aln->seq_type, params.state_freq_set, freq_names);
+
+    if (params.morph_mix_finder) {
+        auto areRatesJustified = [&]() -> bool {
+            return std::all_of(model_names.begin(), model_names.end(),
+                [](const std::string& s) { return s == "MK" || s == "GTRX" || s == "GTR"; });
+        };
+        if (!areRatesJustified()) {
+            outError("Error! Specifying models other than MK and GTRX for MixtureFinder for morphological data is not justified");
+        }
+
+        auto areFreqsJustified = [&]() -> bool {
+            return std::all_of(freq_names.begin(), freq_names.end(),
+                [](const std::string& s) { return s == "+FQ" || s == "+FO"; });
+        };
+        if (!areFreqsJustified()) {
+            outError("Error! Specifying frequencies other than +FQ and +FO for MixtureFinder for morphological data is not justified");
+        }
+
+        auto isMKIncluded = [&]() -> bool {
+            return std::any_of(model_names.begin(), model_names.end(),
+                [](const std::string& s) { return s == "MK"; });
+        };
+        if (!isMKIncluded()) {
+            outError("Error! Not specifying the MK model for MixtureFinder for morphological data is not justified");
+        }
+
+        auto isFQIncluded = [&]() -> bool {
+            return std::any_of(freq_names.begin(), freq_names.end(),
+                [](const std::string& s) { return s == "+FQ"; });
+        };
+        if (!isFQIncluded()) {
+            outError("Error! Not specifying the FQ frequency for MixtureFinder for morphological data is not justified");
+        }
+    }
+
+    auto isOnlyMKAndFQ = [&]() -> bool {
+        bool onlyMK = std::all_of(model_names.begin(), model_names.end(),
+            [](const std::string& s) { return s == "MK"; });
+        bool onlyFQ = std::all_of(freq_names.begin(), freq_names.end(),
+            [](const std::string& s) { return s == "+FQ"; });
+        return onlyMK && onlyFQ;
+    };
+    if (isOnlyMKAndFQ()) {
+        outError("Error! Running MixtureFinder only with the MK model and the FQ frequency is completely meaningless.\nPlease provide additional models and/or frequencies, such as GTRX, +F, and +FO, using -mset and/or -mfreq, if you really want to use MixtureFinder for your multistate data.");
+    }
 
+    if (iqtree->aln->seq_type == SEQ_DNA) {
         nest_network = generateNestNetwork(model_names, freq_names);
     }
-
+    
     // Step 1: run ModelFinder
     params.model_name = "";
     bool under_mix_finder = true;
+
+    string orig_model_set = params.model_set;
+    char *orig_state_freq_set = params.state_freq_set;
+    if (params.morph_mix_finder) {
+        params.model_set = "MK";
+        params.state_freq_set = "+FQ";
+    }
+
     runModelFinder(params, *iqtree, model_info, best_subst_name, best_rate_name, nest_network, under_mix_finder);
 
     // (cancel) Step 2: do tree search for this single-class model
@@ -6571,6 +6682,11 @@ void runMixtureFinderMain(Params &params, IQTree* &iqtree, ModelCheckpoint &mode
     model_info.getString("best_model_list_" + criteria_str, best_model_pre_list);
 
     // Step 3: keep adding a new class until no further improvement
+    if (params.morph_mix_finder) {
+        params.model_set = orig_model_set;
+        params.state_freq_set = orig_state_freq_set;
+    }
+
     if (params.opt_qmix_criteria == 1) {
         cout << endl << "Keep adding an additional class until the p-value from the likelihood ratio test > " << params.opt_qmix_pthres << endl;
     } else {
@@ -6649,13 +6765,29 @@ void runMixtureFinder(Params &params, IQTree* &iqtree, ModelCheckpoint &model_in
 
     if (iqtree->isSuperTree())
         outError("Error! The option -m '" + params.model_name + "' cannot work on data set with partitions");
-
+        
     if (iqtree->aln->seq_type != SEQ_DNA)
-        outError("Error! The option -m '" + params.model_name + "' can only work on DNA data set");
-
-    cout << "--------------------------------------------------------------------" << endl;
-    cout << "|                Running MixtureFinder                             |" << endl;
-    cout << "--------------------------------------------------------------------" << endl;
+        outWarning("MixtureFinder has not been tested for non-DNA data types. Be cautious about interpreting the results");
+
+    if (iqtree->aln->getMaxNumStates() > 6)
+        outWarning("Running MixtureFinder for the given data type can take much time. Please consider restricting the set of the models to test as much as possible");
+
+    if (iqtree->aln->seq_type == SEQ_PROTEIN && !params.force_aa_mix_finder)
+        outError("Error! We already have the mixture frequency vectors C10–C60 which are effective for modeling amino acid data.\nPlease make sure running MixtureFinder for your amino acid data makes sense.\nIf you are determined to do that, please add an option --force-aa-mix-finder to the command line.");
+
+    if (params.morph_mix_finder && iqtree->aln->seq_type != SEQ_MORPH) {
+        outError("Error! The option -morph-mix can only work on morphological data set");
+    }
+
+    if (params.morph_mix_finder) {
+        cout << "--------------------------------------------------------------------" << endl;
+        cout << "|            Running MixtureFinder for morphological data          |" << endl;
+        cout << "--------------------------------------------------------------------" << endl;
+    } else {
+        cout << "--------------------------------------------------------------------" << endl;
+        cout << "|                Running MixtureFinder                             |" << endl;
+        cout << "--------------------------------------------------------------------" << endl;
+    }
 
     // disable the bootstrapping
     int orig_gbo_replicates = params.gbo_replicates;

model/modelmixture.cpp

@@ -1403,6 +1403,8 @@ void ModelMixture::initMixture(string orig_model_name, string model_name, string
                     } else {
                         outError("The user defined frequency model is incorrect");
                     }
+                } else if (fstr == "+FU" || fstr == "+Fu") {
+                    model_freq = FREQ_USER_DEFINED;
                 } else if (fstr == "+FO" || fstr == "+Fo") {
                     model_freq = FREQ_ESTIMATE;
                 } else if (fstr == "+FQ" || fstr == "+Fq") {

model/modelmorphology.cpp

@@ -30,9 +30,11 @@ void ModelMorphology::init(const char *model_name, string model_params, StateFre
 		}
         num_params = 0;
     } else if (name == "GTR" || name == "GTRX") {
-        outWarning("GTRX multistate model will estimate " + convertIntToString(getNumRateEntries()-1) + " substitution rates that might be overfitting!");
-        outWarning("Please only use GTRX with very large data and always test for model fit!");
-        name = "GTRX";
+			if (num_states > 6 || phylo_tree->aln->getNPattern() < 100) {
+				outWarning("GTRX multistate model will estimate " + convertIntToString(getNumRateEntries()-1) + " substitution rates that might be overfitting for the length of your alignment/partition!");
+				outWarning("Please only use GTRX with large data and always test for model fit!");
+			}
+			name = "GTRX";
 	} else {
 		// if name does not match, read the user-defined model
 		readParameters(model_name);
@@ -122,20 +124,22 @@ string ModelMorphology::getName() {
 }
 
 string ModelMorphology::getNameParams(bool show_fixed_params) {
-    if (num_params == 0) return name;
     ostringstream retname;
     size_t pos_plus = name.find('+');
     if (pos_plus != string::npos) {
-        retname << name.substr(0, pos_plus) << '{';
+        retname << name.substr(0, pos_plus);
     } else {
-        retname << name << '{';
+        retname << name;
     }
-    int nrates = getNumRateEntries();
-    for (int i = 0; i < nrates; i++) {
-        if (i>0) retname << ',';
-        retname << rates[i];
+    if (num_params > 0 || show_fixed_params) {
+        retname << '{';
+        int nrates = getNumRateEntries();
+        for (int i = 0; i < nrates; i++) {
+            if (i > 0) retname << ',';
+            retname << rates[i];
+        }
+        retname << '}';
     }
-    retname << '}';
     if (pos_plus != string::npos) {
         retname << name.substr(pos_plus);
     } else {

utils/tools.cpp

@@ -1487,6 +1487,7 @@ void parseArg(int argc, char *argv[], Params &params) {
     params.ignore_checkpoint = false;
     params.checkpoint_dump_interval = 60;
     params.force_unfinished = false;
+    params.force_aa_mix_finder = false;
     params.print_all_checkpoints = false;
     params.suppress_output_flags = 0;
     params.ufboot2corr = false;
@@ -3729,6 +3730,17 @@ void parseArg(int argc, char *argv[], Params &params) {
                     params.opt_qmix_criteria = 2; // using information critera instead of likelihood-ratio test for estimation of number of classes for Q-Mixture model
                 */
                 continue;
+            }
+            if (strcmp(argv[cnt], "-morph-mix") == 0 || strcmp(argv[cnt], "--morph-mix") == 0) {
+                cnt++;
+                if (cnt >= argc)
+                    throw "Use -morph-mix <0|1>";
+                int in_option = convert_int(argv[cnt]);
+                if (in_option < 0 || in_option > 1)
+                    throw "Wrong option for -morpho-mix. Only 0 or 1 is allowed.";
+                if (in_option == 1)
+                    params.morph_mix_finder = true;
+                continue;
             }
 			if (strcmp(argv[cnt], "-a") == 0) {
 				cnt++;
@@ -5475,6 +5487,11 @@ void parseArg(int argc, char *argv[], Params &params) {
 				params.force_unfinished = true;
 				continue;
 			}
+
+      if (strcmp(argv[cnt], "-force-aa-mix-finder") == 0 || strcmp(argv[cnt], "--force-aa-mix-finder") == 0) {
+        params.force_aa_mix_finder = true;
+        continue;
+      }
 
 			if (strcmp(argv[cnt], "-cptime") == 0 || strcmp(argv[cnt], "--cptime") == 0) {
 				cnt++;

utils/tools.h

@@ -1707,6 +1707,11 @@ class Params {
      */
     bool opt_rhas_again;
 
+    /**
+            whether MixtureFinder optimizes the Q-mixture model for morphological data
+     */
+    bool morph_mix_finder;
+
     /**
             The method to optimize (and estimating the number of classes in) the Q-mixture model
             Method 1 (Old method)
@@ -2453,6 +2458,9 @@ class Params {
     /** true if ignoring the "finished" flag in checkpoint file */
     bool force_unfinished;
 
+    /** true if forcing IQ-TREE to run MixtureFinder for amino acid data */
+    bool force_aa_mix_finder;
+
     /** TRUE to print checkpoints to 1.ckp.gz, 2.ckp.gz,... */
     bool print_all_checkpoints;