Imported from majiayu000/claude-skill-registry (skills/data/methylation-aggregation-ammawla-encode-toolkit/SKILL.md). Install upstream with npx skills add majiayu000/claude-skill-registry --skill methylation-aggregation-ammawla-encode-toolkit. Copyright stays with the author.

Aggregate DNA Methylation Data Across Studies

When to Use

• User wants to build a tissue-level DNA methylation landscape from multiple WGBS experiments
• User asks "where is DNA methylated in brain?" or "find hypomethylated regions across donors"
• User needs to identify HMRs (hypomethylated regions), UMRs, or PMDs from aggregated WGBS data
• User wants per-CpG weighted methylation averages from multiple experiments
• Example queries: "aggregate WGBS data for liver", "build methylation map across donors", "find unmethylated CpG islands in pancreas"

Build a comprehensive methylation landscape for a tissue/cell type by merging WGBS bedMethyl files from multiple ENCODE experiments.

Scientific Rationale

The question: "What is the DNA methylation state across the genome in my tissue?"

DNA methylation is fundamentally different from histone marks and accessibility:

The key insight: Unlike histone ChIP-seq where we want the union of all peaks, for methylation we want the average methylation level per CpG site across individuals. Methylation is a quantitative, continuous signal measured at every CpG dinucleotide.

However, for identifying regulatory regions, we focus on hypomethylated regions (HMRs) — stretches of low methylation that mark active regulatory elements. HMRs can be treated more like peaks for union-style aggregation.

Literature Support

• Roadmap Epigenomics (Schultz et al. 2015, Nature, 2,900+ citations): Established that tissue-specific HMRs mark active regulatory elements; demonstrated per-CpG averaging across biological replicates as standard approach
• DMRcate (Peters et al. 2021, Nucleic Acids Research, 65 citations): Method for calling differentially methylated regions from multiple WGBS samples; uses kernel smoothing across CpG sites
• ENCODE Phase 3 (Gorkin et al. 2020, Nature, 301 citations): Integrated methylation data with histone marks and accessibility to define chromatin states
• ENCODE Blacklist (Amemiya et al. 2019, Scientific Reports, 1,372 citations): Problematic genomic regions to filter. DOI
• Zhou et al. 2020 (Nature Genetics): Tissue-specific methylation patterns: ~80% of CpGs are constitutively methylated, ~10% constitutively unmethylated (CpG islands/promoters), ~10% tissue-variable
• Liu et al. 2024 (Briefings in Bioinformatics): Cross-platform comparison (NovaSeq vs DNBSEQ) showing WGBS is gold standard; coverage depth critically affects accuracy; platform differences exist in GC-rich regions
• Ortega-Recalde et al. 2021 (Methods in Molecular Biology): Demonstrated that even low-coverage WGBS can accurately estimate global methylation levels, with bootstrap methods to quantify uncertainty

Two-Level Analysis

1. Per-CpG level: Average methylation at each CpG site across samples (quantitative map)
2. Region level: Identify HMRs, PMDs, and UMRs from the averaged profile (union of regulatory regions)

Step 1: Find All Available WGBS Data

encode_search_experiments(
    assay_title="WGBS",
    organ="pancreas",           # user's tissue of interest
    biosample_type="tissue",
    limit=100
)

Present a summary to the user:

• Total WGBS experiments
• Labs represented
• Unique donors/biosamples
• Genome coverage per experiment

Use encode_get_facets to check availability:

encode_get_facets(assay_title="WGBS", organ="pancreas")

Note: WGBS is expensive to generate. Typical tissues have 2-5 experiments. Even 2 biological replicates are valuable for identifying consistent methylation patterns.

Step 2: Quality-Gate Each Experiment

encode_get_experiment(accession="ENCSR...")

WGBS Quality Checks

• Audit status: no ERROR flags
• Bisulfite conversion rate: >=98% (measured by lambda spike-in or non-CpG methylation)
• Genome coverage: >=10x mean CpG coverage for reliable per-site estimates
• Mapping rate: >=50% (bisulfite-converted reads are harder to map)
• Duplication rate: <30%
• Has methylation state at CpG output files (bedMethyl format)

Include if:

• Bisulfite conversion >=98%
• Mean CpG coverage >=10x
• Has bedMethyl output files for GRCh38

Exclude if:

• ERROR audit flags
• Conversion rate <98% (unconverted reads create false methylation calls)
• Very low coverage (<5x mean) — individual CpG estimates unreliable

Track all included experiments:

encode_track_experiment(accession="ENCSR...")

Step 3: Download bedMethyl Files

For each experiment:

encode_list_files(
    experiment_accession="ENCSR...",
    output_type="methylation state at CpG",
    assembly="GRCh38"
)

bedMethyl format (ENCODE standard):

chr  start  end  name  score  strand  thickStart  thickEnd  color  coverage  percentMethylated

• Column 10: read coverage at this CpG
• Column 11: percent methylation (0-100)

Prefer preferred_default=True files:

encode_download_files(
    file_accessions=["ENCFF...", ...],
    download_dir="/path/to/data/methylation",
    organize_by="flat"
)

Step 4: Per-Sample Quality Filtering

4a. Coverage Filtering (CRITICAL)

Low-coverage CpGs have unreliable methylation estimates. Filter per-sample:

# Keep only CpGs with >= 5x coverage (column 10)
# More stringent: >= 10x for quantitative analysis
awk '$10 >= 5' sample.bedMethyl > sample.covfiltered.bedMethyl

Coverage thresholds by use case:

4b. ENCODE Blocklist Filtering (Amemiya et al. 2019)

# Download from: https://github.com/Boyle-Lab/Blacklist/blob/master/lists/hg38-blacklist.v2.bed.gz
bedtools intersect -a sample.covfiltered.bedMethyl -b hg38-blacklist.v2.bed -v > sample.filtered.bedMethyl

4c. Strand Merging (Optional but Recommended)

CpG methylation is typically symmetric (same on both strands). Merge strand-specific calls to increase per-CpG coverage. Caveat: Some ENCODE bedMethyl files may already be strand-merged — check if both strands are present before applying this step:

# Group CpGs by position (forward and reverse strand of same CpG)
# Sum coverage, calculate weighted average methylation
awk 'BEGIN{OFS="\t"} {
    # CpG position (use the C position as canonical)
    if ($6 == "+") pos = $2
    else pos = $2 - 1
    key = $1"\t"pos
    cov[key] += $10
    meth[key] += ($11/100) * $10
} END {
    for (k in cov) {
        split(k, a, "\t")
        avg_meth = (meth[k] / cov[k]) * 100
        print a[1], a[2], a[2]+2, "CpG", 0, ".", a[2], a[2]+2, "0,0,0", cov[k], avg_meth
    }
}' sample.filtered.bedMethyl | sort -k1,1 -k2,2n > sample.merged_strands.bedMethyl

Step 5: Cross-Sample Aggregation (Per-CpG Averaging)

5a. Create a Unified CpG Matrix

# Step 1: Find CpGs covered in at least M of N samples
# Extract positions from each sample
for f in sample*.merged_strands.bedMethyl; do
    awk 'BEGIN{OFS="\t"} {print $1, $2, $3}' "$f"
done | sort -k1,1 -k2,2n | uniq -c | \
awk -v M=2 '$1 >= M {print $2, $3, $4}' OFS="\t" > shared_cpgs.bed

# Step 2: For each sample, extract methylation at shared CpGs
for f in sample*.merged_strands.bedMethyl; do
    bedtools intersect -a shared_cpgs.bed -b "$f" -wa -wb | \
    awk 'BEGIN{OFS="\t"} {print $1, $2, $3, $NF, $(NF-1)}' > "${f%.bedMethyl}.shared.txt"
    # Columns: chr, start, end, percentMeth, coverage
done

5b. Calculate Average Methylation Per CpG

Weighted average (recommended — accounts for coverage differences):

# Combine all samples, calculate coverage-weighted mean methylation per CpG
cat sample*.shared.txt | \
sort -k1,1 -k2,2n | \
awk 'BEGIN{OFS="\t"} {
    key = $1"\t"$2"\t"$3
    if (key != prev_key && NR > 1) {
        avg = total_weighted_meth / total_cov
        print prev_key, n_samples, total_cov, avg
        n_samples = 0; total_cov = 0; total_weighted_meth = 0
    }
    prev_key = key
    n_samples++
    total_cov += $5
    total_weighted_meth += ($4/100) * $5
} END {
    avg = total_weighted_meth / total_cov
    print prev_key, n_samples, total_cov, avg
}' > tissue_methylation_profile.bed
# Columns: chr, start, end, n_samples, total_coverage, mean_methylation_fraction

Simple average (alternative — equal weight per sample):

# Unweighted mean across samples
awk 'BEGIN{OFS="\t"} {
    key = $1"\t"$2"\t"$3
    meth[key] += $4
    n[key]++
} END {
    for (k in meth) {
        print k, n[k], meth[k]/n[k]
    }
}' <(cat sample*.shared.txt) | sort -k1,1 -k2,2n > tissue_methylation_simple.bed

5c. Calculate Methylation Variability

Track inter-individual variation to identify tissue-variable CpGs:

# Add standard deviation column
# (compute in R or Python for large datasets)

Step 6: Identify Regulatory Methylation Features

6a. Hypomethylated Regions (HMRs)

HMRs mark active regulatory elements. Identify runs of low methylation. The 30% threshold is a commonly used cutoff (Schultz et al. 2015 used similar ranges), but the optimal threshold depends on your tissue and question — some studies use 20%, others 40%. Consider visualizing the methylation distribution first to identify a natural breakpoint:

# Find CpGs with average methylation < 30% (adjust threshold as needed)
awk '$6 < 0.30' tissue_methylation_profile.bed > hypo_cpgs.bed

# Merge adjacent hypomethylated CpGs into regions
# Require minimum 3 CpGs within 1kb of each other
bedtools merge -i hypo_cpgs.bed -d 1000 -c 1 -o count | \
awk '$4 >= 3' > tissue_HMRs.bed
# Columns: chr, start, end, n_hypomethylated_CpGs

6b. Unmethylated Regions (UMRs) — CpG Islands

Very low methylation (<10%) at CpG-dense regions:

awk '$6 < 0.10' tissue_methylation_profile.bed > unmeth_cpgs.bed
bedtools merge -i unmeth_cpgs.bed -d 500 -c 1 -o count | \
awk '$4 >= 5' > tissue_UMRs.bed

6c. Partially Methylated Domains (PMDs)

Large (>10kb) regions of intermediate methylation, often marking repressed regions:

# Find CpGs with methylation 30-70% (partially methylated)
awk '$6 >= 0.30 && $6 <= 0.70' tissue_methylation_profile.bed > partial_cpgs.bed

# Merge with large gap tolerance to find domains
bedtools merge -i partial_cpgs.bed -d 5000 -c 1 -o count | \
awk '$4 >= 20 && ($3-$2) >= 10000' > tissue_PMDs.bed

6d. Tissue-Specific Differentially Methylated Regions

If comparing to another tissue, use DMRcate or similar:

# In R with DMRcate
library(DMRcate)
# Requires a methylation matrix (CpGs x samples with tissue labels)
# Identifies regions where methylation differs between tissues

Step 7: Confidence Annotation

For HMRs/UMRs (region-level features), annotate by sample support:

# Annotate HMRs with sample support
# Intersect each HMR with per-sample hypomethylated CpGs to count support
awk -v N=4 '{
    # Using n_samples from the aggregation
    if ($4 >= N*0.5) conf="HIGH";
    else if ($4 >= 2) conf="SUPPORTED";
    else conf="VARIABLE";
    print $0"\t"conf"\t"$4"/"N
}' tissue_HMRs.bed > tissue_HMRs.annotated.bed

For the per-CpG profile, annotate by coverage confidence:

awk -v N=4 '{
    if ($4 >= N) conf="ALL_SAMPLES";
    else if ($4 >= N*0.5) conf="MAJORITY";
    else conf="PARTIAL";
    print $0"\t"conf
}' tissue_methylation_profile.bed > tissue_methylation.annotated.bed

Step 7b: Summary Statistics

Report to the user:

• Total input experiments: N
• Experiments passing QC: M (bisulfite conversion, coverage)
• Total CpGs per sample (before/after coverage filter)
• Shared CpGs across M+ samples: X
• Mean genome-wide methylation: Y%
• Number of HMRs: Z (with size distribution)
• Number of UMRs: W
• Number of PMDs: V (if applicable)
• High-confidence HMRs: how many in ≥50% of samples
• CpGs with high inter-individual variability

Step 8: Integration with Other ENCODE Data

Methylation data is most powerful when integrated:

1. HMRs + H3K27ac peaks = Active enhancers (use histone-aggregation skill)
2. HMRs + ATAC-seq peaks = Open regulatory elements (use accessibility-aggregation skill)
3. UMRs + H3K4me3 peaks = Active promoters
4. PMDs = Often overlap H3K9me3 (heterochromatin)

# Example: Find HMRs that overlap H3K27ac peaks (active enhancers)
bedtools intersect -a tissue_HMRs.bed -b union_H3K27ac_peaks.bed -wa -u > active_enhancer_HMRs.bed

Step 9: Log Provenance

encode_log_derived_file(
    file_path="/path/to/tissue_methylation.annotated.bed",
    source_accessions=["ENCSR...", "ENCSR...", ...],
    description="Aggregated per-CpG methylation profile across N pancreas WGBS experiments",
    file_type="aggregated_methylation",
    tool_used="bedtools + custom aggregation",
    parameters="coverage >= 5x per sample, shared CpGs in >= 2 samples, coverage-weighted mean, strand-merged"
)

encode_log_derived_file(
    file_path="/path/to/tissue_HMRs.annotated.bed",
    source_accessions=["ENCSR...", "ENCSR...", ...],
    description="Hypomethylated regions from aggregated pancreas methylation profile",
    file_type="aggregated_HMRs",
    tool_used="bedtools merge",
    parameters="mean methylation < 30%, >= 3 CpGs within 1kb, confidence annotated"
)

Pitfalls Specific to Methylation Data

1. Bisulfite conversion rate is critical: Even 1% incomplete conversion creates false methylation at unmethylated CpGs. Always verify >=98% conversion. ENCODE reports this in QC metrics.

2. Coverage drives accuracy: A CpG with 3x coverage has wide confidence intervals (0-100% could easily be 0% or 30%). At 10x, estimates stabilize. At 30x, they are reliable. Always filter by coverage.

3. Non-CpG methylation: Present in some cell types (especially embryonic). ENCODE bedMethyl files typically report CpG context only. If non-CpG methylation is relevant, check experiment metadata.

4. Strand asymmetry: While CpG methylation is typically symmetric, it can be asymmetric at some sites. Strand merging loses this information. For most analyses, merging is appropriate.

5. Cell-type heterogeneity: Bulk WGBS from tissue captures methylation across ALL cell types. A CpG at 50% methylation could mean: (a) all cells are 50% methylated, or (b) half the cells are 0% and half are 100%. These are biologically different. Single-cell methylation data (if available) resolves this.

6. CpG islands vs. open sea: CpG-dense regions (islands) have very different methylation dynamics than CpG-sparse regions. Consider analyzing separately.

7. Do NOT mix assemblies: All files must be GRCh38 or all hg19. CpG positions are exact — even a 1bp offset from liftOver misaligns CpGs.

8. X chromosome: Males have one X (hemimethylation), females have two (one inactivated with different methylation). Handle sex chromosomes separately or filter them.

9. Imprinted regions: Some regions show ~50% methylation in all individuals due to genomic imprinting (one allele methylated, one not). These are normal, not noise.

10. Do NOT use union logic for per-CpG methylation: Unlike histone peaks where union is correct, methylation levels should be AVERAGED. Union logic only applies to the derived HMR/UMR/PMD regions.

11. RRBS is NOT the same as WGBS: Reduced Representation Bisulfite Sequencing (RRBS) covers only CpG-rich regions (~10% of CpGs). Do NOT mix RRBS and WGBS in per-CpG averaging — the CpG universe is different.

12. Sequencing platform matters: Liu et al. 2024 showed that NovaSeq and DNBSEQ-T7 give slightly different methylation estimates, especially in GC-rich regions. Note the platform in provenance if mixing experiments from different sequencers.

Walkthrough: Cross-Tissue CpG Methylation Atlas for Imprinted Gene Regions

Goal: Aggregate whole-genome bisulfite sequencing (WGBS) data across tissues to identify tissue-invariant vs. tissue-specific methylation patterns at imprinted gene loci. Context: Imprinted genes show parent-of-origin-specific methylation. Comparing across tissues reveals which imprinting control regions (ICRs) maintain methylation universally.

Step 1: Find WGBS experiments across tissues

encode_search_experiments(assay_title="WGBS", organism="Homo sapiens", limit=50)

Expected output:

{
  "total": 147,
  "results": [
    {"accession": "ENCSR765JPC", "assay_title": "WGBS", "biosample_summary": "liver", "status": "released"},
    {"accession": "ENCSR832HMR", "assay_title": "WGBS", "biosample_summary": "brain", "status": "released"},
    {"accession": "ENCSR091ENJ", "assay_title": "WGBS", "biosample_summary": "lung", "status": "released"}
  ]
}

Interpretation: 147 WGBS experiments available. Select tissues with ≥2 replicates for reliable per-CpG averaging.

Step 2: List methylation bedGraph files

encode_list_files(accession="ENCSR765JPC", file_format="bed", output_type="methylation state at CpG", assembly="GRCh38")

Expected output:

{
  "files": [
    {"accession": "ENCFF123BED", "output_type": "methylation state at CpG", "file_format": "bed bedMethyl", "file_size_mb": 245.0}
  ]
}

Step 3: Download methylation files

encode_download_files(accessions=["ENCFF123BED", "ENCFF456MET", "ENCFF789CPG"], download_dir="/data/wgbs")

Step 4: Per-CpG weighted averaging across replicates

For each tissue:

1. Merge replicates using weighted average: β = Σ(methylated reads) / Σ(total reads)
2. Filter CpGs with <10× combined coverage
3. Output: per-tissue methylation BED with columns: chr, start, end, β-value, coverage

Step 5: Identify tissue-invariant ICRs

# H19/IGF2 ICR: chr11:2,016,000-2,022,000
bedtools intersect -a merged_methylation.bed -b imprinted_icrs.bed -wa -wb | \
  awk '{sum+=$4; n++} END {print sum/n}'

Interpretation: ICRs showing ~50% methylation across ALL tissues confirm maintained imprinting. Tissue-variable ICRs (range >20%) suggest tissue-specific imprinting loss.

Integration with downstream skills

• Feed differentially methylated regions into → peak-annotation for nearest gene assignment
• Overlay with → histone-aggregation H3K4me3 to find promoter methylation–expression anticorrelation
• Cross-reference CpG variants via → clinvar-annotation for methylation-disrupting mutations
• Compare methylation at regulatory elements from → regulatory-elements

Code Examples

1. Survey WGBS data availability by organ

encode_get_facets(assay_title="WGBS", facet_field="organ", organism="Homo sapiens")

Expected output:

{
  "facets": {
    "organ": {"brain": 32, "liver": 18, "heart": 12, "lung": 10, "blood": 8, "kidney": 6}
  }
}

2. Check experiment quality before aggregation

encode_get_experiment(accession="ENCSR765JPC")

Expected output:

{
  "accession": "ENCSR765JPC",
  "assay_title": "WGBS",
  "biosample_summary": "liver",
  "replicates": 2,
  "status": "released",
  "audit": {"WARNING": 0, "ERROR": 0},
  "pipeline": "WGBS (GRCh38)"
}

3. Track aggregated experiments

encode_track_experiment(accession="ENCSR765JPC", notes="Liver WGBS for cross-tissue methylation atlas")

Expected output:

{
  "status": "tracked",
  "accession": "ENCSR765JPC",
  "notes": "Liver WGBS for cross-tissue methylation atlas",
  "tracked_at": "2025-03-08T12:00:00Z"
}

Integration

Related Skills

• histone-aggregation: HMRs + H3K27ac union peaks identify active enhancers; HMRs + H3K4me3 identify active promoters
• accessibility-aggregation: HMRs typically overlap open chromatin; concordance between HMRs and ATAC/DNase peaks validates both
• hic-aggregation: Hypomethylated enhancers often anchor chromatin loops to target genes
• regulatory-elements: Combine methylation with histone and accessibility data to classify regulatory element types
• epigenome-profiling: Methylation adds a critical layer to chromatin state annotation
• pipeline-wgbs: Process raw WGBS data through the full ENCODE-aligned pipeline
• batch-analysis: Batch processing workflows for systematic methylation aggregation
• publication-trust: Verify literature claims backing analytical decisions

Presenting Results

• Present methylation summary as: total CpGs analyzed | mean coverage | methylation distribution (UMR/LMR/PMD). Show per-sample contribution. Suggest: "Would you like to correlate with histone marks?"

For the request: "$ARGUMENTS"