在R中按Phylum出现次数筛选大型phyloseq熔解表的问题

我正在用R处理基因组数据，需要对一个大型熔解后的phyloseq表做子集筛选——移除表中Phylum（门）出现次数少于100000次对应的所有行。自己写了自定义函数但运行报错，不清楚有没有更简便的实现方式，我以为没法直接用filter()。

自定义函数与报错 #

自定义函数代码 #

phylum_subset <- function(x = melt.ALKSS_few, #melted physeq object 
                          Count = melt.ALKSS$Phylum, #counting phyla
                          Value = 1000   #minimum number of OTUs
                          ){
  phyla.table <- table(x$Count) 
  
  for(Count in x){if(phyla.table[Count]<=100000)
                        subset(x,Phylum != Count)
  }
}

调用代码 #

melt.ALKSS_few.count <- phylum_subset(x = melt.ALKSS_few,Count = melt.ALKSS_few$Phylum,Value = 100000)

报错信息 #

Error in if (phyla.table[Count] <= 1e+05) subset(x, Count != Phylum_col) : 
  the condition has length > 1

数据集片段 #

> dput(droplevels(head(melt.ALKSS_few)))
structure(list(OTU = c("44c21e29adae97a53247abbd73978395", "0f18144d308ada95632ab5193d92073f", 
"d829bee4984f82ffc2453212157caf96", "0f18144d308ada95632ab5193d92073f", 
"0ddcd311e02f742e2e0e61ce02cf9c29", "120eba657e42a11a5c29f97b90f02035"
), Sample = c("S438", "S680", "S437", "S345", "S454", "S513"), 
    Abundance = c(10755, 9568, 8186, 7621, 7506, 7501), BarcodeSequence = c("CATTTTAGGACT", 
    "CGGAATAGAGTA", "CATTTTAGAGTA", "TATAATGGACCA", "CGGAATTGGCAT", 
    "GACGACGGACCA"), PrimerDesc = c("16S", 
    "16S", "16S", "16S", "16S", "16S"), SampleName = c("06222021KC-2-R", 
    "09292021KC-2-R", "06222021KC-1-R", "06032021KC-1-R", "06292921KC-3-R", 
    "06302021KC-3-R"), Project = c("16SLBSKR1-", "16SLBSKR2-", 
    "16SLBSKR1-", "16SLBSKR1-", "16SLBSKR1-", "16SLBSKR2-"), 
    Number = c("456", "694", "455", "363", "471", "491"), Date = c("6_22_2021", "9_29_2021", "6_22_2021", "6_3_2021", 
    "6_29_2021", "6_30_2021"), Year = c(2021L, 2021L, 2021L, 
    2021L, 2021L, 2021L), Season = c("Summer", "Fall", "Summer", 
    "Summer", "Summer", "Summer"), sample_Species = c("Little_Bluestem", 
    "Little_Bluestem", "Little_Bluestem", "Little_Bluestem", 
    "Little_Bluestem", "Little_Bluestem"), SoloOrMixed = c("Solo", 
    "Mixed", "Solo", "Mixed", "Mixed", "Solo"), Location = c("Tyler_SP", "Hy_180", "Tyler_SP", "Roadside_Hy67", 
    "Copper_Breaks_SP", "Caprock_Canyons_SP"), Ecoregion = c("South_Central_Plains", 
    "South_Central_Plains", "South_Central_Plains", "Edwards_Plateau", 
    "Southwestern_Tablelands", "Southwestern_Tablelands"), Habitat = c("Forest", 
    "Roadside", "Forest", "Roadside", "Roadside", "AridRock"), 
    Source = c("Root", "Root", "Root", "Root", "Root", "Root"
    ),  PrecipMonth = c(96.65, 
    37.45, 96.65, 125.94, 125.01, 153.94), PrecipDaysSince = c(1L, 
    1L, 1L, 1L, 1L, 0L), pH = c(6.8, 6.7, 6.8, 8, 8, 7.8), EC = c(139L, 
    182L, 139L, 161L, 125L, 2370L), NO3 = c(0, 4.4, 0, 0.2, 2.2, 
    3.4), P = c(16L, 17L, 16L, 14L, 5L, 6L), K = c(145L, 84L, 
    145L, 114L, 160L, 65L), Ca = c(3918L, 2159L, 3918L, 27256L, 
    6609L, 16508L), Mg = c(166L, 130L, 166L, 188L, 148L, 95L), 
    S = c(10L, 16L, 10L, 24L, 24L, 14299L), Na = c(4L, 3L, 4L, 
    4L, 4L, 4L), Fe = c(19.76, 17, 19.76, 2.31, 1, 0), Zn = c(2.28, 
    15.1, 2.28, 7.01, 0.8, 0.1), Mn = c(64.16, 19, 64.16, 27.01, 
    15, 6), Cu = c(0.16, 0.2, 0.16, 0.16, 0.2, 0.2), Kingdom = c("d__Bacteria", 
    "d__Bacteria", "d__Bacteria", "d__Bacteria", "d__Bacteria", 
    "d__Bacteria"), Phylum = c("Proteobacteria", "Proteobacteria", 
    "Proteobacteria", "Proteobacteria", "Proteobacteria", "Actinobacteriota"
    ), Class = c("Gammaproteobacteria", "Gammaproteobacteria", 
    "Alphaproteobacteria", "Gammaproteobacteria", "Gammaproteobacteria", 
    "Actinobacteria"), Order = c("Xanthomonadales", "Pseudomonadales", 
    "Rhizobiales", "Pseudomonadales", "Pseudomonadales", "Streptomycetales"
    ), Family = c("Rhodanobacteraceae", "Pseudomonadaceae", "Xanthobacteraceae", 
    "Pseudomonadaceae", "Pseudomonadaceae", "Streptomycetaceae"
    ), Genus = c("Rhodanobacter", "Pseudomonas", "Bradyrhizobium", 
    "Pseudomonas", "Pseudomonas", "Streptomyces"), Species = c(NA_character_, 
    NA_character_, NA_character_, NA_character_, NA_character_, 
    NA_character_)), row.names = c(2352002L, 511171L, 7348565L, 
510815L, 468295L, 621043L), class = "data.frame")

问题分析与解决方案 #

原函数问题点 #

• 参数逻辑混乱：Count参数传入的是melt.ALKSS_few$Phylum，但函数里写x$Count，而x是数据框，根本没有名为Count的列，导致table(x$Count)生成空表。
• for循环错误：for(Count in x)会遍历数据框的每一列，Count变成整列向量，导致phyla.table[Count]返回长度大于1的结果，if条件判断不接受这种情况。
• 无返回值：函数里的subset操作没有赋值给变量，循环结束后也没有返回处理后的数据集，调用函数会得到NULL。

简便实现方法 #

不需要写复杂的自定义函数，用base R或者dplyr都能快速实现：

方法1：base R 实现

# 统计每个Phylum的出现次数
phylum_counts <- table(melt.ALKSS_few$Phylum)
# 筛选出次数≥100000的Phylum
keep_phyla <- names(phylum_counts[phylum_counts >= 100000])
# 提取符合条件的行
melt.ALKSS_few_filtered <- melt.ALKSS_few[melt.ALKSS_few$Phylum %in% keep_phyla, ]

方法2：dplyr 实现（更直观）

library(dplyr)

# 方式1：添加临时计数列后筛选
melt.ALKSS_few_filtered <- melt.ALKSS_few %>%
  group_by(Phylum) %>%
  mutate(phylum_count = n()) %>%
  ungroup() %>%
  filter(phylum_count >= 100000) %>%
  select(-phylum_count) # 可选：移除临时计数列

# 方式2：先统计再筛选，更高效
keep_phyla <- melt.ALKSS_few %>%
  count(Phylum) %>%
  filter(n >= 100000) %>%
  pull(Phylum)

melt.ALKSS_few_filtered <- melt.ALKSS_few %>%
  filter(Phylum %in% keep_phyla)

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