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 phastCons against the neutral model to call conserved regions.
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).

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.

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

We provide a fully commented template, config-template.yaml. Copy it from the link above or the internal path below from your local copy of the repository and start filling in the paths and chromosome groups for your project:

cp config-template.yaml my-config.yaml

Open my-config.yaml in an editor. The required inputs from above are all near the top of the file, under a section marked YOU MUST FILL THESE IN, and map onto 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_chr_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 workflow expects the IDs listed in ref_chromosome_groups to match those in ref_gff exactly, and it derives the expected MAF src label for each chromosome from three settings:

• maf_ref_id: the reference species name as it appears in the MAF.
• maf_ref_chr_joiner: the character joining the reference ID and the chromosome name in the MAF (usually ".").
• maf_chr_prefix: an optional prefix on the chromosome name in the MAF that isn't present in the GFF.

For example, if the MAF's src field looks like Homo_sapiens.chr1, and the GFF also call that chromosome chr1, you'd set:

maf_ref_id: "Homo_sapiens"
maf_chr_prefix: ""
maf_ref_chr_joiner: "."

But if the MAF instead labels it Homo_sapiens.chr1 while the GFF just calls it 1, you'd set:

maf_ref_id: "Homo_sapiens"
maf_chr_prefix: "chr"
maf_ref_chr_joiner: "."

and list "1" (not "chr1") under ref_chromosome_groups.

Splitting the alignment into chunks

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.

Note - these defaults aren't a tuned recommendation

The num_seqs_* defaults are a starting point from exploratory analysis on one real alignment, not a tuned recommendation. Review the number and size of chunks they produce on your own data before trusting them.

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.

GC content correction of neutral models

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.

Estimating rho, or using a global value

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.

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:

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.