在R语言中基于两个数据框的唯一匹配填充0/1矩阵

问题背景 #

现有两个R数据框：

1. go.d5g：存储GO注释信息，包含GO ID、关联基因和功能描述，示例数据：

# go.d5g示例
go.d5g <- data.frame(
  ID = c("GO:0001922", "GO:0001922", "GO:0001922", "GO:0001922", "GO:0002901", "GO:0001777"),
  Gene = c("ABL1", "HIF1A", "TNFAIP3", "SH2B2", "ADA", "BAX"),
  Term = c("B-1 B cell homeostasis", "B-1 B cell homeostasis", "B-1 B cell homeostasis", 
           "B-1 B cell homeostasis", "mature B cell apoptotic process", "T cell homeostatic proliferation")
)

2. deg：存储差异表达基因信息，包含log2倍数变化、基因名、表达变化类型和实验对比组，示例数据：

# deg示例
deg <- data.frame(
  L2FC = c(-2.754236, 3.161623, -2.821350, -1.798022, -1.293536, -1.011016),
  Gene = c("SLC13A2", "SNAI2", "STYK1", "CD84", "TLE6", "P2RX1"),
  diffexp = c("Downregulated", "Upregulated", "Downregulated", "Downregulated", "Downregulated", "Downregulated"),
  comp = rep("NS.CB.A,S.ED.A", 6)
)

需求目标 #

生成一个0/1关联矩阵：

• 行：deg中的唯一基因
• 列：go.d5g中的唯一GO ID
• 矩阵值：基因与GO ID关联则为1，否则为0

示例输出格式：

GO:0001922 GO:0002901 GO:0001777 GO:0006924 GO:0033153 GO:0002204
SLC13A2          1          1          0          0          0          0
SNAI2            0          0          0          0          0          0
STYK1            0          1          1          0          1          0
CD84             0          0          0          0          0          0
TLE6             0          1          1          0          0          0
P2RX1            0          0          0          0          0          1

当前低效实现（循环法） #

原代码使用循环遍历基因，数据量大时效率极低：

g.u <- unique(deg$Gene)
goid.u <- unique(go.d5g$ID)
cmat <- matrix(0,nrow=length(g.u),ncol=length(goid.u))
rownames(cmat) <- g.u
colnames(cmat) <- goid.u

for (i in 1:length(g.u)) {
  go.match <- unlist(lapply(g.u[i], function(x) which(go.d5g$Gene %in% x)))
  go.match2 <- go.d5g$ID[go.match]
  cmat[i,which(goid.u %in% go.match2)] <- 1
}

高效替代方案 #

以下三种方法均采用向量化操作，避免循环，大幅提升处理效率：

方法1：tidyverse 流水线实现 #

代码简洁易读，适合习惯tidy语法的用户：

library(tidyverse)

# 1. 提取deg的唯一基因列表
deg_genes <- unique(deg$Gene)
# 2. 去重GO注释（避免同一基因-GO组合重复）
go_unique <- distinct(go.d5g, Gene, ID)
# 3. 构建全基因列表与GO注释的关联，转宽格式并填充0
go_matrix <- tibble(Gene = deg_genes) %>%
  left_join(go_unique, by = "Gene") %>%
  mutate(value = 1) %>%
  pivot_wider(names_from = ID, values_from = value, values_fill = 0) %>%
  column_to_rownames("Gene") %>%
  as.matrix()

方法2：Base R 原生实现 #

无需额外安装包，适合轻量场景：

# 1. 限定因子水平，确保覆盖所有目标基因和GO ID
go_fact <- transform(go.d5g,
                     Gene = factor(Gene, levels = unique(deg$Gene)),
                     ID = factor(ID, levels = unique(go.d5g$ID)))
# 2. 生成频数表并转换为0/1矩阵
go_matrix <- as.matrix(table(go_fact$Gene, go_fact$ID))
go_matrix[go_matrix > 1] <- 1  # 处理同一基因-GO的重复条目
# 3. 补充deg中无GO注释的基因
missing_genes <- setdiff(unique(deg$Gene), rownames(go_matrix))
if (length(missing_genes) > 0) {
  missing_mat <- matrix(0, nrow = length(missing_genes), ncol = ncol(go_matrix),
                        dimnames = list(missing_genes, colnames(go_matrix)))
  go_matrix <- rbind(go_matrix, missing_mat)
}

方法3：稀疏矩阵实现（超大数据集） #

当基因和GO ID数量极多时，使用稀疏矩阵大幅节省内存：

library(Matrix)

# 1. 获取基因和GO ID的索引位置
gene_idx <- match(go.d5g$Gene, unique(deg$Gene))
go_idx <- match(go.d5g$ID, unique(go.d5g$ID))
# 2. 创建稀疏矩阵
sparse_mat <- sparseMatrix(i = gene_idx, j = go_idx, x = 1,
                           dimnames = list(unique(deg$Gene), unique(go.d5g$ID)))
# 3. 按需转换为普通矩阵
go_matrix <- as.matrix(sparse_mat)

结果验证 #

可通过以下代码验证新方法与原循环法结果一致：

# 对比矩阵元素和
sum(go_matrix) == sum(cmat)  # 应返回TRUE

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