# 04.6 methylKit DMR analysis plan

Plan for finding differentially methylated regions (DMRs) with methylKit in the Haws ploidy × pH experiment. It builds on the 04.2/04.3 per-CpG parameter tests and is designed to run alongside the DSS plan in [04.4-DSS-plan.md](04.4-DSS-plan.md).

## What 04.2/04.3 tell us

04.3 per-CpG results (lo.count = 5, hi.perc = 99.9, q < 0.01), from `analyses/Haws_04.3-methylKit-slurm/DML-hiperc-99.9-overdispersion-*/DML-counts-long.csv`:

| samples | min.per.group | CpGs tested | ploidy 25%, `none` | ploidy 25%, `MN` | pH 25%, `none` | pH 25%, `MN` |
|---|---|---|---|---|---|---|
| All | All | 1.96M | 730 | 0 | 1096 | 0 |
| All | 11 | 3.14M | 1784 | 2 | 2416 | 2 |
| All | 8 | 5.09M | 5114 | 7 | 6896 | 5 |
| Outliers removed | All | 2.21M | 1280 | 1 | 1440 | 1 |

`shrinkMN` is still running (array 41067297, as of 2026-10-02).

What this means for DMRs:

1. **The `none` DML are mostly variation between oysters, not treatment effects.** Without overdispersion, methylKit treats the reads as independent binomial draws and pools them across the 12 oysters in each group, so differences between individual oysters look like a treatment effect. Once overdispersion is modeled (`MN`), almost nothing is left. Any DMR approach that **combines `none` DML into regions** (merging nearby significant CpGs) would just group those false positives. That approach is ruled out.
2. **The `none` DML count rises with the number of CpGs tested** (730 → 5114 as min.per.group goes from All to 8). The `MN` count stays near 0. That's the pattern you'd expect from false positives: they scale with the number of tests.
3. **Per-CpG power with `MN` is very low.** With 12 oysters per group and ~5–20x per CpG, a single CpG can't separate a real treatment effect from differences between oysters. Regions help in two ways. Combining CpGs raises coverage per unit, and C. gigas gene-body methylation tends to change across whole regions rather than at single sites, so a region-level effect is biologically reasonable.
4. **Regions also need overdispersion modeling.** Summing reads across CpGs in a tile inflates apparent coverage, because CpGs on the same read aren't independent. Under the binomial model that makes things worse, not better. **Every region test uses `MN` (or `shrinkMN`). `none` is run only as a negative-control comparison.**
5. **04.3 drops all mito CpGs** through `hi.perc` (mito coverage is ~4× the nuclear 99.9th percentile). Fix this here the same way as in the DSS plan.

So the question for this analysis: **with overdispersion modeled, does combining CpGs into regions give enough power to find ploidy or pH effects? And do those regions hold up under label permutation?**

## Design (same as 04.4)

- 2 × 2 factorial, ploidy (2N, 3N) × pH (high, low), n = 6 per cell, 24 oysters.
- **All 24 oysters is the primary analysis** (decided 2026-10-02 from the 04.4 sample QC, `analyses/Haws_04.4-DSS/sample-QC/`). All 24 samples pass the technical checks. 3H-2 is the real outlier: it makes up most of PC1, has the lowest correlation with other samples (0.869) and the highest global methylation (11.7%). 2H-1, 2H-2 and 3H-5 separate on PC2, which tracks global methylation. 2H-3 is not an outlier; it sits in the main cluster and fails no QC flag.
- **Sample-set sensitivity checks** (same sets as 04.4 D4):
  1. `Drop3H2`: drop 3H-2 (sample 14)
  2. `DropPC2`: drop 2H-1, 2H-2, 3H-5 (samples 1, 2, 17)
  3. `DropGeno4`: drop all four genetic outliers, 2H-1, 2H-2, 3H-2, 3H-5 (samples 1, 2, 14, 17)
  4. `OutRM`: drop 2H-3 + 3H-2 (samples 3, 14), the pair removed in 04.2/04.3, kept only for comparison with those runs
- **Genetic background is confounded with pH in the all-24 analysis** (04.4 SNP check, `analyses/Haws_04.4-DSS/SNP-check/`). The four methylation outliers (2H-1, 2H-2, 3H-5, 3H-2) are also genetic outliers, and all four are high-pH (about a 5% chance if spread at random). Genotype similarity tracks methylation similarity (rho = 0.86). Ploidy is balanced (2 diploid, 2 triploid). So an all-24 pH test can't separate pH from genetic background, and the drop sets above are the main pH tests, not optional extras. A genotype-PC covariate doesn't help for pH: in DSS it removed pH signal along with the genetic effect (65 observed vs. a permuted median of 144). See [04.6-methylKit-DMR-results.md](04.6-methylKit-DMR-results.md) for how this affects the results.
- **Model:** `calculateDiffMeth(treatment = <factor tested>, covariates = <other factor>, overdispersion = "MN")`. This is the additive `~ ploidy + pH` model, compared by a likelihood-ratio / F-test against the covariate-only model, as in 04.3.
- **Limitation:** methylKit can't include the ploidy × pH interaction as a nuisance term, because covariates are added to both the full and reduced models and treatment is a single term. The design is balanced, so the additive estimates are unbiased averages over the other factor. But regions with a strong interaction (effect in only one cell, or opposite directions) will lose power or be missed. The DSS factorial model (04.4 A) covers that case. Concordance (step F) shows how much it matters.

## Analyses

### A. Tiled windows (primary)

Fixed, non-overlapping windows from `tileMethylCounts`.

Pipeline for each setting:

```r
raw      <- methRead(... , mincov = 1)               # keep low-coverage CpGs; tile coverage is what gets filtered
raw      <- filterHiCov(raw, perc = 99.9, exempt = "NC_001276.1")   # per-sample nuclear percentile, mito exempt (helper, see Implementation)
tiles    <- tileMethylCounts(raw, win.size = W, step.size = W, cov.bases = K, mc.cores = n)
tiles    <- filterByCoverage(tiles, lo.count = 10) %>% normalizeCoverage()
united   <- unite(tiles, destrand = FALSE, min.per.group = M)
dm       <- calculateDiffMeth(united, covariates = cov, overdispersion = "MN", test = "F", mc.cores = n)
```

- **Order matters:** remove CpGs with very high coverage (likely PCR duplicates or repeats) **before** tiling, so they don't dominate a tile. Filter and normalize **at the tile level** after tiling. 04.3's prep objects are already normalized per CpG, so don't tile them. Tile from the raw coverage files instead.
- Use `step.size = win.size`. Overlapping windows test the same CpGs twice and make BH FDR anti-conservative.
- Strand: the `.cov` files are already merged (`merged_CpG_evidence`), so `destrand = FALSE` as in 04.3.

### B. Annotation-based regions (secondary, has a clear biological interpretation)

`regionCounts(raw, regions = GRanges, cov.bases = K)` over predefined features, then the same filter → unite → `calculateDiffMeth(MN)` steps:

| Region set | Source |
|---|---|
| Gene bodies | `genome-feature-files/cgigas_uk_roslin_v1_gene.gff` (38,208 genes) |
| Exons / UTR | `cgigas_uk_roslin_v1_exonUTR.gff` |
| Introns | gene minus exonUTR (`bedtools subtract`) |
| Promoters / upstream 1 kb | `cgigas_uk_roslin_v1_upstream.gff` |
| Downstream flanks | `cgigas_uk_roslin_v1_downstream.gff` |
| lncRNA | `cgigas_uk_roslin_v1_lncRNA.gff` |
| Transposable elements | Build from the NCBI RepeatMasker output, [GCF_902806645.1_cgigas_uk_roslin_v1_rm.out.gz](https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/29159/102/GCF_902806645.1_cgigas_uk_roslin_v1/GCF_902806645.1_cgigas_uk_roslin_v1_rm.out.gz) (annotation release 102; RepeatMasker open-4.0.8 with Dfam/RepBase 20181026). See the note on TE annotation below the table. |
| Mito genome | NC_001276.1 as one region, plus per-gene if there's an annotation |

**TE annotation (checked 2026-10-02).** The lab's [cgigas_uk_roslin_v1_rm.te.bed](http://owl.fish.washington.edu/halfshell/genomic-databank/cgigas_uk_roslin_v1_rm.te.bed) has the same 781,312 rows as `rm.out.gz`, so it is that file with the class column dropped. That causes two problems:

1. **It isn't only TEs.** It includes Simple_repeat (165,597), Low_complexity (21,815), Satellite (1,461), rRNA (63) and Unknown (9,973) repeats, about 25% of the rows.
2. **Starts are off by one.** BED starts are copied from RepeatMasker's 1-based `begin` (e.g. 4109) instead of `begin − 1`.

Instead, make `cgigas_uk_roslin_v1_TE.bed` in prep from `rm.out.gz`: start = begin − 1, with columns for class/family and strand. Keep DNA (380,683 rows, 116.8 Mb), RC/Helitron (151,061, 84.3 Mb), LINE (29,494, 13.5 Mb), LTR (19,774, 13.6 Mb) and SINE (1,391, 0.2 Mb). Keep Unknown as a separate set. Non-TE repeats can be a separate "simple repeats" region set if they're wanted.

Test regions at two levels: all TEs together, and by class (DNA, RC, LINE, LTR). `bedtools merge` within each set before `regionCounts`, so a CpG isn't counted twice. Many TEs are short (median 141 bp in the old bed), so the `cov.bases` ≥ 3 filter will remove many of them. Report how many are left in each class. There are no repeats on the mito genome. Download the file in prep instead of committing it. (Don't use the *C. virginica* `C_virginica-3.0_TE-all.gff`, which is for a different assembly.)

There are far fewer tests here (~38k genes against millions of CpGs), so multiple-testing correction costs much less. And the results feed straight into 08/09 enrichment.

### C. Data-driven segments (exploratory, group-blind only)

`methSeg()` divides the genome into segments of similar methylation. To avoid double-dipping:

- **Do segment** on group-blind data: methylation pooled across all 24 oysters (`pool()` → `methSeg` on the pooled `methylRaw`). This finds methylated and unmethylated domains (gene-body "mosaic" blocks) without using treatment labels. Then test those segments with `regionCounts` + `MN`, as in B.
- **Do not segment** on a `methylDiff` (treatment differences). Segments chosen because they differ between groups will look significant when tested on the same data, so the p-values can't be trusted.

### D. Parameter grid

| Parameter | Levels | Notes |
|---|---|---|
| Window size `W` (A) | 100, 250, 500, 1000 bp | 100 bp ≈ "clusters of CpGs"; 1000 bp ≈ methylKit default; C. gigas methylated blocks are typically kb-scale |
| Min CpGs per region `cov.bases` (`K`) | 3, 5 | Below 3, a "region" is really a single CpG |
| Tile coverage `lo.count` | 10, 20 (summed across CpGs) | |
| Upper coverage | 99.9th percentile per sample, nuclear only, mito exempt | Fixes the 04.3 mito drop |
| Sample presence `min.per.group` | All, 10, 8 | Drop 11 and 9; 04.3 shows a steady trend across them |
| Samples | All 24 (primary); `Drop3H2`, `DropPC2`, `DropGeno4` (required for pH, see Design); `OutRM` (comparison only) | `SAMPLES` in 04.6-methylKit-DMR-slurm.R |
| Overdispersion | `MN` (primary), `shrinkMN` (after 04.3 finishes), `none` (negative control only) | |
| Thresholds | **Fixed: q < 0.05 and \|Δ\| ≥ 10%** (decided 2026-10-02, before any region results) | Averaging across a region shrinks the difference, so the 25% per-CpG cutoff is too strict for regions. 10% is common in invertebrate DMR studies. Counts at q < 0.01 and \|Δ\| ≥ 15/25% are reported only as supplementary sensitivity checks and do not define the DMR set. **q-only tier (added 2026-10-02, after the first results):** counts at q < 0.05 with any \|Δ\|, compared with permutations, to check whether small but precisely measured differences (median \|Δ\| ≈ 1% among windows with p < 0.001) pass more often than chance. Exploratory; the DMR definition doesn't change. |

Full A grid: 4 W × 2 K × 2 lo.count × 3 presence × 2 sample sets × 2 tests = 192 fits for `MN`. Each fit has far fewer rows than the per-CpG runs (≈ 50k–1M tiles against 2–5M CpGs), so each one is fast. Run the core grid first (W ∈ {250, 1000}, K = 3, lo.count = 10, presence = All and 8, all 24 samples), then expand only along whichever dimensions change the result.

### E. Calibration (decide which counts to believe before reading biology into them)

1. **p-value histograms** for each window size and overdispersion setting. With `MN`, expect a roughly flat histogram with a spike near 0 if there's signal. With `none`, expect a pile-up near 0 even in permuted data.
2. **Label permutations:** use the same scheme as DSS 04.4 D2, sharing the same permutation seeds so the two methods can be compared. Permute ploidy within pH (6/6 per cell kept) and pH within ploidy. Run 20–50 permutations of the chosen setting(s). Report the median and 95th percentile of the number of "DMRs" in permuted data next to the observed count. **A DMR set is only reported if the observed count is clearly above the permutation distribution** (e.g. above its 95th percentile).
3. **Per-CpG → region consistency:** for each DMR, the share of its CpGs whose per-CpG `meth.diff` (from the existing 04.3 `MN` `.rds` files) has the same sign. A real regional effect should be mostly one direction.
4. **Sample-set sensitivity:** for each of `Drop3H2`, `DropPC2`, `DropGeno4` and `OutRM`, count DMRs and how many all-24 DMRs remain. Judge each set against its own permutation null, since individual oysters seem to drive most per-CpG "DML" in 04.4 (100 DSS permutations gave as many DML as the real labels, up to 1,129 in a single permutation). For pH with all 24 oysters, also use `PERM_SCHEME = fixGeno4`, which keeps the four genetic outliers in high pH in every permutation, matching the real confounding.

### F. Concordance and downstream

- **Against DSS (04.4):** overlap of methylKit tiles/segments with DSS `callDMR` regions. Correlation of region-level Δ across shared regions. The interaction diagnostic from 04.4 A shows whether methylKit misses DSS DMRs because of an interaction.
- **Against per-CpG results:** what fraction of the 04.3 `none` DML fall inside `MN` DMRs. A low fraction supports point 1 above (those DML were variation between oysters).
- **Across region types:** do significant tiles (A) concentrate in significant genes (B)?
- **Downstream:** the final DMR set (annotated with gene/exon/intron/flank/TE) and, separately, gene-level results from B go to 07–09. A ranked gene list from B (by signed −log10 p) can be used with GSEA-style enrichment, which works even if few regions pass FDR.

## Implementation

Uses the same staged design as 04.3, and the file numbering used by the DSS pipeline (scripts numbered after their plan):

```
code/Haws/04.6-methylKit-DMR-slurm.R           stages: prep | regions | segments | count | dmr | summary
code/Haws/04.6-methylKit-DMR-slurm.job         same container + R_LIBS_USER as 04.3
code/Haws/04.6-methylKit-DMR-slurm-submit.sh   dependency chain; settings via env vars; skips finished steps
analyses/Haws_04.6-methylKit-DMR/              outputs; .gitignore prep-*/, count-*/, DMR-*/rds/, regions/*.rds, downloads/, logs/
```

Settings are environment variables, as in 04.4-DSS-slurm: `REGIONS`, `COV_BASES`, `LO_COUNT`, `HI_PERC`, `SAMPLES`, `OVERDISPERSION`, `PERM`. The submit script takes space-separated lists for `REGIONS`, `OVERDISPERSION` and `PERM`.

- **prep:** `methRead(mincov = 1)` for all 24 files. For each sample, compute the 99.9th percentile of coverage on nuclear CpGs only and keep CpGs with `coverage <= cutoff` or on `NC_001276.1` (`methylKit::select()`). This replaces `filterByCoverage(hi.perc)`, which can't exempt one chromosome. Writes per-sample cutoffs and mito CpG counts to `high-coverage-filter.tsv`. One object for all 24 samples; the outliers-removed set is made later with `reorganize`.
- **regions:** builds the annotation region sets (B) from `genome-feature-files/` and the NCBI RepeatMasker output (downloaded by the submit script), and writes a 0-based, TE-only `regions/cgigas_uk_roslin_v1_TE.bed` with class/family for 07.
- **segments:** group-blind `methSeg` on counts summed over all 24 oysters (C).
- **count:** one job per (region set, `COV_BASES`): `tileMethylCounts` (non-overlapping windows) or `regionCounts`.
- **dmr:** a 6-task array per run: (ploidy, pH) × `min.per.group` (All, 10, 8). Each task: `reorganize` → `filterByCoverage(lo.count)` → `normalizeCoverage` → `unite` → `calculateDiffMeth(covariates, overdispersion)`. Writes counts at each q/Δ threshold (methylKit SLIM q, plus BH q for comparison with DSS), p-value histogram bins, and for observed labels only, the full results `.rds` and a BED of DMRs at the primary threshold. Gene-level results get gene IDs.
- **Permutations** reuse the dmr stage with `PERM=<seed>`. For annotation region sets (genes, exons, TEs, ...) they also save per-region `meth.diff` and p-values, so gene-set enrichment can be tested against the same permutations. Window runs save counts only. The shuffle is copied from `makeDesign()` in 04.4-DSS-slurm.R, so a seed gives the same permuted labels in both pipelines.
- **summary:** combines every run into `DMR-counts-all.csv`, and writes `DMR-permutation-summary.csv` (observed count, permutation median and 95th percentile, empirical p).
- **Resources:** prep gets 120G (04.3 prep peaked at 34G with mincov = 2) and count 200G (tiling over 4 cores). Both are generous starting points. dmr tasks start at 16 CPUs / 64G; adjust from the first logs.
- The readable `04.6-methylKit-DMR.Rmd` walkthrough (one setting end to end, with figures) comes after the core grid has run.

## Order of work

0. Let the 04.3 `shrinkMN` array finish (41067297) and add it to the per-CpG table. Also check per-CpG `MN` p-value histograms from the existing `.rds` files. If they're far from uniform, look into that before going further.
1. Write `filterHiCov()` and the prep stage. Check per-sample cutoffs and mito retention against 04.4's numbers (mito mean ~317x, nuclear cutoff ~79x for sample 1).
2. Core tile grid (A): W ∈ {250, 1000}, K = 3, lo.count = 10, all 24 samples, presence All and 8, `MN` and `none`. First DMR counts.
3. Permutations (E2) on those settings. **Decision point:** if observed `MN` DMR counts are within the permutation distribution, report "no detectable DMRs at this sample size" as the result. Then concentrate on gene-level (B), ranked-list enrichment, and DSS, instead of increasing the grid until something passes.
4. Run the pH tests without the genetic outliers (`DropPC2`, `Drop3H2`, `DropGeno4`, each with its own permutations) regardless of the all-24 result, because of the genotype–pH confounding. If there's signal: run the rest of the grid (D), shrinkMN, and `OutRM` for comparison with 04.2/04.3.
5. Annotation-based regions (B) and group-blind segments (C).
6. Concordance (F) with DSS and per-CpG results, figures (Δ distribution, chromosome distribution with mito shown separately, heatmaps of top DMRs), then pass the result to 07–09.

## Open questions

- Did 07's `DML-*-TE.bed` overlaps use the *C. virginica* GFF (wrong assembly) or `cgigas_uk_roslin_v1_rm.te.bed`? The Roslin bed would include simple repeats and is off by one. Either way, 07's TE results should be redone with the TE bed built from `rm.out.gz`.
- If DSS and methylKit `MN` both find almost nothing, the paper's conclusion becomes "ploidy and pH have minimal effects on methylation at specific loci." Global and PCA-level differences (06) are still worth reporting.