RNA-seq数据DESeq对象按处理组重复样本过滤低计数基因

需求说明 #

需要过滤掉没有任何一个处理组的两个重复样本计数均≥阈值N的基因，仅保留至少存在一个处理组的两个重复样本计数都满足≥N的基因。以示例数据为例，当N=1时，geneE因所有处理组的双重复本都无法同时满足≥1，会被剔除。

解决方案 #

核心思路是先按处理组拆分样本，再对每个基因逐一检查是否存在符合条件的处理组，具体实现如下：

1. 处理示例数据

# 示例计数矩阵
X1A1 <- c(117, 24, 45, 146, 1)
X1A2 <- c(129, 31, 58, 159, 0)
X1B1 <- c(136, 25, 50, 1293, 0)
X1B2 <- c(131, 24, 50, 1073, 4)
X1C1 <- c(113, 23, 43, 132, 0)
X1C2 <- c(117, 18, 43, 126, 0)
X1D1 <- c(101, 20, 0, 875, 1)
X1D2 <- c(99, 21, 38 , 844, 0)
X24A1 <- c(109, 17, 60, 95, 0)
X24A2 <- c(122, 14, 611, 90, 0)

df <- data.frame(X1A1, X1A2, X1B1, X1B2, X1C1, X1C2, X1D1, X1D2, X24A1, X24A2)
rownames(df) <- c("geneA", "geneB", "geneC", "geneD", "geneE")

# 设置阈值N
N <- 1

# 提取每个样本对应的处理组（去除末尾的重复编号1/2）
groups <- substr(colnames(df), 1, nchar(colnames(df)) - 1)

# 筛选符合条件的基因
keep <- apply(df, 1, function(gene_counts) {
  # 按处理组拆分当前基因的计数
  split_counts <- split(gene_counts, groups)
  # 检查是否存在至少一个处理组的双重复本均≥N
  any(sapply(split_counts, function(group_counts) all(group_counts >= N)))
})

# 过滤后的数据集
df_filtered <- df[keep, ]
df_filtered

运行后会得到仅包含geneA-geneD的数据集，符合预期。

2. 适配DESeq对象

如果你的数据存储在DESeq2的DESeq对象中，只需先提取计数矩阵，再用同样逻辑处理：

library(DESeq2)

# 假设你的DESeq对象名为dds
count_matrix <- counts(dds)

# 提取处理组（根据样本命名规则调整，此处同示例逻辑）
groups_dds <- substr(colnames(count_matrix), 1, nchar(colnames(count_matrix)) - 1)

# 筛选符合条件的基因
keep_dds <- apply(count_matrix, 1, function(gene_counts) {
  split_counts <- split(gene_counts, groups_dds)
  any(sapply(split_counts, function(group_counts) all(group_counts >= N)))
})

# 过滤DESeq对象
dds_filtered <- dds[keep_dds, ]

代码解释 #

• 分组逻辑：通过substr去除样本名末尾的重复编号（1/2），将同处理组的两个样本归为一组；
• 基因筛选：使用apply逐行遍历基因，对每个基因的计数按处理组拆分后，检查是否存在至少一个组的所有样本计数均≥N；
• DESeq适配：从DESeq对象中提取原始计数矩阵后，筛选逻辑与普通数据框完全一致，最后将筛选结果应用回DESeq对象即可。

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