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) #核心步骤,耗时较长

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容