R/ZicoSeq.R
In chloelulu/ZicoSeq: Has not been decided
Defines functions
Ztransform
corHerit
# 2019_02_20
# Beta-mixture to address sequence depth confounding
# Permutation-based FDR control
# Multiple-stage normalization
# Variance modulation by LIMMA - hasn't been finished
# Expectation of F statistics based on posterior sampling
# Omnibus test for combining diffrent effects
# Version 0.9 [PermuteDAA3.2]
#' Computes the GMPR size factor
#' @param comm a matrix of counts, row - features (OTUs, genes, etc) , column - sample
#' @param intersect.no the minimum number of shared features between sample pair, where the ratio is calculated
#' @param ct.min the minimum number of counts required to calculate ratios ct.min = 5 has better results
#' @param verbose suppress output
#' @return A list that contains:
#' \item{gmpr}{the GMPR size factors for all samples; Samples with distinct sets of features will be output as NA.}
#' \item{nss}{number of samples with significant sharing (> intersect.no) including itself}
#' @rdname GMPR
#' @export
GMPR <- function (comm, intersect.no = 4, ct.min = 2, verbose = FALSE) {
# mask counts < ct.min
comm[comm < ct.min] <- 0
if (is.null(colnames(comm))) {
colnames(comm) <- paste0('S', 1:ncol(comm))
}
if (verbose == TRUE)
cat('Begin GMPR size factor calculation ...\n')
comm.no <- numeric(ncol(comm))
gmpr <- sapply(1:ncol(comm), function(i) {
if (i %% 50 == 0) {
if (verbose == TRUE)
cat(i, '\n')
}
x <- comm[, i]
# Compute the pairwise ratio
pr <- x / comm
# Handling of the NA, NaN, Inf
pr[is.nan(pr) | !is.finite(pr) | pr == 0] <- NA
# Counting the number of non-NA, NaN, Inf
incl.no <- colSums(!is.na(pr))
# Calculate the median of PR
pr.median <- colMedians(pr, na.rm=TRUE)
# Record the number of samples used for calculating the GMPR
comm.no[i] <<- sum(incl.no >= intersect.no)
# Geometric mean of PR median
if (comm.no[i] > 1) {
return(exp(mean(log(pr.median[incl.no >= intersect.no]))))
} else {
return(NA)
}
}
)
if (sum(is.na(gmpr))) {
warning(paste0('The following samples\n ', paste(colnames(comm)[is.na(gmpr)], collapse='\n'),
'\ndo not share at least ', intersect.no, ' common taxa with the rest samples! ',
'For these samples, their size factors are set to be NA! \n',
'You may consider removing these samples since they are potentially outliers or negative controls!\n',
'You may also consider decreasing the minimum number of intersecting taxa and rerun the procedure!\n'))
}
if (verbose == TRUE) {
cat('Completed!\n')
cat('Please watch for the samples with limited sharing with other samples based on NSS! They may be outliers! \n')
}
attr(gmpr, 'NSS') <- comm.no
names(gmpr) <- colnames(comm)
return(gmpr * median(colSums(comm)))}
####################################################################################
# Copied from "vegan" package getPermuteMatrix.R to reduce dependency
getPermuteMatrix <- function (perm, N, strata = NULL) {
if (length(perm) == 1) {
perm <- how(nperm = perm)
}
if (!missing(strata) && !is.null(strata)) {
if (inherits(perm, "how") && is.null(getBlocks(perm)))
setBlocks(perm) <- strata
}
if (inherits(perm, "how"))
perm <- shuffleSet(N, control = perm)
if (is.null(attr(perm, "control")))
attr(perm, "control") <- structure(list(within = list(type = "supplied matrix"),
nperm = nrow(perm)), class = "how")
perm
}
# Permutation-based FDR control
perm_fdr_adj <- function (F0, Fp) {
ord <- order(F0, decreasing = T)
F0 <- F0[ord]
perm.no <- ncol(Fp)
Fp <- as.vector(Fp)
Fp <- Fp[!is.na(Fp)]
Fp <- sort(c(Fp, F0), decreasing = F)
n <- length(Fp)
m <- length(F0)
FPN <- (n + 1) - match(F0, Fp) - 1:m
p.adj.fdr <- FPN / perm.no / (1:m)
# p.adj.fdr <- sapply(F0, function(x) sum(Fp >=
# x, na.rm=TRUE) / perm.no)/(1:length(F0))
p.adj.fdr <- pmin(1, rev(cummin(rev(p.adj.fdr))))[order(ord)]
}
# Permutation-based FWER control
# Conservative
perm_fwer_adj <- function (F0, Fp) {
ord <- order(F0, decreasing = T)
F0 <- F0[ord]
col.max <- colMaxs(Fp, na.rm=TRUE)
p.adj.fwer <- sapply(F0, function(x) mean(col.max >= x))[order(ord)]
}
# Westfall young too slow
perm_fwer_adj2 <- function (F0, Fp) {
ord <- order(F0, decreasing = T)
m <- length(F0)
F0 <- F0[ord]
Fp <- Fp[ord, , drop=FALSE]
col.max <- Fp[m, ]
p.adj.fwer <- sapply(m:1, function(i) {
x <- F0[i]
y <- Fp[i, ]
col.max <<- ifelse(y > col.max, y, col.max)
mean(col.max >= x)
})
p.adj.fwer <- rev(p.adj.fwer)
p.adj.fwer <- pmin(1, rev(cummin(rev(p.adj.fwer))))[order(ord)]
}
# Pad back the NA values
na.pad <- function (vec, ind) {
vec0 <- numeric(length(ind))
vec0[!ind] <- vec
vec0[ind] <- NA
vec0
}
# # Pad back the NA values
one.pad <- function (vec, ind) {
vec0 <- numeric(length(ind))
vec0[!ind] <- vec
vec0[ind] <- 1
vec0
}
# Transform to Z score
Ztransform <- function(p.value, e.sign, eff.sign=TRUE, tol=1E-15) {
p.value[p.value <= tol] <- tol
p.value[p.value >= 1 - tol] <- 1 - tol
if (eff.sign == TRUE) {
e.sign[e.sign == 0] <- sample(c(-1, 1), sum(e.sign == 0), replace=T)
z1 <- qnorm(p.value / 2)
z2 <- qnorm(1 - p.value / 2)
z <- ifelse(e.sign > 0, z2, z1)
} else {
z <- qnorm(1 - p.value)
}
return(z)
}
# Internal functions
# For gls with generic correlation structure
corHerit <- function(value, paras, form = ~1, fixed = TRUE) {
# Place holder - check the validity of the parameter
object <- value
attr(object, "formula") <- form
attr(object, "fixed") <- fixed
attr(object, "paras") <- paras
class(object) <- c("corHerit", "corStruct")
object
}
Initialize.corHerit <- function (object, data, ...) {
# Place holder - check the validity of the parameter
form <- formula(object)
if (!is.null(getGroupsFormula(form))) {
attr(object, "groups") <- getGroups(object, form, data = data)
attr(object, "Dim") <- Dim(object, attr(object, "groups"))
} else {
attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
}
attr(object, "covariate") <- getCovariate(object, data = data)
object
}
corMatrix.corHerit <- function (object, covariate = getCovariate(object), ...) {
paras <- attr(object, "paras")
p <- paras[['p']]
I <- diag(p)
hyper <- as.vector(object)
lambda <- 1 / (1 + exp(hyper[1]))
V <- exp(-((paras[['D']])) * (exp(hyper[2])))
cor.m <- (1 - lambda) * V + lambda * I
cor.m
}
# Extract the coefficient
coef.corHerit <- function (object, ...) {
paras <- attr(object, "paras")
coefs <- as.vector(object)
coef1 <- 1 / (1 + exp(coefs[1]))
coef1 <- coef1 / (1 - coef1)
coef2 <- exp(coefs[2])
coefs <- c(coef1, coef2)
names(coefs) <- paste("Hyper", 1:length(coefs), sep="")
coefs
}
# Estimating the hyper parameters
EstHyper <- function (y, D, init.val=c(0, 0)) {
obj.gls <- gls(model = y ~ 1,
correlation = corHerit(value=init.val, paras = list(p=length(y), D = D)))
cc <- c(coef(obj.gls$modelStruct$corStruct), obj.gls$coefficients, obj.gls$logLik)
cc
}
AdjStats <- function (y, V, k, mu, fudge=0.005) {
p <- nrow(V)
# Add small fudge
V.inv <- solve(V + fudge * diag(p))
I <- diag(p)
y.adj <- solve(I + k * V.inv) %*% (k * mu * rowSums(V.inv) + y)
y.adj
}
# Compute the NFP on the null.ind
PermFDR <- function (F0, Fp, null.ind) {
Fp <- as.matrix(Fp)
pct <- sum(null.ind) / length(F0)
Fp <- Fp[null.ind, ]
ord <- order(F0, decreasing = T)
F0 <- F0[ord]
perm.no <- ncol(Fp)
Fp <- as.vector(Fp)
pct.non.na <- mean(!is.na(Fp))
Fp <- Fp[!is.na(Fp)]
Fp <- sort(c(Fp, F0), decreasing = F)
n <- length(Fp)
m <- length(F0)
FPN <- (n + 1) - match(F0, Fp) - 1:m
# Handle identical F0 situations
FPN <- cummax(FPN)
p.adj.fdr <- FPN / pct.non.na / pct / perm.no / (1:m)
p.adj.fdr <- pmin(1, rev(cummin(rev(p.adj.fdr))))[order(ord)]
}
####################################################################################
# Compute tree-based FDR control
TreeFDR <- function (X = NULL, Y = NULL, test.func = NULL, perm.func = NULL, test.obs = NULL, test.perm = NULL, tree = NULL, D,
eff.sign = TRUE, B = 20, q.cutoff = 0.5, alpha = 1,
adaptive = c('Fisher', 'Overlap'), alt.FDR = c('BH', 'Permutation'), verbose = TRUE, ...) {
if (is.null(X) | is.null(Y) | is.null(test.func) | is.null(perm.func)) {
if (is.null(test.obs) | is.null(test.perm)) {
stop('Please specify either "X, Y, test.func, perm.fuc", or "test.obs", "test.perm"!\n')
} else {
X <- matrix(test.obs$p.value, ncol = 1)
rownames(X) <- names(test.obs$p.value)
}
}
adaptive <- match.arg(adaptive)
alt.FDR <- match.arg(alt.FDR)
# Make sure the rows of X and tree tips are labeled to avoid error
if ((is.null(rownames(X)) | is.null(tree$tip.label))) {
warning('Both the data matrix and the tree should have labels (rownames, tip.label) to avoid potential errors!\n')
} else {
if (sum(!(rownames(X) %in% tree$tip.label))){
stop('Some features in the data matrix are not in the tree! Please check!\n')
} else {
if (sum(!(tree$tip.label %in% rownames(X)))) {
warning('The tree have more features than the data matrix! \n')
}
}
}
# Patristic distance
if (!is.null(tree)) {
D <- (cophenetic(tree)) ^ alpha
} else {
D <- D ^ alpha
}
if (!is.null(rownames(X)) & !is.null(tree$tip.label)) {
D <- D[rownames(X), rownames(X)]
}
if (is.null(test.obs)) {
if (verbose)
cat('Test on original data sets ...\n')
test.obs <- test.func(X, Y, ...)
}
if (!is.list(test.obs) | !all(c("e.sign", "p.value") %in%
names(test.obs))) {
stop("test.func should return a list with names e.sign and p.valueif z.transform=TRUE! Please check!\n")
}
null.ind <- test.obs$p.value >= quantile(test.obs$p.value, 1 - q.cutoff)
z.obs <- Ztransform(test.obs$p.value, test.obs$e.sign, eff.sign)
if (verbose)
cat('Test on permuted data sets ...\n')
if (is.null(test.perm)) {
z.perm <- z.perm2 <- matrix(NA, nrow(X), B)
for (i in 1:B) {
perm.obj <- perm.func(X, Y, ...)
if (!is.list(perm.obj) | !all(c("X", "Y") %in% names(perm.obj))) {
stop("perm.func should return a list with names X and Y! Please check!\n")
}
X.perm <- perm.obj$X
Y.perm <- perm.obj$Y
test.perm <- test.func(X.perm, Y.perm, ...)
z.perm[, i] <- z.perm2[, i] <- Ztransform(test.perm$p.value, test.perm$e.sign, eff.sign)
z.perm2[!null.ind, i] <- z.obs[!null.ind]
}
} else {
B <- ncol(test.perm$p.value)
z.perm <- z.perm2 <- matrix(NA, length(test.obs$p.value), B)
for (i in 1:B) {
z.perm[, i] <- z.perm2[, i] <- Ztransform(test.perm$p.value[, i], test.perm$e.sign[, i], eff.sign)
z.perm2[!null.ind, i] <- z.obs[!null.ind]
}
}
if (alt.FDR == 'Permutation') {
if (verbose)
cat('Perform ordinary permutation-based FDR control ...\n')
if (eff.sign == TRUE) {
p.adj0 <- PermFDR(abs(z.obs), abs(z.perm), rep(TRUE, length(z.obs)))
} else {
p.adj0 <- PermFDR(z.obs, z.perm, rep(TRUE, length(z.obs)))
}
}
if (alt.FDR == 'BH') {
if (verbose)
cat('Perform ordinary BH-based FDR control ...\n')
p.adj0 <- p.adjust(test.obs$p.value, 'fdr')
}
if (verbose)
cat("Estimating hyperparameter ... \n")
error <- try(obj <- EstHyper(y = z.obs, D = D))
if (inherits(error, "try-error")) {
if (verbose)
cat('Hyperparameter estimation failed! Ordinary permutation-based FDR control will be used!\n')
# p.adj <- p.adjust(test.obs$p.value,'fdr')
p.adj <- p.adj0
k <- NULL
rho <- NULL
z.adj <- NULL
} else {
if (verbose)
cat("Structure-based adjustment ...\n")
k <- obj[1]
rho <- obj[2]
mu <- obj[3]
V <- exp(-1 * rho * D)
z.adj <- stat.o <- AdjStats(y = z.obs, V = V, k = k, mu = mu)[, 1]
stat.p <- AdjStats(y = z.perm2, V = V, k = k, mu = mu)
if (eff.sign == TRUE) {
stat.o <- abs(stat.o)
stat.p <- abs(stat.p)
}
p.adj <- PermFDR(stat.o, stat.p, null.ind)
# Power loss check
if (adaptive == 'Overlap') {
# Magic numbers
fdr.cutoff <- 0.2
pct.cutoff <- 0.5
ind0 <- p.adj0 <= fdr.cutoff
ind <- p.adj <= fdr.cutoff
if (sum(p.adj[ind0] <= fdr.cutoff) < pct.cutoff * sum(ind0)) {
# Over-adjustment checking
if (verbose)
cat('Potential over-adjustment! Alternative FDR control will be used!\n')
p.adj <- p.adj0
k <- NULL
rho <- NULL
z.adj <- NULL
}
}
if (adaptive == 'Fisher') {
# These cutoffs are used emprically
fdr.cutoff <- 0.2
ind0 <- p.adj0 <= fdr.cutoff
ind <- p.adj <= fdr.cutoff
n <- nrow(X)
test.p <- fisher.test(matrix(c(sum(ind), sum(ind0), n - sum(ind), n - sum(ind0)), 2, 2), alternative = 'less')$p.value
if (test.p <= 0.05) {
# Over-adjustment checking
if (verbose)
cat('Potential over-adjustment! Alternative FDR control will be used!\n')
p.adj <- p.adj0
k <- NULL
rho <- NULL
z.adj <- NULL
}
}
}
if (verbose)
cat("Done!\n")
return(list(p.adj = p.adj, p.unadj = test.obs$p.value, z.adj = z.adj, z.unadj = z.obs, k = k, rho = rho))
}
###################################################################################
# add BBmix
###############################################################################################
bbmix.fit.MM <- function (ct, dep, nIter = 10, winsor.qt = 1.0) {
if (mean(ct == 0) >= 0.95 | sum(ct != 0) < 5) {
stop('The number of nonzero is too small to fit the model! Consider removing it from testing!\n')
}
# Initialization
prop0 <- ct / dep
qt <- quantile(prop0, winsor.qt)
prop0[prop0 >= qt] <- qt
var1 <- var(prop0)
mean1 <- mean(prop0) / 2
var2 <- var1
mean2 <- 3 * mean(prop0) / 2
pi <- 0.5
for (i in 1 : nIter) {
shape1.1 <- ((1 - mean1) / var1 - 1 / mean1) * mean1 ^ 2
shape1.2 <- shape1.1 * (1 / mean1 - 1)
shape2.1 <- ((1 - mean2) / var2 - 1 / mean2) * mean2 ^ 2
shape2.2 <- shape2.1 * (1 / mean2 - 1)
m1 <- shape1.1 / (shape1.1 + shape1.2)
s1 <- shape1.1 + shape1.2
m2 <- shape2.1 / (shape2.1 + shape2.2)
s2 <- shape2.1 + shape2.2
f1 <- dbetabinom(ct, dep, m1, s1)
f2 <- dbetabinom(ct, dep, m2, s2)
q1 <- pi * f1 / (pi * f1 + (1 - pi) * f2)
q2 <- 1 - q1
pi <- mean(q1)
# Rough estimation
mean1 <- sum(prop0 * q1) / sum(q1)
var1 <- sum((prop0 - mean1)^2 * q1) / sum(q1)
mean2 <- sum(prop0 * q2) / sum(q2)
var2 <- sum((prop0 - mean2)^2 * q2) / sum(q2)
}
return(list(shape1.1 = shape1.1, shape1.2 = shape1.2, shape2.1 = shape2.1, shape2.2 = shape2.2, pi = pi, q1 = q1))
}
###############################################################################################
zibb.fit.MM <- function (ct, dep, nIter = 10, winsor.qt = 1.0) {
# Initialization
if (mean(ct == 0) >= 0.95 | sum(ct != 0) < 5) {
stop('The number of nonzero is too small to fit the model! Consider removing it from testing!\n')
}
ct1 <- ct[ct != 0]
dep1 <- dep[ct != 0]
# Initialization
prop0 <- ct1 / dep1
qt <- quantile(prop0, winsor.qt)
prop0[prop0 >= qt] <- qt
var0 <- var(prop0)
mean0 <- mean(prop0)
shape1 <- ((1 - mean0) / var0 - 1 / mean0) * mean0 ^ 2
shape2 <- shape1 * (1 / mean0 - 1)
m0 <- shape1 / (shape1 + shape2)
s0 <- shape1 + shape2
if (s0 < 0) {
s0 <- 1
}
m2 <- 0.75 * m0
s2 <- s0
pi <- 0.75 * mean(ct == 0)
##########################
prop0 <- ct / dep
qt <- quantile(prop0, winsor.qt)
prop0[prop0 >= qt] <- qt
##########################
for (i in 1 : nIter) {
# Expectation step
f1 <- ifelse(ct == 0, 1, 0)
f2 <- dbetabinom(ct, dep, m2, s2)
q1 <- pi * f1 / (pi * f1 + (1 - pi) * f2)
q2 <- 1 - q1
# Maximization step
pi <- mean(q1)
mean2 <- sum(prop0 * q2) / sum(q2)
var2 <- sum((prop0 - mean2)^2 * q2) / sum(q2)
shape2.1 <- ((1 - mean2) / var2 - 1 / mean2) * mean2 ^ 2
shape2.2 <- shape2.1 * (1 / mean2 - 1)
m2 <- shape2.1 / (shape2.1 + shape2.2)
s2 <- shape2.1 + shape2.2
}
shape2.1 <- m2 * s2
shape2.2 <- s2 - shape2.1
return(list(shape2.1 = shape2.1, shape2.2 = shape2.2, pi = pi, q1 = q1))
}
####################################################################################
#' Permutation-based differential abundance analysis
#' @param meta.dat a data frame containing the sample information
#' @param comm a matrix of counts, row - features (OTUs, genes, etc) , column - sample
#' @param grp.name a character, variable of interest; it could be numeric or categorical; should be in "meta.dat"
#' @param adj.name a character vector, variable(s) to be adjusted; they could be numeric or categorical; should be in "meta.dat"
#' @param prev.filter features with prevalence (i.e., nonzero proportion) less than "prev.cutoff" or be filtered
#' @param abund.filter features with a total counts less than "abund.cutoff" or be filtered
#' @param is.winsor a logical value indicating whether winsorization should be performed to replace outliers. The default is TRUE.
#' @param winsor.qt the winsorization quantile, above which the counts will be replaced
#' @param is.prior a logical value indicating whether to perform posterior inference based on some prior distribution on the proportion data
#' @param prior.dist prior distribution, either two-component beta-binomial mixture ("BetaMix") or zeroinflated beta-binomial ("ZIBB")
#' @param post.method method for posterior inference, either based on posterior sampling ("sample") or approximate posterior mean ("mean")
#' @param post.sample.no the number of posterior samples if posterior sampling is used
#' @param link.func a list of functions that connects the ratios to the covariates
#' @param link.d.func a list of the derivative function of "link.func"; only need to specifiy when "post.method" is "mean"
#' @param variance.EB a logical value indicating whehter to perform empirical Bayes based variance shrinkage
#' @param df.prior the degree of freedom of the prior inverse gamma distribution for variance shrinkage
#' @param perm.no the number of permutations; If the raw p values are of the major interest, set "perm.no" to at least 999
#' @param strata a factor indicating the permutation strata; permutation will be confined to each stratum
#' @param stats.combine.func function to combine the F-statistic for the omnibus test
#' @param stage.no the number of stages if multiple-stage ratio stategy is used
#' @param topK the number of dominant features that will be excluded in the initial stage ratio calculation
#' @param stage.fdr the fdr cutoff below which the features will be excluded for calculating the ratio
#' @param stage.max.pct the maximum percentage of features that will be excluded
#' @param is.fwer a logical value indicating whether the family-wise error rate control (West-Young) should be performed
#' @param is.tree.fdr a logical value indicating whether tree-based false discovery rate shuold be carried out
#' @param tree a class of "phylo", the tree relats all the OTUs, and should have the same names in "comm"
#' @param verbose a logical value indicating whether the trace information should be printed out
#' @param return.comm a logical value indicating whether the wisorized, filtered "comm" matrix should be returned
#' @param min.prop Undetermined
#' @param return.perm.F Undetermined
#' @param ... arguments passing to tree-based fdr control
#' @return A list with the elements
#' \item{call}{the call}
#' \item{comm}{the wisorized, filtered "comm" matrix}
#' \item{filter.ind}{a vector of logical values indicating which features are tested}
#' \item{R2}{a matrix of percent explained variance (number of features by number of functions)}
#' \item{F0}{a matrix of F-statistics (number of features by number of functions)}
#' \item{RSS}{a matrix of residual sum squares (number of features by number of functions)}
#' \item{df.model, df.residual}{degree of freedoms for the model and residual space}
#' \item{p.raw}{the raw p-values based on permutations (not accurate if "perm.no" is small)}
#' \item{p.adj.fdr}{permutation-based FDR-adjusted p-values}
#' \item{p.adj.tree.fdr}{permutation-based tree FDR-adjusted p-values}
#' \item{p.adj.fwer}{permutation-based FWER-adjusted (West-Young) p-values}
#' \item{tree.fdr.obj}{the object returned by the "TreeFDR"}
#' @rdname ZicoSeq
#' @import stats
#' @import permute
#' @import nlme
#' @import matrixStats
#' @import vegan
#' @import GUniFrac
#' @importFrom rmutil dbetabinom
#' @export
ZicoSeq <- function (
meta.dat, comm, grp.name, adj.name = NULL,
prev.filter = 0.1, abund.filter = 10, min.prop = 0,
is.winsor = TRUE, winsor.qt = 0.97,
is.prior = TRUE, prior.dist = c('BetaMix', 'ZIBB'),
post.method = c('sample', 'mean'), post.sample.no = 25,
link.func = list(function (x) x^0.25, function (x) x^0.5, function (x) x^0.75),
link.d.func = list(function (x) 0.25 * x^(-0.75), function (x) 0.5 * x^(-0.5), function (x) 0.75 * x^(-0.25)),
variance.EB = FALSE, df.prior = 10,
perm.no = 99, strata = NULL, stats.combine.func = max,
stage.no = 6, topK = NULL, stage.fdr = 0.75, stage.max.pct = 0.50,
is.fwer = FALSE, is.tree.fdr = FALSE, tree = NULL,
verbose = TRUE, return.comm = FALSE, return.perm.F = FALSE, ...) {
# Args:
# meta.dat: a data frame containing the sample information
# comm: a matrix of counts, row - features (OTUs, genes, etc) , column - sample
# grp.name: a character, variable of interest; it could be numeric or categorical; should be in "meta.dat"
# adj.name: a character vector, variable(s) to be adjusted; they could be numeric or categorical; should be in "meta.dat"
# prev.filter: features with prevalence (i.e., nonzero proportion) less than "prev.cutoff" or be filtered
# abund.filter: features with a total counts less than "abund.cutoff" or be filtered
# is.winsor: a logical value indicating whether winsorization should be performed to replace outliers. The default is TRUE.
# winsor.qt: the winsorization quantile, above which the counts will be replaced
# is.prior: a logical value indicating whether to perform posterior inference based on some prior distribution on the proportion data
# prior.dist: prior distribution, either two-component beta-binomial mixture ("BetaMix") or zeroinflated beta-binomial ("ZIBB")
# post.method: method for posterior inference, either based on posterior sampling ("sample") or approximate posterior mean ("mean")
# post.sample.no: the number of posterior samples if posterior sampling is used
# link.func: a list of functions that connects the ratios to the covariates
# link.d.func: a list of the derivative function of "link.func"; only need to specifiy when "post.method" is "mean"
# variance.EB: a logical value indicating whehter to perform empirical Bayes based variance shrinkage
# df.prior: the degree of freedom of the prior inverse gamma distribution for variance shrinkage
# perm.no: the number of permutations; If the raw p values are of the major interest, set "perm.no" to at least 999
# strata: a factor indicating the permutation strata; permutation will be confined to each stratum
# stats.combine.func: function to combine the F-statistic for the omnibus test
# stage.no: the number of stages if multiple-stage ratio stategy is used
# topK: the number of dominant features that will be excluded in the initial stage ratio calculation
# stage.fdr: the fdr cutoff below which the features will be excluded for calculating the ratio
# stage.max.pct: the maximum percentage of features that will be excluded
# is.fwer: a logical value indicating whether the family-wise error rate control (West-Young) should be performed
# is.tree.fdr: a logical value indicating whether tree-based false discovery rate shuold be carried out
# tree: a class of "phylo", the tree relats all the OTUs, and should have the same names in "comm"
# verbose: a logical value indicating whether the trace information should be printed out
# return.comm: a logical value indicating whether the wisorized, filtered "comm" matrix should be returned
# ...: arguments passing to tree-based fdr control
#
# Returns:
# a list that contains:
# call: the call
# comm: the wisorized, filtered "comm" matrix
# filter.ind: a vector of logical values indicating which features are tested
# R2: a matrix of percent explained variance (number of features by number of functions)
# F0: a matrix of F-statistics (number of features by number of functions)
# RSS: a matrix of residual sum squares (number of features by number of functions)
# df.model, df.residual: degree of freedoms for the model and residual space
# p.raw : the raw p-values based on permutations (not accurate if "perm.no" is small)
# p.adj.fdr: permutation-based FDR-adjusted p-values
# p.adj.tree.fdr: permutation-based tree FDR-adjusted p-values
# p.adj.fwer: permutation-based FWER-adjusted (West-Young) p-values
# tree.fdr.obj: the object returned by the "TreeFDR"
this.call <- match.call()
prior.dist <- match.arg(prior.dist)
post.method <- match.arg(post.method)
###############################################################################
# Winsorization to reduce the influence of outlier counts
if (is.winsor == TRUE) {
depth <- colSums(comm)
comm <- apply(comm, 1, function (x) {
x <- x / depth
qt <- quantile(x, winsor.qt)
x[x >= qt] <- qt
round(x * depth)
})
comm <- t(comm)
}
###############################################################################
# Filter to remove very rare taxa
filter.ind <- rowMeans(comm != 0) >= prev.filter & rowSums(comm) >= abund.filter
names(filter.ind) <- rownames(comm)
if (verbose) cat(sum(!filter.ind), ' features are filtered!\n')
comm <- comm[filter.ind, ]
sample.no <- ncol(comm)
otu.no <- nrow(comm)
row.names <- rownames(comm)
depth <- colSums(comm)
if (verbose) cat('The data has ', sample.no, ' samples and ', otu.no, ' features will be tested!\n' )
rabund <- colMeans(t(comm) / colSums(comm))
if (is.null(topK)) {
topK <- round(length(rabund) * 0.4) # Very aggressive
}
size.factor <- colSums(comm[order(rabund)[1 : (length(rabund) - topK)], ])
###############################################################################
# Generate samples from posterior distribution (stacking it)
if (is.prior == TRUE) {
if (verbose) cat('Fitting beta mixture ...\n')
if (post.method == 'sample') {
comm.p <- apply(comm, 1, function (x) {
if (prior.dist == 'BetaMix') {
err1 <- try(res <- bbmix.fit.MM(x, depth))
}
if (prior.dist == 'ZIBB') {
err1 <- try(res <- zibb.fit.MM(x, depth))
}
# Handle error
if (class(err1) != 'try-error') {
if (prior.dist == 'BetaMix') {
prop1.1 <- rbeta(sample.no * post.sample.no, shape1 = x + res$shape1.1, shape2 = res$shape1.2 + depth - x)
prop1.2 <- rbeta(sample.no * post.sample.no, shape1 = x + res$shape2.1, shape2 = res$shape2.2 + depth - x)
prop <- ifelse(runif(sample.no * post.sample.no) <= res$q1, prop1.1, prop1.2)
}
if (prior.dist == 'ZIBB') {
prop1.2 <- rbeta(sample.no * post.sample.no, shape1 = x + res$shape2.1, shape2 = res$shape2.2 + depth - x)
prop <- ifelse(runif(sample.no * post.sample.no) <= res$q1, 0, prop1.2)
}
} else {
prop <- x / depth
v <- var(prop)
m <- mean(prop)
a1 <- ((1 - m) / v - 1 / m) * m ^ 2
a2 <- a1 * (1 / m - 1)
if (is.na(a1) | a1 < 0) {
# uniform prior
prop <- rbeta(sample.no * post.sample.no, shape1 = x + 1, shape2 = otu.no + depth - x)
} else {
prop <- rbeta(sample.no * post.sample.no, shape1 = x + a1, shape2 = a2 + depth - x)
}
}
return(prop)
})
comm.p <- t(comm.p)
comm.p.list <- list()
st <- 1
end <- sample.no
for (i in 1 : post.sample.no) {
comm.p.list[[i]] <- comm.p[, st : end]
st <- st + sample.no
end <- end + sample.no
}
} else {
EV <- apply(comm, 1, function (x) {
if (prior.dist == 'BetaMix') {
err1 <- try(res <- bbmix.fit.MM(x, depth))
}
if (prior.dist == 'ZIBB') {
err1 <- try(res <- zibb.fit.MM(x, depth))
}
# Handle error
if (class(err1) != 'try-error') {
if (prior.dist == 'BetaMix') {
Ep <- res$q1 * (x + res$shape1.1) / (depth + res$shape1.1 + res$shape1.2) +
(1 - res$q1) * (x + res$shape2.1) / (depth + res$shape2.1 + res$shape2.2)
Vp <- res$q1 * ((x + res$shape1.1)^2 / (depth + res$shape1.1 + res$shape1.2)^2 +
(x + res$shape1.1) * (res$shape1.2 + depth - x) / (depth + res$shape1.1 + res$shape1.2)^2 /
(depth + res$shape1.1 + res$shape1.2 + 1)) +
(1 - res$q1) * ((x + res$shape2.1)^2 / (depth + res$shape2.1 + res$shape2.2)^2 +
(x + res$shape2.1) * ( res$shape2.2 + depth - x) / (depth + res$shape2.1 + res$shape2.2)^2 /
(depth + res$shape2.1 + res$shape2.2 + 1)) - Ep^2
}
if (prior.dist == 'ZIBB') {
Ep <- (1 - res$q1) * (x + res$shape2.1) / (depth + res$shape2.1 + res$shape2.2)
Vp <- (1 - res$q1) * ((x + res$shape2.1)^2 / (depth + res$shape2.1 + res$shape2.2)^2 +
(x + res$shape2.1) * ( res$shape2.2 + depth - x) / (depth + res$shape2.1 + res$shape2.2)^2 /
(depth + res$shape2.1 + res$shape2.2 + 1)) - Ep^2
}
} else {
prop <- x / depth
v <- var(prop)
m <- mean(prop)
a1 <- ((1 - m) / v - 1 / m) * m ^ 2
a2 <- a1 * (1 / m - 1)
if (is.na(a1) | a1 < 0) {
# uniform prior
Ep <- (x + 1) / (depth + otu.no)
Vp <- (x + 1) * (depth - x + otu.no - 1) / (otu.no + depth)^2 /
(otu.no + depth + 1)
} else {
Ep <- (x + a1) / (depth + a1 + a2)
Vp <- (x + a1) * (depth - x + a2) / (depth + a1 + a2)^2 /
(depth + a1 + a2 + 1)
}
}
return(c(Ep, Vp))
})
EV <- t(EV)
comm.p.list <- list()
comm.v.list <- list()
comm.p.list[[1]] <- EV[, 1 : sample.no]
comm.v.list[[1]] <- EV[, (sample.no + 1) : (2 * sample.no)]
post.sample.no <- 1
}
} else {
comm.p.list <- list()
comm.p.list[[1]] <- t(t(comm) / depth)
post.sample.no <- 1
}
# Replace zeros or extremely small values for log calculation
for (i in 1 : post.sample.no) {
temp <- comm.p.list[[i]]
temp[temp <= min.prop] <- min.prop
comm.p.list[[i]] <- temp
}
###############################################################################
# Covariate space (including intercept)
if (!is.null(strata)) {
strata <- factor(strata)
}
if (is.null(adj.name)) {
M0 <- model.matrix(~ 1, meta.dat)
} else {
data0 <- meta.dat[, c(adj.name), drop = FALSE]
M0 <- model.matrix( ~ ., data0)
}
data1 <- meta.dat[, c(grp.name), drop = FALSE]
M1 <- model.matrix( ~ ., data1)[, -1, drop = FALSE] # No intercept
M01 <- cbind(M0, M1)
# QR decompostion
qrX0 <- qr(M0, tol = 1e-07)
Q0 <- qr.Q(qrX0)
Q0 <- Q0[, 1:qrX0$rank, drop = FALSE]
H0 <- (Q0 %*% t(Q0))
qrX1 <- qr(M1, tol = 1e-07)
Q1 <- qr.Q(qrX1)
Q1 <- Q1[, 1:qrX1$rank, drop = FALSE]
qrX01 <- qr(M01, tol = 1e-07)
Q01 <- qr.Q(qrX01)
Q01 <- Q01[, 1:qrX01$rank, drop = FALSE]
R0 <- as.matrix(resid(lm(Q1 ~ Q0 - 1)))
pX0 <- ncol(Q0)
pX1 <- ncol(Q1)
pX01 <- ncol(Q01)
df.model <- pX01 - pX0
df.residual <- sample.no - pX01
func.no <- length(link.func)
###############################################################################
# Perform multiple stage normalization
norm.ind <- NULL
for (i in 1:stage.no) {
if (verbose == TRUE)
cat('Stage ', i, ' ')
# Reference proportion
divisor <- size.factor / depth
# Create the giant Y matrix
Y <- matrix(NA, sample.no, func.no * otu.no * post.sample.no)
if (post.method == 'mean') {
D <- matrix(NA, sample.no, func.no * otu.no * post.sample.no)
}
# Change order - func.no * otu.no
for (k in 1 : post.sample.no) {
for (j in 1 : func.no) {
func <- link.func[[j]]
comm.p <- comm.p.list[[k]]
if (post.method == 'mean') {
dfunc <- link.d.func[[j]]
comm.v <- comm.v.list[[k]]
}
Y[, (k - 1) * func.no * otu.no + func.no * (0 : (otu.no - 1)) + j] <-
func(t(comm.p) / divisor) # No scaling first
if (post.method == 'mean') {
# row - sample, column - features
D[, (k - 1) * func.no * otu.no + func.no * (0 : (otu.no - 1)) + j] <-
(dfunc(t(comm.p) / divisor) / divisor)^2 * t(comm.v)
}
}
}
Y <- t(Y)
if (post.method == 'sample') {
TSS <- rowSums(Y^2)
MSS01 <- rowSums((Y %*% Q01)^2)
MSS0 <- rowSums((Y %*% Q0)^2)
} else {
w0 <- colSums((Q0 %*% t(Q0))^2)
w01 <- colSums((Q01 %*% t(Q01))^2)
TSS <- rowSums(Y^2) + colSums(D)
MSS01 <- rowSums((Y %*% Q01)^2) + colSums(w01 * D)
MSS0 <- rowSums((Y %*% Q0)^2) + colSums(w0 * D)
}
MSS <- (MSS01 - MSS0)
RSS <- (TSS - MSS01)
perm.ind <- getPermuteMatrix(perm.no, sample.no, strata = strata)
perm.no <- nrow(perm.ind)
MRSSp <- sapply(1 : perm.no, function(ii) {
if (verbose) {
if (ii %% 10 == 0) cat('.')
}
Rp <- R0[perm.ind[ii, ], , drop = FALSE]
# Project to the reisdual space
Rp <- Rp - H0 %*% Rp
qrRp <- qr(Rp, tol = 1e-07)
Q1p <- qr.Q(qrRp)
Q1p <- Q1p[, 1:qrRp$rank, drop = FALSE]
if (post.method == 'sample') {
MSS01p <- MSS0 + rowSums((Y %*% Q1p)^2)
} else {
w1p <- colSums((Q1p %*% t(Q1p))^2)
MSS01p <- MSS0 + rowSums((Y %*% Q1p)^2) + colSums(w1p * D)
}
MSSp <- (MSS01p - MSS0)
RSSp <- (TSS - MSS01p)
c(MSSp, RSSp)
})
unit <- func.no * otu.no * post.sample.no
MSSp <- MRSSp[1 : unit, ]
RSSp <- MRSSp[(unit + 1) : (2 * unit), ]
# EB is based on the aggregated RSS
RSS.m <- array(RSS, c(func.no, otu.no, post.sample.no))
RSS.m <- t(apply(RSS.m, c(1, 2), mean)) # otu.no by func.no
F0.m <- array((MSS / df.model) / (RSS / df.residual), c(func.no, otu.no, post.sample.no))
F0.m <- t(apply(F0.m, c(1, 2), mean)) # otu.no by func.no
R2.m <- array(MSS / TSS, c(func.no, otu.no, post.sample.no))
R2.m <- t(apply(R2.m, c(1, 2), mean)) # otu.no by func.no
###############################################################################
# Variance moderation - to be improved
if (variance.EB) {
# hist(RSS.m / df.residual)
# out <- apply(RSS.m, 2, function (x) {
# out <- limma::squeezeVar(x / df.residual, df.residual)
# c(out$var.prior, out$df.prior)
# })
# var.prior <- out[1, ]
# df.prior <- out[2, ]
#
# df.prior[df.prior > 1E6] <- 1E6
# df.prior[df.prior < 1E-6] <- 1E-6
var.prior <- colMeans(RSS.m / df.residual, na.rm = TRUE)
F0 <- (MSS / df.model) / ((RSS + df.prior * var.prior) / (df.prior + df.residual))
Fp <- (MSSp / df.model) / ((RSSp + df.prior * var.prior) / (df.prior + df.residual))
# Expectation of F0 and Fp
F0 <- array(F0, c(func.no, otu.no, post.sample.no))
Fp <- array(Fp, c(func.no, otu.no, post.sample.no, perm.no))
F0 <- apply(F0, c(1, 2), mean) # func.no * otu.no
Fp <- apply(Fp, c(1, 2, 4), mean) # func.no * otu.no * perm.no
} else {
F0 <- (MSS / df.model) / (RSS / df.residual)
Fp <- (MSSp / df.model) / (RSSp / df.residual)
# Expectation of F0 and Fp
F0 <- array(F0, c(func.no, otu.no, post.sample.no))
Fp <- array(Fp, c(func.no, otu.no, post.sample.no, perm.no))
F0 <- apply(F0, c(1, 2), mean) # func.no * otu.no
Fp <- apply(Fp, c(1, 2, 4), mean) # func.no * otu.no * perm.no
}
###############################################################################
# Omnibus test by taking maximum
F0 <- apply(F0, 2, stats.combine.func)
Fp <- apply(Fp, c(2, 3), stats.combine.func) # otu.no by perm.no
if (verbose) cat('\n')
if (mean(is.na(F0)) >= 0.1) {
warning('More than 10% observed F stats have NA! Please check! \n')
}
if (mean(is.na(Fp)) >= 0.1) {
warning('More than 10% permuted F stats have NA! Please check! \n')
}
na.ind <- is.na(F0)
F0 <- F0[!na.ind]
Fp <- Fp[!na.ind, ]
rabund0 <- rabund[!na.ind]
which.nan.ind <- which(!na.ind)
p.adj.fdr <- perm_fdr_adj(F0, Fp)
if (i == stage.no) {
break
} else {
# recalculating the size factor
if (mean(p.adj.fdr <= stage.fdr) > stage.max.pct) {
ind <- order(p.adj.fdr, -rabund0)[1 : round(length(p.adj.fdr) * stage.max.pct)]
} else {
ind <- which(p.adj.fdr < stage.fdr)
}
size.factor <- colSums(comm[setdiff(1:nrow(comm), which.nan.ind[ind]), ])
norm.ind <- cbind(norm.ind, !(1:nrow(comm) %in% which.nan.ind[ind]))
}
}
###############################################################################
# Permutation-based false discovery rate control
# Tree-based power improvement
if (is.tree.fdr) {
if (verbose) {
cat("Perform Tree-based FDR control ...\n")
}
# Convert pseudo-F to pseudo-P
if (is.null(tree)) {
stop('Tree FDR requires a tree of class "phylo"!\n')
}
test.obs <- list()
test.perm <- list()
test.obs$p.value <- 1 - pf(F0, df.model, df.residual)
test.perm$p.value <- 1 - pf(Fp, df.model, df.residual)
# Ignore the sign
test.obs$e.sign <- sign(test.obs$p.value)
test.perm$e.sign <- sign(test.perm$p.value)
D <- cophenetic(tree)
D <- D[rownames(comm), rownames(comm)]
tree.fdr.obj <- TreeFDR(X = NULL, Y = NULL, test.obs = test.obs, test.perm = test.perm, D = D,
eff.sign = FALSE, alt.FDR = 'Permutation', verbose = verbose, ...)
p.adj.tree.fdr <- tree.fdr.obj$p.adj
if (verbose) {
cat('FWER control has not used the phylogenetic information!\n')
}
} else {
# scaled F stat
tree.fdr.obj <- NULL
p.adj.tree.fdr <- NULL
}
if (is.fwer) {
p.adj.fwer <- perm_fwer_adj2(F0, Fp)
}
p.raw <- rowMeans(cbind(Fp, F0) >= F0)
p.raw <- na.pad(p.raw, na.ind)
p.adj.fdr <- na.pad(p.adj.fdr, na.ind)
names(p.raw) <- names(p.adj.fdr) <- rownames(R2.m) <- rownames(RSS.m) <- rownames(F0.m) <- row.names
colnames(R2.m) <- colnames(F0.m) <- colnames(RSS.m) <- paste0('Func', 1 : func.no)
if (is.tree.fdr) {
p.adj.tree.fdr <- na.pad(p.adj.tree.fdr, na.ind)
names(p.adj.tree.fdr) <- row.names
}
if (is.fwer) {
p.adj.fwer <- na.pad(p.adj.fwer, na.ind)
names(p.adj.fwer) <- row.names
} else {
p.adj.fwer <- NULL
}
if (!return.comm) {
comm <- NULL
}
if (!return.perm.F) {
Fp <- NULL
}
if (verbose) cat('Completed!\n')
return(list(call = this.call, comm = comm, filter.ind = filter.ind,
R2 = R2.m, F0 = F0.m, Fp = Fp, RSS = RSS.m, df.model = df.model, df.residual = df.residual,
p.raw = p.raw, p.adj.fdr = p.adj.fdr, p.adj.tree.fdr = p.adj.tree.fdr, p.adj.fwer = p.adj.fwer,
tree.fdr.obj = tree.fdr.obj, norm.ind = norm.ind))
}
chloelulu/ZicoSeq documentation built on Nov. 4, 2019, 8:50 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Ztransform corHerit
# 2019_02_20
# Beta-mixture to address sequence depth confounding
# Permutation-based FDR control
# Multiple-stage normalization
# Variance modulation by LIMMA - hasn't been finished
# Expectation of F statistics based on posterior sampling
# Omnibus test for combining diffrent effects
# Version 0.9 [PermuteDAA3.2]
#' Computes the GMPR size factor
#' @param comm a matrix of counts, row - features (OTUs, genes, etc) , column - sample
#' @param intersect.no the minimum number of shared features between sample pair, where the ratio is calculated
#' @param ct.min the minimum number of counts required to calculate ratios ct.min = 5 has better results
#' @param verbose suppress output
#' @return A list that contains:
#' \item{gmpr}{the GMPR size factors for all samples; Samples with distinct sets of features will be output as NA.}
#' \item{nss}{number of samples with significant sharing (> intersect.no) including itself}
#' @rdname GMPR
#' @export
GMPR <- function (comm, intersect.no = 4, ct.min = 2, verbose = FALSE) {
# mask counts < ct.min
comm[comm < ct.min] <- 0
if (is.null(colnames(comm))) {
colnames(comm) <- paste0('S', 1:ncol(comm))
}
if (verbose == TRUE)
cat('Begin GMPR size factor calculation ...\n')
comm.no <- numeric(ncol(comm))
gmpr <- sapply(1:ncol(comm), function(i) {
if (i %% 50 == 0) {
if (verbose == TRUE)
cat(i, '\n')
}
x <- comm[, i]
# Compute the pairwise ratio
pr <- x / comm
# Handling of the NA, NaN, Inf
pr[is.nan(pr) | !is.finite(pr) | pr == 0] <- NA
# Counting the number of non-NA, NaN, Inf
incl.no <- colSums(!is.na(pr))
# Calculate the median of PR
pr.median <- colMedians(pr, na.rm=TRUE)
# Record the number of samples used for calculating the GMPR
comm.no[i] <<- sum(incl.no >= intersect.no)
# Geometric mean of PR median
if (comm.no[i] > 1) {
return(exp(mean(log(pr.median[incl.no >= intersect.no]))))
} else {
return(NA)
}
}
)
if (sum(is.na(gmpr))) {
warning(paste0('The following samples\n ', paste(colnames(comm)[is.na(gmpr)], collapse='\n'),
'\ndo not share at least ', intersect.no, ' common taxa with the rest samples! ',
'For these samples, their size factors are set to be NA! \n',
'You may consider removing these samples since they are potentially outliers or negative controls!\n',
'You may also consider decreasing the minimum number of intersecting taxa and rerun the procedure!\n'))
}
if (verbose == TRUE) {
cat('Completed!\n')
cat('Please watch for the samples with limited sharing with other samples based on NSS! They may be outliers! \n')
}
attr(gmpr, 'NSS') <- comm.no
names(gmpr) <- colnames(comm)
return(gmpr * median(colSums(comm)))}
####################################################################################
# Copied from "vegan" package getPermuteMatrix.R to reduce dependency
getPermuteMatrix <- function (perm, N, strata = NULL) {
if (length(perm) == 1) {
perm <- how(nperm = perm)
}
if (!missing(strata) && !is.null(strata)) {
if (inherits(perm, "how") && is.null(getBlocks(perm)))
setBlocks(perm) <- strata
}
if (inherits(perm, "how"))
perm <- shuffleSet(N, control = perm)
if (is.null(attr(perm, "control")))
attr(perm, "control") <- structure(list(within = list(type = "supplied matrix"),
nperm = nrow(perm)), class = "how")
perm
}
# Permutation-based FDR control
perm_fdr_adj <- function (F0, Fp) {
ord <- order(F0, decreasing = T)
F0 <- F0[ord]
perm.no <- ncol(Fp)
Fp <- as.vector(Fp)
Fp <- Fp[!is.na(Fp)]
Fp <- sort(c(Fp, F0), decreasing = F)
n <- length(Fp)
m <- length(F0)
FPN <- (n + 1) - match(F0, Fp) - 1:m
p.adj.fdr <- FPN / perm.no / (1:m)
# p.adj.fdr <- sapply(F0, function(x) sum(Fp >=
# x, na.rm=TRUE) / perm.no)/(1:length(F0))
p.adj.fdr <- pmin(1, rev(cummin(rev(p.adj.fdr))))[order(ord)]
}
# Permutation-based FWER control
# Conservative
perm_fwer_adj <- function (F0, Fp) {
ord <- order(F0, decreasing = T)
F0 <- F0[ord]
col.max <- colMaxs(Fp, na.rm=TRUE)
p.adj.fwer <- sapply(F0, function(x) mean(col.max >= x))[order(ord)]
}
# Westfall young too slow
perm_fwer_adj2 <- function (F0, Fp) {
ord <- order(F0, decreasing = T)
m <- length(F0)
F0 <- F0[ord]
Fp <- Fp[ord, , drop=FALSE]
col.max <- Fp[m, ]
p.adj.fwer <- sapply(m:1, function(i) {
x <- F0[i]
y <- Fp[i, ]
col.max <<- ifelse(y > col.max, y, col.max)
mean(col.max >= x)
})
p.adj.fwer <- rev(p.adj.fwer)
p.adj.fwer <- pmin(1, rev(cummin(rev(p.adj.fwer))))[order(ord)]
}
# Pad back the NA values
na.pad <- function (vec, ind) {
vec0 <- numeric(length(ind))
vec0[!ind] <- vec
vec0[ind] <- NA
vec0
}
# # Pad back the NA values
one.pad <- function (vec, ind) {
vec0 <- numeric(length(ind))
vec0[!ind] <- vec
vec0[ind] <- 1
vec0
}
# Transform to Z score
Ztransform <- function(p.value, e.sign, eff.sign=TRUE, tol=1E-15) {
p.value[p.value <= tol] <- tol
p.value[p.value >= 1 - tol] <- 1 - tol
if (eff.sign == TRUE) {
e.sign[e.sign == 0] <- sample(c(-1, 1), sum(e.sign == 0), replace=T)
z1 <- qnorm(p.value / 2)
z2 <- qnorm(1 - p.value / 2)
z <- ifelse(e.sign > 0, z2, z1)
} else {
z <- qnorm(1 - p.value)
}
return(z)
}
# Internal functions
# For gls with generic correlation structure
corHerit <- function(value, paras, form = ~1, fixed = TRUE) {
# Place holder - check the validity of the parameter
object <- value
attr(object, "formula") <- form
attr(object, "fixed") <- fixed
attr(object, "paras") <- paras
class(object) <- c("corHerit", "corStruct")
object
}
Initialize.corHerit <- function (object, data, ...) {
# Place holder - check the validity of the parameter
form <- formula(object)
if (!is.null(getGroupsFormula(form))) {
attr(object, "groups") <- getGroups(object, form, data = data)
attr(object, "Dim") <- Dim(object, attr(object, "groups"))
} else {
attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
}
attr(object, "covariate") <- getCovariate(object, data = data)
object
}
corMatrix.corHerit <- function (object, covariate = getCovariate(object), ...) {
paras <- attr(object, "paras")
p <- paras[['p']]
I <- diag(p)
hyper <- as.vector(object)
lambda <- 1 / (1 + exp(hyper[1]))
V <- exp(-((paras[['D']])) * (exp(hyper[2])))
cor.m <- (1 - lambda) * V + lambda * I
cor.m
}
# Extract the coefficient
coef.corHerit <- function (object, ...) {
paras <- attr(object, "paras")
coefs <- as.vector(object)
coef1 <- 1 / (1 + exp(coefs[1]))
coef1 <- coef1 / (1 - coef1)
coef2 <- exp(coefs[2])
coefs <- c(coef1, coef2)
names(coefs) <- paste("Hyper", 1:length(coefs), sep="")
coefs
}
# Estimating the hyper parameters
EstHyper <- function (y, D, init.val=c(0, 0)) {
obj.gls <- gls(model = y ~ 1,
correlation = corHerit(value=init.val, paras = list(p=length(y), D = D)))
cc <- c(coef(obj.gls$modelStruct$corStruct), obj.gls$coefficients, obj.gls$logLik)
cc
}
AdjStats <- function (y, V, k, mu, fudge=0.005) {
p <- nrow(V)
# Add small fudge
V.inv <- solve(V + fudge * diag(p))
I <- diag(p)
y.adj <- solve(I + k * V.inv) %*% (k * mu * rowSums(V.inv) + y)
y.adj
}
# Compute the NFP on the null.ind
PermFDR <- function (F0, Fp, null.ind) {
Fp <- as.matrix(Fp)
pct <- sum(null.ind) / length(F0)
Fp <- Fp[null.ind, ]
ord <- order(F0, decreasing = T)
F0 <- F0[ord]
perm.no <- ncol(Fp)
Fp <- as.vector(Fp)
pct.non.na <- mean(!is.na(Fp))
Fp <- Fp[!is.na(Fp)]
Fp <- sort(c(Fp, F0), decreasing = F)
n <- length(Fp)
m <- length(F0)
FPN <- (n + 1) - match(F0, Fp) - 1:m
# Handle identical F0 situations
FPN <- cummax(FPN)
p.adj.fdr <- FPN / pct.non.na / pct / perm.no / (1:m)
p.adj.fdr <- pmin(1, rev(cummin(rev(p.adj.fdr))))[order(ord)]
}
####################################################################################
# Compute tree-based FDR control
TreeFDR <- function (X = NULL, Y = NULL, test.func = NULL, perm.func = NULL, test.obs = NULL, test.perm = NULL, tree = NULL, D,
eff.sign = TRUE, B = 20, q.cutoff = 0.5, alpha = 1,
adaptive = c('Fisher', 'Overlap'), alt.FDR = c('BH', 'Permutation'), verbose = TRUE, ...) {
if (is.null(X) | is.null(Y) | is.null(test.func) | is.null(perm.func)) {
if (is.null(test.obs) | is.null(test.perm)) {
stop('Please specify either "X, Y, test.func, perm.fuc", or "test.obs", "test.perm"!\n')
} else {
X <- matrix(test.obs$p.value, ncol = 1)
rownames(X) <- names(test.obs$p.value)
}
}
adaptive <- match.arg(adaptive)
alt.FDR <- match.arg(alt.FDR)
# Make sure the rows of X and tree tips are labeled to avoid error
if ((is.null(rownames(X)) | is.null(tree$tip.label))) {
warning('Both the data matrix and the tree should have labels (rownames, tip.label) to avoid potential errors!\n')
} else {
if (sum(!(rownames(X) %in% tree$tip.label))){
stop('Some features in the data matrix are not in the tree! Please check!\n')
} else {
if (sum(!(tree$tip.label %in% rownames(X)))) {
warning('The tree have more features than the data matrix! \n')
}
}
}
# Patristic distance
if (!is.null(tree)) {
D <- (cophenetic(tree)) ^ alpha
} else {
D <- D ^ alpha
}
if (!is.null(rownames(X)) & !is.null(tree$tip.label)) {
D <- D[rownames(X), rownames(X)]
}
if (is.null(test.obs)) {
if (verbose)
cat('Test on original data sets ...\n')
test.obs <- test.func(X, Y, ...)
}
if (!is.list(test.obs) | !all(c("e.sign", "p.value") %in%
names(test.obs))) {
stop("test.func should return a list with names e.sign and p.valueif z.transform=TRUE! Please check!\n")
}
null.ind <- test.obs$p.value >= quantile(test.obs$p.value, 1 - q.cutoff)
z.obs <- Ztransform(test.obs$p.value, test.obs$e.sign, eff.sign)
if (verbose)
cat('Test on permuted data sets ...\n')
if (is.null(test.perm)) {
z.perm <- z.perm2 <- matrix(NA, nrow(X), B)
for (i in 1:B) {
perm.obj <- perm.func(X, Y, ...)
if (!is.list(perm.obj) | !all(c("X", "Y") %in% names(perm.obj))) {
stop("perm.func should return a list with names X and Y! Please check!\n")
}
X.perm <- perm.obj$X
Y.perm <- perm.obj$Y
test.perm <- test.func(X.perm, Y.perm, ...)
z.perm[, i] <- z.perm2[, i] <- Ztransform(test.perm$p.value, test.perm$e.sign, eff.sign)
z.perm2[!null.ind, i] <- z.obs[!null.ind]
}
} else {
B <- ncol(test.perm$p.value)
z.perm <- z.perm2 <- matrix(NA, length(test.obs$p.value), B)
for (i in 1:B) {
z.perm[, i] <- z.perm2[, i] <- Ztransform(test.perm$p.value[, i], test.perm$e.sign[, i], eff.sign)
z.perm2[!null.ind, i] <- z.obs[!null.ind]
}
}
if (alt.FDR == 'Permutation') {
if (verbose)
cat('Perform ordinary permutation-based FDR control ...\n')
if (eff.sign == TRUE) {
p.adj0 <- PermFDR(abs(z.obs), abs(z.perm), rep(TRUE, length(z.obs)))
} else {
p.adj0 <- PermFDR(z.obs, z.perm, rep(TRUE, length(z.obs)))
}
}
if (alt.FDR == 'BH') {
if (verbose)
cat('Perform ordinary BH-based FDR control ...\n')
p.adj0 <- p.adjust(test.obs$p.value, 'fdr')
}
if (verbose)
cat("Estimating hyperparameter ... \n")
error <- try(obj <- EstHyper(y = z.obs, D = D))
if (inherits(error, "try-error")) {
if (verbose)
cat('Hyperparameter estimation failed! Ordinary permutation-based FDR control will be used!\n')
# p.adj <- p.adjust(test.obs$p.value,'fdr')
p.adj <- p.adj0
k <- NULL
rho <- NULL
z.adj <- NULL
} else {
if (verbose)
cat("Structure-based adjustment ...\n")
k <- obj[1]
rho <- obj[2]
mu <- obj[3]
V <- exp(-1 * rho * D)
z.adj <- stat.o <- AdjStats(y = z.obs, V = V, k = k, mu = mu)[, 1]
stat.p <- AdjStats(y = z.perm2, V = V, k = k, mu = mu)
if (eff.sign == TRUE) {
stat.o <- abs(stat.o)
stat.p <- abs(stat.p)
}
p.adj <- PermFDR(stat.o, stat.p, null.ind)
# Power loss check
if (adaptive == 'Overlap') {
# Magic numbers
fdr.cutoff <- 0.2
pct.cutoff <- 0.5
ind0 <- p.adj0 <= fdr.cutoff
ind <- p.adj <= fdr.cutoff
if (sum(p.adj[ind0] <= fdr.cutoff) < pct.cutoff * sum(ind0)) {
# Over-adjustment checking
if (verbose)
cat('Potential over-adjustment! Alternative FDR control will be used!\n')
p.adj <- p.adj0
k <- NULL
rho <- NULL
z.adj <- NULL
}
}
if (adaptive == 'Fisher') {
# These cutoffs are used emprically
fdr.cutoff <- 0.2
ind0 <- p.adj0 <= fdr.cutoff
ind <- p.adj <= fdr.cutoff
n <- nrow(X)
test.p <- fisher.test(matrix(c(sum(ind), sum(ind0), n - sum(ind), n - sum(ind0)), 2, 2), alternative = 'less')$p.value
if (test.p <= 0.05) {
# Over-adjustment checking
if (verbose)
cat('Potential over-adjustment! Alternative FDR control will be used!\n')
p.adj <- p.adj0
k <- NULL
rho <- NULL
z.adj <- NULL
}
}
}
if (verbose)
cat("Done!\n")
return(list(p.adj = p.adj, p.unadj = test.obs$p.value, z.adj = z.adj, z.unadj = z.obs, k = k, rho = rho))
}
###################################################################################
# add BBmix
###############################################################################################
bbmix.fit.MM <- function (ct, dep, nIter = 10, winsor.qt = 1.0) {
if (mean(ct == 0) >= 0.95 | sum(ct != 0) < 5) {
stop('The number of nonzero is too small to fit the model! Consider removing it from testing!\n')
}
# Initialization
prop0 <- ct / dep
qt <- quantile(prop0, winsor.qt)
prop0[prop0 >= qt] <- qt
var1 <- var(prop0)
mean1 <- mean(prop0) / 2
var2 <- var1
mean2 <- 3 * mean(prop0) / 2
pi <- 0.5
for (i in 1 : nIter) {
shape1.1 <- ((1 - mean1) / var1 - 1 / mean1) * mean1 ^ 2
shape1.2 <- shape1.1 * (1 / mean1 - 1)
shape2.1 <- ((1 - mean2) / var2 - 1 / mean2) * mean2 ^ 2
shape2.2 <- shape2.1 * (1 / mean2 - 1)
m1 <- shape1.1 / (shape1.1 + shape1.2)
s1 <- shape1.1 + shape1.2
m2 <- shape2.1 / (shape2.1 + shape2.2)
s2 <- shape2.1 + shape2.2
f1 <- dbetabinom(ct, dep, m1, s1)
f2 <- dbetabinom(ct, dep, m2, s2)
q1 <- pi * f1 / (pi * f1 + (1 - pi) * f2)
q2 <- 1 - q1
pi <- mean(q1)
# Rough estimation
mean1 <- sum(prop0 * q1) / sum(q1)
var1 <- sum((prop0 - mean1)^2 * q1) / sum(q1)
mean2 <- sum(prop0 * q2) / sum(q2)
var2 <- sum((prop0 - mean2)^2 * q2) / sum(q2)
}
return(list(shape1.1 = shape1.1, shape1.2 = shape1.2, shape2.1 = shape2.1, shape2.2 = shape2.2, pi = pi, q1 = q1))
}
###############################################################################################
zibb.fit.MM <- function (ct, dep, nIter = 10, winsor.qt = 1.0) {
# Initialization
if (mean(ct == 0) >= 0.95 | sum(ct != 0) < 5) {
stop('The number of nonzero is too small to fit the model! Consider removing it from testing!\n')
}
ct1 <- ct[ct != 0]
dep1 <- dep[ct != 0]
# Initialization
prop0 <- ct1 / dep1
qt <- quantile(prop0, winsor.qt)
prop0[prop0 >= qt] <- qt
var0 <- var(prop0)
mean0 <- mean(prop0)
shape1 <- ((1 - mean0) / var0 - 1 / mean0) * mean0 ^ 2
shape2 <- shape1 * (1 / mean0 - 1)
m0 <- shape1 / (shape1 + shape2)
s0 <- shape1 + shape2
if (s0 < 0) {
s0 <- 1
}
m2 <- 0.75 * m0
s2 <- s0
pi <- 0.75 * mean(ct == 0)
##########################
prop0 <- ct / dep
qt <- quantile(prop0, winsor.qt)
prop0[prop0 >= qt] <- qt
##########################
for (i in 1 : nIter) {
# Expectation step
f1 <- ifelse(ct == 0, 1, 0)
f2 <- dbetabinom(ct, dep, m2, s2)
q1 <- pi * f1 / (pi * f1 + (1 - pi) * f2)
q2 <- 1 - q1
# Maximization step
pi <- mean(q1)
mean2 <- sum(prop0 * q2) / sum(q2)
var2 <- sum((prop0 - mean2)^2 * q2) / sum(q2)
shape2.1 <- ((1 - mean2) / var2 - 1 / mean2) * mean2 ^ 2
shape2.2 <- shape2.1 * (1 / mean2 - 1)
m2 <- shape2.1 / (shape2.1 + shape2.2)
s2 <- shape2.1 + shape2.2
}
shape2.1 <- m2 * s2
shape2.2 <- s2 - shape2.1
return(list(shape2.1 = shape2.1, shape2.2 = shape2.2, pi = pi, q1 = q1))
}
####################################################################################
#' Permutation-based differential abundance analysis
#' @param meta.dat a data frame containing the sample information
#' @param comm a matrix of counts, row - features (OTUs, genes, etc) , column - sample
#' @param grp.name a character, variable of interest; it could be numeric or categorical; should be in "meta.dat"
#' @param adj.name a character vector, variable(s) to be adjusted; they could be numeric or categorical; should be in "meta.dat"
#' @param prev.filter features with prevalence (i.e., nonzero proportion) less than "prev.cutoff" or be filtered
#' @param abund.filter features with a total counts less than "abund.cutoff" or be filtered
#' @param is.winsor a logical value indicating whether winsorization should be performed to replace outliers. The default is TRUE.
#' @param winsor.qt the winsorization quantile, above which the counts will be replaced
#' @param is.prior a logical value indicating whether to perform posterior inference based on some prior distribution on the proportion data
#' @param prior.dist prior distribution, either two-component beta-binomial mixture ("BetaMix") or zeroinflated beta-binomial ("ZIBB")
#' @param post.method method for posterior inference, either based on posterior sampling ("sample") or approximate posterior mean ("mean")
#' @param post.sample.no the number of posterior samples if posterior sampling is used
#' @param link.func a list of functions that connects the ratios to the covariates
#' @param link.d.func a list of the derivative function of "link.func"; only need to specifiy when "post.method" is "mean"
#' @param variance.EB a logical value indicating whehter to perform empirical Bayes based variance shrinkage
#' @param df.prior the degree of freedom of the prior inverse gamma distribution for variance shrinkage
#' @param perm.no the number of permutations; If the raw p values are of the major interest, set "perm.no" to at least 999
#' @param strata a factor indicating the permutation strata; permutation will be confined to each stratum
#' @param stats.combine.func function to combine the F-statistic for the omnibus test
#' @param stage.no the number of stages if multiple-stage ratio stategy is used
#' @param topK the number of dominant features that will be excluded in the initial stage ratio calculation
#' @param stage.fdr the fdr cutoff below which the features will be excluded for calculating the ratio
#' @param stage.max.pct the maximum percentage of features that will be excluded
#' @param is.fwer a logical value indicating whether the family-wise error rate control (West-Young) should be performed
#' @param is.tree.fdr a logical value indicating whether tree-based false discovery rate shuold be carried out
#' @param tree a class of "phylo", the tree relats all the OTUs, and should have the same names in "comm"
#' @param verbose a logical value indicating whether the trace information should be printed out
#' @param return.comm a logical value indicating whether the wisorized, filtered "comm" matrix should be returned
#' @param min.prop Undetermined
#' @param return.perm.F Undetermined
#' @param ... arguments passing to tree-based fdr control
#' @return A list with the elements
#' \item{call}{the call}
#' \item{comm}{the wisorized, filtered "comm" matrix}
#' \item{filter.ind}{a vector of logical values indicating which features are tested}
#' \item{R2}{a matrix of percent explained variance (number of features by number of functions)}
#' \item{F0}{a matrix of F-statistics (number of features by number of functions)}
#' \item{RSS}{a matrix of residual sum squares (number of features by number of functions)}
#' \item{df.model, df.residual}{degree of freedoms for the model and residual space}
#' \item{p.raw}{the raw p-values based on permutations (not accurate if "perm.no" is small)}
#' \item{p.adj.fdr}{permutation-based FDR-adjusted p-values}
#' \item{p.adj.tree.fdr}{permutation-based tree FDR-adjusted p-values}
#' \item{p.adj.fwer}{permutation-based FWER-adjusted (West-Young) p-values}
#' \item{tree.fdr.obj}{the object returned by the "TreeFDR"}
#' @rdname ZicoSeq
#' @import stats
#' @import permute
#' @import nlme
#' @import matrixStats
#' @import vegan
#' @import GUniFrac
#' @importFrom rmutil dbetabinom
#' @export
ZicoSeq <- function (
meta.dat, comm, grp.name, adj.name = NULL,
prev.filter = 0.1, abund.filter = 10, min.prop = 0,
is.winsor = TRUE, winsor.qt = 0.97,
is.prior = TRUE, prior.dist = c('BetaMix', 'ZIBB'),
post.method = c('sample', 'mean'), post.sample.no = 25,
link.func = list(function (x) x^0.25, function (x) x^0.5, function (x) x^0.75),
link.d.func = list(function (x) 0.25 * x^(-0.75), function (x) 0.5 * x^(-0.5), function (x) 0.75 * x^(-0.25)),
variance.EB = FALSE, df.prior = 10,
perm.no = 99, strata = NULL, stats.combine.func = max,
stage.no = 6, topK = NULL, stage.fdr = 0.75, stage.max.pct = 0.50,
is.fwer = FALSE, is.tree.fdr = FALSE, tree = NULL,
verbose = TRUE, return.comm = FALSE, return.perm.F = FALSE, ...) {
# Args:
# meta.dat: a data frame containing the sample information
# comm: a matrix of counts, row - features (OTUs, genes, etc) , column - sample
# grp.name: a character, variable of interest; it could be numeric or categorical; should be in "meta.dat"
# adj.name: a character vector, variable(s) to be adjusted; they could be numeric or categorical; should be in "meta.dat"
# prev.filter: features with prevalence (i.e., nonzero proportion) less than "prev.cutoff" or be filtered
# abund.filter: features with a total counts less than "abund.cutoff" or be filtered
# is.winsor: a logical value indicating whether winsorization should be performed to replace outliers. The default is TRUE.
# winsor.qt: the winsorization quantile, above which the counts will be replaced
# is.prior: a logical value indicating whether to perform posterior inference based on some prior distribution on the proportion data
# prior.dist: prior distribution, either two-component beta-binomial mixture ("BetaMix") or zeroinflated beta-binomial ("ZIBB")
# post.method: method for posterior inference, either based on posterior sampling ("sample") or approximate posterior mean ("mean")
# post.sample.no: the number of posterior samples if posterior sampling is used
# link.func: a list of functions that connects the ratios to the covariates
# link.d.func: a list of the derivative function of "link.func"; only need to specifiy when "post.method" is "mean"
# variance.EB: a logical value indicating whehter to perform empirical Bayes based variance shrinkage
# df.prior: the degree of freedom of the prior inverse gamma distribution for variance shrinkage
# perm.no: the number of permutations; If the raw p values are of the major interest, set "perm.no" to at least 999
# strata: a factor indicating the permutation strata; permutation will be confined to each stratum
# stats.combine.func: function to combine the F-statistic for the omnibus test
# stage.no: the number of stages if multiple-stage ratio stategy is used
# topK: the number of dominant features that will be excluded in the initial stage ratio calculation
# stage.fdr: the fdr cutoff below which the features will be excluded for calculating the ratio
# stage.max.pct: the maximum percentage of features that will be excluded
# is.fwer: a logical value indicating whether the family-wise error rate control (West-Young) should be performed
# is.tree.fdr: a logical value indicating whether tree-based false discovery rate shuold be carried out
# tree: a class of "phylo", the tree relats all the OTUs, and should have the same names in "comm"
# verbose: a logical value indicating whether the trace information should be printed out
# return.comm: a logical value indicating whether the wisorized, filtered "comm" matrix should be returned
# ...: arguments passing to tree-based fdr control
#
# Returns:
# a list that contains:
# call: the call
# comm: the wisorized, filtered "comm" matrix
# filter.ind: a vector of logical values indicating which features are tested
# R2: a matrix of percent explained variance (number of features by number of functions)
# F0: a matrix of F-statistics (number of features by number of functions)
# RSS: a matrix of residual sum squares (number of features by number of functions)
# df.model, df.residual: degree of freedoms for the model and residual space
# p.raw : the raw p-values based on permutations (not accurate if "perm.no" is small)
# p.adj.fdr: permutation-based FDR-adjusted p-values
# p.adj.tree.fdr: permutation-based tree FDR-adjusted p-values
# p.adj.fwer: permutation-based FWER-adjusted (West-Young) p-values
# tree.fdr.obj: the object returned by the "TreeFDR"
this.call <- match.call()
prior.dist <- match.arg(prior.dist)
post.method <- match.arg(post.method)
###############################################################################
# Winsorization to reduce the influence of outlier counts
if (is.winsor == TRUE) {
depth <- colSums(comm)
comm <- apply(comm, 1, function (x) {
x <- x / depth
qt <- quantile(x, winsor.qt)
x[x >= qt] <- qt
round(x * depth)
})
comm <- t(comm)
}
###############################################################################
# Filter to remove very rare taxa
filter.ind <- rowMeans(comm != 0) >= prev.filter & rowSums(comm) >= abund.filter
names(filter.ind) <- rownames(comm)
if (verbose) cat(sum(!filter.ind), ' features are filtered!\n')
comm <- comm[filter.ind, ]
sample.no <- ncol(comm)
otu.no <- nrow(comm)
row.names <- rownames(comm)
depth <- colSums(comm)
if (verbose) cat('The data has ', sample.no, ' samples and ', otu.no, ' features will be tested!\n' )
rabund <- colMeans(t(comm) / colSums(comm))
if (is.null(topK)) {
topK <- round(length(rabund) * 0.4) # Very aggressive
}
size.factor <- colSums(comm[order(rabund)[1 : (length(rabund) - topK)], ])
###############################################################################
# Generate samples from posterior distribution (stacking it)
if (is.prior == TRUE) {
if (verbose) cat('Fitting beta mixture ...\n')
if (post.method == 'sample') {
comm.p <- apply(comm, 1, function (x) {
if (prior.dist == 'BetaMix') {
err1 <- try(res <- bbmix.fit.MM(x, depth))
}
if (prior.dist == 'ZIBB') {
err1 <- try(res <- zibb.fit.MM(x, depth))
}
# Handle error
if (class(err1) != 'try-error') {
if (prior.dist == 'BetaMix') {
prop1.1 <- rbeta(sample.no * post.sample.no, shape1 = x + res$shape1.1, shape2 = res$shape1.2 + depth - x)
prop1.2 <- rbeta(sample.no * post.sample.no, shape1 = x + res$shape2.1, shape2 = res$shape2.2 + depth - x)
prop <- ifelse(runif(sample.no * post.sample.no) <= res$q1, prop1.1, prop1.2)
}
if (prior.dist == 'ZIBB') {
prop1.2 <- rbeta(sample.no * post.sample.no, shape1 = x + res$shape2.1, shape2 = res$shape2.2 + depth - x)
prop <- ifelse(runif(sample.no * post.sample.no) <= res$q1, 0, prop1.2)
}
} else {
prop <- x / depth
v <- var(prop)
m <- mean(prop)
a1 <- ((1 - m) / v - 1 / m) * m ^ 2
a2 <- a1 * (1 / m - 1)
if (is.na(a1) | a1 < 0) {
# uniform prior
prop <- rbeta(sample.no * post.sample.no, shape1 = x + 1, shape2 = otu.no + depth - x)
} else {
prop <- rbeta(sample.no * post.sample.no, shape1 = x + a1, shape2 = a2 + depth - x)
}
}
return(prop)
})
comm.p <- t(comm.p)
comm.p.list <- list()
st <- 1
end <- sample.no
for (i in 1 : post.sample.no) {
comm.p.list[[i]] <- comm.p[, st : end]
st <- st + sample.no
end <- end + sample.no
}
} else {
EV <- apply(comm, 1, function (x) {
if (prior.dist == 'BetaMix') {
err1 <- try(res <- bbmix.fit.MM(x, depth))
}
if (prior.dist == 'ZIBB') {
err1 <- try(res <- zibb.fit.MM(x, depth))
}
# Handle error
if (class(err1) != 'try-error') {
if (prior.dist == 'BetaMix') {
Ep <- res$q1 * (x + res$shape1.1) / (depth + res$shape1.1 + res$shape1.2) +
(1 - res$q1) * (x + res$shape2.1) / (depth + res$shape2.1 + res$shape2.2)
Vp <- res$q1 * ((x + res$shape1.1)^2 / (depth + res$shape1.1 + res$shape1.2)^2 +
(x + res$shape1.1) * (res$shape1.2 + depth - x) / (depth + res$shape1.1 + res$shape1.2)^2 /
(depth + res$shape1.1 + res$shape1.2 + 1)) +
(1 - res$q1) * ((x + res$shape2.1)^2 / (depth + res$shape2.1 + res$shape2.2)^2 +
(x + res$shape2.1) * ( res$shape2.2 + depth - x) / (depth + res$shape2.1 + res$shape2.2)^2 /
(depth + res$shape2.1 + res$shape2.2 + 1)) - Ep^2
}
if (prior.dist == 'ZIBB') {
Ep <- (1 - res$q1) * (x + res$shape2.1) / (depth + res$shape2.1 + res$shape2.2)
Vp <- (1 - res$q1) * ((x + res$shape2.1)^2 / (depth + res$shape2.1 + res$shape2.2)^2 +
(x + res$shape2.1) * ( res$shape2.2 + depth - x) / (depth + res$shape2.1 + res$shape2.2)^2 /
(depth + res$shape2.1 + res$shape2.2 + 1)) - Ep^2
}
} else {
prop <- x / depth
v <- var(prop)
m <- mean(prop)
a1 <- ((1 - m) / v - 1 / m) * m ^ 2
a2 <- a1 * (1 / m - 1)
if (is.na(a1) | a1 < 0) {
# uniform prior
Ep <- (x + 1) / (depth + otu.no)
Vp <- (x + 1) * (depth - x + otu.no - 1) / (otu.no + depth)^2 /
(otu.no + depth + 1)
} else {
Ep <- (x + a1) / (depth + a1 + a2)
Vp <- (x + a1) * (depth - x + a2) / (depth + a1 + a2)^2 /
(depth + a1 + a2 + 1)
}
}
return(c(Ep, Vp))
})
EV <- t(EV)
comm.p.list <- list()
comm.v.list <- list()
comm.p.list[[1]] <- EV[, 1 : sample.no]
comm.v.list[[1]] <- EV[, (sample.no + 1) : (2 * sample.no)]
post.sample.no <- 1
}
} else {
comm.p.list <- list()
comm.p.list[[1]] <- t(t(comm) / depth)
post.sample.no <- 1
}
# Replace zeros or extremely small values for log calculation
for (i in 1 : post.sample.no) {
temp <- comm.p.list[[i]]
temp[temp <= min.prop] <- min.prop
comm.p.list[[i]] <- temp
}
###############################################################################
# Covariate space (including intercept)
if (!is.null(strata)) {
strata <- factor(strata)
}
if (is.null(adj.name)) {
M0 <- model.matrix(~ 1, meta.dat)
} else {
data0 <- meta.dat[, c(adj.name), drop = FALSE]
M0 <- model.matrix( ~ ., data0)
}
data1 <- meta.dat[, c(grp.name), drop = FALSE]
M1 <- model.matrix( ~ ., data1)[, -1, drop = FALSE] # No intercept
M01 <- cbind(M0, M1)
# QR decompostion
qrX0 <- qr(M0, tol = 1e-07)
Q0 <- qr.Q(qrX0)
Q0 <- Q0[, 1:qrX0$rank, drop = FALSE]
H0 <- (Q0 %*% t(Q0))
qrX1 <- qr(M1, tol = 1e-07)
Q1 <- qr.Q(qrX1)
Q1 <- Q1[, 1:qrX1$rank, drop = FALSE]
qrX01 <- qr(M01, tol = 1e-07)
Q01 <- qr.Q(qrX01)
Q01 <- Q01[, 1:qrX01$rank, drop = FALSE]
R0 <- as.matrix(resid(lm(Q1 ~ Q0 - 1)))
pX0 <- ncol(Q0)
pX1 <- ncol(Q1)
pX01 <- ncol(Q01)
df.model <- pX01 - pX0
df.residual <- sample.no - pX01
func.no <- length(link.func)
###############################################################################
# Perform multiple stage normalization
norm.ind <- NULL
for (i in 1:stage.no) {
if (verbose == TRUE)
cat('Stage ', i, ' ')
# Reference proportion
divisor <- size.factor / depth
# Create the giant Y matrix
Y <- matrix(NA, sample.no, func.no * otu.no * post.sample.no)
if (post.method == 'mean') {
D <- matrix(NA, sample.no, func.no * otu.no * post.sample.no)
}
# Change order - func.no * otu.no
for (k in 1 : post.sample.no) {
for (j in 1 : func.no) {
func <- link.func[[j]]
comm.p <- comm.p.list[[k]]
if (post.method == 'mean') {
dfunc <- link.d.func[[j]]
comm.v <- comm.v.list[[k]]
}
Y[, (k - 1) * func.no * otu.no + func.no * (0 : (otu.no - 1)) + j] <-
func(t(comm.p) / divisor) # No scaling first
if (post.method == 'mean') {
# row - sample, column - features
D[, (k - 1) * func.no * otu.no + func.no * (0 : (otu.no - 1)) + j] <-
(dfunc(t(comm.p) / divisor) / divisor)^2 * t(comm.v)
}
}
}
Y <- t(Y)
if (post.method == 'sample') {
TSS <- rowSums(Y^2)
MSS01 <- rowSums((Y %*% Q01)^2)
MSS0 <- rowSums((Y %*% Q0)^2)
} else {
w0 <- colSums((Q0 %*% t(Q0))^2)
w01 <- colSums((Q01 %*% t(Q01))^2)
TSS <- rowSums(Y^2) + colSums(D)
MSS01 <- rowSums((Y %*% Q01)^2) + colSums(w01 * D)
MSS0 <- rowSums((Y %*% Q0)^2) + colSums(w0 * D)
}
MSS <- (MSS01 - MSS0)
RSS <- (TSS - MSS01)
perm.ind <- getPermuteMatrix(perm.no, sample.no, strata = strata)
perm.no <- nrow(perm.ind)
MRSSp <- sapply(1 : perm.no, function(ii) {
if (verbose) {
if (ii %% 10 == 0) cat('.')
}
Rp <- R0[perm.ind[ii, ], , drop = FALSE]
# Project to the reisdual space
Rp <- Rp - H0 %*% Rp
qrRp <- qr(Rp, tol = 1e-07)
Q1p <- qr.Q(qrRp)
Q1p <- Q1p[, 1:qrRp$rank, drop = FALSE]
if (post.method == 'sample') {
MSS01p <- MSS0 + rowSums((Y %*% Q1p)^2)
} else {
w1p <- colSums((Q1p %*% t(Q1p))^2)
MSS01p <- MSS0 + rowSums((Y %*% Q1p)^2) + colSums(w1p * D)
}
MSSp <- (MSS01p - MSS0)
RSSp <- (TSS - MSS01p)
c(MSSp, RSSp)
})
unit <- func.no * otu.no * post.sample.no
MSSp <- MRSSp[1 : unit, ]
RSSp <- MRSSp[(unit + 1) : (2 * unit), ]
# EB is based on the aggregated RSS
RSS.m <- array(RSS, c(func.no, otu.no, post.sample.no))
RSS.m <- t(apply(RSS.m, c(1, 2), mean)) # otu.no by func.no
F0.m <- array((MSS / df.model) / (RSS / df.residual), c(func.no, otu.no, post.sample.no))
F0.m <- t(apply(F0.m, c(1, 2), mean)) # otu.no by func.no
R2.m <- array(MSS / TSS, c(func.no, otu.no, post.sample.no))
R2.m <- t(apply(R2.m, c(1, 2), mean)) # otu.no by func.no
###############################################################################
# Variance moderation - to be improved
if (variance.EB) {
# hist(RSS.m / df.residual)
# out <- apply(RSS.m, 2, function (x) {
# out <- limma::squeezeVar(x / df.residual, df.residual)
# c(out$var.prior, out$df.prior)
# })
# var.prior <- out[1, ]
# df.prior <- out[2, ]
#
# df.prior[df.prior > 1E6] <- 1E6
# df.prior[df.prior < 1E-6] <- 1E-6
var.prior <- colMeans(RSS.m / df.residual, na.rm = TRUE)
F0 <- (MSS / df.model) / ((RSS + df.prior * var.prior) / (df.prior + df.residual))
Fp <- (MSSp / df.model) / ((RSSp + df.prior * var.prior) / (df.prior + df.residual))
# Expectation of F0 and Fp
F0 <- array(F0, c(func.no, otu.no, post.sample.no))
Fp <- array(Fp, c(func.no, otu.no, post.sample.no, perm.no))
F0 <- apply(F0, c(1, 2), mean) # func.no * otu.no
Fp <- apply(Fp, c(1, 2, 4), mean) # func.no * otu.no * perm.no
} else {
F0 <- (MSS / df.model) / (RSS / df.residual)
Fp <- (MSSp / df.model) / (RSSp / df.residual)
# Expectation of F0 and Fp
F0 <- array(F0, c(func.no, otu.no, post.sample.no))
Fp <- array(Fp, c(func.no, otu.no, post.sample.no, perm.no))
F0 <- apply(F0, c(1, 2), mean) # func.no * otu.no
Fp <- apply(Fp, c(1, 2, 4), mean) # func.no * otu.no * perm.no
}
###############################################################################
# Omnibus test by taking maximum
F0 <- apply(F0, 2, stats.combine.func)
Fp <- apply(Fp, c(2, 3), stats.combine.func) # otu.no by perm.no
if (verbose) cat('\n')
if (mean(is.na(F0)) >= 0.1) {
warning('More than 10% observed F stats have NA! Please check! \n')
}
if (mean(is.na(Fp)) >= 0.1) {
warning('More than 10% permuted F stats have NA! Please check! \n')
}
na.ind <- is.na(F0)
F0 <- F0[!na.ind]
Fp <- Fp[!na.ind, ]
rabund0 <- rabund[!na.ind]
which.nan.ind <- which(!na.ind)
p.adj.fdr <- perm_fdr_adj(F0, Fp)
if (i == stage.no) {
break
} else {
# recalculating the size factor
if (mean(p.adj.fdr <= stage.fdr) > stage.max.pct) {
ind <- order(p.adj.fdr, -rabund0)[1 : round(length(p.adj.fdr) * stage.max.pct)]
} else {
ind <- which(p.adj.fdr < stage.fdr)
}
size.factor <- colSums(comm[setdiff(1:nrow(comm), which.nan.ind[ind]), ])
norm.ind <- cbind(norm.ind, !(1:nrow(comm) %in% which.nan.ind[ind]))
}
}
###############################################################################
# Permutation-based false discovery rate control
# Tree-based power improvement
if (is.tree.fdr) {
if (verbose) {
cat("Perform Tree-based FDR control ...\n")
}
# Convert pseudo-F to pseudo-P
if (is.null(tree)) {
stop('Tree FDR requires a tree of class "phylo"!\n')
}
test.obs <- list()
test.perm <- list()
test.obs$p.value <- 1 - pf(F0, df.model, df.residual)
test.perm$p.value <- 1 - pf(Fp, df.model, df.residual)
# Ignore the sign
test.obs$e.sign <- sign(test.obs$p.value)
test.perm$e.sign <- sign(test.perm$p.value)
D <- cophenetic(tree)
D <- D[rownames(comm), rownames(comm)]
tree.fdr.obj <- TreeFDR(X = NULL, Y = NULL, test.obs = test.obs, test.perm = test.perm, D = D,
eff.sign = FALSE, alt.FDR = 'Permutation', verbose = verbose, ...)
p.adj.tree.fdr <- tree.fdr.obj$p.adj
if (verbose) {
cat('FWER control has not used the phylogenetic information!\n')
}
} else {
# scaled F stat
tree.fdr.obj <- NULL
p.adj.tree.fdr <- NULL
}
if (is.fwer) {
p.adj.fwer <- perm_fwer_adj2(F0, Fp)
}
p.raw <- rowMeans(cbind(Fp, F0) >= F0)
p.raw <- na.pad(p.raw, na.ind)
p.adj.fdr <- na.pad(p.adj.fdr, na.ind)
names(p.raw) <- names(p.adj.fdr) <- rownames(R2.m) <- rownames(RSS.m) <- rownames(F0.m) <- row.names
colnames(R2.m) <- colnames(F0.m) <- colnames(RSS.m) <- paste0('Func', 1 : func.no)
if (is.tree.fdr) {
p.adj.tree.fdr <- na.pad(p.adj.tree.fdr, na.ind)
names(p.adj.tree.fdr) <- row.names
}
if (is.fwer) {
p.adj.fwer <- na.pad(p.adj.fwer, na.ind)
names(p.adj.fwer) <- row.names
} else {
p.adj.fwer <- NULL
}
if (!return.comm) {
comm <- NULL
}
if (!return.perm.F) {
Fp <- NULL
}
if (verbose) cat('Completed!\n')
return(list(call = this.call, comm = comm, filter.ind = filter.ind,
R2 = R2.m, F0 = F0.m, Fp = Fp, RSS = RSS.m, df.model = df.model, df.residual = df.residual,
p.raw = p.raw, p.adj.fdr = p.adj.fdr, p.adj.tree.fdr = p.adj.tree.fdr, p.adj.fwer = p.adj.fwer,
tree.fdr.obj = tree.fdr.obj, norm.ind = norm.ind))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.