workflow.md

The ARGtest workflow

← Back to README

One reasonable post-processing workflow for ARG tree sequences is below. The Snakemake pipeline automates exactly these steps; each can also be run by hand via the linked script (--help for full arguments).

1. Find low-recombination regions. Regions of low recombination are more affected by linked selection, so for analyses assuming a neutral ARG it can help to remove them first. scripts/hapmap_low_rec_mask.py turns a HapMap recombination map plus a .fai into per-chromosome BED masks of the lowest-recombination intervals (set by rec_fraction; see Required Inputs for the intervals-vs-base-pairs caveat). If your ARG came from singer-snakemake with paircoal-reweight: true, you may be able to skip this step (set rec_fraction to 0).
2. Find regions of poor alignment. scripts/find_low_access_regions.py inspects a tree sequence's inferred mutation map and writes windows where more than a cutoff fraction of base pairs are masked to a BED file.
3. Find aberrant mutational load. All samples should have similar derived-mutation loads after accounting for local mutation-rate variation. scripts/mutload_masks.py simulates mutations once on the input ARG to get a per-window, per-individual simulation-based expected load, then flags any (window, individual) pair whose observed load deviates from expectation by more than mutload_cutoff, writing one BED row per outlier window with the affected individuals in column 4. The test is per individual: a diploid/polyploid individual's haploid copies are pooled before testing, so for haploid data it is effectively per sample. scripts/mutload_summary.py gives an HTML diagnostic of residual load after pruning, flagging candidates for whole-sample removal.
4. Remove affected regions. Combine the step 1–3 BED masks per chromosome and apply them to each tree sequence with scripts/trim_regions_single.py. This is the pipeline's coordinate-preserving regional filter: removed spans become empty intervals rather than shifting later coordinates. The resolved mutation-rate map is embedded in each filtered tree sequence's metadata.
5. Remove affected samples. Often only a few samples in a window are problematic — e.g. introgressed (cf. TRACE) or with odd derived-mutation patterns (step 3). scripts/trim_samples.py prunes given individuals over given intervals. The pipeline automatically prunes the step-3 outliers; add your own with trim_individuals (genome-wide IDs) and/or trim_remove_bed (per-individual BEDs), or run the script directly with --individuals / --remove.
6. Validate (optional). scripts/validation_plots_from_ts.py produces compact QC plots (mutational load, diversity, Tajima's D, folded/unfolded SFS) for the cleaned ARG; compare against the same plots on the originals. Alongside each set of plots it writes the numbers behind them as three tab-separated tables — windows.tsv (per-window), samples.tsv (per-sample) and sfs.tsv (per-frequency-bin) — so any panel can be re-checked or re-plotted without reloading the ARGs. 6b. Pool across chromosomes (optional). scripts/genomewide_expected_vs_observed.py reads every chromosome's windows.tsv and plots expected versus observed π and Tajima's D for all windows at once, colouring each point by the recombination rate of that window taken from the HapMap map. The per-chromosome scatters can only show whether one chromosome is off; pooling with a recombination-rate colour axis shows whether the deviation from the 1:1 line tracks recombination rate, which is the signature of ARG inference degrading in low-recombination regions. Requires validation_sim_branch — without a simulated expectation there is nothing to plot against. It pools whatever step 6 produced, so with the default validation_first_chrom_only: true the panels cover only the first chromosome; set it to false for a genuinely genome-wide result. Both the figures and the HTML summary state the actual chromosome coverage, and flag it prominently when it is partial.
7. Merge chromosomes per replicate. scripts/merge_treefiles_by_replicate.py concatenates the per-chromosome tree files into one tree sequence per replicate, merging the embedded mutation-rate ratemaps across chromosomes.
8. Summarise (optional). scripts/pipeline_summary.py collects genome-retention statistics across all steps (mean ± SD across replicates), per-individual outlier counts, the step-6b genome-wide expected-vs-observed panels, and the per-chromosome validation plots into one self-contained HTML file. Final retained Mb are measured directly from the filtered ARGs as positive-rate accessible sequence with a non-empty local genealogy, so they exclude both regions masked in the input ARG and regions removed by this pipeline. Retention percentages use the full FASTA-index reference length as the denominator.

Optionally, set emit_vcf: true to export a per-(chromosome, replicate) VCF from the filtered per-chromosome tree sequences; see VCF export.

Summary statistics after sample pruning

Diversity (π), Tajima's D, and SFS in this pipeline are computed with tskit's built-in methods (ts.diversity, ts.Tajimas_D, ts.allele_frequency_spectrum). These normalize by a constant n · (n − 1) / 2 based on the sample set passed in, so when the sample size varies across regions — as it does after trim_samples.py removes individuals over BED intervals — the per-window statistics are not correctly normalized for the locally retained sample count. Treat post-pruning stats with caution, and prefer comparisons on replicates where sample membership is uniform across the genome.