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 and produce a neutral substitution models and trees (one per 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, optionally GC-corrected per chromosome.
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 model and tree (.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.

The reference genome must be hard masked

For scalability, the pipeline splits the alignment into chunks wherever the reference assembly has a run of Ns (see Filtering parameters), and the resulting chunks are scored for conservation. If your assembly is highly contiguous and mostly free of Ns, there will be few natural places to split the alignment, chunks will be large, and repetitive sequence will be included in the phastCons scan. Practically, the pipeline will take prohibitively long to run.

If your assembly doesn't contain Ns, we recommend hard-masking it with RepeatMasker. Hard-masking replaces repetitive sequence with N, which both gives the pipeline's Ns-based splitting real breakpoints to work with and excludes those repeat regions from conservation scoring entirely.

RepeatMasker -pa [threads] -species "[species or clade name]" -dir [output directory] [genome.fasta]

This produces [genome.fasta].masked alongside a few report files. Specify this masked assembly as ref_fasta in the pipeline config and make sure it's the exact same assembly/coordinates already used to build your MAF since masking a different assembly version will shift coordinates and break the Ns-based chunking.

Checking runs of Ns

You can also check whether your assembly has meaningful runs of Ns, both before and after masking, with a few quick bash commands.

First, extract the runs of Ns and write their lengths to a file:

awk '/^>/{if(seq)print seq; seq=""; next}{seq=seq $0}END{if(seq)print seq}' [genome.fasta] | grep -oE 'N+' | awk '{print length($0)}' > n-run-lengths.txt

You can check how many run there are in total:

wc -l < n-run-lengths.txt

Display a distribution of run lengths:

echo "# digits in length of run: 1 = 1-9bp, 2 = 10-99bp, 3 = 100-999bp, etc."; printf "%-8s %s\n" digits count; awk '{print length($1)}' n-run-lengths.txt | sort -n | uniq -c | awk '{printf "%-8s %s\n", $2, $1}'

And also explicitly display how many are at least 100bp, the pipeline's default min_Ns_to_split_by, and so the actual split points it would find:

awk '$1 >= 100' n-run-lengths.txt | wc -l

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, the reference FASTA/GFF. The workflow expects the IDs listed in ref_chromosome_groups to match those in ref_fasta and 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/FASTA index.

For example, if the MAF's src field looks like Homo_sapiens.chr1, and the GFF/FASTA index 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/FASTA index just call 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.

Optional GC content correction

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. The models can be corrected by adjusting for genome-wide GC content.

Set use_gc_corrected_models: true and provide a sample_file to have the pipeline correct each chromosome's neutral model for the GC content of the genomes in the alignment.

If accessions are provided, the pipeline uses ncbi-datasets-cli to look up the GC content. If GC values exist in the gc column of the sample_file, 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.

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, min_Ns_to_split_by, min_keep_region_len, 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:

• 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.
• Alignment chunk splitting (min_Ns_to_split_by, default 100; min_keep_region_len, default 6): the alignment is split wherever the reference has a run of at least min_Ns_to_split_by Ns, and any resulting chunk shorter than min_keep_region_len bp is discarded before scoring.
• 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 — a proxy for chunks with too little real alignment to score meaningfully.
• 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, FASTA index, 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.