R/mean_relative_abundances.R

In mcalgaro93/HotLoadings: HotLoadings: a fancy way to plot mixOmics loadings and heatmaps

#' Compute Mean Relative Abundances
#'
#' @param top_features A data frame containing \code{n_top} features generated by \code{\link{HotLoadings.top_features}}.
#' @param PSOBJ A phyloseq object.
#' @param Y_name A character string indicating which is the variable associated with component specified. Must be a dicothomous variable.
#' @return A data frame containings \code{n_top} feature associated with a specified \code{component} and information about their mean relative abundance in the 2 \code{Y_name} categories.
#' @seealso \code{\link{HotLoadings.top_features}} to create top features data frame.

HotLoadings.mean_relative_abundances <- function(top_features,PSOBJ,Y_name){
  if(!taxa_are_rows(PSOBJ)){
    PSOBJ@otu_table <- t(PSOBJ@otu_table)
  }
  # Subsetting interesting top n features
  non_other <- match(rownames(top_features),rownames(PSOBJ@otu_table@.Data))
  # Create a data.frame with absolute and relative abundances for top_features and other features as aggregation
  # Aggregate non interesting features
  otuTab_other <- colSums(PSOBJ@otu_table@.Data[-non_other,])
  # Add to interesting features, other features
  otuTab <- rbind(PSOBJ@otu_table@.Data[non_other,],otuTab_other)
  # Find relative abundances
  otuTabRel <- apply(otuTab,2,function(col) col/sum(col))
  rownames(otuTabRel) <- rownames(otuTab) <- c(as.character(top_features$feature),"Other taxa")
  # For each level of Y variable we compute mean relative count value
  comparison_aggr_rel <- aggregate(t(otuTabRel), PSOBJ@sam_data[,Y_name]@.Data, mean)
  # We keep only the levels that are present in top_features
  present_index <- which(comparison_aggr_rel[,1] %in%  unique(top_features[,Y_name]))
  comparison_aggr_rel <- comparison_aggr_rel[present_index,]
  meanRelAbundances_df <- data.frame(rbind(t(comparison_aggr_rel[1,-1]), t(comparison_aggr_rel[2,-1])), rep(comparison_aggr_rel[,1],each = nrow(top_features)+1))
  # Add taxa name as factor variable
  meanRelAbundances_df$feature <- rep(rownames(otuTab),2)
  colnames(meanRelAbundances_df) <- c("Relative_Abundance",Y_name,"feature")
  meanRelAbundances_df[,Y_name] <- factor(meanRelAbundances_df[,Y_name],levels = unique(top_features[,Y_name]),ordered = TRUE)
  return(meanRelAbundances_df[order(meanRelAbundances_df[,Y_name],decreasing = TRUE),])
}