如何用R或Bash按簇筛选仅在单一簇中存在的基因？

R 实现 #

前提准备 #

你需要两个文件：

• gene_matrix.csv：基因×样本的二进制矩阵，第一列为基因名，其余列为样本，值为0（缺失）或1（存在）
• cluster_info.csv：样本簇划分，两列分别为sample（样本名）和cluster（簇ID，1-8）

代码实现 #

library(tidyverse)

# 读取并整理基因矩阵：转成长格式方便后续分组计算
gene_mat <- read_csv("gene_matrix.csv", show_col_types = FALSE) %>%
  pivot_longer(-1, names_to = "sample", values_to = "presence")

# 读取簇划分信息
cluster_df <- read_csv("cluster_info.csv", show_col_types = FALSE)

# 合并数据，计算每个基因在各簇的存在情况
gene_cluster_stats <- inner_join(gene_mat, cluster_df, by = "sample") %>%
  group_by(gene, cluster) %>%
  # 标记该簇是否有至少一个样本携带该基因
  summarise(cluster_has_gene = max(presence), .groups = "drop") %>%
  group_by(gene) %>%
  summarise(
    # 统计该基因存在的簇数量
    total_clusters = sum(cluster_has_gene),
    # 提取唯一存在的簇ID（仅当簇数量为1时）
    exclusive_cluster = ifelse(total_clusters == 1, cluster[cluster_has_gene == 1], NA)
  ) %>%
  # 筛选仅在单个簇存在的基因
  filter(total_clusters == 1)

# 导出结果到CSV
write_csv(gene_cluster_stats, "exclusive_genes.csv")

运行后，exclusive_genes.csv会包含所有符合条件的基因及其对应的专属簇ID。

Bash 实现 #

前提准备 #

• gene_matrix.csv：第一行是样本名，第一列是基因名，其余为0/1值
• cluster_info.txt：每行格式为样本名 簇ID（空格分隔）

代码实现 #

# 第一步：将簇信息转换为样本-簇的CSV映射文件
awk '{print $1","$2}' cluster_info.txt > sample_cluster_map.csv

# 第二步：用awk处理基因矩阵，筛选目标基因
awk -F',' '
# 处理第一行，建立样本名到簇ID的映射
NR==1 {
    for (i=2; i<=NF; i++) {
        sample=$i
        # 从映射文件中查找当前样本的簇ID
        while ((getline line < "sample_cluster_map.csv") > 0) {
            split(line, arr, ",")
            if (arr[1] == sample) {
                sample_cluster[sample] = arr[2]
                break
            }
        }
        close("sample_cluster_map.csv")
        colnames[i] = sample
    }
    next
}
# 处理每一行基因数据
{
    gene=$1
    # 初始化所有簇的存在标记为0
    for (c=1; c<=8; c++) cluster_has[c] = 0
    # 遍历每个样本的数值，更新对应簇的标记
    for (i=2; i<=NF; i++) {
        if ($i == 1) {
            cluster=sample_cluster[colnames[i]]
            cluster_has[cluster] = 1
        }
    }
    # 统计存在该基因的簇数量，并记录唯一簇ID
    count=0
    exclusive_c=""
    for (c=1; c<=8; c++) {
        if (cluster_has[c] == 1) {
            count++
            exclusive_c=c
        }
    }
    # 仅输出符合条件的基因
    if (count == 1) print gene "," exclusive_c
}
' gene_matrix.csv > exclusive_genes_bash.csv

执行后，exclusive_genes_bash.csv即为筛选结果。

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