物种相对丰度热图

# 加载必要的包
library(tidyverse)
library(ggplot2)
library(viridis)

# 读取相对丰度表
# 假设你的数据文件名为 abundance_table.csv，第一列是分类单元，其余列是样本丰度
abundance_table <- read.csv("D:/1、博士期间/9、博士-实验以及项目/11、宏基因组-金华武义-凌恩生物/宏基因分析-只有真菌/BZ20250601995931-zhutianyuan-hangtingting/07.Annotation/NR/nr.tax.profile.txt", header = TRUE, sep = "\t", row.names = 1)

# 提取门、科、属信息
abundance_table <- abundance_table %>%
  rownames_to_column(var = "Taxon") %>%
  mutate(
    Phylum = sapply(strsplit(Taxon, ";"), function(x) x[3]),
    Family = sapply(strsplit(Taxon, ";"), function(x) x[6]),
    Genus = sapply(strsplit(Taxon, ";"), function(x) x[7])
  )

# 定义一个函数，用于绘制指定分类水平下前 10 物种的柱状图
plot_top10 <- function(level, data) {
  # 找出数值列
  numeric_cols <- data %>% 
    select(-Taxon, -Phylum, -Family, -Genus) %>% 
    names()
  
  # 按指定分类水平分组，计算每个处理组的总丰度
  top10 <- data %>%
    group_by(across(all_of(level))) %>%
    summarise(across(all_of(numeric_cols), sum)) %>%
    pivot_longer(-all_of(level), names_to = "Treatment", values_to = "Abundance") %>%
    group_by(Treatment) %>%
    top_n(10, Abundance)
  
  # 绘制柱状图
  if (level == "Family") {
    # 对于科水平，使用 viridis 调色板
    p <- ggplot(top10, aes(x = Treatment, y = Abundance, fill = .data[[level]])) +
      geom_col(position = "stack") +
      labs(title = paste("Top 10", level, "in Different Treatments"),
           x = "Treatment",
           y = "Relative Abundance") +
      scale_fill_viridis_d() +  # 使用 viridis 调色板
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
  } else {
    # 对于其他水平，使用 Set3 调色板
    p <- ggplot(top10, aes(x = Treatment, y = Abundance, fill = .data[[level]])) +
      geom_col(position = "stack") +
      labs(title = paste("Top 10", level, "in Different Treatments"),
           x = "Treatment",
           y = "Relative Abundance") +
      scale_fill_brewer(palette = "Set3") +  # 使用不同颜色区分物种
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
  }
  return(p)
}

# 在门水平绘制柱状图
phylum_plot <- plot_top10("Phylum", abundance_table)
print(phylum_plot)

# 在科水平绘制柱状图
family_plot <- plot_top10("Family", abundance_table)
print(family_plot)

# 在属水平绘制柱状图
genus_plot <- plot_top10("Genus", abundance_table)
print(genus_plot)



  #属作图
  # 按指定分类水平分组，计算每个处理组的总丰度
  top10 <- data %>%
    group_by(across(all_of(level))) %>%
    summarise(across(all_of(numeric_cols), sum)) %>%
    pivot_longer(-all_of(level), names_to = "Treatment", values_to = "Abundance") %>%
    group_by(Treatment) %>%
    top_n(10, Abundance)
  
  if (level == "Genus") {
    p <- ggplot(top10, aes(x = Treatment, y = Abundance, fill = .data[[level]])) +
      geom_col(position = "stack") +
      labs(title = paste("Top 10", level, "in Different Treatments"),
           x = "Treatment",
           y = "Relative Abundance") +
      scale_fill_viridis_d() +  # 使用 viridis 调色板
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
  } else {
    p <- ggplot(top10, aes(x = Treatment, y = Abundance, fill = .data[[level]])) +
      geom_col(position = "stack") +
      labs(title = paste("Top 10", level, "in Different Treatments"),
           x = "Treatment",
           y = "Relative Abundance") +
      scale_fill_brewer(palette = "Set3") +  # 使用 Set3 调色板
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 45, hjust = 1))
  }
  return(p)
}

# 在属水平绘制柱状图
genus_plot <- plot_top10("Genus", abundance_table)
print(genus_plot)