GEO数据挖掘实例

2019.08.22

今天学习了r语言关于GEO数据挖掘方面的知识点，做一些整理。

1.首先是相应数据的下载与检查

利用R包进行下载

rm(list = ls())#清空环境变量
options(stringsAsFactors = F)##字符不作为因子读入
#####数据下载
library(GEOquery)
f='GSE5282.Rdata'
####getGPL获得平台的注释信息，但下载速度会慢很多
####而且注释文件格式大多不如bioconductor包好用
if(!file.exists(f)){
  gset<-getGEO('GSE5282',destdir='.',
               AnnotGPL=F,
               getGPL=F)
  save(gset,file=f)
}
#数据提取
load('GSE5282.Rdata')
ex= exprs(gset[[1]])
pd = pData(gset[[1]])
############代码查看矩阵是否需要log（来自GEO2R）
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
###一般来说如果LogC返回值为True，那么需要再做一步log2处理，此处为T，故做了处理。
##且为了防止出现负数的情况，一般会括号里加1.
ex = log2(ex+1)
##########boxplot可视化数据检查
boxplot(ex,las=2,cex.axis=0.6,main='data check')
boxplot(ex,las=2,cex.axis=0.6,main='data check')
group_list=c(c(rep("EGF_4h",3)),c(rep("Control_4h",3)),c(rep("Control_12h",3)),c(rep("EGF_12h",3)))
####PCA看分组情况
library("FactoMineR")
library("factoextra") 
####a data frame with n rows (individuals) 
####and p columns (numeric variables)
df.pca <- PCA(t(ex), graph = FALSE)
fviz_pca_ind(df.pca,
             geom.ind = "point",
             col.ind = group_list, 
             addEllipses = TRUE, 
             legend.title = "Groups"
)
save(ex,pd,group_list,file='ex_pd.Rdata')

2.构建矩阵，数据处理（DEG）

rm(list = ls())  ## 一键清空~
options(stringsAsFactors = F)
load('ex_pd.Rdata')  ##加载步骤一保存的数据文件
###构建实验设计矩阵
ex=ex[,1:6]   ##这里取数据的前六列（两组）进行处理（limma包一次只能对两组数据进行处理分析）
group_list=group_list[1:6]   ##同上，group_list也取前六列的名称

design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(ex)
library(limma)
####构建比较矩阵,设置用来对比的组别
contrast.matrix<-makeContrasts(contrasts=c('EGF_4h-Control_4h'),levels = design)
###此处('EGF_4h-Control_4h')  中二者位置调换也会改变结果。。。注意！！！
######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入 ！！！！！！重要！！！！！！
fit <- lmFit(ex,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
tT <- topTable(fit2, adjust="fdr", number=nrow(fit2))
dim(tT)
####提取出我们需要的指标
tT <- subset(tT, select=c('adj.P.Val',"P.Value","logFC"))
save(tT,file='DEG.Rdata')

3.进行ID转换

​ 因为geo数据一般来源于基因芯片，故需要对其芯片探针号进行转换，变成相应的基因名称

rm(list = ls())  ## 一键清空~
options(stringsAsFactors = F)
library(BiocManager)
##因为用的例子是大鼠的且芯片和rat2302有关，故用对应的rat2302包去转换
BiocManager::install('rat2302.db')  
library('rat2302.db')
load('ex_pd.Rdata')
load('DEG.Rdata')
#####
tT$probe_id <- rownames(tT)
ls("package:rat2302.db")
#####PROBEID_SYMBOL
a<- toTable(rat2302SYMBOL)
tT <- merge(tT,a,by.x='probe_id',by.y='probe_id')
#####PROBEID_ENTREZID
b<- toTable(rat2302ENTREZID)
tT <- merge(tT,b,by.x='probe_id',by.y='probe_id')
############################存在这种现象，多个探针对应一个symbol####
table(!duplicated(tT$symbol))
ex<- ex[tT$probe_id,]
identical(rownames(ex),tT$probe_id)
####只是一种参考，能自圆其说即可
tT$median=apply(ex,1,median)
tT=tT[order(tT$symbol,tT$median,decreasing = T),]
dim(tT)
tT=tT[!duplicated(tT$symbol),]
save(tT,file='DEG_symbol.Rdata')

4.热图绘制

rm(list=ls())
load('ex_pd.Rdata')
load('DEG.Rdata')
ex <- ex[,1:6]
######pheatmap 这是取前1000个gene进行展示，可根据自己的需求进行调整
choose_gene<-rownames(tT)[1:1000]
#####用normalize后的数据进行展示
data_matrix <-ex[choose_gene,]
#####热图更详细的了解https://www.jianshu.com/p/0e5ce59fa79e
#####scale对数据处理，求得的结果为z-score，即数据与均值之间差几个标准差 
####t是对数据做转置处理，为了适应scale的要求
data_matrix=t(scale(t(data_matrix)))
#####查看scale处理后数据的范围
fivenum(data_matrix)
####目的是避免出现极大极小值影响可视化的效果
###2,-2
data_matrix[data_matrix>1.2]=1.2
data_matrix[data_matrix< -1.2] = -1.2
library(pheatmap)
pheatmap(data_matrix)
####调整下颜色，使正负值颜色的对比会更加鲜明
require(RColorBrewer)
bk = c(seq(-1.2,1.2, length=100))
annotation_col = data.frame(Group = rep(c("EGF_4h", "Control_4h"), c(3,3)))
rownames(annotation_col)<-colnames(ex)
####annotation_col和annotation_row的格式需为数据框
####breaks参数用于值匹配颜色
####看下，color和breaks的对应，进行理解
pheatmap(data_matrix,
         breaks=bk,
         show_rownames = F,
         annotation_col = annotation_col,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(100))
####可以调整的内容有是否聚类、是否分组、是否显示行和列的内容等

5.火山图绘制

rm(list = ls())  
options(stringsAsFactors = F)
load(file = "DEG_symbol.Rdata")
####主要用两个值，logFC和p.value
####可参考进行1的计算
logFC_cutoff <- with(tT,mean(abs( logFC)) + 2*sd(abs( logFC)) )
###1也可自行根据需求设置
####依据logFC和adj.P.val为火山图的颜色分组做准备
tT$change = ifelse(tT$adj.P.Val < 0.05 & abs(tT$logFC) > logFC_cutoff,
                   ifelse(tT$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
####记录上调、下调基因数目，作为title的内容
### round（x，y）取x的y位小数，如ruond(1.132442,2),则结果为1.13
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(tT[tT$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(tT[tT$change =='DOWN',]))

library(ggplot2)
g = ggplot(data=tT, 
           aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile ) + 
  theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))+
  annotate('text',x=tT$logFC[tT$logFC>5],
           y=-log10(tT$P.Value[tT$logFC>5]),
           label=tT$symbol[tT$logFC>5])
print(g)
ggsave(g,filename = 'volcano.png')
save(tT,file='DEG_change.Rdata')

6.KEGG和GO分析

rm(list = ls())  
options(stringsAsFactors = F)
load('DEG_change.Rdata')
tT <- na.omit(tT)
####富集前的基因处理
library(clusterProfiler)
gene_up= tT[tT$change == 'UP','gene_id'] 
gene_down=tT[tT$change == 'DOWN','gene_id'] 
gene_diff=c(gene_up,gene_down)
gene_all = tT$gene_id
###############################KEGG富集分析#####################################################
## 1.KEGG pathway analysis
#上调、下调、差异、所有基因
#clusterProfiler作kegg富集分析：
###我们在意的是表格里的pvalue，对应的通路或term是否显著
kk.up<- enrichKEGG(gene         = gene_up,
                   organism     = 'rno',
                   universe     = gene_all,
                   pvalueCutoff = 5)
kk.down <- enrichKEGG(gene         =  gene_down,
                      organism     = 'rno',
                      universe     = gene_all)
load('annotation.Rdata')
#可视化展示-Barplot
down_kegg <- kk.down@result
down_kegg <- down_kegg[down_kegg$pvalue<1,]
barplot(kk.down,drop=T,showCategory = 12,title = 'bar_down')
barplot(kk.up,drop=T,showCategory = 12,title='bar_up')
#可视化展示-Dotplot
dotplot(kk.down,showCategory = 12,title = 'dot_down')
dotplot(kk.up,showCategory = 12,title='dot_up')
###############################GO############################################################
#细胞组分

library(org.Rn.eg.db)
ego_CC <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "CC",
                   pAdjustMethod = "BH")
#生物过程
ego_BP <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   ont = "BP",
                   pAdjustMethod = "BH")
#分子功能：
ego_MF <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   ont = "MF",
                   pAdjustMethod = "BH")
#可视化展示-Barplot
barplot(ego_CC, showCategory=12,title="bar_CC")
barplot(ego_BP, showCategory=12,title="bar_BP")
barplot(ego_MF, showCategory=12,title="bar_MF")
#可视化展示-Dotplot
dotplot(ego_CC,showCategory = 12,title="dot_CC")
dotplot(ego_BP,showCategory = 12,title="dot_BP")
dotplot(ego_MF,showCategory = 12,title="dot_MF")
save(kk.up,kk.down,ego_BP,ego_CC,ego_MF,file='annotation.Rdata')