有的时候我们需要用多组数据(来自不同GSM)联合进行分析,或者批次效应较为明显,应该进行去除批次效应的操作。
一般可以采取两种方式:
1 先取各自差异基因然后取交集
deg <- intersect(deg1,deg2)
2 先合并,后差异分析:
这种方式有两个注意点
原则上选择来自同一芯片平台的GSE (一般不遵循)
不要选择一个全是处理组,一个全是对照组
(这样的话分组信息和批次信息是重合的,在处理批次信息的话也会将分组之间的差别抹除掉)
批次效应处理前后区别如下图
需要先将两组数据先合并再分析,使用下面的函数去除批次效应
limma::removeBatchEffect()
#或者
sva::ComBat()
下面介绍方式二
使用数据集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")
合并两个或者多个表达矩阵
这里注意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) #合并后的表达矩阵
可以看到较为明显的批次效应
需要修改和检查分组信息,注意多个数据的相同分组要用相同的关键词,例如下面的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 Tumour
下面两个方法任选其一
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 删除。