# 04.4 DSS analysis plan

Plan for redoing differential methylation in the Haws ploidy × pH experiment with DSS, built to line up with the 04.2/04.3 methylKit parameter tests.

## Why DSS now

04.3 ran 04.2 with each overdispersion setting (hi.perc = 99.9, 25% difference, q < 0.01):

| samples | min.per.group | ploidy, `none` | ploidy, `MN` | pH, `none` | pH, `MN` |
|---|---|---|---|---|---|
| All | All | 730 | 0 | 1096 | 0 |
| All | 11 | 1784 | 2 | 2416 | 2 |
| All | 8 | 5114 | 7 | 6896 | 5 |
| Outliers removed | All | 1280 | 1 | 1440 | 1 |

Nearly all methylKit DML go away once variation between oysters is modeled. So the methylKit counts mostly reflect the binomial-only assumption, not treatment effects. DSS sits between those two settings. It fits a beta-binomial model with a dispersion estimated at each CpG and shrunk toward a genome-wide prior. That suits 6 oysters per cell. The main question is whether DSS also finds almost nothing (agreeing with MN) or finds a small, reproducible set.

The existing [04-DSS.Rmd](04-DSS.Rmd) (178 ploidy / 154 pH / 53 interaction DML) has issues to fix, not just rerun:

1. **Main effects are simple effects.** With `~ploidy + pH + ploidy:pH` and treatment coding, `ploidyT` is the ploidy effect **at high pH only**, and `pHL` is the pH effect **in diploids only**. They are not main effects.
2. No coverage filter (1x CpGs included) and no upper coverage filter.
3. `pvals < 0.05 & fdrs < 0.01` means p < 0.05 adds nothing. No effect-size threshold, so the results can't be compared with methylKit's 25%/50% cutoffs.
4. Hard-coded local paths (`/Users/yaaminivenkataraman/...`, `/Volumes/web-1/...`).

The coverage files are on the same Roslin assembly (NC_047559–NC_047568 + NC_001276.1 mito + 225 unplaced scaffolds, ~9.9M CpGs per sample, ~6.1M at ≥5x). Coordinates match the old DSS output, so the old output can be compared directly.

## Experimental design

2 × 2 factorial, fully crossed and balanced: ploidy (2N, 3N) × pH (high, low), n = 6 per cell, 24 oysters. Removing the outliers (2H-3, 3H-2) gives 5/6/5/6, which is fine for a GLM.

**Decisions (2026-10-02):**
- **The ploidy × pH interaction is a nuisance term.** It stays in the model so the ploidy and pH effects are estimated with it accounted for, but interaction DML are not reported as results.
- **The mitochondrial genome (NC_001276.1) stays in.** Mito CpGs have ~4× the genome-wide 99.9th-percentile coverage (sample 1: mean 317x vs. cutoff 79x), so the percentile filter has to be calculated on nuclear CpGs only and not applied to mito. Otherwise every mito CpG is removed, which is what 04.3 does now. Mito methylation is ~2%, close to the bisulfite non-conversion rate, so mito DML need a conversion-rate check before they are interpreted.
- **All 24 oysters is the primary analysis.** The only stated reason for dropping 2H-3 and 3H-2 ([04-methylKit.Rmd:199](04-methylKit.Rmd)) is that they sit apart on the PCA and have the lowest read counts (36–39M). No QC check flagged them. Low depth alone is not a reason to drop a sample in DSS, which weights each CpG by its coverage. Outliers removed is a sensitivity analysis, and the QC step (D3) decides whether that changes.

```r
design <- sampleMetadata %>%
  transmute(ploidy = factor(ploidy, levels = c("2N", "3N")),
            pH     = factor(pH,     levels = c("high", "low")))
```

## Analyses

### A. Factorial model (primary)

`DMLfit.multiFactor` + `DMLtest.multiFactor`, Wald tests per CpG.

- **Primary model:** `~ ploidy * pH` with sum-to-zero contrasts (`options(contrasts = c("contr.sum", "contr.poly"))`). Test `ploidy1` (ploidy effect averaged over pH) and `pH1` (pH effect averaged over ploidy). The interaction is controlled for, not tested as a result.
- **Sensitivity model:** additive `~ ploidy + pH`. This is the closest match to methylKit's "treatment + covariate" setup in 04.2. With a balanced design the estimates should be almost the same as the primary model; large differences mark CpGs where the interaction matters.
- The interaction coefficient is still saved, but only used as a diagnostic (how many CpGs show a strong interaction and whether those overlap the main-effect DML). It is not reported as its own DML set.
- **DMRs:** `callDMR()` on the multiFactor test output (p.threshold = 1e-5, minlen = 50, minCG = 3, dis.merge = 100, pct.sig = 0.5). `delta` is not available for multiFactor, so filter on the mean methylation difference computed afterward.
- Effect size: after testing, compute the group mean difference in % methylation for each significant CpG. Then the 25% / 50% filters can be applied as in 04.2.

### B. Two-group tests with smoothing (secondary)

`DMLtest(BSobj, group1, group2, smoothing = TRUE/FALSE)` → `callDML(delta = 0.25 / 0.10, p.threshold = 0.001)` → `callDMR(delta, p.threshold = 1e-5)`.

| Contrast | n | Use |
|---|---|---|
| 3N vs 2N, pooled over pH | 12 v 12 | smoothing and `delta`-based DML/DMR for ploidy |
| low vs high, pooled over ploidy | 12 v 12 | same for pH |
| 3N vs 2N within each pH; low vs high within each ploidy | 6 v 6 | consistency check only |

Only the two-group tests support smoothing and `delta` (posterior probability that |diff| > delta). They can't adjust for the other factor, but the design is balanced, so the other factor is evenly split between the groups. Use them to check that the DML/DMR from model A hold up with smoothing. Because the interaction is a nuisance term, the within-cell tests only check whether an effect has the same direction in both cells.

### C. Parameter grid (mirrors 04.2)

| Parameter | Levels | 04.2 equivalent |
|---|---|---|
| Lower coverage | ≥1x (DSS default), ≥5x | `lo.count = 5` |
| Upper coverage | none; 99.9th percentile per sample, from **nuclear CpGs only**, mito exempt | `hi.perc = 99.9` (drops all mito) |
| Sample presence | all 24; ≥11/10/9/8 per ploidy (or pH) group; **≥5 of 6 per cell** | `min.per.group` |
| Low-variation filter (`MIN_METH`) | none; mean methylation over all samples 10–90%, mito exempt (default since 2026-10-02) | — |
| Samples | all; outliers removed | OutRM |
| Smoothing (B only) | FALSE; TRUE, span 500 bp | — |
| Thresholds | FDR 0.05 / 0.01; \|Δ\| ≥ 25% / 50% (also 10%) | q 0.01, diff 25/50 |

The per-cell presence filter matters for the factorial model. A CpG with ≥8 of 12 per ploidy group can still have 0 samples in one ploidy × pH cell, and that coefficient is not estimable. Use per-cell presence for analysis A, and per-group presence only for comparing with methylKit.

Run the full grid on the main-effect and interaction tests, which gives about 2 coverage × 2 upper × 6 presence × 2 sample sets = 48 fits. Then narrow down before the smoothing and DMR runs.

### D. Calibration and null checks (deciding which counts to believe)

1. **p-value histograms** for each test and method (DSS, methylKit none, methylKit MN). A well-calibrated test has a flat distribution with a spike near 0. methylKit `none` should be left-skewed overall.
2. **Label permutations:** permute ploidy labels within pH (balanced, 6/6 per cell kept), and pH within ploidy. Run 20–50 permutations through the chosen setting, and use the number of "DML" in permuted data as an empirical false-positive rate. This is the clearest way to decide between methylKit `none` (thousands) and `MN` (~0).
3. **Sample QC (decides the outlier question):** make one table with one row per sample: raw/trimmed reads, Bismark mapping efficiency, duplication, non-CpG methylation (bisulfite conversion), CpGs at ≥5x, median coverage, global % CpG methylation, and PC1/PC2 scores. Use DSS-filtered data and redo the PCA with the mito included. Set the exclusion rule **before** looking at DML counts (e.g. conversion < 99%, mapping or CpG count > 3 MAD from the median, or a PCA position explained by one of those metrics). Exclude 2H-3/3H-2 only if they fail the rule. If they pass, all-24 remains primary. Also check 2H-1 and 3H-3, which 04-methylKit noted as possible outliers.
4. **Sample-set sensitivity (decided 2026-10-02 from the QC results):** all 24 is primary. Sensitivity runs:
   - `Drop3H2`: 3H-2 makes up most of PC1, has the lowest correlation with other samples (0.869), the highest global methylation (11.7%), and the lowest depth. It is flagged only on CpGs at ≥5x (robust z = -3.1).
   - `DropPC2`: 2H-1, 2H-2, 3H-5 separate on PC2. Their technical metrics are normal, and PC2 tracks global methylation (rho = -0.48).
   - `OutRM`: 2H-3 + 3H-2, the pair removed in 04.2/04.3, kept only for comparison with methylKit. 2H-3 passes every QC check and sits in the main cluster.
   For each, report how many all-24 DML are still significant.

### E. Concordance

- DSS vs methylKit 04.3 `none` / `MN` (matched filters): overlap counts, Jaccard index, correlation of methylation differences across shared CpGs.
- New DSS vs old 04-DSS DML (same coordinates).
- Agreement between the factorial (A) and simple-effect (B) results.
- Downstream (07/08/09): annotate whichever final DML/DMR set we choose with gene/exon/intron/TE features, then run enrichment. Hold off until D is done.

## Implementation

Built (2026-10-02), following the 04.3 pattern:

```
code/Haws/04.4-DSS-slurm.R             stages: prep-raw | qc | prep | fit | summary
code/Haws/04.4-DSS-slurm.job           srlab-R4.4 container; DSS 2.54.0 + bsseq 1.42.0 in ~/R_libs_apptainer
code/Haws/04.4-DSS-slurm-submit.sh     targets qc | dss | all; settings as env vars
analyses/Haws_04.4-DSS/                outputs (prep-*/, DSS-*/chunks/, DSS-*/rds/, logs/ gitignored)
```

Still to build: two-group smoothing tests (B), concordance (E), and an Rmd for figures.

**prep-raw:** reads the 24 `.cov` files once into integer count matrices (union of CpGs covered in any sample). The `qc` and `prep` stages for every filter setting use these matrices.

**prep:** for each sample, sets CpGs below `LO_COV`, or above the nuclear-only `HI_PERC` cutoff, to missing (mito is exempt from the upper cutoff), then applies the presence filter. Writes `CpG-filter-counts.tsv`, which includes how many mito CpGs pass.

**fit:** effect-coded design (ploidy +1/2 for 3N, -1/2 for 2N; pH +1/2 for low, -1/2 for high), so the ploidy and pH coefficients are averaged over the other factor even with the interaction in the model. A positive value means higher methylation in 3N or at low pH. A test on real data showed the fit is fast (0.5 s for 7k CpGs × 24), so the default is one task (`NCHUNK = 1`), and DSS estimates the dispersion prior from all CpGs. If `NCHUNK > 1`, the FDR is still recomputed with BH over all CpGs in `summary`. `PERM = <seed>` relabels samples (ploidy shuffled within pH, then pH within ploidy, keeping every group at the same size) for the permutation null.

**summary:** counts at FDR 0.05/0.01 × |Δ| ≥ 0/10/25/50% (`DSS-counts-long.csv`, with mito counts listed separately), DML BED files, `callDMR` DMRs (p < 1e-5), p-value histograms, a methylKit-style `DML-summary-table.csv`, and how many main-effect DML also have a strong interaction (diagnostic). Δ is the mean of the group means, the same averaging the model uses.

## Results so far (2026-10-02)

- **Sample QC** (`analyses/Haws_04.4-DSS/sample-QC/`): all 24 samples pass the technical checks (conversion 99.49–99.54%, mapping 60–62%, duplication 10–13%). 2H-3 is not an outlier. 3H-2 makes up most of PC1. 2H-1, 2H-2, and 3H-5 separate on PC2.
- **No low-variation filter** (3.53M CpGs): FDR < 0.01 and ≥25% difference give 8 ploidy and 4 pH DML, and 0 DMRs. The p-values are conservative (most near 1).
- **`MIN_METH = 10`** (559k CpGs): 187 ploidy and 185 pH DML at FDR < 0.05. The p-value histogram looks well calibrated.
- **Permutations** (20, `permutation-summary-*.csv`): shuffled labels give about as many DML as the real labels. Ploidy: 187 observed vs. a permuted mean of 187 (10 of 20 permutations ≥ observed). pH: 185 vs. 166 (5 of 20). So the empirical FDR is about 1, and **CpG-level DML are not distinguishable from noise** with 6 oysters per group. DSS's nominal FDR is too optimistic here. The permuted counts vary a lot (40–576), which suggests individual oysters drive most of the "DML".
- **100 permutations** (queued from the 04.6 session, same setting): ploidy 187 observed vs. a permuted mean of 223 (43 of 100 ≥ observed). This confirms the 20-permutation result.
- **SNP check** (`SNP-check/`; BS-SNPer per-sample VCFs from the earlier 05-BS-SNPer run, PASS calls only):
  - 424k of 12.8M covered CpGs (3.3%) have a SNP at the C or G in at least one oyster, and so do 12.4% of the 559k `MIN_METH` CpGs.
  - **SNP CpGs do not explain the null.** With `SNP_FILTER = any` (489k CpGs), ploidy gives 146 observed vs. a permuted mean of 153, and pH 150 vs. 138 (20 permutations). SNP CpGs are about equally common among observed and permuted DML (ploidy 18% vs. 15%).
  - **Genetic background accounts for much of the methylation difference between oysters.** Genotype PCA (13,840 SNP CpGs ≥10x in all samples) separates 2H-1, 2H-2, 3H-5, and 3H-2, the same oysters that stand out in the methylation PCA. Genotype PC1 tracks methylation PC2 (rho = -0.72). Genotype similarity vs. methylation similarity across sample pairs: rho = 0.86.
  - **The circularity check passes.** Using only SNPs bisulfite cannot mimic (no C>T at the C, no G>A at the G; 13,127 SNP CpGs), rho = 0.858 and the same oysters separate. Among the other 20 oysters, rho = 0.665, so the link is not only the four outliers.
  - **Genotype is not confounded with treatment.** Among the 20, genotype correlation is the same within and between treatment groups (0.380 vs. 0.380). The four outlier oysters are not close relatives of each other (pairwise genotype correlations at or below the median pair).
- **Genotype PC1 as a covariate** (`GENO_PCS = 1`; bisulfite-safe genotype PCA; 50 permutations that keep each oyster's genotype PC):
  - Ploidy: 265 observed vs. a permuted median of 212 (18 of 50 ≥ observed). Still no effect.
  - pH: 65 observed vs. a permuted median of 144 (47 of 50 > observed). The observed count is lower than chance, a sign of **confounding**.
  - **Correction to the SNP result above:** across all 24 oysters, genotype *is* confounded with pH. All four genetic outliers (2H-1, 2H-2, 3H-5, 3H-2) are high-pH oysters; if they were spread at random, that would happen about 5% of the time (495/10,626). Adjusting for genotype PC1 removes real-label pH signal that the permutations don't have, and the pH DML found without the covariate can't be separated from genetic background. The "not confounded" result holds only among the other 20.
  - So the cleanest pH test is `DropPC2` / `Drop3H2` with their own permutations (running).
- **Genetic outliers removed** (`MIN_METH = 10`, 50 permutations each; observed vs. permuted median, FDR < 0.05):

  | Sample set | CpGs | ploidy | pH | interaction |
  |---|---|---|---|---|
  | `DropPC2` (no 2H-1, 2H-2, 3H-5) | 398k | 205 vs. 145 (14/50 ≥) | 119 vs. 108 (23/50 ≥) | 187 vs. 128 (18/50 ≥) |
  | `Drop3H2` | 501k | 155 vs. 153 (24/50 ≥) | 110 vs. 166 (37/50 ≥) | 142 vs. 159 (26/50 ≥) |

  - **No ploidy or pH effect** at the CpG level once the genetic outliers are removed, so the pH confounding wasn't hiding a real effect.
  - The interaction count that was above most permutations with all 24 oysters (3/20 ≥ without the covariate, 6/50 with genotype PC1) is at the middle of the permutation distribution here, so it came from the outlier oysters.
  - **Conclusion for 04.4:** under every model and sample set tried (no covariate, SNP CpGs removed, genotype PC1, outliers removed), CpG-level DML counts are within the label-permutation distribution. No CpG-level ploidy or pH effect is detectable in this experiment.
- **Next:** check hatchery records for the four genetic outliers (all high pH); feature-level tests (04.6), also run without the outliers (gene body, promoter, 1 kb windows) with permutation-based FDR; 100 permutations for any CpG-level cutoff; SNP check (BS-SNPer).

## Order of work

0. ~~Install DSS/bsseq in the container library~~ Done: DSS 2.54.0 (`DMLfit.multiFactor` supports `smoothing`), bsseq 1.42.0.
1. Sample QC table (D3) from existing trimming/Bismark reports plus the coverage files. Set the exclusion rule and record the result for 2H-3/3H-2.
2. Run prep + the factorial model (A) at the default setting (≥5x, hi.perc 99.9 nuclear-only, all 24, per-cell presence). Get a first count.
3. First run (no low-variation filter): p-values piled up near 1 (conservative), because of the many CpGs at ~0%/~100% in every sample. Added `MIN_METH = 10`, which cuts 3.53M → 559k CpGs; the CpGs it removes can show at most ~20% differences.
4. Calibration (D1, D2) on that setting, against methylKit `none`/`MN`. **Decision point:** choose the method and threshold that hold up under permutation.
5. Full parameter grid (C) and two-group tests (B).
6. DMRs, concordance (E), and figures (heatmaps, chromosome distribution, as in 04-DSS; add mito as its own bar).
7. Pass the final DML/DMR set to 07–09.

## Open questions

- Unplaced scaffolds (225, ~1–2% of CpGs): keep for now, since they are cheap to include and can be dropped at the figure stage.
- 04.3 methylKit currently drops all mito CpGs through `hi.perc`. For the concordance step (E) to compare like with like, either fix 04.3 the same way or leave mito out of the comparison.