Phyloseq PCoA Loadings

get_loadings.R

get_loadings <- function(physeq, ord, method = c("pearson", "spearman"), scaled=TRUE, cor = FALSE, axes=1:2, taxranks=6:7, top = NULL){
	if(!all(rownames(ord$vectors) == sample_names(physeq)))
		stop("Taxa names do not match")
	counts <- t(as(otu_table(physeq),"matrix"))
	scores <- ord$vectors
	cov <- cov(counts, scores, method=method[1])
	if(cor){
		cov <- cor(counts, scores, method=method[1])
	}
	if(scaled){
		for(i in 1:min(ncol(cov), ncol(scores))){
			cov[,i] <- cov[,i]/ (max(range(cov[,i]))/max(range(scores[,i])))
		}
	}
	res <- data.frame(	cov[,axes],
						counts=taxa_sums(physeq), 
						meanrel = taxa_sums(transform_sample_counts(physeq, function(x) x/sum(x)))/nsamples(physeq), 
						length=apply(cov[,axes],1,function(x) sqrt(sum(x^2))),
						as(tax_table(physeq),"matrix")[,taxranks] )
	if(is.null(top))
		return(res)
	res[res$length >= sort(res$length, decreasing=TRUE)[top],]
}