Nothing
#' @title sparcc wrapper
#'
#' @description SpiecEasi R package, says in his package that
#' this is "a reimplementation of SparCC algorithm (Friedman et Alm, PLoS Comp Bio, 2012)."
#' Installation of SpiecEasi can sometimes generate errors, so I have included Dr. Huaying Fang's sparcc wrapper
#' as one of the functions in this package for the estimation of co-abundance networks. His code was acquired from
#' CCLasso (Fang et al, Bioinformatics, 2015), provided in GitHub (https://github.com/huayingfang/CCLasso).
#'
#' @param x count data matrix (OTU table)
#' @param imax Number of iterations in the outer loop
#' @param kmax max iteration steps for SparCC
#' @param alpha threshold for strong correlation
#' @param Vmin absolute value of correlations below this threshold are considered zero by the inner SparCC loop.
#'
#' @return This will estimate an association matrix (network) for observed OTU table.
#' @importFrom stats median
sparcc <- function(x, imax = 20, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
# dimension for w (latent variables)
p <- ncol(x);
n <- nrow(x);
# posterior distribution (alpha)
x <- x + 1;
# store generate data
y <- matrix(0, n, p);
# store covariance/correlation matrix
cov.w <- cor.w <- matrix(0, p, p);
indLow <- lower.tri(cov.w, diag = T);
# store covariance/correlation for several posterior samples
covs <- cors <- matrix(0, p * (p + 1) / 2, imax);
for(i in 1:imax) {
# generate fractions from posterior distribution
y <- t(apply(x, 1, function(x)
gtools::rdirichlet(n = 1, alpha = x)));
# estimate covariance/correlation
cov_cor <- SparCC.frac(x = y, kmax = kmax, alpha = alpha, Vmin = Vmin);
# store variance/correlation only low triangle
covs[, i] <- cov_cor$cov.w[indLow];
cors[, i] <- cov_cor$cor.w[indLow];
}
# calculate median for several posterior samples
cov.w[indLow] <- apply(covs, 1, median);
cor.w[indLow] <- apply(cors, 1, median);
#
cov.w <- cov.w + t(cov.w);
diag(cov.w) <- diag(cov.w) / 2;
cor.w <- cor.w + t(cor.w);
diag(cor.w) <- 1;
#
return(list(cov.w = cov.w, cor.w = cor.w));
}
#-------------------------------------------------------------------------------
# SparCC for fractions known
# function: SparCC.frac
# input:
# x ------ nxp fraction data matrix, row is sample, col is variable
# kmax ------ max iteration steps for SparCC. default is 10
# alpha ------ the threshold for strong correlation. default is 0.1
# Vmin ------ minimal variance if negative variance appears. default is 1e-4
# output: a list structure
# cov.w ------ covariance estimation
# cor.w ------ correlation estimation
SparCC.frac <- function(x, kmax = 10, alpha = 0.1, Vmin = 1e-4) {
# Log transformation
x <- log(x);
p <- ncol(x);
# T0 = var(log(xi/xj)) variation matrix
TT <- stats::var(x);
T0 <- diag(TT) + rep(diag(TT), each = p) - 2 * TT;
# Variance and correlation coefficients for Basic SparCC
rowT0 <- rowSums(T0);
var.w <- (rowT0 - sum(rowT0) / (2 * p - 2))/(p - 2);
var.w[var.w < Vmin] <- Vmin;
#cor.w <- (outer(var.w, var.w, "+") - T0 ) /
# sqrt(outer(var.w, var.w, "*")) / 2;
Is <- sqrt(1/var.w);
cor.w <- (var.w + rep(var.w, each = p) - T0) * Is * rep(Is, each = p) * 0.5;
# Truncated correlation in [-1, 1]
cor.w[cor.w <= - 1] <- - 1;
cor.w[cor.w >= 1] <- 1;
# Left matrix of estimation equation
Lmat <- diag(rep(p - 2, p)) + 1;
# Remove pairs
rp <- NULL;
# Left components
cp <- rep(TRUE, p);
# Do loops until max iteration or only 3 components left
k <- 0;
while(k < kmax && sum(cp) > 3) {
# Left T0 = var(log(xi/xj)) after removing pairs
T02 <- T0;
# Store current correlation to find the strongest pair
curr_cor.w <- cor.w;
# Remove diagonal
diag(curr_cor.w) <- 0;
# Remove removed pairs
if(!is.null(rp)) {
curr_cor.w[rp] <- 0;
}
# Find the strongest pair in vector form
n_rp <- which.max(abs(curr_cor.w));
# Remove the pair if geater than alpha
if(abs(curr_cor.w[n_rp]) >= alpha) {
# Which pair in matrix form
t_id <- c(arrayInd(n_rp, .dim = c(p, p)));
Lmat[t_id, t_id] <- Lmat[t_id, t_id] - 1;
# Update remove pairs
n_rp <- c(n_rp, (p + 1) * sum(t_id) - 2 * p - n_rp);
rp <- c(rp, n_rp);
# Update T02
T02[rp] <- 0;
# Which component left
cp <- (diag(Lmat) > 0);
# Update variance and truncated lower by Vmin
var.w[cp] <- solve(Lmat[cp, cp], rowSums(T02[cp, cp]));
var.w[var.w <= Vmin] <- Vmin;
# Update correlation matrix and truncated by [-1, 1]
#cor.w <- (outer(var.w, var.w, "+") - T0 ) /
# sqrt(outer(var.w, var.w, "*")) / 2;
Is <- sqrt(1/var.w);
cor.w <- (var.w + rep(var.w, each = p) - T0) *
Is * rep(Is, each = p) * 0.5;
# Truncated correlation in [-1, 1]
cor.w[cor.w <= - 1] <- - 1;
cor.w[cor.w >= 1] <- 1;
}
else {
break;
}
#
k <- k + 1;
}
# Covariance
Is <- sqrt(var.w);
cov.w <- cor.w * Is * rep(Is, each = p);
#
return(list(cov.w = cov.w, cor.w = cor.w));
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.