# ChIP-qPCR, drug vs vehicle at one target region: percent input and fold enrichment over IgG, analysed four ways.
# R 4.3.3, base R.
d <- read.csv("chip.csv")                            # one row per qPCR well: experiment, condition, sample (input, IP, IgG), well, ct
f2 <- function(p) formatC(signif(p, 2), format = "g")
DF <- log2(100)                                      # the input is 1% of the chromatin: subtract log2(100) = 6.644 cycles
m <- aggregate(ct ~ experiment + condition + sample, data = d, FUN = mean)
m <- reshape(m, idvar = c("experiment", "condition"), timevar = "sample", direction = "wide")
m$pct_input <- 100 * 2^((m$ct.input - DF) - m$ct.IP); m$pct_igg <- 100 * 2^((m$ct.input - DF) - m$ct.IgG); m$fold_igg <- 2^(m$ct.IgG - m$ct.IP)
ex <- sort(unique(d$experiment)); V <- m[m$condition == "vehicle", ]; G <- m[m$condition == "drug", ]
V <- V[order(V$experiment), ]; G <- G[order(G$experiment), ]
for (i in seq_along(ex)) cat(sprintf("Experiment %d: %% input vehicle %.2f, drug %.2f (ratio %.2f); IgG %% input %.3f and %.3f; fold over IgG vehicle %.0f, drug %.0f\n",
                                     ex[i], V$pct_input[i], G$pct_input[i], G$pct_input[i] / V$pct_input[i], V$pct_igg[i], G$pct_igg[i], V$fold_igg[i], G$fold_igg[i]))
w <- merge(d[d$sample == "IP", ], m[, c("experiment", "condition", "ct.input")], by = c("experiment", "condition"))
w$pct <- 100 * 2^((w$ct.input - DF) - w$ct)          # one % input per IP well
a <- w$pct[w$condition == "drug"]; b <- w$pct[w$condition == "vehicle"]; t <- t.test(a, b)
cat(sprintf("A. qPCR wells as n (%d vs %d), Welch t-test on %% input: drug / vehicle %.2f, p = %s\n", length(a), length(b), mean(a) / mean(b), f2(t$p.value)))
t <- t.test(log2(G$pct_input), log2(V$pct_input))
cat(sprintf("B. Experiments as n but unpaired (4 vs 4), Welch t-test on log2 %% input: drug / vehicle %.2f, p = %s\n", 2^(mean(log2(G$pct_input)) - mean(log2(V$pct_input))), f2(t$p.value)))
r <- log2(G$pct_input / V$pct_input); t <- t.test(r)
cat(sprintf("C. Log2 ratio of %% input per experiment (paired, n = %d): drug / vehicle %.2f (95%% CI %.2f to %.2f), p = %s\n", length(r), 2^mean(r), 2^t$conf.int[1], 2^t$conf.int[2], f2(t$p.value)))
r <- log2(G$fold_igg / V$fold_igg); t <- t.test(r)
cat(sprintf("D. Log2 ratio of fold enrichment over IgG per experiment (paired, n = %d): drug / vehicle %.2f (95%% CI %.2f to %.2f), p = %s\n", length(r), 2^mean(r), 2^t$conf.int[1], 2^t$conf.int[2], f2(t$p.value)))