如何在R中对BLASTn结果按多指标阈值进行无重复基因聚类？

针对你遇到的循环聚类冗余、报错问题，推荐用图论中的连通分量分析来实现，这是处理此类“相似性关联聚类”的标准方法，能避免循环的缺陷，准确把所有满足条件的基因归为同一聚类。

步骤1：筛选符合条件的基因配对 #

先从合并后的SAMPLES数据框中，过滤出满足pident≥90、qcovs≥90、lrqs≥0.9的基因对，仅保留用于关联分析的qseqid和sseqid：

library(dplyr)

# 过滤符合相似性条件的基因对
filtered_pairs <- SAMPLES %>%
  filter(pident >= 90, qcovs >= 90, lrqs >= 0.9) %>%
  select(qseqid, sseqid)

步骤2：构建无向图并计算连通分量 #

使用igraph包处理图结构，把每个基因视为节点，满足条件的基因对视为节点间的无向边，连通的节点集合就是一个聚类：

install.packages("igraph") # 首次使用需安装
library(igraph)

# 将基因对转换为无向图
gene_graph <- graph_from_data_frame(filtered_pairs, directed = FALSE)

# 计算每个基因所属的连通分量（聚类）
cluster_result <- components(gene_graph)

# 转换为数据框，方便后续合并
cluster_df <- data.frame(
  gene_id = names(cluster_result$membership),
  cluster = paste0("Cluster_", cluster_result$membership)
)

步骤3：合并聚类结果到原始数据 #

将聚类ID映射回原始的SAMPLES数据框，同时处理不满足条件的行（可选择设为独立聚类或NA）：

# 合并聚类ID到原始数据
SAMPLES_with_cluster <- SAMPLES %>%
  # 先匹配qseqid的聚类
  left_join(cluster_df, by = c("qseqid" = "gene_id")) %>%
  # 对满足条件的行，确保sseqid的聚类与qseqid一致（可选验证）
  mutate(
    cluster = case_when(
      # 满足条件的行，直接使用已分配的聚类ID
      pident >= 90 & qcovs >= 90 & lrqs >= 0.9 ~ cluster,
      # 不满足条件的基因，单独分配唯一聚类（或改为NA_character_）
      TRUE ~ paste0("Cluster_", nrow(cluster_df) + row_number())
    )
  )

优势说明 #

• 天然解决传递性关联：若A与B匹配、B与C匹配，A/B/C会自动归为同一聚类，不会出现循环方法的冗余分配
• 无遗漏、无重复：每个基因只会被分配到一个聚类
• 代码简洁高效，处理大规模数据的性能远优于循环

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