基于基因组距离与p值修剪R数据框中的SNP对冗余

问题说明 #

现有一个包含SNP对信息的R数据框，字段包括：

• chr_snp1：第一个SNP的染色体编号
• loc_snp1：第一个SNP的碱基对位置
• chr_snp2：第二个SNP的染色体编号
• loc_snp2：第二个SNP的碱基对位置
• pval：该SNP对的p值

数据示例：

# 构造示例数据框
df <- data.frame(
  chr_snp1 = c(1, 1, 1, 1, 2, 2),
  loc_snp1 = c(278474050, 275620856, 144075771, 144075771, 103701229, 102990566),
  chr_snp2 = c(2, 2, 3, 3, 7, 7),
  loc_snp2 = c(57386477, 57386662, 109909704, 111344453, 56369738, 56407691),
  pval = c(7.43e-08, 1.08e-07, 1.02e-06, 2.06e-06, 3.83e-08, 1.07e-07)
)

需要移除冗余SNP对，冗余定义为：

• 染色体组合（chr_snp1与chr_snp2的组合）完全相同
• 两行的loc_snp1差值≤10,000,000
• 两行的loc_snp2差值≤10,000,000

每组冗余对中仅保留p值最低的行，预期结果如下：

# 预期结果
expected_df <- data.frame(
  chr_snp1 = c(1, 1, 2),
  loc_snp1 = c(278474050, 144075771, 103701229),
  chr_snp2 = c(2, 3, 7),
  loc_snp2 = c(57386477, 109909704, 56369738),
  pval = c(7.43e-08, 1.02e-06, 3.83e-08)
)

实现方案 #

我们可以利用聚类算法对同染色体组合内的SNP对进行分组，将位置相近的SNP对归为一类，再从每个类中筛选出p值最小的行。这种方法避免了逐行对比，效率更高。

步骤1：加载所需包 #

使用dplyr进行数据处理，dbscan进行空间聚类：

library(dplyr)
library(dbscan)

步骤2：按染色体组合分组并聚类 #

对每个染色体组合的子集，以loc_snp1和loc_snp2作为二维坐标，用dbscan聚类（距离阈值设为1e7），然后为每个聚类组保留p值最小的行：

cleaned_df <- df %>%
  # 按染色体组合分组
  group_by(chr_snp1, chr_snp2) %>%
  # 对每组进行聚类处理
  mutate(
    # 构造二维坐标矩阵，用于dbscan
    coords = list(cbind(loc_snp1, loc_snp2)),
    # 执行dbscan，eps=1e7，minPts=1（允许单个样本成簇）
    cluster = dbscan(unlist(coords, recursive = FALSE), eps = 1e7, minPts = 1)$cluster
  ) %>%
  ungroup() %>%
  # 按染色体组合+聚类组分组，保留pval最小的行
  group_by(chr_snp1, chr_snp2, cluster) %>%
  filter(pval == min(pval)) %>%
  ungroup() %>%
  # 移除临时字段
  select(-coords, -cluster)

步骤3：验证结果 #

查看处理后的数据框，与预期结果一致：

print(cleaned_df)

方案说明 #

• dbscan聚类：通过将SNP对的两个位置视为二维空间点，距离阈值设为10,000,000，自动将满足位置相近条件的SNP对归为同一簇。
• 分组筛选：在每个染色体组合的每个聚类簇中，筛选出p值最小的行，确保只保留每组冗余对中最显著的结果。
• 效率优势：聚类算法的时间复杂度远低于逐行对比，适合处理大规模数据集。

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