【原】GEO数据库挖掘之多个芯片数据集的合并

健明 2022-04-18

展开全文

众所周知，我们【生信技能树】团队这些年一直在身体力行的推广马拉松学习理念，目前唯一可行的生物信息学入门策略！是60个小时的线上全网直播视频课程，连续4个星期，每个星期5天，每天的晚上8~11点的3个小时互动授课（周三、周日休息）

与十万人一起学生信，你值得拥有下面的学习班：

数据挖掘（GEO,TCGA,单细胞）2022年5月场，快速了解一些生物信息学应用图表
生信入门课-2022年5月场，你的生物信息学第一课

同期也安排了助教和实习生一起协助整理讲师的直播授课时间线和授课内容相关笔记，提炼内容梗概方便学员们查漏补缺，精准复现！

下面是（ GEO数据挖掘）直播配套笔记

举例：GSE83521和GSE89143数据合并

1.下载数据

rm(list = ls())
library(GEOquery)
library(stringr)
gse = "GSE83521"
eSet1 <- getGEO("GSE83521", 
                destdir = '.', 
                getGPL = F)
eSet2 <- getGEO("GSE89143", 
                destdir = '.', 
                getGPL = F)
#(1)提取表达矩阵exp
exp1 <- exprs(eSet1[[1]])
exp1[1:4,1:4]
exp2 <- exprs(eSet2[[1]])
exp2[1:4,1:4]
exp2 = log2(exp2+1)
table(rownames(exp1) %in% rownames(exp2))
length(intersect(rownames(exp1),rownames(exp2)))
exp1 <- exp1[intersect(rownames(exp1),rownames(exp2)),]
exp2 <- exp2[intersect(rownames(exp1),rownames(exp2)),]
boxplot(exp1)
boxplot(exp2)

#(2)提取临床信息
pd1 <- pData(eSet1[[1]])
pd2 <- pData(eSet2[[1]])
if(!identical(rownames(pd1),colnames(exp1))) exp1 = exp1[,match(rownames(pd1),colnames(exp1))]
if(!identical(rownames(pd2),colnames(exp2))) exp2 = exp2[,match(rownames(pd2),colnames(exp2))]

#(3)提取芯片平台编号
gpl <- eSet2[[1]]@annotation

#(4)合并表达矩阵
# exp2的第三个样本有些异常，可以去掉或者用normalizeBetweenArrays标准化，把它拉回正常水平。

exp2 = exp2[,-3]

exp = cbind(exp1,exp2)
boxplot(exp)
Group1 = ifelse(str_detect(pd1$title,"Tumour"),"Tumour","Normal")
Group2 = ifelse(str_detect(pd2$source_name_ch1,"Paracancerous"),"Normal","Tumour")[-3]

Group = c(Group1,Group2)
table(Group)
Group = factor(Group,levels = c("Normal","Tumour"))
save(gse,Group,exp,gpl,file = "exp.Rdata")

两个数据集样本的情况

合并后的数据

2.针对不同数据集数据的差异，需要处理批次效应

2.1 使用limma包里的removeBatchEffect()函数

rm(list = ls())
load("exp.Rdata")
#处理批次效应
library(limma)
#?removeBatchEffect()
batch <- c(rep("A",12),rep("B",5))
exp2 <- removeBatchEffect(exp, batch)
par(mfrow=c(1,2))  # 展示的图片为一行两列
boxplot(as.data.frame(exp),main="Original")
boxplot(as.data.frame(exp2),main="Batch corrected")

2.2 使用sva包中的combat() 函数

rm(list = ls())
load("exp.Rdata")
#处理批次效应(combat)
library(sva)
#?ComBat

batch <- c(rep("A",12),rep("B",5))
mod = model.matrix(~Group)
exp2 = ComBat(dat=exp, batch=batch, 
              mod=mod, par.prior=TRUE, ref.batch="A")
par(mfrow=c(1,2))
boxplot(as.data.frame(exp),main="Original")
boxplot(as.data.frame(exp2),main="Batch corrected")