GEO_多组数据联合分析（去除批次效应）
原创
GEO_多组数据联合分析（去除批次效应）
1 批次效应介绍
有的时候我们需要用多组数据（来自不同GSM）联合进行分析，或者批次效应较为明显，应该进行去除批次效应的操作。
一般可以采取两种方式：
1 先取各自差异基因然后取交集
deg <- intersect(deg1,deg2)
2 先合并，后差异分析：
这种方式有两个注意点
原则上选择来自同一芯片平台的GSE (一般不遵循)
不要选择一个全是处理组，一个全是对照组
（这样的话分组信息和批次信息是重合的，在处理批次信息的话也会将分组之间的差别抹除掉）
批次效应处理前后区别如下图
需要先将两组数据先合并再分析，使用下面的函数去除批次效应
limma::removeBatchEffect()
#或者
sva::ComBat()下面介绍方式二
2 表达矩阵与探针注释
使用数据集GSE83521与GSE89143进行操作
library(tinyarray)
geo1 = geo_download("GSE83521")
geo2 = geo_download("GSE89143")
#(1)提取表达矩阵exp
exp1 <- geo1$exp
exp2 <- geo2$exp
exp2 = log2(exp2+1)
boxplot(exp1)
boxplot(exp2)#exp2的第三个样本有些异常，可以去掉或者用normalizeBetweenArrays标准化，把它拉回正常水平。
#把第三个样本从表达矩阵里去掉
exp2 = exp2[,-3] 
分别获取探针注释ids1,ids2，不同的平台，都要自行修改获得探针注释。
这里的两个数据集恰好来自同一个平台，也分别获取模拟下
ids1
geo1$gpl
#没有注释包(所有的非编码RNA芯片都没有注释R包)，需要读取GPL页面的表格
get_gpl_txt(geo1$gpl)
a = read.delim("GPL19978.txt",skip = 8,comment.char = "!")
colnames(a)
ids = a[,1:2]
head(ids)
#有空行怎么去掉
ids[1,2]
k = ids$circRNA!=""
ids1 = ids[k,]
colnames(ids1)=c("probe_id","symbol")
head(ids1)ids2
geo2$gpl
get_gpl_txt(geo2$gpl)
a = read.delim("GPL19978.txt",skip = 8,comment.char = "!")
colnames(a)
ids = a[,1:2]
head(ids)
ids[1,2]
k = ids$circRNA!=""
ids2 = ids[k,]
colnames(ids2)=c("probe_id","symbol")3 合并表达矩阵
合并两个或者多个表达矩阵
这里注意trans_array()函数的使用，trans_array()可以处理exp和ids,将探针矩阵转换为基因矩阵
exp1 = trans_array(exp1,ids1)
exp2 = trans_array(exp2,ids2)
处理完之后取交集，后合并
s = intersect(rownames(exp1),rownames(exp2));length(s)
exp1 <- exp1[s,]
exp2 <- exp2[s,]
#合并表达矩阵，cbind是按列拼接的函数
exp = cbind(exp1,exp2)
boxplot(exp) #合并后的表达矩阵可以看到较为明显的批次效应
4 分组信息
需要修改和检查分组信息，注意多个数据的相同分组要用相同的关键词，例如下面的Group1和Group2关键词都是"Tumour","Normal"，一模一样。
pd1 <- geo1$pd
#把第三个样本从临床信息表格去掉
pd2 <- geo2$pd[-3,]
library(stringr)
Group1 = ifelse(str_detect(pd1$title,"Tumour"),"Tumour","Normal")
#检查分组是否正确
data.frame(pd1$title,Group1)
Group2 = ifelse(str_detect(pd2$source_name_ch1,"Paracancerous"),"Normal","Tumour")
#检查分组是否正确
data.frame(pd2$source_name_ch1,Group2)
Group = c(Group1,Group2) #两个分组信息合并
table(Group)
Group = factor(Group,levels = c("Normal","Tumour"))
[1] Tumour Tumour Tumour Tumour Tumour Tumour Normal Normal Normal Normal Normal Normal Tumour Normal Normal Tumour
[17] Normal
Levels: Normal Tumour5 处理批次效应
下面两个方法任选其一
library(limma)
#?removeBatchEffect()
#批次信息,原来是第一个数据里的就是A，原来是第二个数据里的就是B。两个数据合并不用改动，如果有多个数据再往后加rep("C",ncol(exp3))...
batch <- c(rep("A",ncol(exp1)),rep("B",ncol(exp2)));batch
[1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B"
exp_re <- removeBatchEffect(exp, batch,group = Group)
#去除批次效应前后的箱线图
par(mfrow=c(1,2))
boxplot(as.data.frame(exp),main="Original")
boxplot(as.data.frame(exp_re),main="Batch corrected")
##去除批次效应前后的PCA图
#draw_pca()函数是简化画pca的包，注意group参数应该是因子
##去除批次效应前后的PCA图
draw_pca(exp,Group)+draw_pca(exp_re,Group)
batch <- factor(batch)
draw_pca(exp,batch)+draw_pca(exp_re,batch)
#对比分析下
(draw_pca(exp,Group)+draw_pca(exp_re,Group))/(draw_pca(exp,batch)+draw_pca(exp_re,batch))
或者采用第二种方法
library(sva)
#?ComBat
batch <- c(rep("A",ncol(exp1)),rep("B",ncol(exp2)))
mod = model.matrix(~Group)
exp_cb = 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(exp_cb),main="Batch corrected")
##去除批次效应前后的PCA图
draw_pca(exp,Group)+draw_pca(exp_cb,Group)PS:
从对比图可以看出，未做批次处理前（左上，左下）的PCA图批次之间的差别（左下）是远远大于分组（左上）的，这样分析出来的差异基因是批次间的。去除批次效应后，批次之间的差异（右下）消除，差异主要是组之间的（右上），这样才是我们所关注的点。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
