GEO數(shù)據(jù)挖掘基本流程與代碼

寫在前面：本文內容出自生信技能樹的生信入門系列課程筆記，感謝小潔老師、Jimmy老師的分享。

GEO數(shù)據(jù)挖掘分析思路：
1.找數(shù)據(jù)，找到GSE編號
2.下載數(shù)據(jù)（表達矩陣、臨床信息、分組信息）
3.數(shù)據(jù)探索（分組之間是否有差異，PCA、整個數(shù)據(jù)的熱圖）
4.limma差異分析及可視化（P值、logFC，火山圖、差異基因的熱圖）
5.富集分析KEGG、GO
注意：該標準流程只適用于表達芯片分析，甲基化、SNP等芯片的流程詳見生信技能樹專題。

GEO表達芯片分析的要點:
解決probe_id與gene symbol、樣本編號GSM與分組之間的對應關系。

GO富集的3個方面：
1.分子功能（Molecular Function，MF ）
2.細胞組分（Cellular Component ，CC）
3.生物過程（Biological Process ，BP）

GO富集的意義：
通過生物屬性的從屬關系，達到比通路富集更深層次的挖掘。比如差異基因都富集到了某些通路上，而這些通路的相關性并不強烈，但是GO富集通過它們從屬關系，可以發(fā)現(xiàn)它們都指向或從屬于某個更深層次的通路。

不同GSE合并要處理批次效應。

以下是步驟及對應代碼，每一個大標題代表實戰(zhàn)中的一個腳本文件，長腳本管理也可參照該方式：

00_pre_install.R

options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

cran_packages <- c('tidyr',
                   'tibble',
                   'dplyr',
                   'stringr',
                   'ggplot2',
                   'ggpubr',
                   'factoextra',
                   'FactoMineR',
                   'devtools',
                   'cowplot',
                   'patchwork',
                   'basetheme',
                   'ggthemes',
                   'paletteer',
                   'AnnoProbe',
                   'tinyarray') 
Biocductor_packages <- c('GEOquery',
                         'hgu133plus2.db',
                         'ggnewscale',
                         "limma",
                         "impute",
                         "GSEABase",
                         "GSVA",
                         "clusterProfiler",
                         "org.Hs.eg.db",
                         "preprocessCore",
                         "enrichplot",
                         "ggplotify")

for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}


for (pkg in c(Biocductor_packages,cran_packages)){
  require(pkg,character.only=T) 
#character.only=T指的是require函數(shù)加載的是pkg所代表的字符串表示的包，而不是pkg這個包

01_start_GEO.R

rm(list = ls())
library(GEOquery)
gse_number = "GSE56649"
eSet <- getGEO(gse_number, 
               destdir = '.', 
               getGPL = F)
class(eSet)
length(eSet)
eSet = eSet[[1]] 
#1.提取表達矩陣exp
exp <- exprs(eSet)
exp[1:4,1:4]
exp = log2(exp+1) #看數(shù)據(jù)是否已經(jīng)取過log值——數(shù)值是否跨數(shù)量級
boxplot(exp) #看一下數(shù)據(jù)是否基本在一條直線上
#若數(shù)據(jù)差異較大，則需標準化一下：exp2=limma::normalizeBetweenArrays(exp)
#(2)提取臨床信息
pd <- pData(eSet)
#(3)調整pd的行名順序與exp列名完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平臺編號
gpl_number <- eSet@annotation
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

02_group_ids.R

# Group(實驗分組)和ids(探針注釋)
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 1.Group----實驗分組要去閱讀臨床信息的表格來獲取，每個GSE的都不一樣
# 第一類，有現(xiàn)成的可以用來分組的列
Group = pd$`disease state:ch1` 

#第二類，自己生成
Group = c(rep("RA",times=13),
          rep("control",times=9))
Group = rep(c("RA","control"),times = c(13,9))

#第三類，匹配關鍵詞，自行分類
Group=ifelse(str_detect(pd$source_name_ch1,"control"),
             "control",
             "RA")

#設置參考水平，指定levels，對照組在前，處理組在后
Group = factor(Group,
               levels = c("control","RA"))
Group

#2.探針注釋的獲取-----------------
#捷徑
find_anno(gpl_number)
ids <- AnnoProbe::idmap('GPL570')
#方法1 BioconductorR包(最常用)
gpl_number 
#http://www.bio-info-trainee.com/1399.html
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
#用toTable把內置數(shù)據(jù)轉換為數(shù)據(jù)框
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
# 方法2 讀取GPL平臺的soft文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：soft文件列名不統(tǒng)一，活學活用，有的表格里沒有symbol列，也有的GPL平臺沒有提供注釋
  a = getGEO(gpl_number,destdir = ".")
  b = a@dataTable@table
  colnames(b)
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol")
  ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"http:///"),]
}

# 方法3 官網(wǎng)下載注釋文件并讀取
##http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus
# 方法4 自主注釋 
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,gse_number,file = "step2output.Rdata")

03_pca_heatmap.R

rm(list = ls())  
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")
#輸入數(shù)據(jù)：exp和Group
#Principal Component Analysis
#http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials

# 1.PCA 圖----
dat=as.data.frame(t(exp))
library(FactoMineR)
library(factoextra) 
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = Group, # color by groups
                         palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot
ggsave(plot = pca_plot,filename = paste0(gse_number,"_PCA.png"))
save(pca_plot,file = "pca_plot.Rdata")

# 2.top 1000 sd 熱圖---- 
#挑出表達矩陣里方差最大的1000個探針的名字！
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,]

# 直接畫熱圖，對比不鮮明——因為有些熱圖顏色對應的取值范圍默認是從數(shù)據(jù)最小值到最大值，離群值會造成大部分數(shù)據(jù)的熱圖顏色相近
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         #cluster_cols = F, #不按分組聚類
         annotation_col=annotation_col
)

# 按行標準化
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row", #標準化的依據(jù)
         breaks = seq(-3,3,length.out = 100) #顏色分割
         ) 
dev.off()

04_DEG.R

rm(list = ls()) 
load(file = "step2output.Rdata")
#差異分析，用limma包來做
#需要表達矩陣和Group，不需要改
library(limma)
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf) #以上操作為用limma包對表達矩陣進行一系列基于統(tǒng)計學原理的處理分析。

#為deg數(shù)據(jù)框添加幾列
#1.加probe_id列，把行名變成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
#2.加上探針注釋
ids = ids[!duplicated(ids$symbol),]
deg <- inner_join(deg,ids,by="probe_id") #隨機去掉對應同一基因的多個探針，只留下一個探針
head(deg)
nrow(deg)
#3.加change列,標記上下調基因
logFC_t=1
P.Value_t = 0.05 #logFC和p-value可以自己定義
k1 = (deg$P.Value < P.Value_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < P.Value_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)
#4.加ENTREZID列，用于富集分析（symbol轉entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol, 
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db)#人類,不同物種的OrgDb不同?。?！
#其他物種http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
dim(deg)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
dim(deg)
length(unique(deg$symbol))
save(Group,deg,logFC_t,P.Value_t,gse_number,file = "step4output.Rdata")

05_plots.R

rm(list = ls()) 
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#1.火山圖----
library(dplyr)
library(ggplot2)
dat  = deg[!duplicated(deg$symbol),] #去掉一個探針對應多個gene symbol的情況

p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, 
             aes(color=change)) +
  ylab("-log10(Pvalue)")+
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()
p

for_label <- dat%>% 
  filter(symbol %in% c("HADHA","LRRFIP1"))

volcano_plot <- p +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(
    aes(label = symbol),
    data = for_label,
    color="black"
  )
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse_number,"_volcano.png"))

#2.差異基因熱圖----

load(file = 'step2output.Rdata')
# 表達矩陣行名替換
exp = exp[dat$probe_id,]
rownames(exp) = dat$symbol
if(F){
  #全部差異基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}else{
  #取前10上調和前10下調
  x=dat$logFC[dat$change !="stable"] 
  names(x)=dat$symbol[dat$change !="stable"] 
  cg=names(c(head(sort(x),10),tail(sort(x),10)))
  length(cg)
}
n=exp[cg,]
dim(n)

#差異基因熱圖
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) 
heatmap_plot <- pheatmap(n,show_colnames =F,
                         scale = "row",
                         #cluster_cols = F, 
                         annotation_col=annotation_col,
                         breaks = seq(-3,3,length.out = 100)
) 
heatmap_plot
ggsave(heatmap_plot,filename = paste0(gse_number,"_heatmap.png"))

# 3.感興趣基因的箱線圖----
g = c(head(cg,3),tail(cg,3))
library(tidyr)
library(tibble)
library(dplyr)
dat = t(exp[g,]) %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  mutate(group = Group)

pdat = dat%>% 
  pivot_longer(cols = 2:(ncol(dat)-1),
               names_to = "gene",
               values_to = "count")

pdat$gene = factor(pdat$gene,levels = cg,ordered = T)
pdat$change = ifelse(pdat$gene %in% head(cg,10),"down","up")
library(ggplot2)
library(paletteer)
box_plot = ggplot(pdat,aes(gene,count))+
  geom_boxplot(aes(fill = group))+
  #scale_fill_manual(values = c("blue","red"))+
  scale_fill_paletteer_d("basetheme::minimal")+
  geom_jitter()+
  theme_bw()+
  facet_wrap(~change,scales = "free")
box_plot
ggsave(box_plot,filename = paste0(gse_number,"_boxplot.png"))

# 4.感興趣基因的相關性----
library(corrplot)
M = cor(t(exp[g,]))
pheatmap(M)

my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
my_color = colorRampPalette(my_color)(10)
corrplot(M, type="upper",
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)
library(cowplot)
cor_plot <- recordPlot() 

# 拼圖
load("pca_plot.Rdata")
library(patchwork)
library(ggplotify)
(pca_plot + volcano_plot +as.ggplot(heatmap_plot))/box_plot

plot_grid(cor_plot,heatmap_plot$gtable)

06_anno.R

rm(list = ls())  
load(file = 'step4output.Rdata')
library(clusterProfiler)
library(ggthemes)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(stringr)
library(enrichplot)

# 1.GO 富集分析----

#(1)輸入數(shù)據(jù)
gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down)

#(2)富集
#以下步驟耗時很長，設置了存在即跳過
if(!file.exists(paste0(gse_number,"_GO.Rdata"))){
  ego <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "ALL",
                  readable = TRUE)
  ego_BP <- enrichGO(gene = gene_diff,
                  OrgDb= org.Hs.eg.db,
                  ont = "BP",
                  readable = TRUE)
  #ont參數(shù)：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,ego_BP,file = paste0(gse_number,"_GO.Rdata"))
}
load(paste0(gse_number,"_GO.Rdata"))

#(3)可視化
#條帶圖
barplot(ego)
#氣泡圖
dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") + 
  scale_y_discrete(labels = function(x) str_wrap(x, width = 45))

#geneList 用于設置下面圖的顏色
geneList = deg$logFC
names(geneList)=deg$ENTREZID

#(3)展示top通路的共同基因，要放大看。
#Gene-Concept Network
cnetplot(ego,categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego, showCategory = 3,foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map,這個函數(shù)最近更新過，版本不同代碼會不同

Biobase::package.version("enrichplot")

if(T){
  emapplot(pairwise_termsim(ego)) #新版本
}else{
  emapplot(ego)#舊版本
}
#(4)展示通路關系 https://zhuanlan.zhihu.com/p/99789859
#goplot(ego)
goplot(ego_BP)

#(5)Heatmap-like functional classification
heatplot(ego,foldChange = geneList,showCategory = 8)

# 2.KEGG pathway analysis----
#上調、下調、差異、所有基因
#（1）輸入數(shù)據(jù)
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)

#（2）對上調/下調/所有差異基因進行富集分析

if(!file.exists(paste0(gse_number,"_KEGG.Rdata"))){
  kk.up <- enrichKEGG(gene         = gene_up,
                      organism     = 'hsa')
  kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa')
  kk.diff <- enrichKEGG(gene         = gene_diff,
                        organism     = 'hsa')
  save(kk.diff,kk.down,kk.up,file = paste0(gse_number,"_KEGG.Rdata"))
}
load(paste0(gse_number,"_KEGG.Rdata"))

#(3)看看富集到了嗎？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)
table(kk.up@result$p.adjust<0.05)
table(kk.down@result$p.adjust<0.05)

#(4)雙向圖
# 富集分析所有圖表默認都是用p.adjust,富集不到可以退而求其次用p值，在文中說明即可
down_kegg <- kk.down@result %>%
  filter(pvalue<0.05) %>% #篩選行
  mutate(group=-1) #新增列

up_kegg <- kk.up@result %>%
  filter(pvalue<0.05) %>%
  mutate(group=1)

source("kegg_plot_function.R") #加載某個函數(shù)
g_kegg <- kegg_plot(up_kegg,down_kegg)
g_kegg
#g_kegg +scale_y_continuous(labels = c(2,0,2,4,6))
ggsave(g_kegg,filename = 'kegg_up_down.png')

# 3.GSEA作kegg和GO富集分析----
#https://www.yuque.com/xiaojiewanglezenmofenshen/dbwkg1/ytawgg

#(1)查看示例數(shù)據(jù)
data(geneList, package="DOSE")
#(2)將我們的數(shù)據(jù)轉換成示例數(shù)據(jù)的格式
geneList=deg$logFC
names(geneList)=deg$ENTREZID
geneList=sort(geneList,decreasing = T)
#(3)GSEA富集
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',
                  verbose      = FALSE)
down_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore < 0,];down_kegg$group=-1
up_kegg<-kk_gse[kk_gse$pvalue<0.05 & kk_gse$enrichmentScore > 0,];up_kegg$group=1
#(4)可視化
g2 = kegg_plot(up_kegg,down_kegg)
g2

07_string.R

ppi網(wǎng)絡分析通過string網(wǎng)站和Cytoscape完成，這里顯示如何用代碼生成輸入數(shù)據(jù)。

rm(list = ls())
load("step4output.Rdata")
gene_up= deg[deg$change == 'up','symbol'] 
gene_down=deg[deg$change == 'down','symbol']
gene_diff = c(gene_up,gene_down)

# 1.制作string的輸入數(shù)據(jù)
write.table(gene_diff,
            file="diffgene.txt",
            row.names = F,
            col.names = F,
            quote = F)
# 從string網(wǎng)頁獲得string_interactions.tsv
# 2.準備cytoscape的輸入文件
p = deg[deg$change != "stable",
        c("symbol","logFC")]
head(p)
write.table(p,
            file = "deg.txt",
            sep = "\t",
            quote = F,
            row.names = F)
# string_interactions.tsv是網(wǎng)絡文件
# deg.txt是屬性表格

一些相關學習資料
GSEA：
https://www.yuque.com/docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?#
富集分析：
http://yulab-smu.top/clusterProfiler-book/index.html
弦圖：https://www.yuque.com/xiaojiewanglezenmofenshen/dbwkg1/dgs65p
GOplot：
https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew

?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容