GEO2R分析结果是否可靠呢？

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是《数据挖掘第五期》学员的分享

听完了生信技能树的数据挖掘课程，想起来了去年导师安排我的一个文献阅读学习笔记，当时并不会使用R代码，就使用了官网的GEO2R工具进行分析，而且文献本身也是那样做的的。正好现在R已经踏入正轨，就把文章的数据分析再复现一遍，R代码版本

文献

  • 题目：Identification of potential core genes in triple negative breast cancer using bioinformatics analysis
  • 数据集：GSE38959，GSE45827，GSE65194
  • 因现如今比较认可的求差异方法是R语言LIMMA包，但此文献直接用官网的GEO2R工具进行分析，现来探究两者结果是否有差异呢？？

R语言LIMMA包对表达芯片做差异分析

这里直接使用jimmy老师的表达芯片数据分析标准代码，我就不过多介绍代码细节了，基本上换一个GEO数据集，需要更改的地方很少很少。

第一步：数据下载

rm(list = ls()) ## 魔幻操作，一键清空~options(stringsAsFactors = F)#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型# 注意查看下载文件的大小，检查数据 f='GSE38959_eSet.Rdata'# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64634library(GEOquery)# 这个包需要注意两个配置，一般来说自动化的配置是足够的。#Setting options('download.file.method.GEOquery'='auto')#Setting options('GEOquery.inmemory.gpl'=FALSE)if(!file.exists(f)){ gset <- getGEO('GSE38959', destdir=".", AnnotGPL = F, ## 注释文件 getGPL = F) ## 平台文件 save(gset,file=f) ## 保存到本地}load('GSE38959_eSet.Rdata') ## 载入数据class(gset) #查看数据类型length(gset) #class(gset[[1]])gset# assayData: 45015 features, 47 samples
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的lista=gset[[1]] #dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数dim(dat)#看一下dat这个矩阵的维度 # GPL4133dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列boxplot(dat,las=2)dat=dat[apply(dat,1,sd)>0,]dim(dat) # 45015 47dat[dat<0]=1boxplot(dat,las=2)
range(dat) # 1.53e+00 4.04e+05 数字太大，需要logdat=log2(dat+1)boxplot(dat,las=2)library(limma)#标准化dat=normalizeBetweenArrays(dat)boxplot(dat,las=2)

第二步：分组&注释

pd=pData(a) #47个样本#通过查看说明书知道取对象a里的临床信息用pData#查看pd，删除后4行table(pd$source_name_ch1)pd=pd[-(44:47),]dat=dat[,-(44:47)]table(pd$source_name_ch1)##分组group_list=ifelse(grepl('normal',pd$title),'normal','TNBC')table(group_list)#normal TNBC 13 30 #设置参考水平，对照在前，处理在后group_list = factor(group_list, levels = c('normal','TNBC'))
# GPL4133 注释ann=read.delim(file='GPL4133.txt',comment.char = '#',stringsAsFactors = F)colnames(ann)dat=as.data.frame(dat)dat$ID=rownames(dat)id_symbol_expr<-na.omit(merge(x=dat,y=ann[c('ID','GENE_SYMBOL')],by='ID',all.x=T))#这里处理一个探针对应不同gene名的情况，随便取一个就行symbol<-lapply(id_symbol_expr$GENE_SYMBOL,FUN = function(x){strsplit(x,'///')[[1]][1]})id_symbol_expr$GENE_SYMBOL<-as.character(symbol)#去除掉NA值，就是说有些探针对应不到已知的基因上，至少在GPL文件中没有对应关系ids=id_symbol_expr[id_symbol_expr$GENE_SYMBOL != 'NA',]ids=ids[ids$ID %in% rownames(dat),]dat=dat[ids$ID,] #处理一个gene名对应多个探针情况，中位数排序取最大#去掉dat最后一列的IDb=dat[,-44]ids$median=apply(b,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行ids=ids[order(ids$GENE_SYMBOL,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的idsids=ids[!duplicated(ids$GENE_SYMBOL),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果sdat=dat[ids$ID,] #新的ids取出ID这一列，将dat按照取出的这一列中的每一行组成一个新的datrownames(dat)=ids$GENE_SYMBOL#把ids的symbol这一列中的每一行给dat作为dat的行名dat[1:4,1:4] #保留每个基因ID第一次出现的信息dat=dat[,-44]dim(dat) #19749 43save(dat,group_list,pd,file = 'GSE38959-output.Rdata')

第三步：差异分析

rm(list = ls()) options(stringsAsFactors = F)load(file = 'GSE38959-output.Rdata')# 每次都要检测数据dat[1:4,1:4] table(group_list) #table函数，查看group_list中的分组个数#非匹配分析library(limma)design=model.matrix(~group_list)fit=lmFit(dat,design)fit=eBayes(fit)deg=topTable(fit,coef=2,number = Inf)## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。logFC_t=2deg$g=ifelse(deg$P.Value>0.05,'stable', ifelse( deg$logFC > logFC_t,'UP', ifelse( deg$logFC < -logFC_t,'DOWN','stable') ))table(deg$g)#DOWN stable UP #308 19102 339save(deg,file = 'GSE38959_deg.Rdata')

就是下载数据，然后根据临床信息分组， 然后做差异分析，选择阈值，确定上下调的基因数量。

GEO2R工具

步骤

这里网页工具的使用，基本上就是鼠标点击，更无须讲解了，关键步骤截图如下：

第一步：分组

 

第二步：点击 “Top 250” 进行差异分析

 

第三步：读入R语言进行筛选

b=read.table(file='GSE38959.txt',header = TRUE)dim(b) #45015 8#去除掉NA值，就是说有些探针对应不到已知的基因上，至少在GPL文件中没有对应关系b=b[order(b$Gene.symbol,decreasing = T),]b[b==""] = NAb=na.omit(b)dim(b) #32995 8#这里处理一个探针对应不同gene名的情况，随便取一个就行symbol<-lapply(b$Gene.symbol,FUN = function(x){strsplit(x,'///')[[1]][1]})b$Gene.symbol<-as.character(symbol)## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。colnames(b)logFC_t=2b$g=ifelse(b$P.Value>0.05,'stable', ifelse(b$logFC > logFC_t,'UP', ifelse( b$logFC < -logFC_t,'DOWN','stable') ))table(b$g)#分组，筛选重复基因up=b[grep('UP',b$g),]down=b[grep('DOWN',b$g),]up=up[!duplicated(up$Gene.symbol),]down=down[!duplicated(down$Gene.symbol),]dim(up) #515 9dim(down) #336 9save(b,up,down,file = 'GEO2R_GSE38959.Rdata')
  • 求出上调515个，下调336个，很大接近文中结果
  •  

结果对比

  • R语言结果上调基因339个
 
  • GEO2R工具上调基因515个
 
  • 可以发现GEO2R工具的结果包含了R语言里面的，并且还多出了很多，可以推出工具用的注释平台不一样。（PS；其实这个时候应该是需要一个韦恩图展现）
  • 查资料发现，GPL文件有两个版本：annot版和soft版。
  • annot版是我们平常在GEO搜索GPL下载的，一般情况下用它进行注释。soft版是用户所提供，也是GEO2R所用到，里面注释的基因也比annot版多很多，文件也挺大，一般网络很难下载。
  • 不过，感觉大众会比较接受annot版的GPL。因为不确定用户版的soft是否有错误，并且里面的symbol如今是否仍在使用。

总结

  1. 通过查阅文献，发现很多高分的文献求差异都是用R语言的LIMMA包，很少直接使用工具分析的。所以在同等条件基础上，更推荐使用前者求差异。并且现在代码也被提炼出来非常套路化，很容易上手。
  2. 其实，GEO2R也是非常受限制，只能做一些基础的分析。例如多组分析，配对或其它批次效应的情况下，就不能用GEO2R来解决我们的问题了。
赞 (0)

相关推荐