RNA-seq 转录组批次效应校正 combat代码

在合并分析不同来源或不同数据集的时候，需要做一个合并前的批次校正，在校正前后，可用PCA聚类查看批次校正的效果，以评估校正的可靠性。

这里我们用到的是sva包的combat函数

代码如下：

library(sva)

expFile <- read.csv("rna_CountMatrix.txt",sep="\t",row.names = 1) #导入表达矩阵

bd <- read.csv("batchdisease.txt",sep="\t",row.names = 1) #表达矩阵样本的疾病信息

batchInfo <- as.vector(expFile["batch",])#批次信息

example <- t(batchInfo)[,1]

dt <- data.frame(t(bd))

mod=model.matrix(~disease,dt)#校正模型参数

expFile <- expFile[-55277,]

library("FactoMineR")

library("factoextra")

pca.plot = function(dat,col){

  df.pca <- PCA(t(dat), graph = FALSE)

  fviz_pca_ind(df.pca,

              geom.ind = "point",

              col.ind = col ,

              addEllipses = TRUE,

              legend.title = "Groups"

  )

}

pca.plot(expFile,factor(batchInfo))

combat_expFile <- ComBat(dat =as.matrix(expFile), batch = example,mod=mod, par.prior = F) #核心步骤，耗时较长