PhastCons HOWTO

ABSTRACT

PhastCons is part of the PHAST (PHylogenetic Analysis with Space/Time models) package, which can be downloaded from the PHAST Website. The latest version of the phastCons help page can be obtained by running the program with the --help option.

CONTENTS

1. Introduction
2. Getting Started
3. State-Transition Parameters
4. Working with Large Data Sets
5. Advanced Topics [NOT YET WRITTEN]
6. References
7. Feedback
8. Credits

PhastCons and the rest of the PHAST package are written in C. The PHAST package has been compiled and installed under various UNIX and linux platforms and under Mac OS X. It should be relatively easy to port the software to other platforms.

2. GETTING STARTED

2.1 Overview of Usage

As input, PhastCons takes a multiple alignment, a phylogenetic model for conserved regions, and a phylogenetic model for nonconserved regions. The model for conserved regions is optional. If it is not given, then this model is defined indirectly by "scaling" the model for nonconserved regions by a factor called rho (see the --rho option to phastCons). The basic idea of the program is to scan along the alignment for regions that are better "explained" by the conserved model than by the nonconserved model; such regions will be output as conserved elements, and the probability that each base is in such a region will be output as the conservation score for that base.

The program also needs a specification of the state-transition probabilities for the HMM. These probabilities determine how freely the model can switch between states, and as a result, they define important features such as the expected length of conserved elements, the expected coverage by conserved elements, and the "smoothness" of the conservation plot (see the phastCons paper for details). The transition probabilities can be set directly with the --transitions option or with the --expected-length and --target-coverage options (see Section 3). Alternatively, they may be left undefined, and will automatically be estimated by maximum likelihood. Probably the trickiest part of getting good results with phastCons is figuring out how to set these parameters. Let's gloss over this potentially complicated issue for now and focus on the basics of using the program. For the moment, assume the transition probabilities are set by the tuning parameters, say, by using --target-coverage 0.25 and --expected-length 12. We'll return to this issue in Section 3.

The program can be run with a command of the form:

    phastCons [OPTIONS] alignment [cons.mod,]noncons.mod > scores.wig

For example,

    phastCons --target-coverage 0.25 --expected-length 12 \
        --msa-format MAF alignment.maf cons.mod,noncons.mod \
        > scores.wig

    phastCons --target-coverage 0.25 --expected-length 12 \
        --rho 0.4 --msa-format MAF alignment.maf noncons.mod \
        > scores.wig

A set of discrete conserved elements can be predicted by adding the --most-conserved option, e.g.,

    phastCons --target-coverage 0.25 --expected-length 12 --rho 0.4 \
        --most-conserved most-cons.bed --msa-format MAF alignment.maf \
        noncons.mod > scores.wig

If you don't care about the conservation scores, you can avoid computing them with the --no-post-probs option. (This can save some time.) If phastCons is run with --no-post-probs and without --most-conserved, then it will neither predict a set of discrete elements nor produce the conservation scores. However, it will still estimate free parameters, if directed to do so (see below).

In summary, the default behavior of the program is to estimate any unspecified free parameters, then to compute conservation scores and output them to stdout. However, it can be made to perform any combination of its three main tasks -- predicting elements, computing scores, and estimating parameters -- by running it with or without --most-conserved, with or without --no-post-probs, and with or without a specification of the values of its free parameters.

2.2 Obtaining Phylogenetic Models

1. Let phastCons itself estimate the models directly from the data, using an "unsupervised" learning algorithm (the EM algorithm described in the phastCons paper).
2. Separately estimate the nonconserved model from labeled training data, e.g., fourfold degenerate (4d) sites or ancestral repeats (ARs), and let phastCons estimate the scaling parameter rho.
3. Separately estimate both the conserved and nonconserved models from labeled training data and don't do any estimation with phastCons.

It is worth noting that in experiments we've done with vertebrate data (including human, mouse, rat, chicken, and Fugu sequences), estimating the nonconserved model from 4d sites rather than letting phastCons estimate them did not make too much of a difference in the predicted conserved elements, even though it did result in quite different branch-length estimates (see Supplementary Material to phastCons paper). The EM algorithm (option 1) seems to be able to learn a separation between conserved and nonconserved regions even if the model it learns for nonconserved regions is not an entirely accurate representation of the neutral substitution process. On the other hand, when the number of species is large (say, greater than a dozen), option 1 can be quite computationally intensive, and options 2 or 3 may be preferable purely for practical reasons (see the ENCODE example, below). As a rule of thumb, I recommend option 1 unless you have (a) more than 12 species or (b) many distant species with low alignment coverage and you are concerned about distortions in branch lengths. I recommend option 2 in these cases.

For options 2 or 3, you can estimate your models using the program phyloFit, also in the PHAST package. Start with an alignment and an annotation file indicating the sites of interest in your alignment. For example, suppose you have an alignment of human, mouse, and rat sequences in a file called data.fa (FASTA format), and suppose you have a file called AR.gff (GFF format) defining the locations of ancestral repeats (ARs) in the coordinate system of the first sequence (human). For example,

human   RepeatMasker    AR      14451829        14451925        .       .       .       type "L1MC4a.LINE.L1"
human   RepeatMasker    AR      14452228        14452443        .       .       .       type "L1MC4a.LINE.L1.2"
human   RepeatMasker    AR      14458153        14458278        .       .       .       type "L2.LINE.L2"
human   RepeatMasker    AR      14458258        14458386        .       .       .       type "L2.LINE.L2.2"
human   RepeatMasker    AR      14458781        14459349        .       .       .       type "L2.LINE.L2.3"
human   RepeatMasker    AR      14460129        14460311        .       .       .       type "L1ME.LINE.L1"
human   RepeatMasker    AR      14460948        14461016        .       .       .       type "L2.LINE.L2.4"
human   RepeatMasker    AR      14463023        14463356        .       .       .       type "MER89-int.LTR.ERV1"
...

        phyloFit --tree "(human,(mouse,rat))" --features AR.gff \
            --do-cats AR --out-root nonconserved data.fa

With 4d sites, the program msa_view may be helpful. Suppose instead of annotations of ARs, you have a file called genes.gff describing coding regions in the human sequence, e.g.,

human    UCSC        CDS     4562    4692    .       -       0       transcript_id "NM_198943" ; gene_id "NM_198943"
human    UCSC        CDS     7469    7605    .       -       0       transcript_id "NM_198943" ; gene_id "NM_198943"
human    UCSC        CDS     7778    7924    .       -       0       transcript_id "NM_198943" ; gene_id "NM_198943"
human    UCSC        CDS     8131    8242    .       -       1       transcript_id "NM_198943" ; gene_id "NM_198943"
human    UCSC        CDS     14601   14749   .       -       0       transcript_id "NM_198943" ; gene_id "NM_198943"
...

    msa_view data.fa --4d --features genes.gff > 4d-codons.ss 

    msa_view 4d-codons.ss --in-format SS --out-format SS \
         --tuple-size 1 > 4d-sites.ss

    phyloFit --tree "(human,(mouse,rat))" --msa-format SS \
        --out-root nonconserved-4d 4d-sites.ss

    for set in set1 set2 set3 ; do 
        msa_view $set.maf --in-format MAF --4d --features $set.gff \
            > $set.codons.ss 
        msa_view $set.codons.ss --in-format SS --out-format SS \
            --tuple-size 1 > $set.sites.ss
    done
    msa_view --aggregate human,mouse,rat *.sites.ss > all-4d.sites.ss
    phyloFit --tree "(human,(mouse,rat))" --msa-format SS \
        --out-root nonconserved-all-4d all-4d.sites.ss

Similar methods could be used to estimate a model for conserved regions, say, from 1st codon positions. This gets slightly more complicated, e.g.:

    msa_view data.maf --in-format MAF --features genes.gff \
        --catmap "NCATS = 3; CDS 1-3" --out-format SS --unordered-ss \
        --reverse-groups transcript_id > codons.ss
    phyloFit --tree "(human,(mouse,rat))" --msa-format SS --do-cats 1 \
        --out-root conserved-codon1 codons.ss

It should now be clear how option 3 would be accomplished: estimate conserved and nonconserved models using appropriate annotations files and phyloFit, then run phastCons as shown in Section 2.1. Option 2 is similar, but requires the use of the --estimate-rho option, e.g.,

    phastCons --target-coverage 0.25 --expected-length 12 \
        --rho 0.4 --estimate-rho mytrees --msa-format MAF \
        alignment.maf nonconserved-4d.mod > scores.wig

Option 1 requires the use of --estimate-trees, which is similar to --estimate-rho:

    phastCons --target-coverage 0.25 --expected-length 12 \
        --estimate-trees mytrees --msa-format MAF alignment.maf \
        init.mod > scores.wig

    phyloFit --tree "(human,(mouse,rat))" --msa-format MAF \
        --out-root init alignment.maf

    phastCons --target-coverage 0.25 --expected-length 12 \
        --estimate-trees mytrees --msa-format MAF alignment.maf \
        init.mod --no-post-probs
    phastCons --target-coverage 0.25 --expected-length 12 \
        --msa-format MAF alignment.maf --most-conserved most-cons.bed \
        mytrees.cons.mod,mytrees.noncons.mod > scores.wig

Note that options 1, 2, and 3 can all be used with or without --target-coverage and --expected-length (see Section 3).

3. STATE-TRANSITION PARAMETERS

3.1 Basics

Here is what is important to understand in order to use the program:

• The probabilities of transitioning between the two states of the HMM are controlled by two parameters: either mu and nu, or if you prefer, gamma and omega. The values of these parameters will influence the coverage by conserved elements, the lengths of conserved elements, and the smoothness of the conservation scores.
• The program allows these parameters to be treated in two ways: as tuning parameters to be set by the user to achieve desired results, or as free parameters to be estimated from the data by maximum likelihood.
• To treat them as tuning parameters, use either the --transitions option (to specify mu and nu), or the --target-coverage and --expected-length options (to specify gamma and omega, respectively).
• To estimate the parameters from the data, simply avoid using the --transitions, --target-coverage, and --expected-length options. The parameters will automatically be estimated if not specified. Alternatively, you can use the --transitions option but precede its arguments with a "~" character (e.g., --transitions ~0.01,0.01); this will cause the parameters to be estimated but initialized to the specified values in the iterative estimation algorithm.
• A hybrid solution is also possible: if you specify --target-coverage only, the program will estimate only the expected-length parameter from the data. In this case, a "~" prefix can be used to specify a starting value for --expected-length.

3.2 Tuning Parameters

The question then becomes what criteria should be used for adjusting the tuning parameters. There are two natural features to focus on. The first, called "coverage," is the fraction of bases in a data set that fall within predicted conserved elements. You might have a reason to aim for a particular level of coverage; for example, with mammalian data, a target of 5% might be selected, since this is considered a reasonable ballpark estimate of the fraction of sites under purifying selection in mammalian genomes. If you don't have a target in mind for the overall level of coverage, you can aim for a particular level of coverage of known functional elements such as protein-coding exons. (This was the method we used in the phastCons paper to calibrate the model in similar ways for vertebrate, insect, worm, and yeast data sets; we chose a target of 65% coverage of coding exons by conserved elements.) The target-coverage parameter (gamma) is a "knob" you can turn to reach a particular level of actual coverage. As you increase this parameter, more bases will be predicted to be in conserved elements, and as you decrease it fewer bases will be predicted to be in conserved elements. (Gamma can be thought of in a Bayesian sense as a prior, which has the effect of imposing a bias on the actual coverage.)

The second feature is called "smoothing." Smoothing describes the degree to which the conservation scores tend to be similar from one base to the next. A high degree of smoothing occurs when the probabilities mu and nu are small -- i.e., when changing states is unlikely -- and a lower degree of smoothing occurs when mu and nu are larger. Higher smoothing means the predicted conserved elements tend to be longer on average, and are more likely to include short sequences of nonconserved bases. Lower smoothing implies shorter, more fragmented conserved elements, containing fewer nonconserved sites. Increasing the expected length parameter (omega) tends to increase the level of smoothing and the average length of predicted conserved elements.

What is the right level of smoothing? This is perhaps a more difficult question than what is the right level of coverage. In many cases, it may be best just to set the smoothing so that the conservation plot "looks right." By plotting the scores next to a multiple alignment, as in the UCSC browser, it is usually not too hard to see if the plot is dramatically "over-smoothed" (in the extreme, it resembles a square wave function, with alternating regions of scores near 0 and scores near 1) or "under-smoothed" (in the extreme, the scores at adjacent sites appear essentially uncorrelated). It may also be useful to look at the shortest predicted conserved elements and see whether they are believable -- are they long enough and conserved enough to reasonably be called "conserved elements"? If not, perhaps the degree of smoothing should be increased. Similarly, one might look at sequences of nonconserved sites within conserved elements. If they are quite long and poorly conserved, yet the elements are not being broken, then perhaps the smoothing should be reduced. One might also look at the breaks between conserved elements (if they are too short, the smoothing should be increased) and sequences of conserved sites that have not been predicted as conserved elements (if they are too long, the smoothing should be decreased).

An alternative, perhaps more principled, way of setting the level of smoothing is to use the "phylogenetic information threshold" (PIT) method described in the phastCons paper. The paper should be consulted for details, but the PIT is essentially a threshold for the amount of "information" (in an information theoretic sense) that is required to predict a conserved element, on average. It is an alternative measure of the degree of smoothing (a higher PIT means more smoothing), which can be meaningfully compared across different data sets and phylogenies. For example, to achieve a given PIT -- say, 10 bits -- a larger expected length might be needed for a two-species data set, for which there is less information about conservation at each site, than for a ten-species data set. (The PIT actually considers the shape of the phylogeny and its branch lengths, not just the number of species.) A program in PHAST called consEntropy can be used to compute the PIT for a given set of parameters. The inputs to the program are the target-coverage and expected-length parameters and the conserved and nonconserved phylogenetic models (see the consEntropy help page and examples below). One can simply run consEntropy after adjusting the tuning parameters and re-estimating any free parameters, and compare the reported PIT to the target value; if necessary, the parameters can be adjusted and the procedure can be repeated (see below). There are cases in which the PIT method will not be appropriate. For example, when an alignment contains a lot of missing data, the PIT will be misleading, because there is actually less information per site than the phylogenetic models would imply. Also, consEntropy will not work with large trees, because, as currently implemented, its running time is exponential in the number of species.

Coverage and smoothing are to some degree separate, with coverage being primarily influenced by the --target-coverage option and smoothing being primarily influenced by the --expected-length option, but they are also interrelated. The final set of predictions depends on all parameters and on the data in a complicated way, and will be affected, sometimes unpredictably, by a change to any parameter. For example, changing the --expected-length option while holding --target-coverage fixed generally will alter coverage somewhat as well as smoothing. For this reason, the --target-coverage and --expected-length parameters should be tuned simultaneously, until constraints on both coverage and smoothing are satisfied.

3.3 An Example

We begin with a guess at the values of the tuning parameters. In this case, a starting target-coverage of 0.05 seems appropriate, since our prior belief is that about 5% of bases are conserved, but we have to consider the fact that only about 40% of the human bases align to mouse or rat, and phastCons will not predict conserved elements in regions of no alignment. Therefore, we start instead with a target coverage of 0.05/0.4 = 0.125, the fraction of aligned bases we expect to be conserved. We'll guess at a value of 20bp for the expected length. The average length of conserved elements is likely to be considerably longer -- probably closer to 100bp -- but we want to influence the model toward the detection of shorter conserved elements, e.g., possible transcription factor binding sites. One way of thinking about it is that we want the conservation score at each base to reflect an average "window" of about 20bp centered at that site. Remember, these are just rough starting values; we're going to adjust them as we go.

A starting estimate of the tree models can be obtained with a command such as:

    phastCons --target-coverage 0.125 --expected-length 20 \
        --estimate-trees v1 --most-conserved most-cons.bed \
        --log optim.log --no-post-probs --msa-format MAF \
        alignment.maf init.mod 

    featureBits -enrichment hg17 CDS.bed most-cons.bed

    CDS.bed 0.159%, most-cons.bed 0.233%, both 0.102%, cover 64.25%, enrich 275.73x

    consEntropy 0.125 20 v1.cons.mod v1.noncons.mod

    Transition parameters: gamma=0.125000, omega=20.000000, mu=0.050000, nu=0.007143
    Relative entropy: H=0.295299 bits/site
    Expected min. length: L_min=49.071119 sites
    Expected max. length: L_max=24.350034 sites
    Phylogenetic information threshold: PIT=L_min*H=14.490669 bits

Now we rerun the program with new values for the tuning parameters. We might try reducing the expected length to 15 and leaving the target-coverage the same, because we don't know how the new expected length will affect the coverage.

    phastCons --target-coverage 0.125 --expected-length 12 \
        --estimate-trees v2 --most-conserved most-cons.bed \
        --log optim.log --no-post-probs --msa-format MAF \
        alignment.maf v1.noncons.mod 

This time we measure a coverage of 58.7% and a PIT of 14.3 bits. Apparently, the reduced smoothing (i.e., lower expected length) has resulted in a greater separation between the conserved and nonconserved models: there are fewer nonconserved sites mixed in with the conserved ones, making the conserved model more conserved, and there are fewer conserved sites mixed in with the nonconserved ones, making the nonconserved model less conserved. As a result, the phylogenetic information per base has increased, and despite the reduced smoothing, the PIT has stayed roughly the same. In addition, the coverage has decreased a bit, possibly because fewer nonconserved sites flanked by conserved sites are being pulled into conserved elements. It seems as though we've lost ground, but we just need to be persistent. Let's try reducing the expected length a bit more and cranking up the target coverage to compensate for it:

    phastCons --target-coverage 0.2 --expected-length 8 \
        --estimate-trees v3 --most-conserved most-cons.bed \
        --log optim.log --no-post-probs --msa-format MAF \
        alignment.maf v2.noncons.mod 

Tuning the program is a bit of an art: the parameters are all interrelated, and it takes a little time to get a feeling for what effect each change will have. As you get better at it, you'll need fewer iterations to reach convergence.

4. WORKING WITH LARGE DATA SETS

The techniques I recommend are, basically, (1) to "divide and conquer" and make as much use as possible of parallel processing, and (2) to make the estimation problem as easy as possible. In this section, I will discuss two different examples of large data sets, one consisting of a moderate number of species and a large number of sites, and the other consisting of a large number of species and a moderate number of sites. A "divide and conquer" strategy is the key to addressing the first problem and making the estimation problem easier is the key to addressing the second problem. Note that the two strategies I will describe are essentially orthogonal, and could be used in combination with a data set having many sequences *and* many sites.

4.1 Many Sites, Not Too Many Sequences

Here's the recommended "divide and conquer" approach for this case:

1. Split the alignments into small fragments (each spanning about 1Mb of the human genome)
2. Estimate parameters for each fragment, separately and in parallel.
3. Combine the separately estimated parameters by averaging.
4. Predict conserved elements and conservation scores globally using the combined estimates
5. Adjust tuning parameters and return to step 2, if necessary. Repeat until coverage and smoothing targets are met.

Let's discuss the steps one by one.

1. Split the alignments into small fragments

   This can be accomplished by running a program called msa_split on each MAF. Use a command like the following:

           mkdir -p CHUNKS            # put fragments here
           for file in chr*.maf ; do
               root=`basename file .maf`
               msa_split $file --in-format MAF --refseq $root.fa \
                   --windows 1000000,0 --out-root CHUNKS/$root --out-format SS \
                   --min-informative 1000 --between-blocks 5000 
           done
   

   This will cause each MAF file ($root.maf) and the corresponding reference sequence ($root.fa), to be read in, assembled in memory into a kind of pseudo-global alignment (with local alignments tiled beneath the human sequence), and then split into non-overlapping chunks of about 1,000,000 sites each; these chunks will be written in SS format to files with a prefix of $root in the CHUNKS directory. The --min-informative option causes fragments with fewer than 1000 "informative" sites (i.e., with at least two aligned bases) to be discarded, and the --between-blocks option tells the program to try to split between alignment blocks if possible, e.g., in recent transposon insertions in the human genome. (It will adjust the split point by at most 5000 sites.) The file $root.fa must be in FASTA format and must use the same coordinate system as the multiple alignment, as will be true if the FASTA files and MAFs from the same version (assembly) of the UCSC Genome Browser are used. See the msa_split help page for additional details.

   If desired, this step can be done in parallel, with one process per MAF file. Beware of the fact that each process will be quite I/O intensive.

2. Estimate parameters for each fragment.

   In this step, we will estimate free parameters using phastCons for each alignment fragment. Because the number of species is relatively small, we can get away with full estimation of the conserved and nonconserved models (option 3 from Section 2.2). Suppose we already have a starting model, init.mod. We might use a command like:

           mkdir -p TREES     # put estimated tree models here
           rm -f TREES/*      # in case old versions left over
           for file in chr*.*.ss ; do 
               root=`basename $file .ss` 
               phastCons --target-coverage 0.125 --expected-length 20 \
                   --gc 0.4 --estimate-trees TREES/$root \
                   $file init.mod --no-post-probs
           done
   

   Here, the --gc option is used to ensure that the same G+C content is assumed for all fragments, since the same G+C content will be assumed during prediction. Note that the --log and --quiet options may be useful here as well.

   If you have a compute cluster at your disposal, you will probably want to execute these individual processes in parallel.

3. Combine the separately estimated parameters by averaging A program in PHAST called phyloBoot, which is primarily designed for bootstrapping, supports averaging of model parameters. It can be used to obtain average conserved and nonconserved models as follows:

           ls TREES/*.cons.mod > cons.txt
           phyloBoot --read-mods '*cons.txt' --output-average ave.cons.mod 
           ls TREES/*.noncons.mod > noncons.txt
           phyloBoot --read-mods '*noncons.txt' --output-average ave.noncons.mod 
   

   PhyloBoot will output a summary of the means, standard deviations, medians, etc. of each parameter. You may want to scan this summary to convince yourself that the standard deviations are not too large and that no outliers have severely skewed the means. The means and medians should be roughly equal and the standard deviations should be fairly small compared to the means. (Some larger standard deviations for branches near the root of the tree are normal.) There are sometimes a few zeroes in the 'min' column, especially for short branches; this is okay, but beware of very large values in the 'max' column, which could result from a failure to converge properly, and might skew the mean. If you see any bad outliers, you should re-estimate with different starting parameters or track down and discard the responsible .mod files. I've never had to do this; the estimates generally are well behaved.

4. Predict conserved elements and conservation scores globally using the combined estimates.

   This is simply a matter of running phastCons on each fragment, using the averaged tree models.

           mkdir -p ELEMENTS SCORES
           rm -f ELEMENTS/* SCORES/*
           for file in chr*.*.ss ; do 
               root=`basename $file .ss` 
               phastCons --target-coverage 0.125 --expected-length 20 \
                   --most-conserved ELEMENTS/$root.bed --score \
                   $file ave.cons.mod,ave.noncons.mod > SCORES/$root.wig
           done
   

   Note that the same --target-coverage and --expected-length arguments are used as in step 2, that this time we have asked the program to predict both discrete elements and base-by-base conservation scores, and that the tree models from step 3 have been used. The --score option tells phastCons to compute log odds scores for predicted elements. The --gc option is no longer necessary (ave.cons.mod and ave.noncons.mod contain the background distribution), but the --seqname and --idpref options could be useful in this step. (The default names and ids are obtained from the filenames, and in this case should be okay).

   The *.bed and *.wig files can now be combined across the whole genome or across each chromosomes, if desired, e.g.,

           cat ELEMENTS/*.bed | sort -k1,1 -k2,2n > most-conserved.bed
   

5. Adjust tuning parameters and return to step 2, if necessary. Repeat until coverage and smoothing targets are met.

   Suppose our targets are 65% coverage of coding exons and a PIT of 10 bits. At UCSC, we would test the first constraint using the featureBits program, e.g.,

           featureBits -enrichment hg17 knownGene:cds most-conserved.bed
   

   The second constraint could be tested with consEntropy, e.g.,

           consEntropy .125 20 ave.cons.mod ave.noncons.mod
   

   As in Section 3, we adjust the tuning parameters as needed, only now we repeat steps 2, 3, and 4 after each adjustment, until convergence.

4.1.1 Random Samples and Local Parameter Estimates

    ls chr*.*.ss | chooseLines -k 100 - > random-sample.txt

    mkdir -p TREES     # put estimated tree models here
    rm -f TREES/*
    for file in `cat random-sample.txt` ; do 
        root=`basename $file .ss` 
        phastCons --target-coverage 0.125 --expected-length 20 \
            --gc 0.4 --estimate-trees TREES/$root \
            $file init.mod --no-post-probs
    done

    phastCons --target-coverage 0.125 --expected-length 20 \
        --most-conserved ELEMENTS/$root.bed --score \
        $file TREES/$root.cons.mod,TREES/$root.noncons.mod \
        > SCORES/$root.wig

4.2 Many Sequences, Not Too Many Sites

In this case, I recommend making the estimation problem easier. The best way to do that here is probably to estimate the nonconserved model from fourfold degenerate sites, as described in Section 2.2, and then to let phastCons estimate only the parameter rho (option 2 from section 2.2). We'll work with an aggregate alignment constructed from the original alignments for the purposes of parameter estimation. Suppose we start with a set of FASTA-formatted files named for ENCODE regions, ENm001.fa, ..., ENr334.fa. We can concatenate these into a single file using msa_view with a command like:

    msa_view --aggregate {species-list} EN*.fa > all-encode.fa

    msa_view all-encode.fa --out-format SS > all-encode.ss

For parameter estimation, we can use a command like:

    phastCons --target-coverage 0.05 --expected-length 20 \
        --gc 0.4 --most-conserved all-encode.bed --no-post-probs \
        --estimate-rho newtree all-encode.ss nonconserved-4d.mod 

    awk '{total += $3 - $2} END {print total/29909540}' all-encode.bed

ConsEntropy cannot be used in this case to test the smoothness constraint, because of the large number of sequences. Instead, I recommend plotting the conservation scores and judging the smoothness by eye, as discussed in Section 3.2.

When the coverage and smoothness constraints have been met, predictions for the individual targets can be obtained with commands such as:

    for reg in EN*.fa ; do
        root=`basename reg .fa`
        phastCons --target-coverage  --expected-length  \
            --most-conserved $root.bed all-encode.ss \
            newtree.cons.mod newtree.noncons.mod > $root.wig
    done