数据挖掘 | Count数据去除批次效应后不是整数甚至还出现负值导致无法进行差异分析怎么办？

之前咱们介绍过数据挖掘 | 批次效应的鉴定与处理 | 附完整代码 + 注释 | 看完不会来揍我，但是很多小伙伴遇到了Count数据批次处理后不是整数甚至还出现负值的问题，这就导致无法使用某些包包进行差异分析（对差异分析感兴趣的小伙伴可以查看：看完还不会来揍/找我 | 差异分析三巨头 —— DESeq2、edgeR 和 limma 包 | 附完整代码 + 注释），如下所示：

那该怎么办呢？ComBat-seq帮你搞定！走！一起看看去！

---

• 如果小伙伴们有需求的话，可以加入我们的交流群：一定要知道 | 永久免费的环境友好型生信学习交流群又双叒叕来啦！| 伴随不定期群友好物分享！在这里，你可以稍有克制地畅所欲言！

• 超级建议大家在入群前或入群后可以看一下这个：干货满满 | 给生信小白的入门小建议 | 掏心掏肺版！绝对干货满满！让你不虚此看！

• 如果有需要个性化定制分析服务的小伙伴，可以看看这里：你要的个性化生信分析服务今天正式开启啦！定制你的专属解决方案！全程1v1答疑！！绝对包你满意！

---

ComBat-seq介绍

其实咱们之前有浅浅提过一嘴，Count数据可以用ComBat_seq：

但是没有细讲！咱们今天展开说说！

ComBat-seq是一种专门针对Count数据的批次效应处理工具，它是基于ComBat的改进版。

ComBat-seq发表于NAR Genomics and Bioinformatics，有兴趣的小伙伴们可以查看原文：ComBat-seq: batch effect adjustment for RNA-seq count data | NAR Genomics and Bioinformatics | Oxford Academic

原文链接：https://doi.org/10.1093/nargab/lqaa078

ComBat-seq接受未经转换的原始计数矩阵作为输入，与ComBat相同，它也需要批次信息（也就是你有几批数据，也就是不同的数据集，也就是不同的实验批次等等）。

ComBat-seq使用负二项回归模型来去除批次效应，通过将原始数据映射到预期分布来提供调整后的数据。与ComBat假设的高斯分布相比，这种方法更能捕捉RNA-Seq计数数据的特性。ComBat-seq去除批次效应后得到的数据保留了计数数据的整数特性，使其可以与常见差异分析工具（例如edgeR、DESeq2等，这些工具要求的输入数据为原始计数数据）的数据需求相兼容。

方法原理

Combat-seq的原理如下图所示：

有兴趣的小伙伴可以去看看原文！这里我就不展开说了哈哈哈哈哈哈！

代码实操

包的安装与加载

# 安装并加载包包
# devtools::install_github("zhangyuqing/sva-devel")
library(sva)
1
2
3

如果GitHub上的包安装出现问题，噢不，如果包的安装过程中出现任何问题，都可以查看：看完不会来揍我 | R包的下载与安装 | 再也没有一个包可以逃出你的手掌心啦

基本用法

ComBat_seq要求至少输入两个参数 —— 来自RNA-Seq的原始计数矩阵（不进行任何归一化或转换）以及批次信息。

# 基本用法

# 生成一个50行8列的负二项分布矩阵
set.seed(0425)
count_matrix <- matrix(rnbinom(400, size = 10, prob = 0.1), nrow = 50, ncol = 8)
head(count_matrix)
#      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
# [1,]   93   73   50   54  100  107   90   81
# [2,]  104   85  142  112  146   99   89   71
# [3,]   79  121   83  110  126   66   60   95
# [4,]   61   71   99   85   85   53   84   55
# [5,]  117  143   55  101   76  115  103  129
# [6,]  106   71   47   84   71   85  107   75

# 都是整数嗷！

# 创建一个批次(batch)向量，一共8个样本，其中前4个样本属于批次1，后4个样本属于批次2
batch <- c(rep(1, 4), rep(2, 4))
head(batch)
# [1] 1 1 1 1 2 2

# 使用ComBat_seq函数进行批次效应的调整
adjusted <- ComBat_seq(count_matrix, batch = batch, group = NULL)
# count_matrix: 待调整的数据矩阵
# batch: 批次信息向量，指定每个样本所属的批次
# group: 不指定分组信息，因此默认使用批次作为调整的依据
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

指定单个原有生物学分组

样本中的有些分组，本身就是有差别的，不可以矫枉过正，所以我们在去除批次的同时，也需要保留它们原本生物学上的差异哟！数据挖掘 | 批次效应的鉴定与处理 | 附完整代码 + 注释 | 看完不会来揍我中有具体介绍过，ComBat函数的mod选项就是用来指定原有分组的，咱们这里的ComBat_seq函数用的是group选项。我们指定分组后，原有的生物学差异就会在调整后的数据中保留下来啦！

# 指定单个原有生物学分组

# 创建一个长度为8的分组(group)向量
group <- rep(c(0,1), 4)
head(group)
# [1] 0 1 0 1 0 1

# 使用ComBat_seq函数进行批次效应的调整，并加入分组信息作为调整的依据
adjusted_counts <- ComBat_seq(count_matrix, batch = batch, group = group)
# count_matrix: 待调整的数据矩阵
# batch: 批次信息向量，指定每个样本所属的批次
# group: 分组信息向量，用于更精细地调整数据
1
2
3
4
5
6
7
8
9
10
11
12

指定多个原有生物学分组

如果你希望指定多个生物学分组变量，可以将它们作为矩阵或数据框传递给参数covar_mod。

# 指定多个原有生物学分组

# 创建两个长度为4的协变量向量
cov1 <- rep(c(0,1), 4)
cov2 <- c(0,0,1,1,0,0,1,1)

# 将这两个协变量合并成一个协变量矩阵
covar_mat <- cbind(cov1, cov2)
head(covar_mat)
# cov1 cov2
# [1,]    0    0
# [2,]    1    0
# [3,]    0    1
# [4,]    1    1
# [5,]    0    0
# [6,]    1    0

# 使用ComBat_seq函数进行批次效应的调整，并加入协变量信息作为调整的依据
adjusted_counts <- ComBat_seq(count_matrix, batch = batch, group = NULL, covar_mod = covar_mat)
# count_matrix: 待调整的数据矩阵
# batch: 批次信息向量，指定每个样本所属的批次
# group: 不指定分组信息，因为你已经把想要指定的分组信息放在covar_mat中啦！所以这里就不需要了！
# covar_mod: 协变量矩阵，用于更精细地调整数据
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

文末碎碎念

那今天的分享就到这里啦！我们下期再见哟！

最后顺便给自己推荐一下嘿嘿嘿！

如果我的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！

蟹蟹你们的喜欢和支持！！！

参考资料

1. https://github.com/zhangyuqing/ComBat-seq
2. https://doi.org/10.1093/nargab/lqaa078
3. https://www.jianshu.com/p/ed8a1b1ab7b7