Provided by: phast_1.5+dfsg-1_amd64 
NAME
phastCons - Identify conserved elements or produce conservation scores, given
SYNOPSIS
The alignment file can be in any of several file formats (see --msa-format). The
phylogenetic models must be in the .mod format produced by the phyloFit program.
DESCRIPTION
Identify conserved elements or produce conservation scores, given a multiple alignment and
a phylo-HMM. By default, a phylo-HMM consisting of two states is assumed: a "conserved"
state and a "non-conserved" state. Separate phylogenetic models can be specified for
these two states, e.g.,
phastCons myfile.ss cons.mod,noncons.mod > scores.wig
or a single model can be given for the non-conserved state, e.g.,
phastCons myfile.ss --rho 0.5 noncons.mod > scores.wig
in which case the model for the conserved state will be obtained by multiplying all branch
lengths by the scaling parameter rho (0 < rho < 1). If the --rho option is not used, rho
will be set to its default value of 0.3.
By default, the phylogenetic models will be left unaltered, but if the --estimate-trees
option is used, e.g.,
phastCons myfile.ss init.mod --estimate-trees newtree > scores.wig
then the phylogenetic models for the two states will be estimated from the data, and the
given tree model (there must be only one in this case) will be used for initialization
only. It is also possible to estimate only the scale factor --rho, using the
--estimate-rho option. The transition probabilities for the HMM can either be specified
at the command line or estimated from the data using an EM algorithm. To specify them at
the command line, use either the --transitions option or the --target-coverage and
--expected-length options. The recommended method is to use --target-coverage and
--expected-length, e.g.,
phastCons --target-coverage 0.25 --expected-length 12 myfile.ss
cons.mod,noncons.mod > scores.wig
The program produces two main types of output.
The primary output, sent to stdout in fixed-step WIG format
(http://genome.ucsc.edu/goldenPath/help/wiggle.html), is 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. The scores are
reported in the coordinate frame of a designated reference sequence (see --refidx), which
is by default the first sequence in the alignment. They can be suppressed with the
--no-post-probs option. The secondary type of output, activated with the --most-conserved
(aka --viterbi) option, is a set of discrete conserved elements. These elements are
output 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 the --score
option.
Other uses are also supported, but will not be described in detail here. For example, it
is 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 is
possible to produce a "coding potential" score instead of a conservation score (see
--coding-potential). It is 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 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 used with --estimate-trees or --estimate-rho, the specified value
will be used for initialization only (rho will be estimated). This option is
ignored if two tree models are given.
--estimate-trees, -T <fname_root> Estimate free parameters of tree models and write new
models to <fname_root>.cons.mod and <fname_root>.noncons.mod.
--estimate-rho, -O <fname_root>
Like --estimate-trees, but estimate only the parameter rho.
--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, for 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 first character of
argument is '~', estimate parameters, but initialize to 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 transition parameters such that the
expected fraction of sites in conserved elements is <gamma> (0 < <gamma> < 1).
This is a *prior* rather than *posterior* expectation and 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 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 format based on
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 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 total log likelihood using the forward algorithm and write to 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 log of optimization
procedure to specified file.
--refidx, -r <refseq_idx> Use coordinate frame of specified sequence in output. Default
value is 1, first sequence in alignment; 0 indicates coordinate frame of entire
multiple alignment.
--seqname, -N <name> (Optionally use with --viterbi) Use specified string for 'seqname'
(GFF) or 'chrom' field in output file. Default is obtained from input file name
(double filename root, e.g., "chr22" if input file is "chr22.35.ss").
--idpref, -P <name>
(Optionally use with --viterbi) Use specified string as prefix of generated ids in
output file. Can be used to ensure ids are unique. Default is obtained from input
file name (single filename root, e.g., "chr22.35" if input file is "chr22.35.ss").
--quiet, -q Proceed quietly (without updates to stderr).
--help, -h
Print this help message. (Indels) [experimental]
--indels, -I
Expand HMM state space to model indels as described in Siepel & Haussler (2004).
--max-micro-indel, -Y <length> (Optionally use with --indels) Maximum length of an
alignment gap to be considered a "micro-indel" and therefore addressed by the indel
model. Gaps longer than this threshold will be treated as missing data. Default
value is 20.
--indel-params, -D [~]<alpha_0,beta_0,tau_0,alpha_1,beta_1,tau_1>
(Optionally use with --indels and default two-state HMM) Fix the indel parameters
at (alpha_0, beta_0, tau_0) for the conserved state and at (alpha_1, beta_1, tau_1)
for the non-conserved state, rather than estimating them by maximum likelihood.
Alternatively, if first character of argument is '~', estimate parameters, but
initialize with specified values. Alpha_j is the rate of insertion events per
substitution per site in state j (typically ~0.05), beta_j is the rate of deletion
events per substitution per site in state j (typically ~0.05), and tau_j is
approximately the inverse of the expected indel length in state j (typically
0.2-0.5).
--indels-only, -J Like --indels but force the use of a single-state HMM. This option
allows the effect of the indel model in isolation to be observed. Implies
--no-post-probs. Use with --lnl. (Felsenstein/Churchill model) [rarely used]
--FC, -X
(Alternative to --hmm; specify only one *.mod file with this option) Use an HMM
with a state for every rate category in the given phylogenetic model, and
transition probabilities defined by an autocorrelation parameter lambda (as
described by Felsenstein and Churchill, 1996). A rate constant for each state
(rate category) will be multiplied by the branch lengths of the phylogenetic model,
to create a "scaled" version of the model for that state. If the phylogenetic
model was estimated using Yang's discrete gamma method (-k option to phyloFit),
then the rate constants will be defined according to the estimated shape parameter
'alpha', as described by Yang (1994). Otherwise, a nonparameteric model of rate
variation must have been used (-K option to phyloFit), and the rate constants will
be as defined (explicitly) in the *.mod file. By default, the parameter lambda
will be estimated by maximum likelihood (see --lambda).
--lambda, -l [~]<lambda>
(Optionally use with --FC) Fix lambda at the specified value rather than estimating
it by maximum likelihood. Alternatively, if first character is '~', estimate but
initialize at 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. Thus, lambda = 0 implies no
autocorrelation and lambda = 1 implies perfect autocorrelation. (Coding potential)
[experimental]
--coding-potential, -p
Use parameter settings that cause output to be interpretable as a coding potential
score. By default, a simplified version of exoniphy's phylo-HMM is used, 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 and a set of
default *.mod files (all of which can be overridden). This option can be used with
or without --indels.
--extrapolate, -e <phylog.nh> | default
Extrapolate to a larger set of species based on the given phylogeny
(Newick-format). The trees in the given tree models (*.mod files) must be subtrees
of the larger phylogeny. For each tree model M, a copy will be created of the
larger phylogeny, then scaled such that the total branch length of the subtree
corresponding to M's tree equals the total branch length of M's tree; this new
version will then be used in place of M's tree. (Any species name present in this
tree but not in the data will be ignored.) If the string "default" is given
instead of a filename, then a phylogeny for 25 vertebrate species, estimated from
sequence data for Target 1 (CFTR) of the NISC Comparative Sequencing Program
(Thomas et al., 2003), will be assumed.
--alias, -A <alias_def>
Alias names in input alignment according to given definition, e.g., "hg17=human;
mm5=mouse; rn3=rat". Useful with default *.mod files, e.g., with
--coding-potential. (Default models use generic common names such as "human",
"mouse", and "rat". This option allows a mapping to be established between the
leaves of trees in these files and the sequences of an alignment that uses an
alternative naming convention.)
Custom HMMs [rarely used]
--hmm, -H <hmm_fname>
Name of HMM file explicitly defining the probabilities of all state transitions.
States in the file must correspond in number and order to phylogenetic models in
<mod_fname_list>. Expected file format is as produced by 'hmm_train.'
--catmap, -c <fname>|<string> (Optionally use with --hmm) Mapping of feature types to
category numbers. 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>
States of interest in the phylo-HMM, specified by number (indexing starts with 0),
or if --catmap, by category name. Default value is 1. Choosing --states "0,1,2"
will cause output of the sum of the posterior probabilities for states 0, 1, and 2,
and/or of regions in which the Viterbi path coincides with (any of) states 0, 1, or
2 (see --viterbi).
--reflect-strand, -U <pivot_states>
(Optionally use with --hmm) Given an HMM describing the forward strand, create a
larger HMM that allows for features on both strands by "reflecting" the original
HMM about the specified "pivot" states. The new HMM will be 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> Require "informative" columns (i.e., columns with more
than two non-missing-data characters, excluding sequences specified by
--not-informative) in specified HMM states, to help eliminate false positive
predictions. States can be specified by number (indexing starts with 0) or, if
--catmap is used, by category name. Non-informative columns will be given 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 completely.
--not-informative, -F <list>
Do not consider the specified sequences (listed by name) when deciding whether a
column is informative. This option may be useful when sequences are present that
are very close to the reference sequence and thus do not contribute much in the way
of phylogenetic information. E.g., one might use "--not-informative chimp" with a
human-referenced multiple alignment including chimp sequence, to avoid
false-positive predictions based only on human/chimp alignments (can be a problem,
e.g., with --coding-potential).
--ignore-missing, -z
(For use when estimating transition probabilities) Ignore regions of missing data
in all sequences but the reference sequence (excluding sequences specified by
--not-informative) when estimating transition probabilities. Can help avoid
too-low estimates of <mu> and <nu> or too-high estimates of <lambda>. Warning:
this option should not be used with --viterbi because coordinates in output will be
unrecognizable.
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.