GEO数据挖掘补充(三)——多数据联合分析

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)#注意要判断exp是否需要取log
table(rownames(exp1) %in% rownames(exp2))#x%in%y——x的每个元素都在y中存在吗(一对多)
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
 
if (F) {
  library(tinyarray)
find_anno(gpl)#是否找到GPL平台文件
get_gpl_txt(gpl)#得到GPL平台文件连接，下载GPL文件并放到目录文件下
#检查下载的GPL文件，需处理异常
a = read.delim("GPL19978.txt",
               check.names=F,
               skip = 8,
               comment.char="!",) #skip跳过前几行，去掉几行，点开文件看一下，comment.char是跳过以！开头的行
colnames(a)
ids = a[,1:2]
#有空行怎么去掉
ids[1,2]
k = ids$circRNA!=""   #非空的行的条件
ids1 = ids[k,]  #只保留非空的行
colnames(ids1)=c("probe_id","symbol")
 
###若不是一个平台，需进行另一个平台的注释ids2
ids2 = ids[k,]
colnames(ids2)=c("probe_id","symbol")
 
#探针矩阵转换为基因矩阵
exp1 = trans_array(exp1,ids1)   #exp1本来是探针矩阵
exp2 = trans_array(exp2,ids1)
}#需要将探针注释转变为基因注释时再使用
 
#(4)合并表达矩阵
# exp2的第三个样本有些异常，可以去掉或者用normalizeBetweenArrays标准化，把它拉回正常水平。
 
exp2 = exp2[,-3]#删除第3个样本
 
exp = cbind(exp1,exp2)
boxplot(exp,las=2)#帮助为graphics::boxplot
#las值设置坐标轴数值方向，可取0，1，2，3
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")
 

rm(list = ls())
load("exp.Rdata")
#处理批次效应：removeBatchEffect()函数
library(limma)
#?removeBatchEffect()
batch <- c(rep("A",12),rep("B",5))
    #由之前的Group1和Group2可知两组的构成;1组为12个样本，1组为5个样本，消除他们之间的批次效应
exp_rectify<- removeBatchEffect(exp, batch)
par(mfrow=c(1,2))
boxplot(as.data.frame(exp),main="Original",las=2)
boxplot(as.data.frame(exp_rectify),main="Batch corrected",las=2)#解决批次效应后的数据作图
 
#画PCA图
library(tinyarray)
batch = factor(batch)
draw_pca(exp_rectify,batch)  #批次效应校正后的PCA图
#draw_pca(exp,batch,addEllipses = F)  #不显示圈
draw_pca(exp,Group) #批次效应校正前的PCA图
#draw_pca(exp,Group,addEllipses = F)  #不显示圈

rm(list = ls())
load("exp.Rdata")
#处理批次效应(combat)
library(sva)
#?ComBat
 
batch <- c(rep("A",12),rep("B",5))
mod = model.matrix(~Group)
exp_rectify = 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_rectify),main="Batch corrected")
 
#画PCA图
library(tinyarrary)
batch = factor(batch)
draw_pca(exp_rectify,batch)  #批次效应矫正后PCA图
#draw_pca(exp,batch,addEllipases = F)  #不显示圈
draw_pca(exp,Group) #批次效应矫正前PCA图
#draw_pca(exp,Group,addEllipases = F)  #不显示圈