A typical set of loci for analysis with PhyloAcc are conserved non-exonic elements (CNEEs). The phyloacc-workflows repository contains Snakemake workflows that take a whole-genome alignment (in MAF format) and a reference genome annotation (in GFF format) and produce a set of neutral substitution models and trees (one per scaffold/chromosome), conserved elements, and a final set of CNEE alignments ready to hand to PhyloAcc (see the README for how PhyloAcc uses these as input).

The pipeline works in three broad stages, each of which can be turned on or off independently in the config file:

1. Neutral model estimation: 4-fold degenerate codons are extracted from the alignment and used to fit a neutral substitution model with phyloFit.
2. Conservation scoring: the alignment is split into manageable chunks, and scored with either phastCons or phyloP against the neutral model to call conserved regions. Users can choose one or both methods, depending on how their data scale both in terms of power and computationally. See below for more details.
3. CNEE extraction: conserved regions overlapping coding sequence (from a GFF) are removed, short remaining fragments are dropped, and the surviving elements are extracted as individual FASTA or MAF alignments.

The pipeline splits the alignment by chromosome/scaffold, both for scalability and to allow for chromosome-specific neutral models. This also means you will end up running PhyloAcc for each chromosome/scaffold separately, rather than on the whole genome at once.

Outputs of the pipeline include the neutral models and trees (.mod files), and the final CNEE alignments (FASTA files).

Conservation scoring method and power

The workflow can use either phastCons or phyloP to score conservation. phastCons aggregates information across sites to directly predict conserved regions, while phyloP scores each site independently and requires a separate clustering step to call conserved regions.

In general, phyloP lacks the power to predict conserved sites on small trees, and therefore phastCons is usually the preferred method (and is the default). However, it can be slow on very large trees, and phyloP+clustering may be the only practical option in those cases.

If running phyloP, a power calculation will be performed prior to execution to determine whether the tree is too small to reliably detect conserved sites. If so, the workflow will exit with an error and suggest switching to phastCons, or overriding and running phyloP anyway, which may still result in no conserved sites being detected.

One or both of phyloP and phastCons can be run in the same workflow, and the config file allows you to specify which method to use with the main workflow switches (run_phylop and run_phastcons; see below).

phyloP Power Calculations (click to expand)

Per-site phyloP scores each base for conservation with the following formula: \[ p = \operatorname{erfc}\!\left(\sqrt{\mathrm{LRT}/2}\right), \qquad \mathrm{LRT} = 2\left[\ln L(\hat\rho) - \ln L(1)\right] \] where \(\mathrm{LRT}\) measures how much better a slower-than-neutral (conserved) model fits the base than the neutral model: a larger \(\mathrm{LRT}\) gives a smaller \(p\), i.e. stronger evidence of conservation, and \(\operatorname{erfc}\) is the complementary error function. phyloP reports these as scores, \(-\log_{10} p\).

These scores have a ceiling (\(\mathrm{LRT}_{\max}\)), which is determined mainly by the total neutral tree depth (substitutions/site) and the number of species in the tree, with a smaller contribution from the substitution rate itself (a model whose fastest base turns over more quickly than average raises the ceiling). Concretely, the largest \(\mathrm{LRT}\) a site can reach, which sets that ceiling, is capped by whichever of two limits is smaller: \[ \mathrm{LRT}_{\max} \approx \min\!\Big(\underbrace{2\,T\,\max_b(-Q_{bb})}_{\text{depth}\,\times\,\text{rate}},\ \ \underbrace{2(n-1)\big(-\ln \min_b \pi_b\big)}_{\text{species}\,\times\,\text{composition}}\Big) \] The first term grows with the tree depth \(T\) and the fastest base's exit rate \(\max_b(-Q_{bb})\); the second with the number of species \(n\) and the rarest base's frequency \(\min_b\pi_b\). The ceiling follows whichever term is smaller, so too little depth, or too few species, each cap it.

Then, since millions of sites are being tested, we must correct the resulting scores for multiple tests. This results in a multiple-testing threshold that is set by the number of sites being tested and the desired false-positive rate: \[ \text{threshold} = \log_{10}\!\left(M/\alpha\right) \] where \(M\) is the number of sites tested and \(\alpha\) the desired false-positive rate — the more sites tested, the higher the bar.

If phyloP's score ceiling falls below the multiple-testing threshold, no site is detectable, regardless of how conserved it is.

In other words, conserved sites are detectable only when: \[ -\log_{10}\operatorname{erfc}\!\left(\sqrt{\mathrm{LRT}_{\max}/2}\right)\;\ge\;\log_{10}(M/\alpha) \] This is the condition under which the tree's best-possible score clears the correction for \(M\) tested sites at level \(\alpha\). The ceiling rises with tree depth but eventually saturates, so two things can leave phyloP powerless: too little total tree depth, or too few species. Deep, taxon-rich trees are comfortably detectable; shallow or few-taxon trees need element-based conservation instead.

1. Clone the repository

Use the following command to clone the repository:

git clone https://github.com/phyloacc/phyloacc-workflows.git

If you do not have or wish to use git, download the archive directly from GitHub or with the following command:

wget https://github.com/phyloacc/phyloacc-workflows/archive/main.zip

2. Ensure conda is installed

The workflow uses conda to manage dependencies. You can check if conda is installed with conda --version. If you don't have conda (i.e. if conda command returns a "command not found" error), check out our tutorial to install it:

3. Set up the conda environment

The repository includes a small wrapper script, phyloacc_workflows, that manages a dedicated conda environment for you with all the required dependencies. To create it, run:

./phyloacc_workflows setup

This creates a conda environment named phyloacc-workflows. You can confirm the environment is ready at any time with:

./phyloacc_workflows check

If you ever change envs/environment.yml or pull an update that changes it, re-running ./phyloacc_workflows setup will update the existing environment rather than recreate it from scratch.

The wrapper has a few other subcommands too (run ./phyloacc_workflows --help for the full list).

If ./phyloacc_workflows: Permission denied shows up when you try to run it, make the script executable and try again:

chmod +x phyloacc_workflows

1. Required inputs

This workflow works in two steps: 1) predict neutral models (one per chromosome) and then 2) using the neutral models, predict conserved elements. These steps require the following inputs:

The paths to these files and other pipeline options are specified in a single YAML config file, described in the pipeline configuration section.

2. Pipeline config file

Everything the workflow needs, including the paths to the raw inputs above is specified in a single YAML config file. YAML is a format that works by pairing keys and values as key: value pairs. The keys are provided and represent specific settings the workflow needs, and you fill in the values.

Config template

There are two ways to get a starting config file. Either works; pick whichever suits you.

Quick start (recommended): generate a config with the required fields left blank:

./phyloacc_workflows init -o my-config.yaml

Full reference: alternatively, copy the fully-commented template if you'd rather have every option documented inline as you fill it in (also available at config-template.yaml). Copy it from the link above or the internal path below from your local copy of the repository:

cp config-template.yaml my-config.yaml

Open my-config.yaml in an editor and fill in the following config keys, along with a few other required settings:

There are many other settings in the config file that are commented within it. In the following sections we highlight a few that are important to understand to know if you need to adjust them for your dataset.

Matching chromosome IDs

Relevant config keys: ref_chromosome_groups, maf_ref_id, maf_ref_chr_joiner, maf_prefix, gff_prefix

A common source of early errors is that the reference chromosome/scaffold IDs don't line up between the MAF and the reference GFF.

The easiest solution to this is for the user to match the IDs prior to running the pipeline, most likely by editing the GFF file.

However, the pipeline itself provides a method to handle this automatically, as long as the labels share a core ID (e.g. "chr1" and "1", with "1" being the core), by specifying a prefix for the MAF and GFF chromosome names. The settings for this are confusing, which is why pre-editing the IDs is recommended. Click below if you'd like to read more about how the prefix system works.

maf_ref_id: "Homo_sapiens"
maf_ref_chr_joiner: "."
maf_prefix: "chr"
gff_prefix: ""

maf_ref_id: "Homo_sapiens"
maf_ref_chr_joiner: "."
maf_prefix: ""
gff_prefix: "chr"

GC content correction of neutral models (phyloFit)

Relevant config keys: use_gc_corrected_models, sample_file

Because the neutral models are estimated from 4-fold degenerate sites and subsequently applied to the whole genome, if those sites have different GC content the models may be inaccurate. By default, the models are corrected for this by computing GC content directly from the MAF.

Alternatively, a sample_file can be provided for GC correction. If a column called accessions exists in the sample file, the pipeline uses ncbi-datasets-cli to look up the GC content. If pre-computed GC values exist in a gc column, those values are used instead. Both columns can exist and different samples can use different methods to provide GC content.

With the GC content read, the pipeline uses PHAST's mod_freqs script to adjust neutral models for each chromosome.

Set use_gc_corrected_models: false to disable GC correction.

Splitting the alignment into chunks (phastCons)

Relevant config keys: split_strategy, num_seqs_max_for_gap, num_seqs_min_gap_bp, num_seqs_min_keep_region_len

Before scoring conservation, each chromosome's alignment is split into smaller chunks, both for scalability, and so that long stretches with little or no real alignment data don't get scored at all. split_strategy controls how the pipeline decides where to split.

By default (split_strategy: num_seqs), the pipeline looks directly at the MAF's own alignment blocks: any stretch where too few species are aligned is treated as a gap and used as a split point. Three settings control split behavior:

• num_seqs_max_for_gap (default 3): a MAF block with this many species or fewer counts toward a "gap" run.
• num_seqs_min_gap_bp (default 100): the minimum contiguous length of such a low-coverage run to actually count as a split point. Short blocks with low coverage may not be true split points.
• num_seqs_min_keep_region_len (default 6): chunks shorter than this, after splitting, are dropped.

Other split strategies include ns (split by Ns in a provided reference genome) and fixed_windows (split into fixed-size windows). These are documented in the config template.

Estimating rho, or using a global value (phastCons)

Relevant config keys: rho_mode, fixed_rho, global_rho_stat

phastCons needs a single "rho" parameter describing how conserved the alignment is overall relative to the neutral model, and the pipeline applies one such value per chromosome to every chunk it scores. By default (rho_mode: fixed), that's simply the value you set for fixed_rho (default 0.3), which we've found to be a reasonable value for typical vertebrate datasets.

Alternatively, set rho_mode: estimate to instead have phastCons estimate rho separately for each alignment chunk, then summarize those per-chunk estimates into a single chromosome-wide value using global_rho_stat (p90 (value of the 90th percentile of chunk estimates) by default, or median/mean of the chunk estimates). Any chunk whose own estimated rho exceeds that chromosome-wide value is skipped for conservation calling, rather than scored with an inflated rho.

Conserved site prediction and clustering (phyloP)

Relevant config keys: run_phylop, phylop_alpha, phylop_power_gate, phylop_power_override, phylop_power_num_sites, phylop_power_fallback_num_sites, phylop_cluster_method, hmm_t0_0, hmm_t1_1, hmm_e0_0, hmm_e1_1, hmm_s0, hmm_min_len, hmm_max_len, naive_merge_gap_bp, naive_min_region_sites, naive_min_region_len_bp, windowed_window_bp, windowed_min_sites_per_window

The default pipline mode predicts conserved elements with phastCons, however this may be intractable for large alignments. In that case, you can use phyloP to predict conserved sites and then cluster them into elements.

To enable this, set run_phylop: true to additionally call conserved and accelerated sites with phyloP, which tests each site in the alignment individually against the neutral model, rather than scoring alignment chunks the way phastCons does (see conservation scoring method above). Significant sites are then clustered into regions and fed into the same CNEE-building stage as phastCons, producing a second, independent set of CNEEs.

Calling significant sites. A site is called conserved or accelerated if it passes an FDR threshold, phylop_alpha (default 0.05), after Benjamini-Hochberg correction across all sites on the chromosome.

Clustering sites into regions. Once significant sites are called, they need to be grouped into contiguous regions before they can be turned into CNEEs. phylop_cluster_method controls how:

• hmm (default): a simple 2-state online HMM (in/out of a conserved element) run over the per-position conserved/not-conserved stream.
• gap_merge: Merge significant sites that are close together, then drop regions that are too small or too short.
• windowed: bin the chromosome into fixed-size windows; a window is called conserved if it has enough significant sites, and adjacent conserved windows are merged.

Each method has its own set of tuning parameters, which are documented in the full table below and in the config template.

Filtering parameters

Relevant config keys: filter_threshold_4d, max_gap_pct, cnee_ces_merge_gap_bp, cnee_min_len_bp

Several thresholds control how aggressively data is filtered at different stages of the pipeline (see also splitting the alignment into chunks for the settings that control where chunk boundaries themselves are placed):

• 4-fold degenerate sites (filter_threshold_4d, default 0.5): sites used to fit the neutral model are dropped if more than this fraction of sequences are missing at that site.
• Chunk quality (max_gap_pct, default 0.9): after splitting, a chunk is dropped entirely if more than this fraction of its non-reference alignment columns are gaps.
• Final CNEE filtering (cnee_ces_merge_gap_bp, default 5; cnee_min_len_bp, default 50): conserved regions within cnee_ces_merge_gap_bp bp of each other (after coding sequence is removed) are merged into a single element, and anything shorter than cnee_min_len_bp bp afterward is dropped.

The defaults are reasonable for typical vertebrate-scale alignments, but you may want to loosen them for smaller or more divergent datasets, or tighten them for very large ones.

Specifying resources

Relevant config keys: rule_resources

At the bottom of the config is a list of per-rule cluster resources, which the workflow passes to Snakemake when submitting jobs. The required resources depend on the number of species in the alignment and the size of the genomes. Many rules are fast and light and will use the default resouces. Others can be slow and memory-intensive and have their own resource settings.

Values in the template are based on a benchmark of a 15 species alignment of mammals. The config notes which rules should scale with genome size and which with sample size. If you run out of memory or time on a rule, increase the resources for that rule in your config file and re-run the workflow.

Full config reference

Every key recognized by the config file, in the order it appears in config-template.yaml:

The workflow is executed through the same phyloacc_workflows wrapper used for setup. Its run subcommand activates the conda environment and passes everything you give it straight through to snakemake, defaulting to the Snakefile in the repository unless you specify your own with -s.

1. Dry run

Always start with a dry run to make sure the config file is valid and to see what jobs Snakemake plans to run, before anything is actually submitted or executed:

./phyloacc_workflows run --configfile my-config.yaml -j 20 -e slurm --dryrun

Here, -j 20 is the maximum number of jobs Snakemake will have in flight at once, and -e slurm tells Snakemake to submit jobs to a SLURM cluster using the resources you set per rule under rule_resources in your config file. If you're testing on a single machine instead of a cluster, drop -e slurm and Snakemake will run everything locally using up to -j CPU cores.

Here is an example rulegraph for fitting neutral models and extracting CNEEs

2. Executing the workflow

Once the dry run looks right, drop --dryrun to actually run it:

./phyloacc_workflows run --configfile my-config.yaml -j 20 -e slurm

Depending on the size of your alignment and how many chromosomes/scaffolds you're analyzing, this can take anywhere from minutes to many hours. Snakemake will print progress as jobs are submitted and complete, and each rule also writes its own log under <output_dir>/logs/<rule name>/ for closer inspection.

3. Re-running and troubleshooting

Snakemake only re-runs rules whose outputs are missing or out of date, so if a run is interrupted or errors out, address the cause of the failure and then run the exact same command again and it will pick up where it left off rather than starting over.

If a particular rule keeps failing, check its log file first, both under <output_dir>/logs/<rule name>/ and (for cluster runs) in the SLURM job's own output. Common early culprits are a chromosome/scaffold ID that doesn't match between the MAF and GFF (see Matching chromosome IDs), or a cluster partition/resource in rule_resources that doesn't exist on your system.

All outputs are written under the output_dir you set in your config file, organized into numbered subdirectories reflecting the stage of the pipeline that produced them. The ones you'll care about most are:

If run_phylop is also enabled, CNEE building runs a second time on phyloP's clustered regions, writing an equivalent, independent set of outputs under 05-cnees/phylop/ instead of 05-cnees/phastcons/. With both stages enabled you end up with two CNEE sets, one per method.

From here, the CNEE alignment directory for a chromosome (or all of them pooled together) is ready to hand straight to phyloacc.py along with the neutral model produced above. See the PhyloAcc README for how to set up and run PhyloAcc itself on these inputs.