Weighted pseudo-bulk samples

weighted.R

library(edgeR)
g <- factor(rep(1:4, c(50, 20, 10, 5)))
N <- 20000
mus <- 100
# mus <- 100 * 2^rnorm(N*length(g)) # Uncomment for variable mu
y <- matrix(rnbinom(N*length(g), mu=mus, size=1), nrow=N, byrow=TRUE)

design <- model.matrix(~gl(2,2))

# Summation.
summed <- sumTechReps(y, g)
d <- DGEList(summed)
d <- estimateDisp(d, design)
fit <- glmQLFit(d, design)
res <- glmQLFTest(fit)
findInterval(c(0.001, 0.01, 0.05), sort(res$table$PValue))/N

# Averaging.
n <- as.integer(table(g))
averaged <- sweep(summed, 2, n, "/")
d2 <- DGEList(averaged)
d2$weights <- matrix(n, ncol=length(n), nrow=nrow(y), byrow=TRUE)

d2 <- estimateDisp(d2, design)
fit2 <- glmQLFit(d2, design)
res2 <- glmQLFTest(fit2)
findInterval(c(0.001, 0.01, 0.05), sort(res2$table$PValue))/N