微生物污染筛选中，不同维度数据框的cor.test Spearman检验实现

我正在遵循一篇《Nature》文章的工作流程，通过Spearman相关性检验识别数据中的污染物。前期步骤已计算物种流行度（prevalence）并过滤数据，得到两个数据集：

• FilPv1：已识别的污染物数据集（281行×11列），包含流行度低于阈值的物种
• FilPV1NONCont：潜在污染物数据集（1743行×11列），包含流行度高于阈值的物种

当前需完成的核心步骤：若某潜在污染物物种与同批次内任意已识别污染物物种的Spearman相关系数ρ＞0.7，则判定该潜在污染物为真正的污染物。计算相关性时需基于中心化对数比转换（clr）后的微生物相对丰度，使用compositions包的clr函数做转换，用R原生的cor.test做Spearman检验。

我之前尝试用嵌套循环逐个对比两个数据集里物种的Prev列值，代码如下：

for (i in 1:nrow(FilPv1)){
  for (x in 1:nrow(FilPV1NONCont))
  cor.test(FilPv1$Prev[i], FilPV1NONCont$Prev[x], method = "spearman")
}

但因样本量不足报错，意识到这是在拿单个值对比单个值，完全错误。请问该如何正确处理？

两个数据集的结构如下：

> head(FilPV1NONCont)
        Prev     totAb  Kingdom         Phylum              Class           Order           Family        Genus
101571     1 0.8581138 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae Burkholderia
1503054    1 0.3168894 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae Burkholderia
87883      1 0.2786909 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae Burkholderia
60550      1 0.1883047 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae Burkholderia
152480     1 0.3580607 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae Burkholderia
488447     1 0.1884915 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae Burkholderia
            Species      clr cont
101571    ubonensis 7.404605    0
1503054   stagnalis 6.408421    0
87883   multivorans 6.279972    0
60550    pyrrocinia 5.887930    0
152480    ambifaria 6.530571    0
488447  contaminans 5.888922    0

> head(FilPv1)
             Prev        totAb  Kingdom         Phylum              Class           Order           Family
1636423 0.3636364 0.0040465871 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae
1478019 0.3636364 0.0028885249 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae
96344   0.7272727 0.0047268082 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae
82541   0.3636364 0.0014331855 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae
1855616 0.4545455 0.0008442892 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae
2878150 0.6363636 0.0072747175 Bacteria Pseudomonadota Betaproteobacteria Burkholderiales Burkholderiaceae
                   Genus       Species        clr cont
1636423     Burkholderia  sp. MSMB1588  0.1411178    1
1478019      Cupriavidus      sp. KK10 -0.1960100    1
96344        Cupriavidus    oxalaticus  0.2964942    1
82541        Cupriavidus      gilardii -0.8968564    1
1855616 Polynucleobacter sp. MWH-UH25E -1.4260161    1
2878150     Caballeronia    sp. Lep1P3  0.7276490    1

核心逻辑：相关性检验的是同一批次样本中，两个物种的丰度变化趋势的关联，需要用每个物种在所有样本中的丰度向量来计算，而非单个流行度数值。

步骤1：整理正确的丰度数据结构 #

假设你有原始的样本×物种丰度矩阵otu_table（每行是一个样本，每列是一个物种的相对丰度），先筛选出目标物种的丰度：

# 提取两类数据集的物种名
contaminant_species <- FilPv1$Species
potential_species <- FilPV1NONCont$Species

# 从原始丰度表中筛选对应物种的丰度列
contaminant_abundances <- otu_table[, colnames(otu_table) %in% contaminant_species]
potential_abundances <- otu_table[, colnames(otu_table) %in% potential_species]

步骤2：执行clr转换 #

compositions::clr不允许数据中有0，需先添加极小值替换0：

library(compositions)

# 处理0值并做clr转换
contaminant_clr <- clr(contaminant_abundances + 1e-6)
potential_clr <- clr(potential_abundances + 1e-6)

步骤3：批量计算Spearman相关性 #

用cor函数直接生成相关系数矩阵，避免低效的嵌套循环：

# 计算相关性矩阵（行：潜在污染物，列：已识别污染物）
# 转置是因为cor默认计算列之间的相关，我们需要物种（原矩阵的列）之间的相关
cor_matrix <- cor(t(potential_clr), t(contaminant_clr), method = "spearman")

# 筛选出任意一个相关系数>0.7的潜在污染物
potential_contaminants <- rownames(cor_matrix)[apply(cor_matrix, 1, function(x) any(x > 0.7))]

特殊情况处理：若clr列是物种的丰度向量 #

如果你的数据框中clr列存储的是该物种在所有样本中的clr丰度向量（每个行对应一个物种的向量），可以直接提取向量计算：

# 提取两类物种的clr丰度向量列表
cont_clr_list <- FilPv1$clr
pot_clr_list <- FilPV1NONCont$clr

# 初始化判定结果向量
is_contaminant <- logical(length(pot_clr_list))

# 遍历每个潜在污染物，检查与已识别污染物的相关性
for (j in seq_along(pot_clr_list)) {
  pot_vec <- pot_clr_list[[j]]
  for (k in seq_along(cont_clr_list)) {
    cont_vec <- cont_clr_list[[k]]
    # 确保两个向量对应同批次样本（长度一致）
    if (length(pot_vec) != length(cont_vec)) next
    # 计算Spearman相关系数
    rho <- cor(pot_vec, cont_vec, method = "spearman")
    if (rho > 0.7) {
      is_contaminant[j] <- TRUE
      break  # 只要有一个符合条件就停止当前物种的检查
    }
  }
}

# 将判定结果添加到潜在污染物数据框中
FilPV1NONCont$is_contaminant <- is_contaminant

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