QL.results: Obtain p- and q-values using results from 'QL.fit'

In emilygoren/BinQuasi: Analyzing Replicated ChIP Sequencing Data Using Quasi-Likelihood

Description

Obtain significance results for quasi-likelihood models fit to ChIP-seq data partitioned into counts.

Usage

QL.results(fit, Dispersion = "Deviance", one.sided = TRUE,
  spline.df = NULL, Plot = FALSE, padj = TRUE)

Arguments

fit	The list returned by the function QL.fit.
Dispersion	Must be one of "Deviance" or "Pearson", specifying which type of estimator should be used for estimating the quasi-likelihood dispersion parameter.
one.sided	logical. If TRUE, a one-sided test for the ChIP coefficient is reported. Otherwise, if FALSE, a two-sided test is reported.
spline.df	Optional. User may specify the degrees of freedom to use when fitting a cubic spline to log-scale(estimated dispersion) versus the log(average count). Default uses cross-validation in sreg function to pick appropriate degrees of freedom.
Plot	logical. If TRUE, the estimated dispersion versus the average count are plotted on a log-scale with the corresponding cubic spline fit overlaid.
padj	logical. If TRUE, Benjamini & Hochberg's adjustment for multiple comparisons is applied.

Details

Obtain significance results from an object fitted using QL.fit. Used within the main peak calling function, BQ.

Value

A list containing:

P.values	List of matrices providing p-values for the QL, QLShrink and QLSpline methods, respectively. The i-th column of each element of pvals corresponds to the hypothesis test assigned in the i-th window.
Q.values	List of matrices providing q-values for the QL, QLShrink and QLSpline methods, respectively. The i-th column of each element of qvals corresponds to the hypothesis test assigned in the i-th window.
F.stat	List of matrices providing F-statistics for the QL, QLShrink and QLSpline methods, respectively. The i-th column of each element of F.stat corresponds to the hypothesis test assigned in the i-th window.
d0	Vector containing estimated additional denominator degrees of freedom gained from shrinking dispersion estimates in the QLShrink and QLSpline procedures, respectively.

Author(s)

Emily Goren ([e-mail]) based on modifications of code by Steve Lund.

References

Goren, Liu, Wang and Wang (2018) "BinQuasi: a peak detection method for ChIP-sequencing data with biological replicates" Bioinformatics.

Benjamini and Hochberg (1995) "Controlling the false discovery rate: a practical and powerful approach to multiple testing" Journal of the Royal Statistical Society Series B, 57: 289-300.

Lund, Nettleton, McCarthy and Smyth (2012) "Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates" SAGMB, 11(5).

Examples

 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

emilygoren/BinQuasi documentation built on May 17, 2019, 12:08 p.m.