1. Configuration
   1. YAML config file
   2. Wrapper Snakefile
   3. Resource management
   4. Reference assemblies
   5. Experiments (differential methylation)
   6. Conda environments
   7. Analysis module setup
      1. External tool repositories → repos.* config keys
      2. Analysis reference files → haplotype.* and deconv.* config keys
2. emseq.smk — Core Processing
   1. Wrapper variables
   2. Processing steps
   3. Outputs
3. emseq_analysis.smk — Downstream Analysis
   1. Wrapper variables
   2. Processing steps
   3. Outputs
4. Testing
   1. Test wrappers
   2. Test data limitations
5. Continuous Integration
6. Change Log

The EM-seq repository provides two modular Snakemake workflows for enzymatic methylation sequencing data. Both are designed to be included from a project-specific wrapper Snakefile that defines samples, references, and resource limits.

Configuration is split between a YAML config file and a wrapper Snakefile. The YAML holds all declarative settings: sample lists, reference genome definitions, tool parameters, environment paths, and resource limits. The wrapper Snakefile reads the YAML and assigns every variable that the modules consume — environments, directories, sample lists, tool parameters, and project-level decisions like which references to align against and which spike-ins were included. The wrapper also defines the rule all output targets and any project-specific custom rules (e.g. input symlinks, sample aliases).

This separation serves two purposes. First, the YAML is portable across machines and projects while the wrapper encodes the specific run configuration. Second, when a project composes multiple pipeline modules (e.g. EM-seq + cfDNA CNA + fragmentation analysis), the wrapper is the single place where all variable names from all modules are visible. This makes namespace collisions immediately obvious rather than hidden inside separate module files.

Passed via --configfile. See config/test.yaml for a complete working example.

Core keys:

Analysis-specific keys (required only for emseq_analysis.smk; see Analysis module setup below for how to prepare these):

The wrapper reads the YAML config and assigns every variable that the included modules expect. This is deliberate: by centralizing all variable assignments in one file, you can see at a glance every name in scope across all modules and catch collisions before they become runtime bugs. The wrapper also makes project-level decisions that do not belong in the YAML — which reference genomes to align against and which spike-in controls were included — and defines the rule all output targets.

A minimal wrapper follows this structure:

import os
configfile: "config/my_project.yaml"

def resolve_config_paths(config_dict):
    for k, v in config_dict.items():
        if isinstance(v, str):
            config_dict[k] = os.path.expandvars(os.path.expanduser(v))
        elif isinstance(v, dict):
            resolve_config_paths(v)
        elif isinstance(v, list):
            config_dict[k] = [os.path.expandvars(os.path.expanduser(i))
                              if isinstance(i, str) else i for i in v]
resolve_config_paths(config)

# --- Environments ---
ENV_EMSEQ = config['envs']['emseq']
ENV_METHYLKIT = config['envs']['methylkit']

# --- Repositories ---
R_EMSEQ = config['repos']['emseq']

# --- Data directories ---
D_DATA = config['main-data-dir']
D_EMSEQ = f"{D_DATA}/emseq"
D_REF = f"{D_DATA}/ref"
D_LOGS = f"{D_DATA}/logs"
D_BENCHMARK = f"{D_DATA}/benchmark"
D_INPUTS = f"{D_DATA}/inputs"

# --- Tool parameters ---
MOSDEPTH_QUANT_LEVELS = config.get("mosdepth-quant-levels", "1,5,10,20")
EMSEQ_MINCOV = config.get("emseq-mincov", 2)
FASTP_EXTRA = config.get("fastp", {}).get("extra", "")
EMSEQ_REF_INPUTS = {k: v['input'] for k, v in config['emseq_ref_assemblies'].items()}

# --- Samples and references ---
emseq_library_ids = config["library-ids"]
emseq_ref_names = ["ncbi_decoy_hg38"]   # project decision
emseq_align_methods = ["bwa_meth"]       # project decision
spike_builds = ["puc19", "unmeth_lambda"]  # project decision
KEEP_BED = config["keep-bed"]
EXCL_BED = config["exclude-bed"]
meth_map = config["meth-map"]

rule all:
    input:
        expand(f"{D_EMSEQ}/dmr/tabix/{{lib}}.{{ref}}.bwa_meth.methyldackel.txt.bgz",
               lib=emseq_library_ids, ref=emseq_ref_names),
        # ... additional targets

include: "emseq.smk"
include: "emseq_analysis.smk"   # optional

The test wrappers (workflows/test.smk, workflows/test-analysis.smk) and config/example-config.yaml serve as templates for creating project-specific wrappers.

Resource management uses two mechanisms: threads and concurrency.

Threads is a built-in Snakemake feature. Some rules use threads: workflow.cores to consume all available cores (e.g. alignment), while others use a hardcoded thread count for I/O-bound operations where more cores don't help (e.g. threads: 8 for samtools sort, dedup, BAM filtering). Rules that use workflow.cores are automatically capped at the --cores flag value, so the same workflow runs correctly on any machine.

Concurrency is a custom resource pattern using Snakemake's resources: directive. Each rule declares a concurrency cost, and a global budget is set at runtime via --resources concurrency=N. Snakemake schedules jobs as long as the total concurrency cost of running jobs stays within the budget. The system is calibrated so that a rule with concurrency = 100 runs one at a time when --resources concurrency=100. On a larger machine, setting --resources concurrency=300 allows three such jobs in parallel. Lighter-weight rules declare lower costs — for example, a rule with concurrency = 10 allows 10 simultaneous instances at --resources concurrency=100, or 30 at --resources concurrency=300. This ratio-based approach means the same rule definitions work across machines of different sizes; only the global budget changes.

Not every rule needs concurrency. Rules that are I/O-bound and individually lightweight (BAM filtering, samtools stats, M-bias, spike-in methylation calling) have no concurrency declaration — Snakemake schedules them freely based on available cores. Concurrency is reserved for rules that are either resource-intensive (alignment: concurrency = 100), memory-constrained (methylKit R jobs: concurrency = 50), or benefit from controlled parallelism (methylation calling, dedup: concurrency = 25).

To run the pipeline, pass both --cores and --resources concurrency=N:

# Small machine (8 cores)
snakemake -s workflows/test.smk --cores 8 --resources concurrency=100

# Large workstation (48 cores)
snakemake -s workflows/test.smk --cores 48 --resources concurrency=100

# Very large machine (96+ cores, run more jobs in parallel)
snakemake -s workflows/test.smk --cores 96 --resources concurrency=300

Reference genomes are specified as a nested map in emseq_ref_assemblies. Each entry includes a download URL, short name, and the expected input filename. The core module handles indexing automatically. This map should include both the primary alignment reference and the spike-in control references used for conversion rate estimation.

emseq_ref_assemblies:
  # Primary alignment reference
  ncbi_decoy_hg38:
    url: https://ftp.ncbi.nlm.nih.gov/...
    name: ncbi_decoy_hg38
    input: GCA_000001405.15_GRCh38_no_alt_plus_hs38d1_analysis_set.fna.gz
  # Spike-in controls (EM-seq kit includes these for conversion rate QC)
  unmeth_lambda:
    url: https://www.neb.com/...
    name: unmeth_lambda
    input: lambda.fa.gz
  puc19:
    url: https://www.neb.com/...
    name: puc19
    input: pUC19.fa.gz

The meth-map config key defines one or more differential methylation experiments. Each experiment specifies the samples, reference, aligner, treatment vector, and methylKit parameters.

meth-map:
  tumor_vs_normal:
    libs: ["sample_A", "sample_B", "sample_C", "sample_D"]
    emseq_ref_name: ["ncbi_decoy_hg38"]
    align_method: ["bwa_meth"]
    tx: [1, 1, 0, 0]
    mincov: 10
    mingroup: 1
    chunksize: "1e9"
    win_size: 10000

Four conda environment YAMLs are provided in config/. Wrappers reference them by relative path. When running with --use-conda, Snakemake creates isolated prefix-based environments from these files.

The core pipeline (emseq.smk) has no external prerequisites — all tools are installed via conda and reference indexing is automated.

The analysis pipeline (emseq_analysis.smk) requires external tool repositories and reference files that must be prepared once before running. These are one-time setup steps that depend on the genome build. The results feed into the analysis-specific YAML config keys listed above.

Clone the three external repos:

./tools/setup_repos.sh

Then build each tool:

# mHapTools — ships macOS binary; must be compiled from source on Linux
cd ~/repos/mHapTools
cd htslib && ./configure && make -j4 && cd ..
sed -i 's/libhts.dylib/libhts.so/' CMakeLists.txt
cmake -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j4
cp build/mhaptools .

# wgbs_tools — C extensions must be built inside the deconv conda env
cd ~/repos/wgbs_tools
/path/to/deconv-conda-prefix/bin/python setup.py

# UXM_deconv — pure Python, no build needed

Set the repos config keys to the resulting paths:

repos:
  mhaptools: ~/repos/mHapTools
  wgbs_tools: ~/repos/wgbs_tools
  uxm_deconv: ~/repos/UXM_deconv

Three reference files are needed. Detailed preparation steps are documented in emseq.org under ** Reference preparation.

CpG reference → haplotype.cpg-ref — a tabixed file of all CpG dinucleotide positions in the reference genome. Generate from the reference FASTA:

awk 'BEGIN{OFS="\t"} /^>/{chr=substr($1,2);pos=0;next}
  {seq=toupper($0); for(i=1;i<length(seq);i++)
    if(substr(seq,i,2)=="CG") print chr,pos+i-1,pos+i;
  pos+=length(seq)}' ref.fa | bgzip > hg38_CpG.gz
tabix -b 2 -e 3 hg38_CpG.gz

MHB reference → haplotype.mhb-bed — methylation haplotype blocks from Guo et al. 2017 Nature Genetics. The published coordinates are hg19; liftOver to hg38:

# 1. Download Supplementary Table 1 from Guo 2017 (Excel)
# 2. Extract chr/start/end to BED (see emseq.org Reference preparation)
# 3. liftOver hg19 → hg38
liftOver guo2017_mhb_hg19.bed hg19ToHg38.over.chain.gz guo2017_mhb_hg38.bed unmapped.bed
sort -k1,1 -k2,2n guo2017_mhb_hg38.bed > guo2017_mhb_hg38_sorted.bed

UXM atlas → deconv.atlas — available in the cloned UXM_deconv repo:

# Use the atlas matching your genome build:
~/repos/UXM_deconv/supplemental/Atlas.U25.l4.hg38.tsv

Paired-end FASTQ files in, per-sample CpG methylation calls and QC metrics out.

The wrapper assigns all variables before include: "emseq.smk" (see the Wrapper Snakefile section for the full pattern):

• Trimming and QC: Adapter removal and quality trimming (fastp), read quality assessment before and after trimming (FastQC)
• Alignment: Bisulfite-aware alignment to one or more reference genomes via BWA-meth and/or Biscuit, coordinate sorting
• Duplicate handling: Duplicate marking with dupsifter, then BAM filtering that removes marked duplicates along with secondary, supplementary, and failed QC reads
• Methylation calling: CpG methylation extraction in methylKit format (MethylDackel), per-sample methylKit object creation, tabix indexing
• Coverage: Depth profiling with threshold summaries and quantization tracks (mosdepth), aggregated coverage plots
• Spike-in controls: Alignment to pUC19 (methylated) and Lambda (unmethylated) spike-in references for conversion rate estimation
• Bias detection: Position-dependent methylation bias assessment (MethylDackel mbias)
• Reporting: Aggregated QC report (MultiQC)

Design decision — duplicate marking vs. removal: Deduplication is deliberately split into two rules. emseq_dedup uses dupsifter to mark duplicates (SAM flag 0x400) without removing them. emseq_filter_bam then removes marked duplicates as part of a combined filter (-F 3840: secondary + supplementary + failed QC + duplicates) while also requiring proper pairs (-f 2), MAPQ ≥ 30, and intersection with keep/exclude BED regions. Separating these steps ensures that the marked-but-unfiltered BAM is available for M-bias assessment, where the full read set is needed to accurately detect position-dependent biases before duplicate removal.

Filtered BAMs and methylation calls in, differential methylation results, haplotype metrics, and tissue deconvolution out.

All core variables (above), plus these analysis-specific variables assigned in the wrapper:

External tool repositories must be cloned before running (see tools/setup_repos.sh).

• Differential methylation: Per-sample methylKit unite across an experiment, differential methylation calling at both per-base and tiled (windowed) resolution, extraction of positional methylation values for significant sites (methylKit)
• CpG annotation: Genomic feature annotation of differentially methylated regions — CpG islands, shores, shelves, inter-CGI regions, gene promoters, exons, introns (annotatr)
• Methylation haplotypes: BAM-to-mhap format conversion (mHapTools), per-MHB haplotype metrics including MHL, PDR, and entropy (mhaptk)
• Tissue deconvolution: Chromosome-prefix BAM reheading, pat/beta file generation (wgbs_tools), cell-type-of-origin deconvolution against a reference atlas (UXM)

The repository includes in-repo test data and wrapper Snakefiles for both modules. Test data consists of real EM-seq reads subsetted to chr22 with matching spike-in references (pUC19, Lambda), blacklist regions, and analysis references (CpG sites, MHB blocks).

• workflows/test.smk — tests the core module (emseq.smk) end-to-end: FASTQ → trimming → alignment (BWA-meth + Biscuit) → dedup → filter → methylation calling → QC
• workflows/test-analysis.smk — tests the analysis module (emseq_analysis.smk): DMR calling → annotation
• config/test.yaml — shared test configuration defining samples, references, experiments, and conda environments

To run locally:

# Core pipeline
snakemake -s workflows/test.smk --configfile config/test.yaml --cores 4 --use-conda --resources concurrency=100

# Analysis pipeline (requires core outputs + external repos)
snakemake -s workflows/test-analysis.smk --configfile config/test.yaml --cores 4 --use-conda --resources concurrency=100

The in-repo test data is subsetted to chr22 to keep the repository small. This is sufficient for validating the full core pipeline and the DMR/annotation portion of the analysis pipeline. Two analysis features are excluded from CI testing and included as commented-out targets in the test wrapper:

• Methylation haplotypes (mHapTools/mhaptk) — mHapTools requires C++ compilation from source with a bundled htslib version; too fragile for CI runners.
• Tissue deconvolution (wgbs_tools/UXM) — requires broader genomic coverage than a single chromosome provides.

Both have been validated on production-scale data.

Core processing (emseq.smk):

Downstream analysis (emseq_analysis.smk):

• Development since last tag
  • None
• wf/emseq/v5.1.0
  • [2026-04-07 Tue] Parameterized align_method in modules:
    • New required wrapper variable emseq_align_methods — list of alignment methods to run. Replaces hardcoded ["biscuit","bwa_meth"] in emseq_multiqc and emseq_analysis_uxm_deconv rules.
    • Wrappers must now define emseq_align_methods before including modules. Updated test.smk, test-analysis.smk, README documentation, and module header comments.
    • Removed obsolete test-wrap.smk wrapper and old untangled development workflow from emseq.org.
• wf/emseq/v5.0.0
  • [2026-03-31 Tue] Two-module architecture with comprehensive documentation:
    • Split README into per-module documentation (intro, wrapper variables, processing steps, DAG, outputs).
    • Added Configuration section: YAML/wrapper architecture, analysis module setup, resource management with basecamp concurrency model.
    • Section banners and rationale comments throughout both modules.
    • Updated Writeup section for two-module split.
    • Biscuit align: added samtools view -u -t ref.fai pipe fix.
    • Relocated orphaned R scripts to active Methylkit section; removed duplicates.
    • Untracked stale diff/ files.
  • [2026-03-31 Tue] CI overhaul:
    • Added analysis CI workflows (smk-analysis-dry, smk-analysis-run). All four CI workflows pass.
    • Analysis CI tests DMR + annotation; haplotypes and deconvolution excluded (validated on production data only).
    • Added --conda-frontend conda to all CI workflows (mamba prefix collision workaround).
    • Added test reference fixtures: chr22 CpG ref, MHB BED, UXM atlas subset.
  • [2026-03-31 Tue] Resource and architecture changes:
    • Reverted to direct thread/concurrency values per rule (basecamp pattern).
    • Removed project-specific AERO wrapper files.
  • [2026-03-25 Wed] Production validation fixes:
    • Added mkdir -p to emseq_align_bwameth_spike shell block for samtools sort tmp dir.
    • Added wildcard_constraints: experiment = "[^.]+" to emseq_analysis_methylkit_unite and emseq_analysis_methylkit_diff to prevent greedy matching of .tiled suffix.
    • Haplotype conda env: added mhaptk, pandas, matplotlib, scipy, seaborn, tqdm as pip deps.
    • Fixed mosdepth quantize string bug (single → global comma replacement).
  • [2026-03-16 Mon] Added emseq_analysis.smk module: DMR (migrated from emseq.smk), methylation haplotypes (mHapTools/mhaptk), tissue deconvolution (wgbs_tools/UXM), CpG annotation. New conda envs, test-analysis.smk wrapper, reference preparation docs.
  • Added CI workflows: core dry-run, core full run, analysis dry-run, analysis full run.
• wf/emseq/v4.0.0
  • Minor testing updates
  • Fixed test.yaml file locations
  • Changed bam filtering step to remove duplicates marked by dupsifter. So rule emseq_dedup MARKS duplicates while rule emseq_filterbam removes them as part of `samtools view -@ {threads} -u -f 2 -q 30 -F 3840 "{input.bam}`
  • Changed mosdepth to take filtered bam and quantify by keep.bed, ignoring filtered regions in depth calculation.
  • Conforms to v1.3 snakemake style guide
• wf/emseq/v3.1.0
  • [2025-09-19 Fri] Added a first github workflow test
  • [2025-09-19 Fri] Robust annotation of methylkit outputs validated as rscript
  • [2025-09-19 Fri] Added mature github tests for building in-repo test data and smk_dry
• [2025-09-19 Fri] Updated EM-seq main pipeline to wf/emseq/v3.0.0.
  • Includes in-repo small test data for a complete run of emseq.smk
  • Includes test.smk wrapper and corresponding test.yaml for in-repo small test run
  • emseq.smk expanded to include differential methylation from nested list map
  • Many small fixes for consistent naming and run condition optimization
• [2025-09-18 Thu] Updated EM-seq main pipeline to wf/emseq/v2.0.0. Mainly improved and simplified variable naming.