PHYLOP(1)

🎀 phyloP(1) — User Commands 🎀

okay bestie, this is the manual for phyloP, and it's literally about computing conservation or acceleration p-values. let's get into it 💅

---

✨ NAME ✨

phyloP — compute conservation or acceleration p-values based on an alignment. The phylogenetic model must be in the .mod format produced by the phyloFit program. The alignment file can be in any of several file formats (see --msa-format). No alignment is required with the --null option.

---

🧬 DESCRIPTION 🧬

Compute conservation or acceleration p-values based on an alignment and a model of neutral evolution. Will also compute p-values of conservation/acceleration in a subtree and in its complementary supertree given the whole tree (see --subtree). P-values can be produced for entire input alignments (the default), pre-specified intervals within an alignment (see --features), or individual sites (see --wig-scores and --base-by-base).

The default behavior is to compute a null distribution for the total number of substitutions from the tree model, an estimate of the number of substitutions that have actually occurred, and the p-value of this estimate wrt the null distribution. These computations are performed as described by Siepel, Pollard, and Haussler (2006). In addition to the SPH method, phyloP can compute p-values or conservation/acceleration scores using a likelihood ratio test (--method LRT), a score-based test (--method SCORE), or a procedure similar to that used by GERP (Cooper et al., 2005) (--method GERP). These alternative methods are currently supported only with --base-by-base, --wig-scores, or --features.

The main advantage of the SPH method is that it can provide a complete and exact description of distributions over numbers of substitutions. However, simulation experiments suggest that the LRT and SCORE methods have somewhat better power than SPH for identifying selection, especially when the expected number of substitutions is small (e.g., with short branch lengths and/or short intervals/individual sites). These two methods are also faster. They are generally similar to one another in power, but in many cases SCORE is considerably faster than LRT. On the other hand, SCORE appears to have slightly less power than LRT at low false positive rates, i.e., for cases of extreme selection. Thus, when using --base-by-base, --wig-scores, or --features, LRT is recommended for most purposes, but SCORE is a good alternative if speed is an issue. When computing p-values with the SPH method, the default is to use the posterior expected number of substitutions as an estimate of the actual number. This is a conservative estimate, because it is biased toward the mean of the null distribution by the prior. These p-values can be made less conservative with --fit-model and more conservative with --confidence-interval (see below).

---

🌸 EXAMPLE 🌸

1. Using the SPH method, compute and report p-values of conservation and acceleration for a given alignment with respect to a neutral model of evolution. Estimated numbers of substitutions are also reported.

phyloP neutral.mod alignment.fa > report.txt

The file neutral.mod could be produced by running phyloFit on data from ancestral repeats or fourfold degenerate sites with an appropriate tree topology and substitution model.

2. Compute and report p-values of conservation and acceleration for a particular subtree of interest (using SPH).

phyloP --subtree human-mouse_lemur neutral.mod alignment.fa > report.txt

Here human-mouse_lemur denote the most recent common ancestor of human and mouse_lemur, which is the node that defines the primate clade in this phylogeny. The tree_doctor program with the --name-ancestors option can be used to assign names to ancestral nodes of the tree.

3. Describe the complete null distribution over the number of substitutions for a 10bp alignment given the specified neutral model (using SPH).

phyloP --null 10 neutral.mod > null.txt

A two-column table is produced with numbers of substitutions and their probabilities, up to an appropriate upper limit.

4. Describe the complete posterior distribution over the number of substitutions in a given alignment (using SPH).

phyloP --posterior neutral.mod alignment.fa > posterior.txt

5. Compute conservation scores (-log10 p-values) for each site in an alignment and output them in the fixed-step wig format (see http://genome.ucsc.edu/goldenPath/help/wiggle.html). Use the likelihood ratio test (LRT) method.

phyloP --wig-scores --method LRT neutral.mod alignment.fa > scores.wig

The --mode option can be used instead to produce acceleration scores (ACC), scores of nonneutrality (NNEUT), or scores that summarize conservation and acceleration (CONACC). The --base-by-base option can be used to output additional statistics of interest (estimated scale factors, log10 likelihood ratios, etc.). As discussed above, several arguments to --method are possible.

6. Similarly, compute scores describing lineage-specific conservation in primates.

phyloP --wig-scores --method LRT --subtree human-mouse_lemur neutral.mod alignment.fa > scores.wig

7. Compute conservation p-values and associated statistics for each element in a BED file. This time use a score test and allow for acceleration as well as conservation, flagging elements under acceleration by making their p-values negative (CONACC mode).

phyloP --features elements.bed --method SCORE --mode CONACC neutral.mod alignment.fa > element-scores.txt

This option can also be used with --subtree. The --gff-scores option can be used to output the original features in GFF format with scores equal to -log10 p. Note that the input file can be in GFF instead of BED format.

✨ OPTIONS ✨

Okay bestie, this is the main event. Every single flag, no skips, we're doing this 💅

• --msa-format, -i FASTA|PHYLIP|MPM|MAF|SS — tells it what format your alignment is in. If you don't say, it'll just sniff the file contents and guess for you (default behaviour).
• --method, -m SPH|LRT|SCORE|GERP — which method computes your p-values or your conservation/acceleration scores. Default is SPH. Here's the tea:
  • The likelihood ratio test (LRT) and the score test (SCORE) compare an alternative model that has a free scale parameter against the given neutral model — or, if you're using --subtree, an alternative model with free scale parameters for the supertree and the subtree against a null model with a single free scale parameter.
  • Those p-values come from comparing test statistics against asymptotic chi-square null distributions.
  • The GERP-like method (GERP) estimates the number of "rejected substitutions" per base by comparing the (per-site) maximum likelihood expected number of substitutions against the expected number under the neutral model.
  • > Heads up: right now LRT, SCORE, and GERP only work with --base-by-base, --wig-scores, or --features. No other combos, sorry babe 😭
• --wig-scores, -w — compute separate p-values per site, then turn those into site-specific conservation (acceleration) scores as -log10(p). Output is base-by-base scores in fixed-step wig format, using the coordinate system of the reference sequence (peep --refidx). In GERP mode, you get rejected substitutions per site instead of -log10 p-values. Cute format, love that for us ✨
• --base-by-base, -b — literally like --wig-scores, but it spits out multiple values per site, and exactly what depends on your method:
  • SPH → mean and variance of the posterior distribution
  • LRT and SCORE → the estimated scale factor(s) and test statistics
  • GERP → the estimated numbers of neutral, observed, and rejected substitutions, plus the number of species available at each site
• --refidx, -r <refseq_idx> — (for use with --wig-scores or --base-by-base) use the coordinate frame of the specified sequence in the output. Default is 1, which is the first sequence in the alignment; 0 means the coordinate frame of the entire multiple alignment.
• --mode, -o CON|ACC|NNEUT|CONACC — (for use with --wig-scores, --base-by-base, or --features) decides the flavour of your p-values:
  • one-sided so that small p (big -log10 p) means unexpected conservation → CON (this is the default 💖) or acceleration → ACC
  • two-sided so that small p means an unexpected departure from neutrality → NNEUT
  • the fourth one, CONACC, uses positive values (p-values or scores) to flag conservation and negative values to flag acceleration.
  • > In GERP mode specifically: CON and CONACC both report the number of rejected substitutions R (which, yes, may be negative), ACC reports -R, and NNEUT reports abs(R).
• --features, -f <file> — read features from <file> (GFF or BED format) and output a table of p-values and related stats, one row per feature. The features are assumed to use the coordinate frame of the first sequence in the alignment. Not for use with --null or --posterior. Also see --gff-scores below 👀
• --gff-scores, -g — (for use with features) instead of a table, output a GFF and give each feature a score equal to its -log10 p-value.
• --subtree, -s <node-name> — (not available in GERP mode) partition the tree into the subtree beneath the node whose name you gave, plus the complementary supertree, and look at conservation/acceleration in the subtree given the supertree. The branch above the specified node goes with the subtree. So like, given the tree "((human,chimp)primate,(mouse,rat)rodent)", the option "--subtree primate" makes one partition of human, chimp, and the branch leading to them, and another partition of the rest of the tree. Meanwhile "--subtree human" makes one partition of just human and the branch leading to it, and another partition of the rest of the tree. In SPH mode, a reversible substitution model is assumed.
• --branch, -B <node-name(s)> — (not available in GERP or SPH mode) like subtree, but it partitions the tree into the set of named branches (each named by its child node) and all the remaining branches. Then it tests for conservation/acceleration in the set of named branches relative to the others. The argument is a comma-delimited list of child nodes.
• --chrom, -N <name> — (optionally use with --wig-scores or --base-by-base) chromosome name for wig output. Default is the root of the multiple alignment filename.
• --log, -l <fname> — write a log to <fname> describing the details of parameter optimization. Super handy for debugging.

Warning: may produce a large file. Like, huge. Watch your disk bestie 🫣

• --seed, -d <seed> — give it a random number seed, should be an integer >=1. Random numbers get used in some cases to generate starting values for optimization. If you don't specify one, it'll use a seed based on the current time.
• --no-prune, -P — don't prune species from the tree that aren't in the alignment. Instead, treat those species as having missing data in the alignment. And fr fr, missing data does affect the results when --method SPH is used.
• --help, -h — produce this help message.

---

🌸 Options for SPH mode only 🌸

• --null, -n <nsites> — compute just the null (prior) distribution of the number of substitutions, as defined by the tree model and the given number of sites, and output it as a table. The 'alignment' argument gets ignored entirely. If you use it with --subtree, you instead get the joint distribution over the number of substitutions in the specified supertree and subtree.
• --posterior, -p — compute just the posterior distribution of the number of substitutions, given the alignment and the model, and output it as a table. If you use it with --subtree, you instead get the joint distribution over the number of substitutions in the specified supertree and subtree.
• --fit-model, -F — fit the model to the data before computing the posterior distribution, by estimating a scale factor for the whole tree, or (if --subtree) separate scale factors for the specified subtree and supertree. This makes p-values less conservative. It has no effect with --null and currently cannot be used with --features. It can be used with --wig-scores and --base-by-base.
• --epsilon, -e <val> — (default 1e-10, or 1e-6 if --wig-scores or --base-by-base) the threshold used when truncating tails of distributions; tail probabilities smaller than this value get discarded. If you want accurate p-values smaller than 1e-10, you'll need this option, and it'll cost you some speed. Note that truncation only affects the right tails, not the left tails, so it should only be an issue with p-values of acceleration.
• --confidence-interval, -c <val> — allow for uncertainty in the estimate of the actual number of substitutions by using a (central) confidence interval about the mean of the specified size (0 < val < 1). To stay conservative, the maximum of this interval is used when computing a p-value of conservation, and the minimum is used when computing a p-value of acceleration. The variance of the posterior is computed exactly, but the confidence interval is based on the assumption that the combined distribution will be approximately normal (which is true for large numbers of sites by the central limit theorem). Math ate that 💅
• --quantiles, -q — (for use with --null or --posterior) report the quantiles of the distribution rather than the whole distribution.

---

👀 SEE ALSO 👀

• Cooper GM, Stone EA, Asimenos G, NISC Comparative Sequencing Program, Green ED, Batzoglou S, Sidow A. Distribution and intensity of constraint in mammalian genomic sequence. Genome Res. 2005 15(7):901-13.
• Siepel A, Pollard KS, and Haussler D. New methods for detecting lineage-specific selection. In Proceedings of the 10th International Conference on Research in Computational Molecular Biology (RECOMB 2006), pp. 190-205.

---

okay that's literally the whole page, you made it 🎀 this was phyloP version 1.4, dated May 2016, and yes it lives at section PHYLOP(1) — go be the main character of your own comparative genomics era, bestie 💖🦋✨