Last active
December 4, 2015 22:21
-
-
Save eclarke/ae1dc3bb5bebb63ffb0f to your computer and use it in GitHub Desktop.
Indicator species (indicator value) functions
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Indicator value functions ----------------------------------------------- | |
#' Returns the indicator value (Dufrene, 1997) for a given row of species counts | |
#' along with a vector of class assignments. | |
#' @param row: vector of counts (usually a row in a counts matrix) | |
#' @param class: which level in the grouping variable to test | |
#' @param classes: factor describing the grouping of the counts vector | |
.indval <- function(row, class, classes) { | |
idxs <- classes == class | |
A.ij <- sum(row[idxs]) / sum(row) | |
B.ij <- sum(row[idxs] > 0) / sum(idxs) | |
return(A.ij*B.ij) | |
} | |
#' Returns the indicator value and the p-value (as determined by permutation) | |
#' for each feature (row) in a species matrix. | |
#' @param mat: a species matrix (features=rows, samples=columns) | |
#' @param class: which level in the grouping vector to test | |
#' @param classes: a factor describing the grouping of the samples in the matrix | |
#' @param n.permutations: number of permutations to do to calculate the p value | |
#' @param p.adjust.method: which method to use to account for multiple testing | |
#' (default: 'fdr') | |
get_indvals <- function(cts, group.factor, test.group, n.permutations=10000, p.adjust.method='fdr') { | |
message("Performing indicator species permutation test") | |
stopifnot(test.group %in% levels(group.factor)) | |
stopifnot(ncol(cts) == length(group.factor)) | |
cl <- makePSOCKcluster(detectCores()) | |
indvals <- adply(cts, 1, function(x) { | |
indval <- .indval(x, test.group, group.factor) | |
rand.indvals <- parSapply(cl=cl, c(1:n.permutations), function(y, test.group, group.factor, .indval){ | |
.indval(x, test.group, sample(group.factor)) | |
}, test.group, group.factor, .indval) | |
p.value <- sum(rand.indvals >= indval) / length(rand.indvals) | |
return(data.frame(indval=indval, p.value=p.value)) | |
}) | |
stopCluster(cl) | |
indvals$adj.p.value <- p.adjust(indvals$p.value, method=p.adjust.method) | |
indvals$otu <- indvals$X1 | |
indvals$X1 <- NULL | |
return(indvals) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment