7.修改注释表做物种组成图_taxonomy.tsv-CSDN博客

本文链接：https://blog.csdn.net/weixin_44282754/article/details/136840883

我们通过qiime2注释出来的很多都是未分类的注释，但是在注释的结果里会显示NA
但是看看高分文章，细心的人就会发现,出现了UC Cryomorphaceae.也就是说如果注释到了属的分类为unclassified或者uncultured，他会直接找上一级科Cryomorphaceae，缩写一下就是UC Cryomorphaceae。

这意味着我们要适当对注释表进行修改，修改成下图的样子：

话不多说，先准备好你的Otu表table.tsv和注释表taxonomy.tsv
这两个表是qiime2的结果，下图所示

table.tsv

taxonomy.tsv

然后话不多说上代码！

data <- read.table('table.tsv',header = T,encoding='UTF-8',check.names = F,stringsAsFactors = F)
colnames(data)[1] <- "Feature ID"
library(readr)
# Read the .tsv file
taxonomy <- read_tsv("taxonomy.tsv")

# 合并相对丰度表和分类表
merged_data <- merge(data, taxonomy, by = "Feature ID") # 假设otu_id是连接两个表的键


# 确保已经加载了必要的包
library(dplyr)
library(tidyr)

# 使用separate()函数分列，这次我们想将数据分为七列
merged_data <- merged_data %>% 
  separate(col = Taxon, into = c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), sep = ";", remove = TRUE, extra = "merge")



fill_na_or_uncultured_with_all_previous_info <- function(row) {
  levels <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
  
  # 遍历每个级别，从“Phylum”开始
  for (i in 2:length(levels)) {
    # 检查是否为NA或包含"uncultured"
    if (is.na(row[i]) || grepl("uncultured", row[i], ignore.case = TRUE)) {
      # 初始化当前级别的新注释
      new_annotation <- paste0(tolower(substr(levels[i], 1, 1)), "__", ifelse(is.na(row[i]), "unclassified", "uncultured"))
      
      # 向上遍历所有上级分类
      for (j in (i-1):1) {
        # 如果上一级也是unclassified或uncultured，也要加上它
        if (is.na(row[j]) || grepl("uncultured", row[j], ignore.case = TRUE)) {
          new_annotation <- paste(new_annotation, paste0(tolower(substr(levels[j], 1, 1)), "__unclassified"), sep = ";")
        } else {
          new_annotation <- paste(new_annotation, row[j], sep = ";")
          break # 找到非unclassified和非uncultured的级别后停止
        }
      }
      
      row[i] <- new_annotation
    }
  }
  
  return(row)
}

# 应用这个函数到每一行数据
merged_data_corrected <- t(apply(merged_data[, c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")], 1, fill_na_or_uncultured_with_all_previous_info))

# 将处理后的数据转换回数据框并恢复列名
merged_data_corrected <- as.data.frame(merged_data_corrected)

merged_data[, 96:102] <- merged_data_corrected

# 初始化一个列表来保存每个分类级别的相对丰度数据框
relative_abundance_list <- list()

# 定义所有的分类级别
levels <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")

# 遍历每个分类级别
for(level in levels) {
  # 为当前级别创建相对丰度数据框
  data_selected <- merged_data %>%
    select(matches(level), 2:95)  # 选择对应分类级别和相对丰度列
  
  # 计算相对丰度
  data_relative_abundance <- data_selected %>%
    mutate(across(2:ncol(data_selected), ~ ./sum(.x, na.rm = TRUE), .names = "rel_{.col}"))
  
  # 将结果保存到列表中
  relative_abundance_list[[level]] <- data_relative_abundance
}

# 初始化一个新列表用于保存加和合并后的数据
summed_abundance_list <- list()

# 遍历每个分类级别
for(level in levels) {
  # 对指定的分类级别进行加和合并
  level_sum <- relative_abundance_list[[level]] %>%
    group_by(.data[[level]]) %>%  # 按当前级别分组
    summarise(across(starts_with("rel_"), sum, na.rm = TRUE))  # 对所有相对丰度列进行加和，忽略NA值
  
  # 将结果保存到新列表中
  summed_abundance_list[[level]] <- level_sum
}

# 将列表中的数据框写入一个Excel文件
write_xlsx(summed_abundance_list, path = "七个分级合并.xlsx")