# NOT RUN {
set.seed(5)
####################################################
# Simulate data three replicates with one chromosome
####################################################
reps <- 3
chr.length <- 1e5
window.width <- 200
K <- chr.length / window.width
start <- seq(1, chr.length, by = window.width)
end <- start + window.width - 1
n.peaks <- 100 # No. of true peaks
peak.idx <- sample.int(K, n.peaks)
samples <- c(paste0('C', 1:reps), paste0('I', 1:reps))
# Set parameters
beta0 <- runif(K, log(10), log(100))
beta1 <- rep(0, K); beta1[peak.idx] <- log(5) / runif(n.peaks)^(1/5)
# Set means
mu.ChIP <- exp(beta0 + beta1)
mu.input <- exp(beta0)
# Negative binomial dispersion parameter
phi <- 1/rchisq(K, df = 5)
# Draw read counts using a negative binomial distribution
C <- lapply(1:reps, function(r) rpois(K, (mu.ChIP * rgamma(K, 1/phi))/(1/phi)))
I <- lapply(1:reps, function(r) rpois(K, (mu.input * rgamma(K, 1/phi))/(1/phi)))
counts <- do.call('cbind', append(C, I))
colnames(counts) <- samples
rownames(counts) <- start
####################################################
# Fit quasi-negative binomial model to each window.
####################################################
design.matrix <- cbind(rep(1, reps*2), # Intercept
rep(c(1,0), each = reps)) # Indicates ChIP sample
chip.col.indicator <- c(0,1) # Second column of design matrix indicates ChIP sample
fit <- QL.fit(counts, design.matrix, chip.col.indicator,
log.offset = rep(1, ncol(counts)), Model = 'NegBin')
window.results <- QL.results(fit)
# Number of significant windows.
sum(window.results$Q.values$QLShrink < 0.05)
# Compare significant windows to truth.
res <- as.numeric(window.results$Q.values$QLShrink < 0.05)
# Number of true positives
TP <- sum(res[peak.idx] == 1)
TP
# Number of false negatives
FN <- n.peaks - TP
FN
# Number of false positives
FP <- sum(res) - TP
FP
# Number of true negatives
TN <- (K - sum(res)) - FN
TN
# }
Run the code above in your browser using DataCamp Workspace