之前咱们介绍过数据挖掘 | 批次效应的鉴定与处理 | 附完整代码 + 注释 | 看完不会来揍我，但是很多小伙伴遇到了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)
如果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: 不指定分组信息，因此默认使用批次作为调整的依据指定单个原有生物学分组
样本中的有些分组，本身就是有差别的，不可以矫枉过正，所以我们在去除批次的同时，也需要保留它们原本生物学上的差异哟！数据挖掘 | 批次效应的鉴定与处理 | 附完整代码 + 注释 | 看完不会来揍我中有具体介绍过，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: 分组信息向量，用于更精细地调整数据
指定多个原有生物学分组
如果你希望指定多个生物学分组变量，可以将它们作为矩阵或数据框传递给参数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: 协变量矩阵，用于更精细地调整数据
文末碎碎念
那今天的分享就到这里啦！我们下期再见哟！
最后顺便给自己推荐一下嘿嘿嘿！
如果我的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！
蟹蟹你们的喜欢和支持！！！
参考资料
https://github.com/zhangyuqing/ComBat-seq https://doi.org/10.1093/nargab/lqaa078 https://www.jianshu.com/p/ed8a1b1ab7b7
