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来解决我们的问题了。
生物信息学

空白云服务器的一些生物信息学设置

2020-8-29 23:53:41

生物信息学

用matplotlib实现画中画

2020-9-1 16:02:42

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索