Skip to content

Instantly share code, notes, and snippets.

@eclarke eclarke/indvals.R
Last active Dec 4, 2015

Embed
What would you like to do?
Indicator species (indicator value) functions
# 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
You can’t perform that action at this time.