#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.