Masel Lab documentation changes

doc/Advanced-Tutorial.md

@@ -523,19 +523,16 @@ This will print an output file `example.phy.mlrate` that looks like:
 	10      0.00001
 
 
-Trimming alignment sites by likelihood
+Trimming alignments
 --------------------------------------
 <div class="hline"></div>
 
-Phylogenetic inference can be highly sensitive to fast-evolving, saturated or
-erroneous sites in a sequence alignment. Many studies used Gblocks 
-([Castresana, 2000]) or trimAl ([Capella-Gutiérrez et al., 2009])
-to trim alignment sites prior to phylogenetic reconstruction.
-Here, we present an alternative approach called  **trimmed log-likelihood**, a 
-robust phylogenetics method, that automatically detects and trims such sites 
-during tree search directly.
+Aligned sites that do not share a common ancestry do not provide useful information for phylogenetic inference. Prequal is designed to remove non-homologous sequences prior to alignment ([Whelan et al., 2018]). 
 
-The trimed log-likelihood method works by dynamically excluding a user-defined proportion of sites
+Even when sequences as a whole are homologous, there can still be sites that are not homologous for all sequences in the alignment. Removing whole alignment columns can risk losing information about closely related species. Three programs that convert predicted misaligned characters to dashes without removing whole columns are Divvier ([Ali et al., 2019]), Taper ([Zhang et al., 2021]), and CLOAK ([Wheeler et al., 2026]). There is a tradeoff between precision and recall in alignmet filtering, with stricter filters removing more errors, but also removing more correctly aligned characters. These three programs performed best out of the tested alignment filtering methods on a precision recall curve, with Divvier being the strictest filter and CLOAK being the gentlest filter ([Wheeler et al., 2026]). 
+
+IQ-Tree contains a method for trimming sites directly during tree search called **trimmed log-likelihood**. 
+The trimmed log-likelihood method works by dynamically excluding a user-defined proportion of sites
 with the lowest log-likelihood values during the tree search. As the search
 progresses,  the likelihood of each site is recalculated at each step using
 current tree and model parameters.  This ensures that site removal is always
@@ -580,8 +577,6 @@ Where to go from here?
 
 See [Command Reference](Command-Reference) for a complete list of all options available in IQ-TREE.
 
-[Capella-Gutiérrez et al., 2009]: https://doi.org/10.1093/bioinformatics/btp348
-[Castresana, 2000]: https://doi.org/10.1093/oxfordjournals.molbev.a026334
 [Gadagkar et al., 2005]: https://doi.org/10.1002/jez.b.21026
 [Kishino et al., 1990]: https://doi.org/10.1007/BF02109483
 [Kishino and Hasegawa, 1989]: https://doi.org/10.1007/BF02100115
@@ -595,4 +590,7 @@ See [Command Reference](Command-Reference) for a complete list of all options av
 [Strimmer and Rambaut, 2002]: https://doi.org/10.1098/rspb.2001.1862
 [Mayrose et al., 2004]: https://doi.org/10.1093/molbev/msh194
 [Yang, 1995]: http://www.genetics.org/content/139/2/993.abstract
-
+[Whelan et al., 2018]: https://doi.org/10.1093/bioinformatics/bty448
+[Ali et al., 2019]: https://doi.org/10.1093/molbev/msz142
+[Zhang et al., 2021]: https://doi.org/10.1111/2041-210X.13696
+[Wheeler et al., 2026]: https://doi.org/10.64898/2025.12.01.691663

doc/Estimating-amino-acid-substitution-models.md

@@ -171,6 +171,17 @@ To estimate a non-reversible model from a folder of alignments:
 	grep -A 22 "can be used as input for IQ-TREE" train_plant.NONREV.iqtree | tail -n21 > NQ.plant
 
 
+Filtering substitution model training alignments
+-----------------------------------
+
+Multiple sequence alignments often contain errors that can inflate the rates of less mutationally accessible amino acid substitutions that require 2-3 nucleotide substitutions ([Wheeler et al., 2026]). To train cleaned QC substitution models, we recommend using the program Divvier to filter each alignment in the training set, using the partial filtering option ([Ali et al., 2019]). Divvier can be found at (https://github.com/simonwhelan/Divvier). 
+
+The filtered alignments can then be used to infer Q matrices as described above. Because the amino acid frequency parameters inferred for a QC matrix will reflect the frequencies of the filtered alignments rather than the originals, QC substitution models should always use emperical amino acid frequencies with the +F setting. Note that the +F option is not valid for nonreversible NQ matrices, which cannot be separated into symmetric exchangeability scores and amino acid frequencies. Therefore filtered alignments should not be used to train NQ models.
+
+If you use the new models trained on cleaned multiple sequence alignments (QC.pfam…), please cite:
+
+>  Andrew L. Wheeler, Chiragdeep Chatur, Peter W Goodman, Robert C Edgar, Gavin A Huttley, Joanna Masel (2026), Improved gene tree inference from removing alignment errors both from focal genes and when training substitution models. in press at Molecular Biology & Evolution. <https://doi.org/10.64898/2025.12.01.691663>
+
 Estimating linked exchangeabilities
 -----------------------------------
 
@@ -240,3 +251,5 @@ If you use the GTRspmix model, please cite the following paper:
 [Ran et al., 2018]: https://doi.org/10.1098/rspb.2018.1012
 [Banos et al., 2024]: https://doi.org/10.1093/molbev/msae174
 [Harada et al., 2026]: https://doi.org/10.64898/2026.06.18.729217
+[Wheeler et al., 2026]: https://doi.org/10.64898/2025.12.01.691663
+[Ali et al., 2019]: https://doi.org/10.1093/molbev/msz142

doc/Quickstart.md

@@ -122,9 +122,13 @@ A few typically analyses are listed in the following. Note that it is assumed th
 
         iqtree -s example.phy
 
-* Infer maximum-likelihood tree using `GTR+I+G` model:
+* Infer maximum-likelihood tree using `GTR+I+R` model:
 
-        iqtree -s example.phy -m GTR+I+G
+        iqtree -s example.phy -m GTR+I+R
+
+* Infer maximum-likelihood tree for protein sequence using `QC.pfam+F+I+R` model:
+
+      iqtree -s aa_example.phy -m QC.pfam+F+I+R
 
 * Perform ModelFinder without subsequent tree inference:
 

doc/Substitution-Models.md

@@ -197,10 +197,21 @@ IQ-TREE supports all common empirical amino-acid exchange rate matrices (alphabe
 | Q.pfam   | nuclear | General Q matrix ([Minh et al., 2021]) estimated from Pfam version 31 database ([El-Gebali et al., 2018]). | 
 | Q.plant  | nuclear | Q matrix ([Minh et al., 2021]) estimated for plants ([Ran et al., 2018]). | 
 | Q.yeast  | nuclear | Q matrix ([Minh et al., 2021]) estimated for yeasts ([Shen et al., 2018]). | 
+| QC.archaea  | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for archaeal from the HAMAP database ([Lima et al., 2009]). | 
+| QC.bacteria   | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for bacteria from the HAMAP database ([Lima et al., 2009]). | 
+| QC.bird | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for birds ([Jarvis et al., 2015]) . | 
+| QC.insect | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for insects ([Misof et al., 2014]) . | 
+| QC.mammal | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for mammals from OrthoMaM database version 12 ([Allio et al., 2024]) . | 
+| QC.pfam   | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) from Pfam version 31 database ([El-Gebali et al., 2018]). | 
+| QC.plant | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for plants ([Ran et al., 2018]) . | 
+| QC.yeast | nuclear | Q matrix ([Minh et al., 2021]) estimated on cleaned alignments ([Wheeler et al., 2026]) for for yeasts ([Shen et al., 2018]) . | 
 | rtREV    | viral | Retrovirus ([Dimmic et al., 2002]). |
 | VT       | nuclear | General 'Variable Time' matrix ([Mueller and Vingron, 2000]). |
 | WAG      | nuclear | General matrix ([Whelan and Goldman, 2001]). |
 
+>**TIP**: Cleaned "QC" substitution models make mammalian single copy ortholog trees on average more similar to the mammalian species tree ([Wheeler, et al. 2026]). Nevertheless, alignment errors in the focal gene alignment will cause ModelFinder to prefer unfiltered Q to QC. If QC is desired, ModelFinder should not be used. QC models should always be used in conjunction with the +F option to use empirical amino acid frequencies. 
+{: .tip}
+
 ### Protein mixture models
 
 IQ-TREE also supports a series of protein mixture models:
@@ -425,10 +436,10 @@ IQ-TREE supports all common rate heterogeneity across sites models:
 | +G       | discrete Gamma model ([Yang, 1994]) with default 4 rate categories. The number of categories can be changed with e.g. `+G8`. |
 | +GC       | continuous Gamma model ([Yang, 1994]) (for AliSim only). |
 | +I+G     | invariable site plus discrete Gamma model ([Gu et al., 1995]). |
-| +R       | FreeRate model ([Yang, 1995]; [Soubrier et al., 2012]) that generalizes the `+G` model by relaxing the assumption of Gamma-distributed rates. The number of categories can be specified with e.g. `+R6` (default 4 categories if not specified). The FreeRate model typically fits data better than the `+G` model and is recommended for analysis of large data sets. |
+| +R       | FreeRate model ([Yang, 1995]; [Soubrier et al., 2012]) that generalizes the `+G` model by relaxing the assumption of Gamma-distributed rates. The number of categories can be specified with e.g. `+R6` (default 4 categories if not specified). Gamma models produce artifacts, especially for large datasets ([Ferretti et al., 2026]). FreeRate models are therefore recommended. |
 | +I+R     | invariable site plus FreeRate model. |
 
->**TIP**: The new ModelFinder (`-m MFP` option) tests the FreeRate model, whereas the standard procedure (`-m TEST`) does not.
+>**TIP**: TIP: ModelFinder (-m MFP option) tests the FreeRate model, whereas the outdated (-m TEST) does not. To test only the FreeRate model, use (-mrate I,R,I+R).
 {: .tip}
 
 Users can fix the parameters of the model. For example, `+I{0.2}` will fix the proportion of invariable sites (pinvar) to 0.2; `+G{0.9}` will fix the Gamma shape parameter (alpha) to 0.9; `+I{0.2}+G{0.9}` will fix both pinvar and alpha. To fix the FreeRate model parameters, use the syntax `+Rk{w1,r1,...,wk,rk}` (replacing `k` with the number of categories). Here, `w1, ..., wk` are the weights and `r1, ..., rk` the rates for each category. 
@@ -501,3 +512,7 @@ Users can fix the parameters of the model. For example, `+I{0.2}` will fix the p
 [davidcerny/GEOS26100-Fall2022]: https://github.com/davidcerny/GEOS26100-Fall2022
 [Černý & Simonoff (2023)]: https://doi.org/10.1038/s41598-023-35784-3
 [Harada et al., 2026]: https://doi.org/10.64898/2026.06.18.729217
+[Wheeler et al., 2026]: https://doi.org/10.64898/2025.12.01.691663
+[Allio et al., 2024]: https://doi.org/10.1093/nar/gkad834
+[Lima et al., 2009]: https://doi.org/10.1093/nar/gkn661
+[Ferretti et al., 2026]: https://doi.org/10.1093/sysbio/syag037

doc/Tutorial.md

@@ -178,7 +178,7 @@ change the prefix with:
 This prevents output files being overwritten when you perform multiple analyses on the same alignment within the same folder.
 
 
-Choosing the right substitution model
+Choosing the best-fitting substitution model
 -------------------------------------
 <div class="hline"></div>
 

doc/recipes/How-to-clean-an-alignment-before-making-a-gene-tree.md

@@ -0,0 +1,42 @@
+---
+layout: workshop
+title: "How to clean an alignment before making a gene tree"
+author: Andrew Wheeler
+date:    2026-06-07
+docid: 100
+---
+
+# Introduction
+
+This recipe explains how to filter out errors from a single gene multiple sequence alignment, prior to gene tree inference, using the "CLOAK" function in Muscle5. 
+
+CLOAK filters based on consensus among an ensemble of alternative alignments. First, an alignment ensemble is inferred by varying the initial guide tree and HMM parameters used for multiple sequence alignment. If characters are found to be homologous across all alignments in the ensemble, they are retained in the filtered alignment. Otherwise they are replaced by dashes. If a whole alignment column does not have consistent homology, but subsets of the column do, the column is split into 2 or columns that are both retained in the filtered alignment. This makes the program act as a gentle filter, since it removes a significant number of non-homologous characters without losing many informative homologies. CLOAK is recommended for filtering single gene alignments, which are often sensitive to any loss of data from their limited number of informative sites ([Wheeler, et al. 2026]).
+
+# What you need
+
+Muscle5: https://www.drive5.com/muscle/
+
+IQ-TREE: http://www.iqtree.org/
+
+# Example Usage
+
+The first step is to generate an ensemble of alignments using the "stratified" option in Muscle5. This creates a set of 16 variant alignments, all output in a single ensemble fasta file, instead of a single alignment.
+
+```
+muscle -align sequences.fasta -stratified -output ensemble.efa
+```
+
+CLOAK can then generate a single filtered alignment, based on consensus among the ensemble, using the following command:
+
+```
+muscle -cloak input_ensemble_file -mincol <integer> -output <output_file_name>
+```
+
+Arguments:
+- -input_ensemble_file : Path to the input alingment ensemble (EFA) file
+- -mincol <integer>   : Minimum number of non-gap characters required per column
+                        for that column to be retained in the output.
+                        Default value of 2 if not specified
+- -output <filename>  : Name of the file where the filtered MSA will be written. By default this will be {input_file_name}.cloak.fa
+
+[Wheeler, et al. 2026]: https://doi.org/10.64898/2025.12.01.691663