# 04.6 methylKit DMR results

Results of the region-level methylKit analysis planned in [04.6-methylKit-DMR-plan.md](04.6-methylKit-DMR-plan.md), run 2026-10-02 with [04.6-methylKit-DMR-slurm.R](04.6-methylKit-DMR-slurm.R). Outputs are in `analyses/Haws_04.6-methylKit-DMR/`.

## Summary

**With overdispersion modeled and label permutations as the null, methylKit finds no region-level effect of ploidy or pH**, whether the regions are 250 bp windows, 1 kb windows, or whole genes. Gene-set enrichment (GSEA with NCBI GO terms) finds none either. The real labels give about as many hits as shuffled labels at every threshold tested.

This agrees with the per-CpG results in 04.3 and the DSS results in 04.4: at n = 6 oysters per ploidy × pH group, neither treatment produces methylation differences that can be told apart from differences between individual oysters.

**The genetic structure found in 04.4 doesn't change this.** That's four genetically distinct oysters, all at high pH. The result is the same with those oysters removed, three different ways, and with a permutation null that keeps them in place. The one weak pH lead (9 genes at q < 0.05 with 3 outliers dropped) disappears when all four are dropped. The top all-24 pH gene set weakens from FDR 0.078 to 0.275 under the genotype-matched null. See [Does the genetic structure found in 04.4 change this?](#does-the-genetic-structure-found-in-044-change-this)

Two methodological findings matter for any further methylation analysis of this data:

1. **Without overdispersion (`none`), methylKit reports thousands of DMRs**, 4,539–9,056 at q < 0.05 and |Δ| ≥ 10%. These reflect differences between individual oysters, not treatment.
2. **With overdispersion (`MN`), methylKit's q-values are still too optimistic.** With shuffled labels, BH q < 0.05 flags 27–185 windows in *every* permutation, where at most about 1 in 20 permutations should flag anything. Significance for this data has to come from permutations, not from methylKit q-values.

## Methods in brief

- **Data:** the 24 Bismark coverage files (Roslin assembly), all CpGs with at least 1 read, about 10M per oyster.
- **High-coverage filter:** CpGs above each oyster's 99.9th-percentile coverage are removed (cutoffs 55–102x), calculated on nuclear CpGs only. Mito CpGs are kept, unlike 04.3, which drops them all. `methRead` was called with `header = FALSE`; the default treats the first CpG of each file as a header and drops it.
- **Regions:**
  - **Windows:** non-overlapping 250 bp and 1 kb windows with at least 3 covered CpGs.
  - **Genes:** whole gene bodies from `cgigas_uk_roslin_v1_gene.gff`.
  - **Coverage:** at least 10 reads per region in each oyster.
  - **Presence:** coverage in all 24 oysters (primary), with at least 10 or 8 per group as checks.
- **Test:** `calculateDiffMeth` with `overdispersion = "MN"` (F-test). The other factor is a covariate, so this is the additive model `~ ploidy + pH`; methylKit can't include the ploidy × pH interaction.
- **Negative control:** the same tests with `overdispersion = "none"`.
- **DMR definition (fixed before any results):** q < 0.05 and |Δ| ≥ 10 percentage points. Counts at q < 0.01, |Δ| ≥ 15/25%, and a q-only tier (any |Δ|) were also recorded. The q-only tier was added after the first results.
- **Permutations:** ploidy is shuffled within pH and pH within ploidy, using the same code and seeds as 04.4 DSS. That gives 20 for windows and 100 for genes.
- **Enrichment:** genes are ranked by signed −log10 p. NCBI GO terms (`gene_ontology.gaf.gz`, with GO.db ancestors added) are kept for sets of 10–500 genes, giving 2,196 sets. Each set gets a GSEA enrichment score. NES, nominal p, and FDR all come from the 100 label permutations. Each permutation was also scored against the others, as a calibration check.

## Results

### DMR counts (q < 0.05, |Δ| ≥ 10%, all 24 oysters)

| Regions | Regions tested | Ploidy, `MN` | pH, `MN` | Ploidy, `none` | pH, `none` |
|---|---|---|---|---|---|
| 250 bp windows | 721,545 (65 mito) | 1 | 2 | 8,704 | 9,052 |
| 1 kb windows | 424,348 (18 mito) | 0 | 2 | 4,539 | 5,164 |
| genes | 34,058 | 0 | 0 | — | — |

No mito region was a DMR in any test.

### Real labels vs. permutations (`MN`)

| Regions | Test | DMRs: real | DMRs: permuted median (range) | q < 0.05, any Δ: real | Permuted median (range) |
|---|---|---|---|---|---|
| 1 kb windows | ploidy | 0 | 1 (0–6) | 37 | 65 (27–104) |
| 1 kb windows | pH | 2 | 1 (0–3) | 70 | 70.5 (≤108) |
| 250 bp windows | ploidy | 1 | 2.5 | 131 | 117 (≤199) |
| 250 bp windows | pH | 2 | 2 | 130 | 117.5 (≤170) |
| genes | ploidy | 0 | 0 | 0 | 0 (0–3) |
| genes | pH | 0 | 0 | 1 | 0 (0–2) |

Across all 72 combinations of region set, test, presence setting, and threshold, the real count is within the permutation distribution, usually near the median. The one borderline case (250 bp, pH, q < 0.01 with any Δ: 25 against a permuted median of 9.5, empirical p = 0.05) is what chance gives at least once in 72 comparisons. Full table: `DMR-permutation-summary.csv`.

### Why there are no DMRs: significance and effect size don't overlap

- **Regions that pass q < 0.05 differ by very little.** The median |Δ| among windows with p < 0.001 is about 1 percentage point. These are well-covered regions where a small shift can be measured precisely.
- **Regions with large differences aren't significant.** Thousands of regions differ by ≥ 10%, but they're noisy, low-coverage regions that the overdispersion model correctly doesn't trust.
- **Conversion efficiency doesn't explain the small shifts.** Bisulfite conversion differs between groups by only about 0.01 percentage points (99.50% in every group).

### P-value calibration

- **Overall shape:** p-value histograms are close to flat. 5.5–6.6% of p-values are below 0.05, and shuffled labels give the same share (5.5–7.0%).
- **The extreme tail:** a small number of regions get extremely small p-values (down to about 1e-12, with |Δ| about 1%) whether the labels are real or shuffled. That's why BH q < 0.05 flags windows in every permutation. The likely cause is well-covered regions where methylKit's overdispersion estimate falls to its floor. This hasn't been confirmed yet.

### Gene-set enrichment (GSEA, GO, 2,196 sets, 100 permutations)

| | Ploidy | pH |
|---|---|---|
| Sets at FDR < 0.05 | 0 | 0 |
| Sets at FDR < 0.25 | 0 | 1 |
| Nominal p < 0.05 | 81 (3.7%) | 40 (1.8%) |
| Shuffled labels with ≥ 1 set at FDR < 0.05 | 17 / 100 | 15 / 100 |
| Shuffled labels with ≥ 1 set at FDR < 0.25 | 50 / 100 | 48 / 100 |

- **The one pH set:** *mitotic spindle assembly checkpoint signaling* (GO:0007094, 12 genes, NES −1.95, FDR 0.078). Its genes are 2.2 points less methylated at low pH on average.
- **It isn't distinguishable from chance:** about half of the shuffled labelings also give at least one set at FDR < 0.25.
- **Its neighbours aren't independent support:** nearby sets (nuclear division, chromosome segregation, sister chromatid separation) lean the same way, but they share its genes.
- **The FDR estimate runs liberal:** most shuffled labelings flag nothing, but a few flag hundreds of sets (up to 652 for ploidy and 1,168 for pH).

Tables: `enrichment-gene-cb3-lo10-hiperc99.9-All-MN/GSEA-GO-{ploidy,pH}-All.csv`, plus `-calibration.csv`.

## Does the genetic structure found in 04.4 change this?

The 04.4 SNP check (`analyses/Haws_04.4-DSS/SNP-check/`) found that **genetic background accounts for much of the methylation difference between oysters**:

- **Genotype tracks methylation:** genotype similarity tracks methylation similarity across oyster pairs (rho = 0.86).
- **Four genetic outliers:** 2H-1, 2H-2, 3H-5 and 3H-2 separate in both the genotype and methylation PCAs.
- **All four are high-pH oysters.** If they were spread at random, that would happen about 5% of the time. Two are diploid and two triploid, so ploidy is balanced.

**Ploidy: unlikely to change.** The outliers are split 2/2 between diploids and triploids, and among the other 20 oysters genotype isn't related to treatment. DSS with genotype PC1 as a covariate still finds no ploidy effect (265 observed vs. a permuted median of 212).

**pH: it could change, and the current analysis doesn't settle it.** The four genetic outliers can affect a pH test in two opposite ways:

1. **False signal.** If genotype shifts their methylation, a high-vs-low comparison is partly "these four oysters vs. the rest". That would inflate the real pH count above the permutations. It didn't happen: real pH counts are at or below the permutation median. So the outliers aren't creating false DMRs here.
2. **Hidden signal.** The outliers add between-oyster variation, all of it in the high-pH group. `MN` treats that variation as noise, so the pH test loses power. A real but modest pH effect could be masked. Dropping the outliers is the way to find out.

The permutations also aren't a perfect null for pH. They spread the four outliers across both pH groups at random, but the real labels have all four in high pH.

**What wouldn't help:**
- **A genotype covariate for pH.** Genotype PC1 is confounded with pH, so adjusting for it removes pH signal along with the genetic effect. In DSS, pH with genotype PC1 gave 65 observed vs. a permuted median of 144, below chance.
- **Masking SNP CpGs.** In DSS, CpGs with SNPs don't explain the null (noSNP: 146 observed vs. 153 permuted for ploidy). At the region level, a few SNP CpGs are diluted further.

**Check run (2026-10-02).** 1 kb windows and genes were rerun without the genetic outliers, each set with its own permutations (20 for windows, 100 for genes). The all-24 analysis was also rerun with a restricted permutation, `fixGeno4`: the four outliers keep their real labels, and only the other 20 oysters are shuffled. That matches the real genotype–pH link exactly.

| Sample set | Dropped | Oysters (high / low pH) |
|---|---|---|
| `DropPC2` | 2H-1, 2H-2, 3H-5 | 21 (9 / 12) |
| `Drop3H2` | 3H-2 | 23 (11 / 12) |
| `DropGeno4` | all four | 20 (8 / 12) |
| All, `fixGeno4` permutations | none | 24 (12 / 12) |

Real counts against each set's own permutations (coverage in all oysters, `MN`):

| Sample set | Regions | Test | DMRs (q < 0.05, \|Δ\| ≥ 10%): real / perm median | q < 0.05, any Δ: real / perm median (95th pct) | empirical p |
|---|---|---|---|---|---|
| `DropPC2` | 1 kb | ploidy | 4 / 2 | 115 / 126 (164) | 0.57 |
| `DropPC2` | 1 kb | pH | 3 / 2.5 | 119 / 121.5 (174) | 0.67 |
| `DropPC2` | genes | ploidy | 0 / 0 | 1 / 1 (4) | 0.57 |
| `DropPC2` | genes | pH | 0 / 0 | **9 / 1 (4)** | **0.02** |
| `Drop3H2` | 1 kb | ploidy | 0 / 1 | 46 / 78.5 (103) | 0.86 |
| `Drop3H2` | 1 kb | pH | 0 / 1 | 50 / 61 (114) | 0.81 |
| `Drop3H2` | genes | ploidy | 0 / 0 | 0 / 1 (4) | 1.00 |
| `Drop3H2` | genes | pH | 0 / 0 | 0 / 0 (4) | 1.00 |
| `DropGeno4` | 1 kb | ploidy | 2 / 3.5 | 120 / 134 (164) | 0.86 |
| `DropGeno4` | 1 kb | pH | 1 / 3.5 | 135 / 119.5 (167) | 0.48 |
| `DropGeno4` | genes | ploidy | 0 / 0 | 0 / 1 (4) | 1.00 |
| `DropGeno4` | genes | pH | 0 / 0 | 0 / 1 (4) | 1.00 |
| All, `fixGeno4` | 1 kb | ploidy | 0 / 0 | 37 / 66.5 (80) | 0.95 |
| All, `fixGeno4` | 1 kb | pH | 2 / 2 | 70 / 67.5 (98) | 0.48 |
| All, `fixGeno4` | genes | ploidy | 0 / 0 | 0 / 0 (2) | 1.00 |
| All, `fixGeno4` | genes | pH | 0 / 0 | 1 / 0 (5) | 0.49 |

The empirical p is for the q-only column.

GO enrichment (GSEA, 2,196 sets, 100 permutations per set):

| Sample set | Ploidy: sets at FDR < 0.25 | pH: sets at FDR < 0.25 | Lowest pH FDR (top set) | Shuffled labels with ≥ 1 set at FDR < 0.25 (pH) |
|---|---|---|---|---|
| All, free permutations | 0 | 1 | 0.078 (*mitotic spindle assembly checkpoint signaling*) | 48 / 100 |
| All, `fixGeno4` | 0 | 0 | 0.275 (same set) | 53 / 100 |
| `DropPC2` | 0 | 0 | 0.862 (*negative regulation of nuclear division*) | 51 / 100 |
| `Drop3H2` | 0 | 1 | 0.173 (*U5 snRNP*, 10 genes, NES +1.91) | 47 / 100 |
| `DropGeno4` | 0 | 0 | 0.561 (*G protein-coupled receptor signaling*) | 53 / 100 |

No set reaches FDR < 0.05 in any run.

**Conclusion: accounting for the genetic structure doesn't change the result.** Neither removing the genetic outliers nor using a null that keeps them in place reveals a ploidy or pH effect.

- **The one lead didn't replicate.** `DropPC2` had 9 genes at q < 0.05 for pH against a permutation median of 1 (p = 0.02). They're small shifts; none reaches |Δ| ≥ 10%. The lead is one of about 70 comparisons in this check, and it's gone in `DropGeno4` (0 genes, which also drops 3H-2). So it's most likely chance, not a pH effect masked by the outliers.
- **The all-24 pH gene set was partly driven by the genetic outliers.** *Mitotic spindle assembly checkpoint signaling* weakens from FDR 0.078 to 0.275 under the null that keeps the four outliers in high pH. It falls out of the top sets when they are removed (`DropPC2`, `DropGeno4`).
- **The restricted null was the more appropriate one for pH, and it changes nothing.** The `fixGeno4` distributions are close to the free ones, and the real counts sit near the middle of both.
- **Removing oysters costs power.** `DropGeno4` keeps only 8 high-pH oysters, so a null there means "no detectable effect", not "no effect".

## Limitations

- **No interaction term:** methylKit can't fit the ploidy × pH interaction, so regions where the effect differs between the two groups of the other factor could be missed. The DSS factorial model covers that case.
- **Whole-gene tests dilute effects:** the median gene differs by only 0.3 points. The exon, intron, upstream and TE region sets are built but not yet tested.
- **Permutation numbers:** windows have 20 permutations and genes 100. Nominal GSEA p-values can't go below about 0.02, because each set is compared only with permutations whose score has the same sign (about half).
- **The overdispersion floor explanation is unconfirmed:** the cause of the extreme p-values in a few regions hasn't been checked directly.