GEO数据挖掘全流程分析

最新推荐文章于 2024-08-21 16:27:03 发布

Eric's blog

最新推荐文章于 2024-08-21 16:27:03 发布

阅读量2.6w

点赞数 38

分类专栏：数据库学习文章标签： GEO

本文链接：https://blog.csdn.net/Eric_blog/article/details/104580174

版权

数据库学习专栏收录该内容

6 篇文章 26 订阅

订阅专栏

声明：以下学习资料根据“生信技能树”网络系列免费教学材料整理而成，代码来自“生信技能树”校长jimmy的github。GEO数据库挖掘系列知识分享课程，于2016年首发于生信菜鸟团博客。配套教学视频在B站，特此声明。

前言：关于GEO数据

我们的目标是要从读懂文献到复刻文献实验，再到掌握GEO数据挖掘的能力。首先便是要广泛阅读，在读文献时，提炼脉络，读懂文献使用了哪个或哪些GSE数据集，对数据做了哪些处理。了解清楚后，便可下载相应的数据集，得到表达矩阵，作差异分析，注释等一系列下游分析。
一篇文章可以有一个或多个GSE数据集，一个GSE里可以有一个或多个GSM样本。多个研究的GSM样本可以根据研究目的整合为一个GDS，每个数据集有着自己对应的芯片平台（GPL），一个GSE里可能有多个平台测出的数据。
本分析是基于R语言的平台，所以需要一些R语言的基础知识。
了解GEO，表达芯片与R

第一部分：GEO芯片数据下载及整理

GEO官网
本例以GSE42872数据集为例，学习GEO数据挖掘分析，分析文献
在这里插入图片描述
通过阅读文献找到相应的GSE数据集，并在官网可以下相应的数据集信息及背景知识
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872

这个数据集使用的是GPL6244这个芯片平台，由6个样本组成，前三个为对照组，后三个为处理组。了解了这个数据集的相关背景后，接下来就需要进行数据下载了。数据下载的方式有很多种，这里我们用R包GEOquery来下载，具体下载代码如下：

rm(list = ls())  
## 魔幻操作，一键清空~当前环境中对象全部删除
options(stringsAsFactors = F)
#在调用as.data.frame的时，将stringsAsFactors设置为FALSE可以避免character类型自动转化为factor类型
f='GSE42872_eSet.Rdata'
#把GSE42872_eSet.Rdata赋值给f，方便后面流程化处理
##根据数据集不同修改相应的GSE号

library(GEOquery)
# 这个包需要注意两个配置，一般来说自动化的配置是足够的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)

if(!file.exists(f)){
  gset <- getGEO('GSE42872', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
##这是一个函数，利用包将数据集的表达信息下载下来，赋值给了gset，而不下载注释信息和平台信息，病保存到本地，文件名为f。

load('GSE42872_eSet.Rdata') 
 ## 载入数据
class(gset) 
 #查看数据类型
length(gset) 
 ##看一下有几个元素
gset[[1]]
#取第一个元素
class(gset[[1]])
 #查看改元素的数据类型
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list

a=gset[[1]] 
##取出第一个元素赋值给一个对象a
dat=exprs(a) 
#a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数，该函数得到表达矩阵
#现在 得到的dat就是一个表达矩阵，只不过基因的ID是探针名
dim(dat)
#看一下dat这个矩阵的维度
dat[1:5,1:5] 
#查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
#这个表达矩阵是已经log之后的，表达量一般是0-10左右，如果是原始芯片表达的信号值一般是几千到一万，则需要log处理。

boxplot(dat,las=2) 
#画个图看一下各样本之间有没有批次效应，一般中位数都差不多，las是将横坐标样本信息竖着排列

pd=pData(a) 
#通过查看说明书知道取对象a里的临床信息用pData
View(pd)
## 查看一下，挑选一些感兴趣的临床表型，这里我们欲得到其分组title信息。
library(stringr)
#运行一个字符分割包
group_list=str_split(pd$title,' ',simplify = T)[,4]
#抽取title一列，按照空格分割，取第四个元素即Control和Vemurafenib
table(group_list)
#看一下两个分组各有几个

也可以使用中国镜像下载：

library(devtools)
install_github("jmzeng1314/AnnoProbe")
library(AnnoProbe)
gset<-geoChina("GSE12345")

GPL6244这个是该数据的芯片平台，首先在网页查一下有没有对应的包，有的话直接下载相应的包，没有用下面的函数直接下载。
##网址：http://www.bio-info-trainee.com/1399.html
或者在博客查找对应关系
##可以查到对应的包为 hugene10sttranscriptcluster
##因此我们优先使用包来寻找探针和基因名的对应关系，在bioconductor下载这个包。
如果下载了相应的芯片平台的包，运行以下程序：

library(hugene10sttranscriptcluster.db)
##运行这个包
ls("package:hugene10sttranscriptcluster.db")
#对这个包进行探索，看一下有多少元素
##找到有SYMBOL的那个元素就是我们需要的对应关系
##例如：[34] "hugene10sttranscriptclusterSYMBOL"    
ids=toTable(hugene10sttranscriptclusterSYMBOL) 
#toTable这个函数：通过看hgu133plus2.db这个包的说明书知道提取probe_id（探针名）和symbol（基因名）的对应关系的表达矩阵的函数为toTable。并赋值给ids
##这时候我们得到19825个探针对应的基因名。
###刚才我们的表达矩阵中是33297个基因探针，这就意味着刚才的表达矩阵中可能存在多个探针重复对应一个基因名。这就需要我们对数据进行进一步筛选、处理。
head(ids) 
#head为查看前六行

如果没有查到芯片平台对应的R包就需要下载平台注释信息，按以下操作

if(F){
  library(GEOquery)
  #Download GPL file, put it in the current directory, and load it:
  gpl <- getGEO('GPL6244', destdir=".")
##需要修改相应的平台号，把平台信息赋值给gpl
  colnames(Table(gpl))  
  head(Table(gpl)[,c(1,15)]) 
## you need to check this , which column do you need
  probe2gene=Table(gpl)[,c(1,15)]
  head(probe2gene)
  library(stringr)  
  save(probe2gene,file='probe2gene.Rdata')
}
load(file='probe2gene.Rdata')
ids=probe2gene

现在我们就通过两种方法得到了GPL探针和基因的对应关系。

head(ids)
#为查看前六行
colnames(ids)=c('probe_id','symbol')  
#将列名统一改为'probe_id','symbol'方便后续统一操作。

length(unique(ids$symbol)) 
#[1] 18832个独特的基因探针，意味着本来19825个里面有一部分是重复的
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
#每个对象出现的个数
plot(table(sort(table(ids$symbol))))
#画图观察

ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in%  rownames(dat),]
##%in%用于判断是否匹配，然后取匹配的几行，去掉无法匹配的信息。

dat[1:5,1:5]   
dat=dat[ids$probe_id,] 
#取表达矩阵中可以与探针名匹配的那些，去掉无法匹配的表达数据，这时只剩下19825个探针及表达信息，其余已被剔除。

ids$median=apply(dat,1,median) 
#ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#对ids$symbol按照ids$median中位数从大到小排列的顺序排序
##即先按symbol排序，相同的symbol再按照中位数从大到小排列，方便后续保留第一个值。
##将对应的行赋值为一个新的ids，这样order()就相当于sort()
ids=ids[!duplicated(ids$symbol),]
#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果.最后得到18832个基因。
dat=dat[ids$probe_id,] 
#新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol
#把ids的symbol这一列中的每一行给dat作为dat的行名
head(dat)
##至此我们就得到了该数据集的表达矩阵，最后将结果保存。
save(dat,group_list,file = 'step1-output.Rdata')
write.csv(dat,file="expressionmetrix_GSE.csv")

第二部分：数据分布判断及分析

检查样本分布

exprSet=dat
exprSet['GAPDH',]
boxplot(exprSet[,1])
boxplot(exprSet['GAPDH',])
exprSet['ACTB',]
library(reshape2)
?melt()
class(exprSet)
exprSet_L=melt(exprSet,)
colnames(exprSet_L)=c('symbol','sample','value')
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)

关于melt函数的用法及介绍可以参考博客
下面进行绘图

### ggplot2 
library(ggplot2)
p=ggplot(exprSet_L,
         aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
pdf('result1.pdf')#width = 5,height = 10)
p
dev.off()

p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
pdf('result2.pdf')#width = 5,height = 10)
p
dev.off()
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
pdf('result3.pdf')#width = 5,height = 10)
p
dev.off()
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
pdf('result4.pdf')#width = 5,height = 10)
p
dev.off()
p=ggplot(exprSet_L,aes(value,col=group))+geom_density() 
pdf('result5.pdf')#width = 5,height = 10)
p
dev.off()
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())

pdf('result6.pdf')#width = 5,height = 10)
p
dev.off()

PCA图

hc=hclust(dist(t(exprSet)))
plot(hc)
colnames(exprSet)=paste(group_list,1:78,sep = " ")
hc=hclust(dist(t(exprSet)))
p<-plot(hc)
pdf('result7.pdf')#width = 5,height = 10)
p
dev.off()
library(ggfortify)
df=as.data.frame(t(exprSet))
df$group<-group_list
p<-autoplot(prcomp(df[,1:(ncol(df)-1)]),
         data=df,
         colour="group")
pdf('result8.pdf')#width = 5,height = 10)
p
dev.off()

第三部分：利用limma包差异表达分析

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = 'step1-output.Rdata')
# 每次都要检测数据
dat[1:4,1:4] 
table(group_list) #table函数，查看group_list中的分组个数
#通过为每个数据集绘制箱形图，比较数据集中的数据分布
boxplot(dat[1,]~group_list) #按照group_list分组画箱线图

bp=function(g){         #定义一个函数g，函数为{}里的内容
  library(ggpubr)
  df=data.frame(gene=g,stage=group_list)
  p <- ggboxplot(df, x = "stage", y = "gene",
                 color = "stage", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}
bp(dat[1,]) ## 调用上面定义好的函数，避免同样的绘图代码重复多次敲。
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])
dim(dat)

用limma包做GEO芯片合并后的差异表达分析：

#source("http://bioconductor.org/biocLite.R")
#biocLite("limma")

logFoldChange=1
adjustP=0.05

library(limma)
setwd("C:\\Users\\lexb4\\Desktop\\geoBatch\\07.diff")
rt=read.csv("namalized.csv")
rt=as.matrix(rt)
###解决基因名重复的问题
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
##删除基因名重复

#differential定义样品类型
modType=c(rep("normal",25),rep("tumor",25),rep("normal",5),rep("tumor",5))
#定义肿瘤类型
design <- model.matrix(~0+factor(modType))
colnames(design) <- c("normal","tumor")

fit <- lmFit(rt,design)
allDiff[1:5,1:5]
cont.matrix<-makeContrasts(tumor-normal,levels=design)
#处理组减去对照组；
#去了log之后的值做减法就是除法，得到的值就是logFC
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

allDiff=topTable(fit2,adjust='fdr',number=200000)
#校正方法FDR，定义200000就保证了所有的基因的差异情况
write.csv(alldiff,file="alldiff.csv")

#write table
diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ]
write.table(diffSig,file="diff.xls",sep="\t",quote=F,row.names=F)
diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ]
write.table(diffUp,file="up.xls",sep="\t",quote=F,row.names=F)
diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ]
write.table(diffDown,file="down.xls",sep="\t",quote=F,row.names=F)
#提取差异基因表达矩阵
hmExp=rt[as.vector(diffSig[,1]),]
diffExp=rbind(id=colnames(hmExp),hmExp)
write.table(diffExp,file="diffExp.txt",sep="\t",quote=F,col.names=F)

##画火山图
pdf(file="vol.pdf")
xMax=max(-log10(allDiff$adj.P.Val))
yMax=max(abs(allDiff$logFC))
##最大值可以手动输入
plot(-log10(allDiff$adj.P.Val), allDiff$logFC, xlab="-log10(adj.P.Val)",ylab="logFC",
     main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.8)
diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange)
points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.8)
diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange))
points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="green",cex=0.8)
##cex为点的大小，可以自己设置
abline(h=0,lty=2,lwd=3)
dev.off()

接着来画个热图（并且定义样本类型）

#install.packages("pheatmap")

setwd("C:\\....pheatmap")      #设置工作目录
rt=read.table("diffExp.txt",sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
##检查是否有重复的基因名
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)

#rt=log2(rt+1)
#rt[rt>15]=15

library(pheatmap)
##注释文件
Geo=c(rep("GSE33335",50),rep("GSE56807",10))
##前面10个来源于GSE33335,后面10个来源于GSE56807
Type=c(rep("N",25),rep("T",25),rep("N",5),rep("T",5))    #修改正常和癌症样品数目

names(Geo)=colnames(rt)
ann=cbind(Geo,Type)
ann=as.data.frame(ann)

tiff(file="heatmap.tiff",
       width = 20,            #图片的宽度
       height =25,            #图片的高度
       units ="cm",
       compression="lzw",
       bg="white",
       res=500)  ##图片分辨率
pheatmap(rt, annotation=ann, 
         color = colorRampPalette(c("green", "black", "red"))(50),  ##地，中，高表达颜色
         cluster_cols =F,
#         scale="row",   校正，视图形而定
         fontsize_row=3,
         fontsize_col=5)
dev.off()

第四部分功能富集分析

在做富集分析之前，要先把基因的symble转换成基因的entrezID
首先准备基因symble文件，一列为symble，一列为logFC

#install.packages("RSQLite")
#source("https://bioconductor.org/biocLite.R")
#biocLite("org.Hs.eg.db")

setwd("C:\\Users\\lexb4\\Desktop\\geoBatch\\09.symbol2id")

library("org.Hs.eg.db")
rt=read.table("symbol.txt",sep="\t",check.names=F,header=T)
genes=as.vector(rt[,1])
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
write.table(out,file="id.txt",sep="\t",quote=F,row.names=F)

运行后就多了一列为基因的entrezID。接下来就可以进行基因富集分析，结果会以柱状图和气泡图的类型呈现。

#install.packages("colorspace")
#install.packages("stringi")
#source("http://bioconductor.org/biocLite.R")
#biocLite("DOSE")
#biocLite("clusterProfiler")
#biocLite("pathview")

setwd("C:\\...GO")
library("clusterProfiler")
library("org.Hs.eg.db")
rt=read.table("id.txt",sep="\t",header=T,check.names=F)
rt=rt[is.na(rt[,"entrezID"])==F,]
##保留ID不是NA的基因列表
geneFC=rt$logFC
gene=rt$entrezID
names(geneFC)=gene

#GO富集分析
kk <- enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =0.05, qvalueCutoff = 0.05)
write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F)

#柱状图
tiff(file="barplot.tiff",width = 20,height = 30,units ="cm",compression="lzw",bg="white",res=300)
barplot(kk, drop = TRUE, showCategory = 47)
dev.off()

#点图
tiff(file="dotplot.tiff",width = 20,height = 30,units ="cm",compression="lzw",bg="white",res=300)
dotplot(kk,showCategory = 47)
dev.off()

KEGG富集分析（柱状图，气泡图，通路图）
通路中有颜色的就是差异基因
输入文件：ID.txt

#install.packages("colorspace")
#install.packages("stringi")
#source("http://bioconductor.org/biocLite.R")
#biocLite("DOSE")
#biocLite("clusterProfiler")
#biocLite("pathview")

setwd("C:\\.....KEGG")
library("clusterProfiler")
rt=read.table("id.txt",sep="\t",header=T,check.names=F)
rt=rt[is.na(rt[,"entrezID"])==F,]

geneFC=rt$logFC
gene=rt$entrezID
names(geneFC)=gene

#kegg富集分析
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05)
##然后看一下有多少个通路有意义
##太少了可以把p值或者q值适当调整
write.table(kk,file="KEGG.txt",sep="\t",quote=F,row.names = F)

#柱状图
tiff(file="barplot.tiff",width = 20,height = 20,units ="cm",compression="lzw",bg="white",res=300)
barplot(kk, drop = TRUE, showCategory = 12)
dev.off()

#点图
tiff(file="dotplot.tiff",width = 20,height = 20,units ="cm",compression="lzw",bg="white",res=300)
dotplot(kk)
dev.off()

#通路图
library("pathview")
keggxls=read.table("KEGG.txt",sep="\t",header=T)
for(i in keggxls$ID){
  pv.out <- pathview(gene.data = geneFC, pathway.id = i, species = "hsa", out.suffix = "pathview")
}