This repository contains Nextflow-based analysis tools for Enzymatic Methylation Sequencing (EM-seq) and Enzymatic 5hmC-seq (E5hmC-seq) data processing.
Complete EM-seq processing pipeline that accepts UBAM inputs:
- Adapter trimming and read alignment with (fastp, bwa-meth)
- Duplicate marking (Picard)
- Methylation calling (MethylDackel)
- Quality control metrics and statistics (Picard, Samtools, FastQC, MultiQC)
- Optional BED file intersection for targeted analysis (bedtools)
If your files are in fastq format, you will need to convert them to uBams prior to running the main pipeline, e.g.:
nextflow run fastq_to_ubam.nf \
--input_glob "tests/fixtures/fastq/emseq-test*{.ds.1,.ds.2}.fastq.gz" \
--read_format 'paired-end'| Parameter | Description | Default |
|---|---|---|
--input_glob |
glob for your gzipped fastq files | ['*.{1,2}.fastq.gz'] |
--read_format |
'paired-end' or 'single-end' | 'paired-end' |
Quantifies read depth and CpG methylation across major genomic feature categories — exons, CDS, genes, mRNA, CpG islands, and EPD promoters — for a set of aligned BAM files.
Outputs
feature_depth_bokeh.html— interactive per-feature depth plot (switchable in-browser between hex-density and violin display, log/linear y-axis, per-feature or shared y-range, and a depth-threshold counter)
combined_feature_counts.tsv— merged featureCounts table across all feature types- Per-sample methylation TSVs under
<output_dir>/features/<sample>/<context>/
Usage
nextflow run feature_cov_meth.nf \
--bam_files_glob '*.md.{bam,bam.bai}' \
--mk_files '*.methylKit.gz' \
--human_t2t2Use --mouse for GRCm39, --human_t2t2 for T2T CHM13v2.0, or supply all genome params manually for a custom assembly.
| Parameter | Description | Default |
|---|---|---|
--mouse |
GRCm39 reference shortcuts | false |
--human_t2t2 |
T2T CHM13v2.0 reference shortcuts | false |
--local_ref_files_path |
Root path to locally cached reference files | (set in config) |
--context |
Methylation context(s), comma-separated (CpG, CHG, CHH) |
CpG |
--count_dup_reads |
Include duplicate reads in feature counts | false |
--output_dir |
Output directory | cov_vs_meth.output |
For custom assemblies set --genome, --ucsc_cpg_islands_gtf, --refseq_gff_url, --ncbi_assembly_report_url, --epd_promoter_bed_url, --old_new_chain_url, --cpg_chr_lookup, and --refseq_chr_lookup directly.
- Install miniforge and bioconda (see Requirements)
- Install Nextflow (e.g. conda install nextflow, or see Nextflow installation guide)
- Clone this repository (
git clone https://github.com/nebiolabs/EM-seq.git). Modifynextflow.configas needed for your environment, e.g. if running locally, change executor block to 'local' and set, e.g.--max_cpus 10 --max_memory 30.GB. - Download or prepare a genome reference FASTA file (see Reference Genomes)
- Create a bwameth index for the fasta and add it to your references in conf/references.config
- Run the pipeline with appropriate parameters (see Basic Usage)
- Examine results in the EM-seq_output directory
EM-seq-Alignment-Summary-<FLOWCELL_ID>_multiqc_report.htmlin em-seq_output for overall QC summary- Mbias files
em-seq_output/methylDackelExtracts/mbias(to identify sample-dependent positional biases) - Methylation output files in
em-seq_output/methylDackelExtracts(suitable for analysis with methylKit) - Aligned reads in
em-seq_output/markduped_bams(methylation coloring is recommended for visualization in IGV)
nextflow run main.nf \
--genome 'test' \
--ubam_dir './' \
--email your.email@example.com \
--flowcell FLOWCELL_ID \
--adapter_set 'Illumina'ubam_dir should be the folder where your ubam files are.
| Parameter | Description | Default |
|---|---|---|
--genome |
reference genome found in conf/references.config | Required |
--email |
Email for notifications | Required |
--flowcell |
Flowcell identifier | Optional |
--outputDir |
Output directory | em-seq_output |
--enable_neb_agg |
Enable NEB aggregation reporting | False |
--adapter_set |
adapter set found in conf/adapter_sequences.config (Illumina, Element, MGI, Nextera-mosaic) |
Required |
--multiqc_gc_group |
contig group whose GC bias curve MultiQC plots (see GC bias curves) | largest group |
Modify the conf/references.config file to specify your genome files
genome_fapath to your genome fasta filegenome_faipath to your genome fasta fai filebwameth_indexpath to your genome fasta file where bwameth indices existtarget_bedBED file for targeted analysis, Optionalgc_groups_dirdirectory of per-organism GC bias references, Optional (see below)
Picard measures GC bias by comparing where reads start against the GC of every window in the
reference it is given, so a composite reference reports one curve blended across every organism in
it. Setting gc_groups_dir for a genome adds a separate curve per organism alongside the
whole-reference one. The directory holds a contig_groups.tsv assigning each contig to a group,
plus a <group>.fa/.fa.fai/.dict per group; the contigs to keep and the region BED are derived
from each group's .fa.fai at run time. The directory is validated against genome_fai and
against contig_groups.tsv at startup.
MultiQC plots a single curve per library. By default that is the curve of the group with the most
windows -- the host organism for the +meth_controls references -- because the whole-reference curve
counts the spike-in controls, and pUC19 (51% GC), lambda (50%) and Xp12 (68%) pull it noticeably
away from a host-only curve. Use --multiqc_gc_group <name> to plot a different group.
Without gc_groups_dir you get the whole-reference curve alone, controls included; the run warns
that it is not directly comparable to a host-only curve.
Group definitions for the published references ship in
assets/contig_groups/; the per-group FASTA subsets they refer to are
built separately and are not part of this repository.
--tmp_dir- Temporary directory (default:/tmp)--workflow- Workflow identifier (default:EM-seq)--enable_neb_agg- Enable NEB aggregation reporting (default:False)
Pre-built reference genomes with methylation spike-in controls, plus matching UCSC CpG island annotations. These are the same references NEB uses internally:
| Species | Assembly | Genome FASTA | CpG islands gtf |
|---|---|---|---|
| Human | T2T CHM13v2.0 | T2T_chm13v2.0+meth_controls | t2t_cpg_islands |
| Human | GRCh38 | grch38_core+meth_controls | grch38_cpg_islands |
| Mouse | GRCm39 | grcm39+meth_controls | grcm39_cpg_islands |
Download and unzip the FASTA, then add its paths to conf/references.config (see
References Config). The pipeline reads the reference from disk and never
fetches it, so it also needs a .fai and a bwameth index alongside; main.nf checks for both at
startup and prints the command to build whichever is missing.
The CpG island files are only used by feature_cov_meth.nf. Pass one with
--ucsc_cpg_islands_gtf, or save the downloads under --local_ref_files_path (default
~/nebnext_projects/em-seq/em-seq_ref_files) keeping their names, and the --mouse and
--human_t2t2 shortcuts will find them. GRCh38 has no shortcut, so it always needs these
parameters passed explicitly.
Per-organism GC bias curves for these references are configured by the group definitions in
assets/contig_groups/ — see GC bias curves.
If you already have one of these references from before September 2026, re-download it. The
phage_lambda,phage_T4andphage_Xp12control sequences were corrected to match RefSeq, so an older copy gives control metrics that are not comparable. You must rebuild the bwameth index and the.faiafter replacing the FASTA.
- Create your own reference by appending the control sequences to your preferred genome fasta (e.g.
cat genome.fa methylation_controls.fa > genome+methylation_controls.fa)Sequence Methylation State Purpose lambda All Cs are unmodified (included in kits) confirmation that APOBEC enzyme is working optimially (unprotected C->T) pUC19c All CpG sites contain 5mC (included EM-seq kits) confirmation that TET2 + T4-BGT step is protecting 5mC (5mC -> 5ghmC/5caC T4 All Cs are 5hmC (included in the 5-hmC Seq kits) confirmation that T4-BGT is protecting 5hmC optimally (5hmC -> 5ghmC)
- Nextflow
- Miniforge, Micromamba, or Conda for dependency management
- Bioconda channel configured
- Sufficient computational resources (memory scales with input size)
These in the "legacy" folder are retained for reference and reproducibility but are not actively maintained and are not compatible with the latest Nextflow versions. Use NXF_VER=22.10.4 nextflow run ... to reproduce the results in the EM-seq paper.
em-seq.nf- Original alignment and methylation calling workflowbins.nf- TSS-centered binned coverage analysiscov_vs_meth.nf- Coverage vs methylation analysis for genomic features
Analysis methods in this repository were used in the following publication:
Vaisvila R, Ponnaluri VKC, Sun Z, et al. Enzymatic methyl sequencing detects DNA methylation at single-base resolution from picograms of DNA. Genome Res. 2021;31(7):1280-1289. doi:10.1101/gr.266551.120
- We have a related project to call variants in EM-seq converted libraries: EM-seq variant calling
- You may also be interested in the nf-core methylseq project
- git tag -f current_production
- git push -f origin current_production
- development workflow will run from master branch
- Tests are run using nf-test and are integrated into github actions
- install nf-test from bioconda using conda/mamba
- To run all tests:
nf-test test- When new tests are added or results change, to update the results snapshot:
nf-test test --updateSnapshotNextflow dropped upstream support for v24 and older in July 2026, but this pipeline still targets
24.10.x -- that is what CI pins and what the reference runs use -- so main.nf ships with
nextflow.preview.topic = true as its top line.
Nextflow 25.x made topic channels stable and removed that directive. To run on 25.x or newer, comment out or delete that top line of main.nf.