10. Epigenomics (ChIP-seq, ATAC-seq, Methyl-seq)

Model-based Analysis of ChIP-Seq (MACS3) captures the dynamic Poisson distribution of read counts across local genomic backgrounds to identify punctate transcription factor binding summits or broad histone mark enrichment domains.

# 1. ATAC-Seq Open Chromatin Peak Calling (shifted for Tn5 insertion)
macs3 callpeak \
  -t atac_aligned.bam \
  -f BAMPE -g hs \
  -n atac_sample \
  --nomodel --shift -100 --extsize 200 \
  --qvalue 0.01 \
  --outdir macs3_atac_peaks/

# 2. Transcription Factor Narrow Peak Calling with Input Control
macs3 callpeak \
  -t chip_tf.bam -c input_control.bam \
  -f BAM -g hs \
  -n tf_narrow \
  --qvalue 0.05 \
  --outdir macs3_tf_peaks/

# 3. Broad Histone Mark Peak Calling (H3K27me3, H3K36me3)
macs3 callpeak \
  -t chip_h3k27me3.bam -c input_control.bam \
  -f BAM -g hs \
  --broad --broad-cutoff 0.1 \
  -n h3k27me3_broad \
  --outdir macs3_broad_peaks/

Hypergeometric Optimization of Motif EnRichment (HOMER) discovers de novo and known transcription factor binding motifs in ChIP-seq and ATAC-seq peak summits and annotates peaks with distance to nearest transcription start sites (TSS).

# 1. Discover enriched transcription factor motifs in peak summits
findMotifsGenome.pl \
  atac_sample_peaks.narrowPeak \
  hg38 \
  homer_motifs_out/ \
  -size 200 \
  -mask \
  -p 12

# 2. Annotate peaks with nearest genes and genomic features
annotatePeaks.pl \
  atac_sample_peaks.narrowPeak \
  hg38 \
  > peak_gene_annotations.txt

Genrich is a peak caller designed specifically for ATAC-seq and DNase-seq data. It analyzes aligned paired-end reads directly from queryname-sorted BAMs, natively handling cut site intervals, unpaired cut ends, and mitochondrial exclusions.

# 1. Sort BAM by queryname (required by Genrich)
samtools sort -n -@ 8 -o atac_qsorted.bam atac_aligned.bam

# 2. Call peaks in ATAC mode with mitochondrial removal
Genrich \
  -t atac_qsorted.bam \
  -o genrich_peaks.narrowPeak \
  -j \
  -r \
  -e chrM \
  -v \
  -q 0.05

ArchR is a high-performance R package for single-cell chromatin accessibility. By storing data in disk-backed Arrow files, it processes hundreds of thousands of scATAC-seq cells with gene score modeling, peak calling, and trajectory inference.

library(ArchR)
addArchRThreads(threads = 16)
addArchRGenome("hg38")

# 1. Create Arrow files from 10x fragments.tsv.gz
inputFiles <- c("sample_fragments.tsv.gz")
ArrowFiles <- createArrowFiles(inputFiles = inputFiles, sampleNames = c("scATAC"),
                                minTSS = 4, minFrags = 1000, addTileMat = TRUE, addGeneScoreMat = TRUE)

# 2. Build ArchR Project and filter doublets
proj <- ArchRProject(ArrowFiles = ArrowFiles, outputDirectory = "ArchR_out")
proj <- filterDoublets(proj)

# 3. LSI Dimensionality reduction and UMAP
proj <- addIterativeLSI(proj, useMatrix = "TileMatrix", name = "IterativeLSI")
proj <- addClusters(proj, reducedDims = "IterativeLSI", resolution = 0.8)
proj <- addUMAP(proj, reducedDims = "IterativeLSI")
plotEmbedding(proj, colorBy = "cellColData", name = "Clusters")

Signac is an extension of Seurat for the analysis of single-cell chromatin data (scATAC-seq and 10x Multiome). It performs peak matrix construction, TF-IDF normalization, motif enrichment, and co-embedding with scRNA-seq.

library(Signac)
library(Seurat)

# 1. Create ChromatinAssay from fragment file
counts <- Read10X_h5("filtered_peak_bc_matrix.h5")
chrom_assay <- CreateChromatinAssay(counts = counts, sep = c(":", "-"),
                                     genome = 'hg38', fragments = 'fragments.tsv.gz')
seurat_obj <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

# 2. Normalization via TF-IDF and SVD (LSI)
seurat_obj <- RunTFIDF(seurat_obj)
seurat_obj <- FindTopFeatures(seurat_obj, min.cutoff = 'q0')
seurat_obj <- RunSVD(seurat_obj)

# 3. UMAP on chromatin components
seurat_obj <- RunUMAP(seurat_obj, reduction = 'lsi', dims = 2:30)
DimPlot(seurat_obj, reduction = "umap")

Bismark is the benchmark pipeline for mapping bisulfite-treated sequencing reads to in silico transformed genomes and extracting cytosine methylation calls across CpG, CHG, and CHH contexts in WGBS and RRBS assays.

# 1. Prepare in silico bisulfite converted reference genome
bismark_genome_preparation --bowtie2 /ref/GRCh38/

# 2. Align paired-end bisulfite reads
bismark --bowtie2 --multicore 8 /ref/GRCh38/ \
  -1 sample_val_1.fq.gz -2 sample_val_2.fq.gz

# 3. Deduplicate alignments and extract methylation calls
deduplicate_bismark -p --bam sample_val_1_bismark_bt2_pe.bam
bismark_methylation_extractor --paired-end --bedGraph --counts --buffer_size 10G \
  --multicore 8 sample_val_1_bismark_bt2_pe.deduplicated.bam

bwameth maps bisulfite-sequenced reads against in silico converted reference genomes using BWA-MEM. It runs substantially faster than legacy bisulfite aligners, scaling easily to whole-genome bisulfite sequencing cohorts.

# 1. Index reference with bwameth
bwameth.py index GRCh38.fa

# 2. Align bisulfite paired-end reads and stream to sorted BAM
bwameth.py -t 16 --reference GRCh38.fa \
  sample_R1.fq.gz sample_R2.fq.gz \
  | samtools sort -@ 4 -o bwameth_aligned.sorted.bam -

# 3. Extract methylation metrics with MethylDackel
MethylDackel extract GRCh38.fa bwameth_aligned.sorted.bam -o sample_methyldackel

methylKit is an R package for analyzing and comparing DNA methylation profiles from WGBS and RRBS. It performs coverage filtering, normalization, and tests for differentially methylated cytosines (DMCs) and regions (DMRs) using logistic regression.

library(methylKit)

# 1. Read methylation coverage calls
file.list <- list("ctrl1.CpG.txt", "ctrl2.CpG.txt", "treat1.CpG.txt", "treat2.CpG.txt")
myobj <- methRead(file.list, sample.id=list("c1","c2","t1","t2"),
                  assembly="hg38", treatment=c(0,0,1,1), context="CpG")

# 2. Filter by coverage and unite samples
filtered.myobj <- filterByCoverage(myobj, lo.count=10, hi.perc=99.9)
meth <- unite(filtered.myobj, destrand=TRUE)

# 3. Calculate differential methylation using logistic regression
myDiff <- calculateDiffMeth(meth, overdispersion="MN", test="Chisq")
myDiff25p <- getMethylDiff(myDiff, difference=25, qvalue=0.01)
write.csv(as.data.frame(myDiff25p), "significant_dmrs.csv")

modkit is Oxford Nanopore’s official tool for analyzing modified base calls directly from raw electrical currents without chemical bisulfite damage. It parses Mm/Ml BAM tags to produce site-level bedMethyl files for 5mC, 5hmC, and 6mA.

# 1. Aggregate native 5mC/5hmC methylation probabilities into bedMethyl
modkit pileup \
  ont_modified.sorted.bam \
  sample_ont_5mc.bed \
  --cpg \
  --ref GRCh38.fa \
  --threads 16 \
  --combine-strands

# 2. Extract summary statistics of genome-wide methylation
modkit summary ont_modified.sorted.bam --ref GRCh38.fa