数据分析：提取微生物不同分类学水平谱

微生物具有7个分类学水平的特征，记录下如何提取每个分类学水平的微生物profile。更多知识分享请到 https://zouhua.top/。

knitr::opts_chunk$set(warning = F, message = F)
library(dplyr)
library(tibble)
library(tidyr)

load data

otu_table <- read.csv("origin_data/otu_table.csv")
tax_name <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species", "OTUID")
head(otu_table %>% select(tax_name, everything()))

get profile

get_profile <- function(tag="Species", occurrence=0.2, ncount=10){
  
  # tag="Species"
  # occurrence=0.2
  # ncount=10
  
  prof <- otu_table %>% select(tag, starts_with("A"))
  colnames(prof)[which(colnames(prof) == tag)] <- "tax"
  prof_sum <- prof %>% group_by(tax) %>%
    summarise(across(everything(), sum)) %>%
    ungroup() %>%
    column_to_rownames("tax")
  
  prof_filter <- prof_sum %>% rownames_to_column("Type") %>% 
    filter(apply(dplyr::select(., -one_of("Type")), 1, 
                 function(x){sum(x > 0)/length(x)}) > occurrence) %>%
            data.frame(.) %>% 
    column_to_rownames("Type")
  prof_filter <- prof_filter[rowSums(prof_filter) > ncount, ]
  
  res <- prof_filter %>% rownames_to_column(tag)
 
  return(res) 
}

output_fun <- function(dat=species_filter, tag="Species"){
  
  dirout <- "./profile/"
  if(!dir.exists(dirout)){
    dir.create(dirout)
  }
  
  filename <- paste0(dirout, tag, ".filterd.csv")
  
  write.csv(dat, file = filename, row.names = F)
}

taxonomy

otu_filter <- get_profile(tag = tax_name[8])
output_fun(dat = otu_filter, tag = tax_name[8])

species_filter <- get_profile(tag = tax_name[7])
output_fun(dat = species_filter, tag = tax_name[7])

genus_filter <- get_profile(tag = tax_name[6])
output_fun(dat = genus_filter, tag = tax_name[6])

参考文章如引起任何侵权问题，可以与我联系，谢谢。