GEO数据库是最常用的生信挖掘数据库之一，其中储存了很多芯片、二代测序和其他高通量测序数据，包括Affymetrix、Agilent和Illumina等多种平台。

今天介绍GPL570芯片平台的GEO原始数据的处理及合并，包括GSE39582 GSE37892 GSE14333三个数据集。这次的三个数据集来源于2021年发表在 Theranostics 上的 m6A分型文章，PMID：33500720，处理方法参照文章方法学部分。

原始cel文件的下载在官网找到对应数据集后点击Download下的http即可

GPL570.png

首先介绍GSE39582的表达矩阵和临床资料的获取，GSE39582是数据量最大、应用最广泛的结肠癌数据集。GSE39582临床资料完善，其中的生存数据包括OS和RFS。

##加载包
library(affy)
library(GEOquery)
affy_data=ReadAffy(celfile.path = 'E:\\GEO原始文件\\GSE39582_RAW')
#RMA算法。
eset.rma=rma(affy_data)
eset=exprs(eset.rma)
colnames(eset)=unlist(lapply(colnames(eset),function(x){
  strsplit(x,'[_]')[[1]][1]
}))
#加载注释包
library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
head(ids)
length(unique(ids$symbol))
table(sort(table(ids$symbol)))
library(tidyverse)
eset=as.data.frame(eset)
eset<- eset %>% mutate(probe_id=rownames(eset)) %>% inner_join(ids, by = 'probe_id') %>% 
  dplyr::select(probe_id,symbol,everything())
#多个基因，保留表达量最高的
eset2<- eset %>% group_by(symbol)%>%
  summarise_all(max) %>%
  dplyr::select(-probe_id) %>%
  column_to_rownames('symbol')

#临床信息的处理
gset=getGEO('GSE39582',destdir = '.',AnnotGPL = F, getGPL = F)
pdata=pData(gset[[1]])
phe=pdata[,12:25]
library(stringr)
colnames(phe)=c('sex','age','stage','T','N','M','location','chemo','regimen','RFS_event','RFS_time','OS_event','OS_time','MMR')
phe$sex=str_split(phe$sex,' ',simplify = T)[,2]
phe$age=str_split(phe$age,' ',simplify = T)[,3]
phe$stage=str_split(phe$stage,' ',simplify = T)[,2]
phe$T=str_split(phe$T,' ',simplify = T)[,2]
phe$N=str_split(phe$N,' ',simplify = T)[,2]
phe$M=str_split(phe$M,' ',simplify = T)[,2]
#这里我根据phe的各列字符分布情况把phe分为两部分tmp1和tmp2,然后cbind
tmp1=phe[,1:6]
tmp2=apply(phe[,7:ncol(phe)],2,function(x){
  str_split(x,' ',simplify = T)[,2]
})
rownames(tmp2)=rownames(phe)
tmp=cbind(tmp1,tmp2)
rownames(tmp)=rownames(phe)
#去掉正常组织的表达
group_list=ifelse(str_detect(pdata$source_name_ch1,'Adenocarcinoma'),'tumor','normal')
group_list=factor(group_list,levels = c('normal','tumor'))
tmp_tumor=tmp[group_list=='tumor',]
eset3=eset2[,rownames(tmp_tumor)]

eset_39582=eset3
phe_39582=tmp_tumor

save(eset_39582,phe_39582,file='GSE39582_RAW.Rdata')

随后是GSE37892的处理，它的生存数据只有RFS而无OS。

library(affy)
library(GEOquery)
affy_data=ReadAffy(celfile.path = 'E:\\GEO原始文件\\GSE37892_RAW')
eset.rma=rma(affy_data)
eset=exprs(eset.rma)
colnames(eset)=unlist(lapply(colnames(eset),function(x){
  strsplit(x,'[_]')[[1]][1]
}))

library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
head(ids)
length(unique(ids$symbol))
table(sort(table(ids$symbol)))
library(tidyverse)
eset=as.data.frame(eset)
eset<- eset %>% mutate(probe_id=rownames(eset)) %>% inner_join(ids, by = 'probe_id') %>% 
  dplyr::select(probe_id,symbol,everything())

eset2<- eset %>% group_by(symbol)%>%
  summarise_all(max) %>%
  dplyr::select(-probe_id) %>%
  column_to_rownames('symbol')
gset7=getGEO('GSE37892',destdir = '.',AnnotGPL = F, getGPL = F)
##这里命名为pdata7是因为自己整理过一些数据集，排序为7
pdata7=pData(gset7[[1]])
phe7=pdata7[,14:16]
colnames(phe7)=c('sugery.date','metastasis.date','last.contact.date')
library(stringr)
phe7$sugery.date= str_split(phe7$sugery.date,' ',simplify = T)[,4]
phe7$metastasis.date= str_split(phe7$metastasis.date,' ',simplify = T)[,5]
phe7$last.contact.date= str_split(phe7$last.contact.date,' ',simplify = T)[,5]
phe7$fustat=ifelse(phe7$metastasis.date == 'NA',0,1)
phe7$metastasis.date=ifelse(phe7$metastasis.date == 'NA',phe7$last.contact.date,phe7$metastasis.date)
#这里用到时间处理函数as.Data
phe7$futime=as.numeric(as.Date(phe7$metastasis.date)-as.Date(phe7$sugery.date))
phe7=phe7[,4:5]
phe7$futime=round(phe7$futime/30,2)

ss=intersect(rownames(phe7),colnames(eset2))
eset2=eset2[,ss]
phe7=phe7[ss,]
eset_37892=eset2
phe_37892=phe7
save(eset_37892,phe_37892,file='GSE37892.Rdata')

GSE14333的生存资料是DFS。

library(affy)
library(GEOquery)
affy_data=ReadAffy(celfile.path = 'E:\\GEO原始文件\\GSE14333_RAW')
eset.rma=rma(affy_data)
eset=exprs(eset.rma)
colnames(eset)=unlist(lapply(colnames(eset),function(x){
  strsplit(x,'[.]')[[1]][1]
}))

library(hgu133plus2.db)
ids=toTable(hgu133plus2SYMBOL)
head(ids)
length(unique(ids$symbol))
table(sort(table(ids$symbol)))
library(tidyverse)
eset=as.data.frame(eset)
eset<- eset %>% mutate(probe_id=rownames(eset)) %>% inner_join(ids, by = 'probe_id') %>% 
  dplyr::select(probe_id,symbol,everything())

eset2<- eset %>% group_by(symbol)%>%
  summarise_all(max) %>%
  dplyr::select(-probe_id) %>%
  column_to_rownames('symbol')
library(GEOquery)
gset=getGEO('GSE14333',destdir = '.',AnnotGPL = F, getGPL = F)
pdata=pData(gset[[1]])
phe=pdata[,c(2,38)]
library(stringr)
phe$fustat=str_split(phe$`Location:ch1`,';',simplify = T)[,6]
phe$futime=str_split(phe$`Location:ch1`,';',simplify = T)[,5]
phe=phe[,c(3,4)]
phe$fustat=str_split(phe$fustat,':',simplify = T)[,2]
phe$futime=str_split(phe$futime,':',simplify = T)[,2]
table(phe$fustat)
phe=as.data.frame(phe)

phe$futime=as.numeric(phe$futime)
phe$fustat=as.numeric(phe$fustat)
phe=na.omit(phe)

ss=intersect(rownames(phe),colnames(eset2))
eset2=eset2[,ss]
phe=phe[ss,]
eset_14333=eset2
phe_14333=phe
save(eset_14333,phe_14333,file='GSE14333.Rdata')

接下来是多个GSE数据集的合并，使用sva包的Combat函数去除批次效应

samegene=intersect(rownames(eset_39582),rownames(eset_37892))

eset_39582=eset_39582[samegene,]
eset_37892=eset_37892[samegene,]
eset_14333=eset_14333[samegene,]

colnames(eset_39582)=paste0('GSE39582_',colnames(eset_39582))
colnames(eset_37892)=paste0('GSE37892_',colnames(eset_37892))
colnames(eset_14333)=paste0('GSE14333_',colnames(eset_14333))

meta=cbind(eset_39582,eset_37892,eset_14333)

library(sva)##加载包
exp_sva=as.matrix(meta)
dimnames=list(rownames(exp_sva),colnames(exp_sva))
data=matrix(as.numeric(as.matrix(exp_sva)),nrow = nrow(exp_sva),dimnames = dimnames)
dim(data)
batchType=c(rep('GSE39582',566),rep('GSE37892',130),rep('GSE14333',226))
out=data.frame(ComBat(data,batchType),par.prior=T)
outTab=out[,-ncol(out)]
outTab=rbind(geneNames=colnames(outTab), outTab)
write.table(outTab, file="exp.merge.txt", sep="\t", quote=F, col.names=F)

后续会尝试跟着文章的分析做下去，并更新推文，希望小伙伴们可以继续关注，欢迎大家关注‘文献与代码学习’，同步更新。谢谢大家，敬请批评指正。