如何用phyloseq的estimate_richness批量计算各门ASV丰富度？

问题背景 #

需要批量获取每个站位下各门（Phylum）的ASV丰富度，目前只能手动逐个对目标门进行子集化处理，处理前20个门时操作繁琐。现有单门类处理代码如下：

ps.Prymnesiophyceae <- subset_taxa(ps, Phylum == "Prymnesiophyceae")  
alpha_diversity<-estimate_richness(ps.Prymnesiophyceae,measure=c("Shannon","Observed"))  
H<-alpha_diversity$Shannon  
S1<-alpha_diversity$Observed  
S<-log(S1)  
evenness<-H/S
alpha<-cbind(Shannon=H,Richness=S1,Evenness=evenness,sample_data(Prymnesiophyceae))

补充问题（按建议尝试后出现样本名匹配错误） #

尝试批量处理后，发现样本名格式不匹配，无法正确合并元数据，相关数据输出如下：

> dput(head(sample_names(ps.transect), n=2)) 
c("2-1-DCM_S21_L001_R1_001.fastq", "2-1-SA_S9_L001_R1_001.fastq" )

> dput(head(alpha, n=2)) 
structure(list(Observed = c(31, 25), Shannon = c(2.84184012598765, 
2.53358345702604), taxon = c("Prymnesiophyceae", "Prymnesiophyceae" ), sample_id = c("X2.1.DCM_S21_L001_R1_001.fastq", "X2.1.SA_S9_L001_R1_001.fastq" ), S = c(3.43398720448515,
3.2188758248682), evenness = c(0.827562817437384, 
0.787101955736294)), row.names = c("X2.1.DCM_S21_L001_R1_001.fastq",  "X2.1.SA_S9_L001_R1_001.fastq"), class = "data.frame")

> dput(head(smpl_data, n=1)) 
new("sample_data", .Data = list("001_DCM", 125L, structure(1L, .Label = "DCM", class = "factor"), structure(1L, .Label = "Transect", class = "factor"), structure(1L, .Label = "STZ", class = "factor"), 
    structure(1L, .Label = "STFW", class = "factor"), "Oligotrophic", 
    16L, -149.9978333, -29.997, 130.634, 17.1252, 35.4443, 1025.835008, 
    1.1968, 1e-12, 5.387, 2.8469, 52.26978546, 98.0505, 0, 0, 
    0.02, 0.9, 0, 0, 2069.47, 8.057, 377.3), names = c("Station_neat",  "Depth_our", "Depth_bin", "Loc", "Front", "Water", "Zone", "Bottle",  "Lon", "Lat", "pressure..db.", "Temperature", "Salinity", "Density_kgm.3",  "Fluorescence_ugL", "PAR", "BottleO2_mLL", "CTDO2._mLL", "OxygenSat_.",  "Beam_Transmission", "N_umolL", "NO3_umolL", "PO4_umolL", "SIL_umolL",  "NO2_umolL", "NH4_umolL", "DIC_uMkg", "pH", "pCO2_matm"), row.names = "2-1-DCM_S21_L001_R1_001.fastq", 
    .S3Class = "data.frame")

解决方案 #

1. 批量计算各门类群的Alpha多样性 #

使用purrr包实现批量循环，自动处理指定门类群的alpha多样性计算，无需手动逐个子集化：

library(phyloseq)
library(dplyr)
library(tibble)
library(purrr)
library(stringr)

# 步骤1：筛选前20个最丰富的门（按总相对丰度排序）
top_phyla <- ps %>%
  tax_glom(taxrank = "Phylum") %>% # 按门合并OTU/ASV
  transform_sample_counts(function(x) x / sum(x)) %>% # 转换为相对丰度
  psmelt() %>% # 转换为长格式数据框
  group_by(Phylum) %>%
  summarise(total_abundance = sum(Abundance)) %>%
  arrange(desc(total_abundance)) %>%
  slice(1:20) %>% # 取前20个
  pull(Phylum) # 提取门名称向量

# 步骤2：定义单门类群alpha多样性计算函数
calc_phylum_alpha <- function(target_phylum, ps_object) {
  # 子集化目标门类群
  ps_subset <- subset_taxa(ps_object, Phylum == target_phylum)
  # 计算指定alpha指数
  alpha_results <- estimate_richness(ps_subset, measure = c("Shannon", "Observed")) %>%
    rownames_to_column("sample_id") %>% # 将行名转为样本ID列
    mutate(
      taxon = target_phylum,
      S = log(Observed),
      evenness = Shannon / S
    ) %>%
    rename(Richness = Observed) # 重命名列
  
  return(alpha_results)
}

# 步骤3：批量计算所有目标门的alpha多样性
all_alpha_data <- purrr::map_dfr(top_phyla, calc_phylum_alpha, ps_object = ps)

2. 修复样本名匹配问题 #

出现匹配错误是因为R自动将含特殊字符的样本名转换为合法变量名（添加前缀X、替换-为.），只需统一样本名格式即可：

# 修正alpha数据中的样本名：移除前缀X，将.替换回-
all_alpha_data <- all_alpha_data %>%
  mutate(
    sample_id = str_remove(sample_id, "^X") %>% 
    str_replace_all("\\.", "-")
  )

# 将phyloseq的样本数据转换为普通数据框，添加样本ID列
sample_metadata <- as.data.frame(sample_data(ps)) %>%
  rownames_to_column("sample_id")

# 合并alpha多样性数据与样本元数据
final_alpha_data <- all_alpha_data %>%
  left_join(sample_metadata, by = "sample_id")

执行上述代码后，final_alpha_data将包含所有目标门的Shannon指数、丰富度、均匀度，以及对应的站位等元数据，可直接用于后续分析。

内容的提问来源于stack exchange，提问作者CoccoLoco