PHYLOFIT(1)

📦 version 1.8, in debian bookworm, debian trixie, ubuntu jammy, ubuntu noble, ubuntu questing, ubuntu resolute, ubuntu stonking

PHYLOFIT(1) — User Commands 💖

okay bestie, sit down, hydrate, because this is the phyloFit manual and it is a lot — but in the best way. It's long, so it's coming to you in three parts. This is part one. Let's go. 🫶

🎀 NAME 🎀

phyloFit — Fits one or more tree models to a multiple alignment of DNA. slay ✨

🌸 DESCRIPTION 🌸

So basically: phyloFit fits one or more tree models to a multiple alignment of DNA sequences by maximum likelihood, using the tree topology and substitution model you tell it to use. 💅

If you define categories of sites via --features and --catmap (more on those below, bestie), then it's gonna estimate a separate model for each category. Each model gets its own little description written out to its own file, with the suffix .mod. 🎀

Those .mod files minimally include:

  • a substitution rate matrix,
  • a tree with branch lengths,
  • and estimates of nucleotide equilibrium frequencies.

They may also include info about parameters for modeling rate variation. Obsessed. 🦋

✨ SYNOPSIS ✨

phyloFit [OPTIONS] <msa_fname>

okay so <msa_fname> needs to be a multiple alignment in FASTA format, or one of several alternative formats (peep --msa-format for that 👀). For backward compatibility, this argument may be preceded by -m or --msa. Also, heads up: --tree is required in most cases. And by default, every output file gets the prefix "phyloFit" (see --out-root if you wanna change that).

🦋 EXAMPLE 🦋

(If you're like me, you want some basic examples first and the list of all the options later. Same, bestie. Same.)

  1. Compute the distance between two aligned sequences (in FASTA file pair.fa) under the REV model.
phyloFit pair.fa

(output is to phyloFit.mod; distance in substitutions per site appears in the TREE line in the output file)

  1. Fit a phylogenetic model to an alignment of human, chimp, mouse, and rat sequences. Use the HKY85 substitution model. Write output to files with prefix "myfile".
phyloFit --tree "((human,chimp),(mouse,rat))" --subst-mod HKY85
--out-root myfile primate-rodent.fa
  1. As above, but use the discrete-gamma model for rate variation, with 4 rate categories.
phyloFit --tree "((human,chimp),(mouse,rat))" --subst-mod HKY85
--out-root myfile --nrates 4 primate-rodent.fa
  1. As above, but use genome-wide data, stored in the compact "sufficient-statistics" format (can be produced with "msa_view -o SS").
phyloFit --tree "((human,chimp),(mouse,rat))" --subst-mod HKY85
--out-root myfile --nrates 4 --msa-format SS primate-rodent.ss
  1. Fit a context-dependent phylogenetic model (U2S) to an alignment of human, mouse, and rat sequences. Use an EM algorithm for parameter optimization and relax the convergence criteria a bit (recommended with context-dependent models). Write a log file for the optimization procedure. Consider only non-overlapping pairs of sites.
phyloFit --tree "(human,(mouse,rat))" --subst-mod U2S --EM
--precision MED --non-overlapping --log u2s.log --out-root
hmr-u2s hmr.fa
  1. As above, but allow overlapping pairs of sites, and compute likelihoods by assuming Markov-dependence of columns (see Siepel & Haussler, 2004). The EM algorithm can no longer be used (optimization will be much slower).
phyloFit --tree "(human,(mouse,rat))" --subst-mod U2S --precision
MED --log u2s-markov.log --markov hmr.fa
  1. Compute a likelihood using parameter estimates obtained in (5) and an assumption of Markov dependence. This provides a lower bound on the likelihood of the Markov-dependent model.
phyloFit --init-model hmr-u2s.mod --lnl --markov hmr.fa
  1. Given an alignment of several mammalian sequences (mammals.fa), a tree topology (tree.nh), and a set of gene annotations in GFF (genes.gff), fit separate models to sites in 1st, 2nd, and 3rd codon positions. Use the REV substitution model. Assume coding regions have feature type 'CDS'.
phyloFit --tree tree.nh --features genes.gff --out-root
mammals-rev --catmap "NCATS = 3; CDS 1-3" --do-cats 1,2,3
mammals.fa

(output will be to mammals-rev.cds-1.mod, mammals-rev.cds-2.mod, and mammals-rev.cds-3.mod)

💅 OPTIONS 💅

  • --tree, -t <tree_fname>|<tree_string> — (Required if more than three species, or more than two species and a non-reversible substitution model, e.g., UNREST, U2, U3.) Name of file or literal string defining tree topology. The tree must be in Newick format, with the label at each leaf equal to the index or name of the corresponding sequence in the alignment (indexing begins with 1). Examples: --tree "(1,(2,3))", --tree "(human,(mouse,rat))". Currently, the topology must be rooted. When a reversible substitution model is used, the root is ignored during the optimization procedure. 🌳
  • --subst-mod, -s — one of JC69|F81|HKY85|HKY85+Gap|REV|SSREV|UNREST|R2|R2S|U2|U2S|R3|R3S|U3|U3S (default REV). Nucleotide substitution model, bestie. JC69, F81, HKY85, REV, and UNREST have the usual meanings (see, e.g., Yang, Goldman, and Friday, 1994). SSREV is a strand-symmetric version of REV. HKY85+Gap is an adaptation of HKY that treats gaps as a fifth character (courtesy of James Taylor 🫶). The others — all considered "context-dependent" — are as defined in Siepel and Haussler, 2004. The options --EM and --precision MED are recommended with context-dependent models (see below).
  • --msa-format, -i — FASTA|PHYLIP|MPM|MAF|SS (default is to guess the format from the file contents). Alignment format. FASTA is as usual. PHYLIP is compatible with the formats used in the PHYLIP and PAML packages. MPM is the format used by the MultiPipMaker aligner and some other of Webb Miller's older tools. MAF ("Multiple Alignment Format") is used by MULTIZ/TBA and the UCSC Genome Browser. SS is a simple format describing the sufficient statistics for phylogenetic inference (distinct columns or tuple of columns and their counts). Note that the program "msa_view" can be used for file conversion. 🔄
  • --out-root, -o <output_fname_root> — (default "phyloFit"). Use the specified string as the root filename for all files created. 📁
  • --min-informative, -I <ninf_sites> — Require at least <ninf_sites> "informative" sites — i.e., sites at which at least two non-gap and non-missing-data ('N' or '*') characters are present. Default is 50. ✅
  • --gaps-as-bases, -G — Treat alignment gap characters ('-') like ordinary bases. By default, they are treated as missing data. 🙈
  • --ignore-branches, -b <branches> — Ignore the specified branches in likelihood computations and parameter estimation, and treat the induced subtrees as independent. Can be useful for likelihood ratio tests. The argument <branches> should be a comma-separated list of nodes in the tree, indicating the branches above these nodes, e.g., human-chimp,cow-dog. (See tree_doctor --name-ancestors regarding names for ancestral nodes.) This option does not currently work with --EM. ⚠️
  • --quiet, -q — Proceed quietly. 🤫
  • --help, -h — Print this help message. (Options for controlling and monitoring the optimization procedure.)
  • --lnl, -L — (for use with --init-model) Simply evaluate the log likelihood of the specified tree model, without performing any further optimization. Can be used with --post-probs, --expected-subs, and --expected-total-subs. 📊
  • --EM, -E — Fit model(s) using EM rather than the BFGS quasi-Newton algorithm (the default). 🔥
  • --precision, -p — HIGH|MED|LOW (default HIGH). Level of precision to use in estimating model parameters. Affects convergence criteria for iterative algorithms: higher precision means more iterations and longer execution time. ⏳
  • --log, -l <log_fname> — Write a log to <log_fname> describing details of the optimization procedure. 📝
  • --init-model, -M <mod_fname> — Initialize with the specified tree model. By choosing good starting values for parameters, it is possible to reduce execution time dramatically. If this option is chosen, --tree is not allowed. The substitution model used in the given model will be used unless --subst-mod is also specified. Note: currently only one mod_fname may be specified; it will be used for all categories. 🧠
  • --init-random, -r — Initialize parameters randomly. Can be used multiple times to test whether the m.l.e. is real. 🎲
  • --seed, -D <seed> — Provide a random number seed for choosing initial parameter values (usually with --init-random, though random values are used in some other cases as well). Should be an integer >=1. If not provided, the seed is chosen based on the current time. ⏰
  • --init-parsimony, -y — Initialize branch lengths using parsimony counts for the given data. Only currently implemented for models with a single character state (ie, not di- or tri-nucleotides). Other --init options such as --init-random or --init-model can be used in conjunction to initialize substitution matrix parameters. 🌿
  • --print-parsimony, -Y <filename> — Print parsimony score to the given file, and quit. (Does not optimize or report likelihoods.) 🖨️
  • --clock, -z — Assume a molecular clock in estimation. Causes the distances to all descendant leaves to be equal for each ancestral node and cuts the number of free branch-length parameters roughly in half. 🕰️
  • --scale-only, -B — (for use with --init-model) Estimate only the scale of the tree, rather than individual branch lengths (branch proportions fixed). Equilibrium frequencies and rate-matrix parameters will still be estimated unless --no-freqs and --no-rates are used. 📏
  • --scale-subtree, -S <node_name> — (for use with --scale-only) Estimate separate scale factors for the subtree beneath the identified node and the rest of the tree. The branch leading to the subtree is included with the subtree. If ":loss" or ":gain" is appended to <node_name>, subtree scale is constrained to be greater than or less than (respectively) scale for the rest of the tree. 🌲
  • --estimate-freqs, -F — Estimate equilibrium frequencies by maximum likelihood, rather than approximating them by the relative frequencies in the data. If using the SSREV model, this option implies --sym-freqs. 💯
  • --sym-freqs, -W — Estimate equilibrium frequencies, assuming freq(A)=freq(T) and freq(C)=freq(G). This only works for an alphabet ACGT (and possibly gap). This option implies --estimate-freqs. ⚖️
  • --no-freqs, -f — (for use with --init-model) Do not estimate equilibrium frequencies; just use the ones from the given tree model. 🙅
  • --no-rates, -n — …and bestie, we are NOT done, this is literally where part 1 cuts off. Part 2 is coming, don't move. 💅
  • (for use with --init-model) Do not estimate rate-matrix parameters; just use the ones from the given tree model. No notes, no edits, just vibes from the model you handed over. 💅
  • --ancestor, -A <seqname> — treat the sequence you name as the root of the tree 🫶. The tree topology has to define this sequence as a child of the root (and fr, in practice the branch from the root down to that sequence is kept, but it gets constrained to length zero).
  • --error, -e <filename> — for every single parameter, report the estimate, the variance, and the 95% confidence interval, all printed to the filename you gave, one parameter per line. Receipts. ✨
  • --no-opt, -O <param_list> — hold the parameters in your comma-separated param_list constant at their initial values. Lowkey important: this applies only to the "main" model, and not to any models you defined with --alt-mod.

Your param list can contain values like branches to hold branch lengths constant, ratematrix, backgd, or ratevar to hold the entire rate matrix, the equilibrium frequencies, or the rate variation parameters constant (respectively, bestie). There are also substitution model-specific parameters like kappa (the transition/transversion rate ratio).

Note: to hold certain branches constant but still optimize others, put an exclamation point in the newick-formatted tree after the branch lengths that should be held constant. This can be useful for enforcing a star-phylogeny. But heads up — the two branches coming from the root of the tree are treated as one. So they should both be held constant, or not held constant at all. No in-between, no halfsies. 🙈
This option does not work with --scale-only or --clock. Not compatible, periodt.
  • --bound <param_name[lower_bound,upper_bound]> — set boundaries for a parameter. lower_bound or upper_bound may be an empty string to keep the default. For example, --bound gc_param[1,] will set the lower bound for gc_param to 1 (keeping the upper bound at infinity), for a GC model. 🌸

This only applies to parameters for the model in the "main" tree, but similar syntax can be used within the --alt-mod arguments. Can be used multiple times to set boundaries for different parameters. Stack 'em up. 🎀

  • --selection <selection_param> — use selection in the model (it's also implied if --init-model is used and contains a selection parameter). Selection scales rate matrix entries by selection_param/(1-exp(-selection-param)); this is done after the rate matrix is scaled to set the expected number of substitutions per unit time to 1. If you're using codon models, selection acts only on non-synonymous mutations.

(Options for modeling rate variation) — yeah, this is where we're headed next. 🦋

  • --nrates, -k <nratecats> (default 1). Number of rate categories to use. Specifying a value greater than one causes the discrete gamma model for rate variation to be used (Yang, 1994). The icon himself. 👑
  • --alpha, -a <alpha> (for use with --nrates). Initial value for alpha, the shape parameter of the gamma distribution. Default is 1. ⭐
  • --rate-constants, -K <rate_consts> — use a non-parametric mixture model for rates, instead of assuming a gamma distribution. The argument <rate_consts> must be a comma-delimited list explicitly defining the rate constants to be used. The "weight" (mixing proportion) associated with each rate constant will be estimated by EM (so yes, this option implies --EM). If --alpha is used together with this option, then the mixing proportions will be initialized to reflect a gamma distribution with the shape parameter you specified. It's giving ✨custom✨.

(Options for separate handling of sites in different annotation categories)

  • --features, -g <fname> — an annotations file (GFF or BED format) describing features on one or more sequences in the alignment 🍓. Together with a category map (see --catmap), it will be taken to define site categories, and a separate model will be estimated for each category. If no category map is specified, a category will be assumed for each type of feature, and they'll be numbered in the order of appearance of the features. Features are assumed to use the coordinate frame of the first sequence in the alignment and should be non-overlapping (see refeature --unique).
  • --catmap, -c <fname>|<string> — (optionally use with --features) mapping of feature types to category numbers. You can either give a filename or an "inline" description of a simple category map, like --catmap "NCATS = 3 ; CDS 1-3" or --catmap "NCATS = 1 ; UTR 1". Note that category 0 is reserved for "background" — i.e. everything that is not described by a defined feature type. 💜
  • --do-cats, -C <cat_list> — (optionally use with --features) estimate models for only the specified categories (a comma-delimited list of categories, by name or number). Default is to fit a model for every category. Go big or go selective, your call. 💅
  • --reverse-groups, -R <tag> — (optionally use with --features) group features by <tag> (e.g. "transcript_id" or "exon_id") and reverse complement segments of the alignment corresponding to groups on the reverse strand. Groups must be non-overlapping (see refeature --unique). Super useful with categories corresponding to strand-specific phenomena (e.g. codon positions). 🧚

(Options for context-dependent substitution models)

  • --markov, -N — (for use with context-dependent substitution models, and not available with --EM.) Assume Markov dependence of alignment columns, and compute the conditional probability of each column given its N-1 predecessors using the two-pass algorithm described by Siepel and Haussler (2004). Here, N is the "order" of the model, as defined by --subst-mod; e.g. N=1 for REV, N=2 for U2S, N=3 for U3S. 🌀

The alternative (which is the default) is simply to work with joint probabilities of tuples of columns. (You can make sure those tuples are non-overlapping with the --non-overlapping option.) Using joint probabilities during parameter estimation lets you use the --EM option and can be much faster; plus it appears to produce nearly equivalent estimates. If you want, you can estimate parameters without --markov, and then evaluate the likelihood using --lnl and --markov together. This gives you a lower bound on the likelihood of the Markov-dependent model. No cap.

  • --non-overlapping, -V — (for use with context-dependent substitution models; not compatible with --markov, --features, or --msa-format SS) avoid using overlapping tuples of sites in parameter estimation. If a dinucleotide model is selected, every other tuple will be considered; if a nucleotide triplet model is selected, every third tuple will be considered. This option cannot be used with an alignment represented only by unordered sufficient statistics. Sorry bestie, rules are rules. 😭

(Option for lineage-specific models)

  • --label-branches branch1,branch2,branch3...:label — create a group of branches by giving a set of branches a single label 🏷️. The label should be a word that doesn't contain special characters (in particular: no spaces, brackets, parentheses, pound signs, commas, or colons). The label is for use with the --alt-model option below, so an alternate model can be defined for a set of branches. A branch is specified by the name of the node which is a descendant of that branch.

For example, --label-branches hg18,chimp,hg18-chimp:HC will apply the label "HC" to the hg18, chimp, and hg18-chimp branches in the following tree:

(((hg18,chimp)hg18-chimp, (mouse,rat)mouse-rat)

The same label could be defined without using --label-branches by specifying the tree (either on the command-line or within init-model) like this:

(((hg18 # HC, chimp #HC)#HC, (mouse,rat))
  • --label-subtree node[+]:label — similar to --label-branches, except it labels the entire subtree of the named node 🌳. If the node name is followed by a "+" sign, then it includes the branch leading up to the node in the subtree.
  • --alt-model, -d <label:(model|param_list)> — create a lineage-specific substitution model on a group of branches. The group is defined by a label, which can be specified inside the tree string (using the # sign notation), or by using the --label-branches or --label-subtree arguments. If the alt-model applies to only a single branch, labelling isn't necessary and you can just use the name of the node descending from that branch instead. (See --label-branches above for more deets on labelling groups of branches. 💖)

If you give the name of a substitution model (HKY85, REV, UNREST, etc.) after the colon, then that model will be used for the group of branches, and parameters relevant to the model will be estimated separately in this group. This model may be different (or the same) as the model used in the rest of the tree, but it must have the same number of states and be of the same order as the model used for the rest of the tree. Non-negotiable. 🔥

Alternately, you can give a list of parameter names after the colon. In that case, the same substitution model is used for the entire tree, but the parameters you listed will be estimated separately in the specified group of branches.

The parameter names are model-specific, and include "kappa" for HKY models, "alpha" for GC models, "ratematrix" to specify all rate-matrix parameters in general models, and "backgd" for the equilibrium background frequencies. The parameter names may optionally be followed by boundaries in square brackets to specify parameter bounds, as described in the --bound option. 🎀

The --alt-model argument may be used multiple times, if you want to (for example) optimize a parameter independently on several different groups of branches. Example:

phyloFit align.fa --subst-mod HKY85 \
 --tree "(human, (mouse#MR, rat#MR)#MR, cow)"\
 --alt-model "MR:kappa[0, 1]"

will estimate the HKY85 parameter kappa separately on the mouse/rat subtree, and constrain kappa between 0 and 1. The quotes are often necessary to prevent the square brackets from being parsed by the shell. 🐚 The same model could be achieved with:

phyloFit align.fa --subst-mod HKY85 \
 --tree "(human, (mouse,rat)mouse-rat, cow)"\
 --label-branches mouse,rat,mouse-rat:MR \
 --alt-model "MR:kappa[0,1]"

or

✨ (Options for posterior probabilities) ✨

…and here's the last lil' bit of that example, in case you're following along at home 🫶:

phyloFit align.fa --subst-mod HKY85 \ --tree "(human,
(mouse,rat)mouse-rat, cow)" \ --label-subtree "mouse-rat+:MR" \
--alt-model "MR:kappa[0,1]"

Okay so now we're in the posterior probabilities corner of the options, bestie. This whole cluster is the "the model is fit, now tell me what actually happened at the ancestral nodes" era. 💖

  • --post-probs, -P — outputs the posterior probabilities of all bases at all ancestral nodes. Output goes to auxiliary file(s) with suffix ".postprobs". That's a file suffix with the quotes literally in it, iconic behavior, we love her 💅
  • --expected-subs, -X — outputs the posterior expected number of substitutions on each branch at each site, summed across all types of substitutions. Output goes to auxiliary file(s) with suffix ".expsub".
  • --expected-subs-col, -J — the deluxe version, no summing across substitution types here: posterior expected number of substitutions of each type on each branch, for each site. Output goes to auxiliary file(s) with suffix .expcolsub. (Yeah, no quotes on that suffix. She's different. Periodt.)
  • --expected-total-subs, -Z — posterior expected number of substitutions of each type on each branch, but summed across all sites instead. Output goes to auxiliary file(s) with suffix ".exptotsub".
  • --column-probs, -U — 🎀 heads up, this one is (for use with -init-model; implies --lnl). It outputs a separate log probability for each type of column in the input. Output goes to a file with suffix ".colprobs". And the values? They're log base 2, not base e. Write that down bestie, it matters.

✨ (Options for estimation in sliding window) ✨

Now for the sliding-window girls, the ones who can't commit to just one tree 😭

  • --windows, -w <size,shift> — applies a sliding window to the alignment and fits a separate tree to each window. The arguments specify the size of the window and the amount by which it shifts on each iteration, both counted in bases of the first sequence in the alignment — which is assumed to be the reference sequence, okay, that's the whole point. And because one set of outputs was never enough, separate versions of all output files get created for each window. Yes, all of them. Good luck bestie 💀
  • --windows-explicit, -v <window_coord_list> — same energy as --windows, except every start and end coordinate must be explicitly specified</> by you. Each successive pair of numbers is interpreted as defining the start and end of a window. It can be used with a two-column file and the '*' operator, e.g., --windows-explicit '*mycoords'. Fancy ✨

✨ SEE ALSO ✨

If you want the receipts (and bestie, you do), go read these:

  1. A. Siepel and D. Haussler. 2004. Phylogenetic estimation of context-dependent substitution rates by maximum likelihood. Mol. Biol. Evol., 21:468-488.
  2. Z. Yang, N. Goldman, and A. Friday. 1994. Comparison of models for nucleotide substution used in maximum likelihood phylogenetic estimation. Mol. Biol. Evol., 11:316-324.
  3. Z. Yang. 1994. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J. Mol. Evol., 39:306-314.

And that's the whole page, bestie, you ate 🍓 You just read phyloFit 1.4, dated May 2016, and this gorgeous thing is the manual page PHYLOFIT(1) 👑 Now go fit a tree and text me how it went 💖✨🦋

📦 where this page comes from