R语言分析1：RNA-seq的批次效应校正

批次效应（batch effect），表示样品在不同批次中处理和测量产生的与试验期间记录的任何生物变异无关的技术差异。其既可能来自实验，也可能是来自分析流程。实验中样品收集、建库、测序的不同批次可能带来系统性的偏差；分析中不同工具的使用也有一定偏差。

注意:批次校正只能降低批次效应的影响，而不能完全消除批次效应，在假定实验设计没有问题时，可以通过探究数据结构的方法去评估批次效应，衡量的方法如: 1. distance measures(距离法)；2. clustering(样品层次聚类法)； 3.spatial methods(空间统计法)。

 ：加载相关数据

# 安装并加载所需的R包
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("GEOquery")
library(GEOquery)


# 用GEOquery包来下载 ExpressionSet 对象
eSet_GSE13861 <- getGEO("GSE13861", destdir = '.')
eSet_GSE26899 <- getGEO("GSE26899", destdir = '.')

# 提取表达矩阵exp
exp_GSE13861 <- exprs(eSet_GSE13861[[1]])
exp_GSE26899 <- exprs(eSet_GSE26899[[1]])

# 提取临床信息
pd_GSE13861 <- pData(eSet_GSE13861[[1]])
pd_GSE26899 <- pData(eSet_GSE26899[[1]])

# 数据合并并记录分组
exp = cbind(exp_GSE13861, exp_GSE26899)
condition <- data.frame(sample = colnames(exp),
                        batch = c(rep('GSE13861', 90), rep('GSE26899', 108)))
Modtype <- c(rep('N-GC', 19), rep('GC', 7), rep('GIST', 2), rep('GC', 14), rep('GIST', 1), rep('GC', 27), 
             rep('GIST', 1), rep('GC', 1), rep('GIST', 1), rep('GC', 12), rep('GC', 1), rep('GC', 4), 
             rep('GC', 1), rep('N-GC', 1), rep('GC', 1), rep('GIST', 1), rep('GC', 1), rep('N-GC', 3), 
             rep('GC', 1), rep('N-GC', 2), rep('GIST', 1), rep('GC', 1), rep('N-GC', 2), rep('GC', 2), 
             rep('N-GC', 3), rep('GC', 8), rep('N-GC', 1), rep('GC', 3), rep('GIST', 1), rep('GC', 75))
mod <- model.matrix(~as.factor(Modtype))

# 使用PCA查看批次效应
library(FactoMineR)
library(factoextra)

pre.pca <- PCA(t(exp), graph = FALSE)
fviz_pca_ind(pre.pca,
             geom = c("point", "text"),
             col.ind = condition$batch,
             addEllipses = TRUE)

# 也通过使用Hierarchical clustering的聚类方法
dist_mat <- dist(t(exp))
clustering <- hclust(dist_mat, method = "complete")
par(mfrow = c(2, 1))
plot(clustering, labels = condition$batch)
plot(clustering, labels = Modtype)

batch-1

★ 整个数据不能被圈到一起，说明批次效应影响较大

1. 使用校正

ComBat 使用参数或非参数经验贝叶斯模型，输入数据为干净的、标准化的表达数据，通常是芯片数据
ComBat_seq 使用负二项回归的ComBat改进模型，专门针对RNA-Seq count数据

# BiocManager::install("sva")
library(sva)

combat_count <- ComBat(as.matrix(exp),
                       batch = condition$batch,
                       mod = mod # 添加生物分组信息
)

combat.pca <- PCA(t(combat_count), graph = FALSE)
fviz_pca_ind(combat.pca, col.ind = condition$batch,
             geom = c("point", "text"),
             addEllipses = TRUE,
             legend.title="Group")

batch-2

2. 使用校正

library(limma)
library(ggplot2)
limma_count <- removeBatchEffect(exp,
                                batch = condition$batch,
                                design = mod)

# 对负值取0（pmax()函数），对小数点四舍五入取整（round()函数）
# limma_count_zero_int <-round(pmax(limma_count, 0)) 

limma.pca <- PCA(t(limma_count), graph = FALSE)
fviz_pca_ind(limma.pca, col.ind = condition$batch,
             geom = c("point", "text"),
             addEllipses = TRUE,
             legend.title="Group")

batch-3

★ 结果与combat差不多