PHASTCONS(1)

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

🌸 PHASTCONS — User Commands 🌸

okay bestie, this is the phastCons manual page, part 1 of 2, and we are getting through it together. hydrate. stretch. let's go 💅

✨ NAME ✨

phastCons — identify conserved elements or produce conservation scores, given…

and yes, the original sentence literally just stops at "given" 😭 the manpage did that, not me. we move.

🦋 SYNOPSIS 🦋

Your alignment file can be in any of several file formats (peep --msa-format for the full menu 👀).

The phylogenetic models must be in the .mod format — the one the phyloFit program produces. that's the only flavor on the shelf, no cap.

💖 DESCRIPTION 💖

so the whole job: identify conserved elements or produce conservation scores, given a multiple alignment and a phylo-HMM.

By default we're working with a phylo-HMM that has two states: a "conserved" state and a "non-conserved" state. You can hand over separate phylogenetic models for those two states, like:

phastCons myfile.ss cons.mod,noncons.mod > scores.wig

or you can hand over just one model, for the non-conserved state:

phastCons myfile.ss --rho 0.5 noncons.mod > scores.wig

and then the conserved-state model gets built by multiplying all the branch lengths by the scaling parameter rho (0 < rho < 1). If you don't pass --rho, rho just sits at its default value of 0.3. no notes. 🎀

By default the phylogenetic models are left unaltered — hands off, we don't touch. But if you use the --estimate-trees option, e.g.

phastCons myfile.ss init.mod --estimate-trees newtree > scores.wig

then the phylogenetic models for the two states get estimated from the data, and the tree model you gave (there must be only one in this case, bestie, don't be greedy 🙈) is used for initialization only. You can also estimate just the scale factor --rho, using --estimate-rho.

The transition probabilities for the HMM? Two paths: you specify them at the command line, or the program estimates them from the data using an EM algorithm. Command-line girlies use either --transitions, or --target-coverage together with --expected-length. The recommended method is --target-coverage + --expected-length, like:

phastCons --target-coverage 0.25 --expected-length 12 myfile.ss cons.mod,noncons.mod > scores.wig

The program produces two main types of output.

(not a typo, that is genuinely the heading, she's a whole sentence 💀)

The primary output is sent to stdout in fixed-step WIG format (http://genome.ucsc.edu/goldenPath/help/wiggle.html), and it's a set of base-by-base conservation scores. The score at each base is equal to the posterior probability that that base was "generated" by the conserved state of the phylo-HMM. Scores get reported in the coordinate frame of a designated reference sequence (see --refidx), which by default is the first sequence in the alignment. You can suppress them with --no-post-probs.

The secondary type of output — activated with --most-conserved (aka --viterbi) — is a set of discrete conserved elements. Those elements come out in either BED or GFF format, also in the coordinate system of the reference sequence (see --most-conserved). They can be assigned log-odds scores using --score. ✨

Other uses are also supported, but they're not getting the full breakdown here. For example, it's possible to produce conservation scores and conserved elements using a k-state phylo-HMM of the kind described by Felsenstein and Churchill (1996) (see --FC), and it's possible to produce a "coding potential" score instead of a conservation score (see --coding-potential). It's also possible to give the program a custom HMM and to specify any subset of its states to use for prediction (see --hmm and --states).

See the phastCons HOWTO for additional details. 🫶

🎀 EXAMPLE 🎀

1. Given phylogenetic models for conserved and nonconserved regions and HMM transition parameters, compute a set of conservation scores.

phastCons --transitions 0.01,0.01 mydata.ss cons.mod,noncons.mod > scores.wig

2. Similar to (1), but define the conserved model as a scaled version of the nonconserved model, with rho=0.4 as the scaling parameter. Also predict conserved elements as well as conservation scores, and assign log-odds scores to predictions.

phastCons --transitions 0.01,0.01 --most-conserved mostcons.bed --score --rho 0.4 mydata.ss noncons.mod > scores.wig

(if the output file were "mostcons.gff," then output would be in GFF instead of BED format 👑

3. This time, estimate the parameter rho from the data. Suppress both the scores and the conserved elements. Specify the transition probabilities using --target-coverage and --expected-length instead of --transitions.

phastCons --target-coverage 0.25 --expected-length 12 --estimate-rho newtree --no-post-probs mydata.ss noncons.mod

4. This time estimate all free parameters of the tree models.

phastCons --target-coverage 0.25 --expected-length 12 --estimate-trees newtree --no-post-probs mydata.ss noncons.mod

5. Estimate the state-transition parameters but not the tree models. Output the conservation scores but not the conserved elements.

phastCons mydata.ss cons.mod,noncons.mod > scores.wig

6. Estimate just the expected-length parameter and also estimate rho.

phastCons --target-coverage 0.25 --estimate-rho newtree mydata.ss noncons.mod > scores.wig

✨ OPTIONS ✨

🌳 Tree models

  • --rho, -R <rho> — set the scale (overall evolutionary rate) of the model for the conserved state to be <rho> times that of the model for the non-conserved state (0 < <rho> < 1; default 0.3). If you use it with --estimate-trees or --estimate-rho, your value is only for initialization (rho will be estimated anyway). This option is ignored if two tree models are given. 🙈
  • --estimate-trees, -T <fname_root> — estimate the free parameters of the tree models and write the new models to <fname_root>.cons.mod and <fname_root>.noncons.mod.
  • --estimate-rho, -O <fname_root> — like --estimate-trees, but it only estimates the parameter rho. minimal effort, honestly iconic.
  • --gc, -G <val> — (optionally use with --estimate-trees or --estimate-rho) assume a background nucleotide distribution consistent with the given average G+C content (0 < <val> < 1) when estimating tree models. (The frequencies of G and C will be set to <val>/2 and the frequencies of A and T will be set to (1-<val>)/2.) This option overrides the default behavior of estimating the background distribution from the data (if --estimate-trees) or obtaining them from the input model (if --estimate-rho).
  • --nrates, -k <nrates> | <nrates_conserved,nrates_nonconserved> — (optionally use with a discrete-gamma model and --estimate-trees) assume the specified number of rate categories, instead of the number given in the *.mod file. The shape parameter 'alpha' will be as given in the *.mod file. In the case of the default two-state HMM, two values can be specified — the numbers of rates for the conserved and the nonconserved states, resp.

🔄 State-transition parameters

  • --transitions, -t [~]<mu>,<nu> — fix the transition probabilities of the two-state HMM as specified, rather than estimating them by maximum likelihood. Alternatively, if the first character of the argument is '~', estimate the parameters but initialize them to the specified values. The argument <mu> is the probability of transitioning from the conserved to the non-conserved state, and <nu> is the probability of the reverse transition. The probabilities of self transitions are thus 1-<mu> and 1-<nu>, and the expected lengths of conserved and nonconserved elements are 1/<mu> and 1/<nu>, respectively.
  • --target-coverage, -C <gamma> — (alternative to --transitions) constrain the transition parameters such that the expected fraction of sites in conserved elements is <gamma> (0 < <gamma> < 1). Heads up bestie: this is a prior rather than posterior expectation, and it assumes stationarity of the state-transition process. Adding this constraint causes the ratio mu/nu to be fixed at (1-<gamma>)/<gamma>. If used with --expected-length, the transition probabilities will be completely fixed; otherwise the expected-length parameter <omega> will be estimated by maximum likelihood.
  • --expected-length, -E [~]<omega> — (--expected-lengths also allowed, for backward compatibility) (for use with --target-coverage, alternative to --transitions) set the transition probabilities such that the expected length of a conserved element is <omega>. Specifically, the parameter mu is set to 1/<omega>. If preceded by '~', <omega> will be estimated, but will be initialized to the specified value.

📁 Input/output

  • --msa-format, -i PHYLIP|FASTA|MPM|SS|MAF — alignment file format. Default is to guess the format based on the file contents. Note that the msa_view program can be used to convert between formats.
  • --viterbi (alternatively --most-conserved), -V <fname> — predict discrete elements using the Viterbi algorithm and write to the specified file. Output is in BED format, unless <fname> has suffix ".gff", in which case output is in GFF.
  • --score, -s — (optionally use with --viterbi) assign a log-odds score to each prediction.
  • --lnl, -L <fname> — compute the total log likelihood using the forward algorithm and write to the specified file.
  • --no-post-probs, -n — suppress output of posterior probabilities. Useful if only discrete elements or likelihood is of interest.
  • --log, -g <log_fname> — (optionally use when estimating free parameters) write a log of the optimization procedure to the specified file.
  • --refidx, -r <refseq_idx> — use the coordinate frame of the specified sequence in the output. Default value is 1, the first sequence in the alignment; 0 indicates the coordinate frame of the entire multiple alignment.
  • --seqname, -N <name> — (optionally use with --viterbi) use the specified string for the 'seqname' (GFF) or 'chrom' field in the output file. Default is obtained from the input file name (double filename root, e.g. "chr22" if the input file is "chr22.35.ss").
  • --idpref, -P <name> — (optionally use with --viterbi) use the specified string as a prefix of the generated ids in the output file. Can be used to ensure ids are unique. Default is obtained from the input file name (single filename root, e.g. "chr22.35" if the input file is "chr22.35.ss").
  • --quiet, -q — proceed quietly (without updates to stderr). 🤫
  • --help, -h — print this help message. (and then the original page has a little (Indels) [experimental] floating right here with no home 🫠 that's just how the source file is laid out, I'm reporting it faithfully)
  • --indels, -I — expand the HMM state space to model indels as described in Siepel & Haussler (2004).
  • --max-micro-indel, -Y <length> — (pair it with --indels if you want) The longest a gap in the alignment can be and still count as a "micro-indel" — aka something the indel model will actually deal with. Anything longer than this threshold gets treated as missing data, periodt. Default is 20.
  • --indel-params, -D [~]<alpha_0,beta_0,tau_0,alpha_1,beta_1,tau_1> — (works with --indels and the default two-state HMM) Locks the indel parameters in place instead of letting maximum likelihood estimate them: (alpha_0, beta_0, tau_0) for the conserved state and (alpha_1, beta_1, tau_1) for the non-conserved state. Or, plot twist: if the first character of the argument is ~, the parameters still get estimated, but they start out at the values you gave. Here's the lore:
    • alpha_j is the rate of insertion events per substitution per site in state j — typically like ~0.05
    • beta_j is the rate of deletion events per substitution per site in state j — also typically ~0.05
    • tau_j is approximately the inverse of the expected indel length in state j — typically 0.2-0.5
  • --indels-only, -J — Basically --indels but it forces a single-state HMM. This lets you stare directly at the indel model working all by itself, no other vibes in the way. It implies --no-post-probs. Use it with --lnl. (Felsenstein/Churchill model) [rarely used]
  • --FC, -X — (this is an alternative to --hmm; with this option you specify only one *.mod file) Uses an HMM that has a state for every rate category in the given phylogenetic model, with transition probabilities defined by an autocorrelation parameter lambda (as described by Felsenstein and Churchill, 1996). A rate constant for each state (rate category) gets multiplied by the branch lengths of the phylogenetic model, which creates a "scaled" version of the model for that state. If the phylogenetic model was estimated using Yang's discrete gamma method (the -k option to phyloFit), the rate constants get defined according to the estimated shape parameter 'alpha', as described by Yang (1994). Otherwise, a nonparametric model of rate variation must have been used (the -K option to phyloFit), and the rate constants will be exactly as defined (explicitly) in the *.mod file. By default, the parameter lambda is estimated by maximum likelihood (see --lambda).
  • --lambda, -l [~]<lambda> — (optionally use with --FC) Pins lambda at the value you give instead of estimating it by maximum likelihood. And if the first character is ~, it gets estimated but initialized at your specified value. Allowable range is 0-1. With k rate categories, the transition probability between state i and state j will be lambda * I(i == j) + (1-lambda)/k, where I is the indicator function. So lambda = 0 means no autocorrelation at all, and lambda = 1 means perfect autocorrelation. (Coding potential) [experimental]
  • --coding-potential, -p — Parameter settings that make the output readable as a coding potential score. Obsessed. By default it uses a simplified version of exoniphy's phylo-HMM, with a noncoding (background) state, a conserved non-coding (CNS) state, and states for the three codon positions. This option implies:
--catmap "NCATS=4; CNS 1; CDS 2-4" --hmm <default-HMM-file> --states CDS --reflect-strand background,CNS

plus a set of default *.mod files — and all of those can be overridden. It can be used with or without --indels.

  • --extrapolate, -e <phylog.nh> | default — Extrapolates out to a bigger set of species based on the phylogeny you hand it (Newick-format). The trees in the given tree models (*.mod files) must be subtrees of that larger phylogeny. Here's how it works: for each tree model M, a copy of the larger phylogeny gets created, then scaled so that the total branch length of the subtree matching M's tree equals the total branch length of M's tree — and that new version gets used in place of M's tree. (Any species name that shows up in this tree but not in the data is just ignored, no stress.) If you give the string "default" instead of a filename, then it assumes a phylogeny for 25 vertebrate species, estimated from sequence data for Target 1 (CFTR) of the NISC Comparative Sequencing Program (Thomas et al., 2003).
  • --alias, -A <alias_def> — Renames things in the input alignment according to your definition, e.g. "hg17=human; mm5=mouse; rn3=rat". Super handy with the default *.mod files, like with --coding-potential. (The default models use generic common names such as "human", "mouse", and "rat". This option lets you build a mapping between the leaves of the trees in those files and the sequences of an alignment that's using a whole different naming convention. 💅)

🎀 Custom HMMs [rarely used] 🎀

  • --hmm, -H <hmm_fname> — Names the HMM file that explicitly defines the probabilities of all the state transitions. The states in the file must correspond in number and order to the phylogenetic models in <mod_fname_list>. Expected file format is whatever 'hmm_train.' produces.
  • --catmap, -c <fname>|<string> — (optionally use with --hmm) Maps feature types to category numbers. You can give either a filename or an "inline" description of a simple category map, e.g. --catmap "NCATS = 3 ; CDS 1-3".
  • --states, -S <state_list> — The states of interest in the phylo-HMM, specified by number (indexing starts with 0) or, if --catmap is in play, by category name. Default value is 1. Choosing --states "0,1,2" makes the output the sum of the posterior probabilities for states 0, 1, and 2, and/or the regions where the Viterbi path coincides with (any of) states 0, 1, or 2 (see --viterbi).
  • --reflect-strand, -U <pivot_states> — (optionally use with --hmm) If you've got an HMM describing the forward strand, this creates a bigger HMM that allows for features on both strands by "reflecting" the original HMM about the specified "pivot" states. That new HMM then gets used for prediction on both strands. States can be specified by number (indexing starts with 0), or if --catmap, by category name.

🌸 Missing data [rarely used] 🌸

  • --require-informative, -M <states> — Requires "informative" columns in the specified HMM states, to help eliminate false positive predictions. (An informative column is one with more than two non-missing-data characters, not counting the sequences specified by --not-informative.) States can be specified by number (indexing starts with 0) or, if --catmap is used, by category name. Non-informative columns get emission probabilities of zero. By default this option is active, with <states> equal to the set of states of interest for prediction (as specified by --states). Use "none" to disable it completely.
  • --not-informative, -F <list> — Don't consider the sequences you list (by name) when deciding whether a column is informative. This is a lifesaver when you've got sequences that are super close to the reference and therefore aren't really adding phylogenetic info. E.g., you might use "--not-informative chimp" with a human-referenced multiple alignment that includes chimp sequence, so you dodge false-positive predictions based only on human/chimp alignments (which can be a problem, e.g., with --coding-potential).
  • --ignore-missing, -z — (use this when estimating transition probabilities) Ignores regions of missing data in all sequences but the reference sequence (excluding sequences specified by --not-informative) when estimating transition probabilities. Can help you avoid estimates of <mu> and <nu> that come out too low, or estimates of <lambda> that come out too high.
⚠️ Warning: do not use this option with --viterbi, because the coordinates in the output will be unrecognizable. Like, genuinely unreadable. Not the vibe.

👑 SEE ALSO 👑

  • J. Felsenstein and G. Churchill. 1996. A hidden Markov model approach to variation among sites in rate of evolution. Mol. Biol. Evol., 13:93-104.
  • A. Siepel, G. Bejerano, J. S. Pedersen, et al. 2005. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. (in press)
  • A. Siepel and D. Haussler. 2004. Computational identification of evolutionarily conserved exons. Proc. 8th Annual Int'l Conf. on Research in Computational Biology (RECOMB '04), pp. 177-186.
  • J. Thomas et al. 2003. Comparative analyses of multi-species sequences from targeted genomic regions. Nature 424:788-793.
  • 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 made it to the end, iconic behavior fr fr. This one's phastCons version 1.4, dated May 2016, and it lives at PHASTCONS(1) — go forth and find some conserved elements ✨🦋💖

📦 where this page comes from