suppressPackageStartupMessages({library(DSS); library(data.table)})
b <- deparse(callDMR); cat(head(grep("class|multiFactor|delta|stat", b, value = TRUE), 15), sep = "\n")
dataDir <- "/mmfs1/gscratch/scrubbed/sr320/github/project-oyster-oa/data/Haws"
f <- file.path(dataDir, paste0("zr3644_", 1:24, "_R1_val_1_val_1_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov"))
d <- lapply(f, function(x) { z <- fread(cmd = paste("awk '$1==\"NC_047559.1\" && $2 < 3000000'", shQuote(x)), select = c(1,2,5,6)); setnames(z, c("chr","pos","X","U")); z[, .(chr, pos, N = X + U, X)] })
BS <- makeBSseqData(d, paste0("s", 1:24))
cv <- getCoverage(BS); keep <- rowSums(cv >= 5) == 24; BS <- BS[keep][1:min(20000, sum(keep))]
cat("CpGs in test:", length(BS), "\n")
ploidy <- rep(c("2N","3N"), each = 12); pH <- rep(rep(c("high","low"), each = 6), 2)
design <- data.frame(ploidyEffect = ifelse(ploidy == "3N", 0.5, -0.5), pHEffect = ifelse(pH == "low", 0.5, -0.5))
t0 <- Sys.time(); fit <- DMLfit.multiFactor(BS, design = design, formula = ~ ploidyEffect * pHEffect); el <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
cat("fit seconds:", el, " per 150k CpGs (min):", el / length(BS) * 150000 / 60, "\n")
print(colnames(fit$X))
r <- DMLtest.multiFactor(fit, coef = "ploidyEffect"); print(class(r)); print(head(r, 3))
cat("p<0.001:", sum(r$pvals < 0.001, na.rm=TRUE), " NA p:", sum(is.na(r$pvals)), "\n")
dmr <- callDMR(r, p.threshold = 1e-3); print(class(dmr)); print(head(dmr, 3))
r2 <- data.frame(chr = r$chr, pos = r$pos, stat = r$stat, pvals = r$pvals, fdrs = r$fdrs); class(r2) <- class(r)
dmr2 <- callDMR(r2, p.threshold = 1e-3); cat("rebuilt input gives same DMRs:", identical(dim(dmr), dim(dmr2)), "\n")