R/reshape_metaphlan2_result.R

In yeguanhuav/visual16S: Visualization functions for 16S and Metagenomics analysis

#' reshape_metaphlan2_result
#'
#' This is a function that reshape metaphlan2 result for generating phyloseq object.
#'
#' @param metaphlan2_merged Input a metaphlan2 merged result.
#'
#' @param remove_archaea Default is FALSE. If TRUE, remove Archea.
#'
#' @param remove_viruses Default is FALSE. If TRUE, remove Viruses.
#'
#' @param level Choose a taxonomy level to extract. The level name should be one of "Kingdom",
#' "Phylum", "Class", "Order", "Family", "Genus", "Species". Default is "Species".
#'
#' @param abundance_filtering Set metaphlan results that lower than abundance_filtering to 0.
#' Default is 0.01, metaphlan results are relative abundance.
#'
#' @param rescale Default is FALSE. If TRUE, rescale the metaphlan results after filtering.
#'
#' @export
#'
#' @examples
#' reshape_metaphlan2_result(metaphlan2_merged, level = "Species")

reshape_metaphlan2_result <- function(
  metaphlan2_merged,
  remove_archaea = FALSE,
  remove_viruses = FALSE,
  level = "Species",
  abundance_filtering = 0.01,
  rescale = FALSE
) {
  # 去掉Archaea
  if(remove_archaea) {
    metaphlan2_merged <- metaphlan2_merged %>% dplyr::filter(!str_detect(ID, "Archaea"))
  }
  # 去掉Viruses
  if(remove_viruses) {
    metaphlan2_merged <- metaphlan2_merged %>% dplyr::filter(!str_detect(ID, "Viruses"))
  }
  if(level == "Species") {
    level_retain <- "s__"
    level_remove <- "t__"
  } else if(level == "Genus") {
    level_retain <- "g__"
    level_remove <- "s__"
  } else if(level == "Family") {
    level_retain <- "f__"
    level_remove <- "g__"
  } else if(level == "Order") {
    level_retain <- "o__"
    level_remove <- "f__"
  } else if(level == "Class") {
    level_retain <- "c__"
    level_remove <- "o__"
  } else if(level == "Phylum") {
    level_retain <- "p__"
    level_remove <- "c__"
  } else if(level == "Kingdom") {
    level_retain <- "k__"
    level_remove <- "p__"
  }
  ### 根据level filter metaphlan2结果
  otu_tbl <- metaphlan2_merged %>%
    # 去掉所有level下一级的结果
    dplyr::filter(., !str_detect(.$ID, level_remove)) %>%
    # 去掉level以上的结果
    dplyr::filter(., str_detect(.$ID, level_retain)) %>%
    # 只保留level的名字
    dplyr::mutate(ID = str_remove(ID, paste0(".*", level_retain)))
  # 相对丰度小于abundance_filtering都可以看作是0
  otu_tbl[otu_tbl <= abundance_filtering] <- 0
  # 去掉都是0的level
  otu_tbl <- otu_tbl %>%
    filter_if(is.numeric, any_vars(. != 0)) %>%
    column_to_rownames("ID")
  # Re-scale relative abundance
  if(rescale) {
    otu_tbl <- (otu_tbl %>% convert_to_percentage(row_sum = F)) * 100
  }
  return(otu_tbl)
}