---
title: "01 - Data download and sample metadata"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(eval = FALSE)
library(data.table)
```

# Coverage files

Bismark CpG-report coverage files (one per sample, 24 total) are the input.
They are hosted at:

`https://gannet.fish.washington.edu/spartina/project-oyster-oa/Haws/bismark-2/`

with names of the form
`zr3644_<i>_R1_val_1_val_1_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov`
for `i = 1..24`. Each is ~400-450 MB and holds ~10 million CpG positions in the
6-column Bismark format: `chrom start end pct_meth count_methylated count_unmethylated`.

```{r download}
dir.create("data/cov", recursive = TRUE, showWarnings = FALSE)
base <- "https://gannet.fish.washington.edu/spartina/project-oyster-oa/Haws/bismark-2/"
for (i in 1:24) {
  fn <- sprintf("zr3644_%d_R1_val_1_val_1_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov", i)
  dest <- file.path("data/cov", fn)
  if (!file.exists(dest))
    download.file(paste0(base, fn), dest, mode = "wb", quiet = TRUE)
}
```

# Sample metadata

The design is a 2x2 factorial: **ploidy** (2N/3N) x **pH** (high/low), n = 6 per
cell, 24 samples. Metadata is reconstructed from the analysis code that
accompanied the data (the original sample sheet was not distributed publicly).
Sample IDs follow `{ploidy}{pH}-{rep}`: 2H = 2N/high, 2L = 2N/low, 3H = 3N/high,
3L = 3N/low.

```{r metadata}
cov_files <- sprintf(
  "data/cov/zr3644_%d_R1_val_1_val_1_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov",
  1:24)

sampleID <- c(paste0("2H-", 1:6), paste0("2L-", 1:6),
              paste0("3H-", 1:6), paste0("3L-", 1:6))

meta <- data.table(
  sample_number = 1:24,
  sampleID      = sampleID,
  cov_file      = cov_files,
  ploidy        = rep(c("2N", "3N"), each = 12),
  pH            = rep(c("high", "low", "high", "low"), each = 6)
)
# methylKit treatment codes: encode each of the four groups as an integer
meta[, group := paste0(ploidy, "_", pH)]
meta[, group_code := as.integer(factor(group,
      levels = c("2N_high", "2N_low", "3N_high", "3N_low"))) - 1L]  # 0,1,2,3

fwrite(meta, "sample_metadata.csv")
print(meta)
```

`group_code` (0 = 2N_high, 1 = 2N_low, 2 = 3N_high, 3 = 3N_low) is the
methylKit `treatment` vector for the full 24-sample object. Two-group contrasts
later subset this and pass their own 0/1 treatment.