GEO数据库中芯片数据分析思路

#实战代码有很多注意事项， 请不要不听课直接跑代码。
rm(list = ls())
library(GEOquery)
#先去网页确定是否是表达芯片数据，不是的话不能用本流程。
gse_number = "GSE42872"
eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
##下载并读取数据，本地有的话会识别，本地不完整的话要删掉
##destdir = '.'下载在当前目录， getGPL = F不下载注释文件。
class(eSet)

## [1] "list"

length(eSet)

## [1] 1

eSet = eSet[[1]]
#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)

## [1] 33297     6

exp[1:4,1:4]

##         GSM1052615 GSM1052616 GSM1052617 GSM1052618
## 7892501    7.24559    6.80686    7.73301    6.18961
## 7892502    6.82711    6.70157    7.02471    6.20493
## 7892503    4.39977    4.50781    4.88250    4.36295
## 7892504    9.48025    9.67952    9.63074    9.69200

#检查矩阵是否正常，如果是空的就会报错，空的和有负值的、有异常值的矩阵需要处理原始数据。
#如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)。
#自行判断是否需要log
#exp = log2(exp+1)
boxplot(exp)

#(2)提取临床信息
pd <- pData(eSet)
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p

## [1] TRUE

if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
# 分组信息来自临床信息，分组信息需要与表达矩阵列名一一对应，
# 临床信息需要和表达矩阵列一一对应
#(4)提取芯片平台编号
gpl_number <- eSet@annotation;gpl_number

## [1] "GPL6244"

save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

# Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
# 生成Group向量的三种常规方法，三选一，选谁就把第几个逻辑值写成T，另外两个为F。如果三种办法都不适用，可以继续往后写else if
if(F){
  # 1.Group----
  # 第一种方法，有现成的可以用来分组的列
  Group = pd$`disease state:ch1` 
}else if(F){
  # 第二种方法，自己生成
  Group = c(rep("RA",times=13),
            rep("control",times=9))
  Group = rep(c("RA","control"),times = c(13,9))
}else if(T){
  # 第三种方法，使用字符串处理的函数获取分组
 Group=ifelse(str_detect(pd$title,"Control"),
               "control",
               "RA")##匹配关键词
}
​
# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("control","RA"))
Group

## [1] control control control RA      RA      RA     
## Levels: control RA

##factor因子有水平，即取值。
##levels水平有顺序，第一个位置是领头羊，是参考水平。
##levels水平可以默认生成，也可以自行指定。
##参考水平的用处：差异分析时自动作为对照组。

#捷径
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
ids <- AnnoProbe::idmap('GPL6244')
##是曾建明老师2020年开发的一款用于下载GEO数据集并注释的R包
?idmap##根据所给的GPL号，返回探针的注释
?geoChina##根据所给的GSE号，下载对应的表达矩阵
?annoGene##根据gencode中的GTF文件注释基因ID
​
#四种方法，方法1里找不到就从方法2找，以此类推。
#方法1 BioconductorR包(最常用)
gpl_number 

## [1] "GPL6244"

#在该网站找到芯片平台编号对应的物种和探针的注释信息R包
##http://www.bio-info-trainee.com/1399.html
if(!require(hugene10sttranscriptcluster.db))BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")

##  [1] "hugene10sttranscriptcluster"              "hugene10sttranscriptcluster.db"          
##  [3] "hugene10sttranscriptcluster_dbconn"       "hugene10sttranscriptcluster_dbfile"      
##  [5] "hugene10sttranscriptcluster_dbInfo"       "hugene10sttranscriptcluster_dbschema"    
##  [7] "hugene10sttranscriptclusterACCNUM"        "hugene10sttranscriptclusterALIAS2PROBE"  
##  [9] "hugene10sttranscriptclusterCHR"           "hugene10sttranscriptclusterCHRLENGTHS"   
## [11] "hugene10sttranscriptclusterCHRLOC"        "hugene10sttranscriptclusterCHRLOCEND"    
## [13] "hugene10sttranscriptclusterENSEMBL"       "hugene10sttranscriptclusterENSEMBL2PROBE"
## [15] "hugene10sttranscriptclusterENTREZID"      "hugene10sttranscriptclusterENZYME"       
## [17] "hugene10sttranscriptclusterENZYME2PROBE"  "hugene10sttranscriptclusterGENENAME"     
## [19] "hugene10sttranscriptclusterGO"            "hugene10sttranscriptclusterGO2ALLPROBES" 
## [21] "hugene10sttranscriptclusterGO2PROBE"      "hugene10sttranscriptclusterMAP"          
## [23] "hugene10sttranscriptclusterMAPCOUNTS"     "hugene10sttranscriptclusterOMIM"         
## [25] "hugene10sttranscriptclusterORGANISM"      "hugene10sttranscriptclusterORGPKG"       
## [27] "hugene10sttranscriptclusterPATH"          "hugene10sttranscriptclusterPATH2PROBE"   
## [29] "hugene10sttranscriptclusterPFAM"          "hugene10sttranscriptclusterPMID"         
## [31] "hugene10sttranscriptclusterPMID2PROBE"    "hugene10sttranscriptclusterPROSITE"      
## [33] "hugene10sttranscriptclusterREFSEQ"        "hugene10sttranscriptclusterSYMBOL"       
## [35] "hugene10sttranscriptclusterUNIPROT"

ids <- toTable(hugene10sttranscriptclusterSYMBOL)##以data.frame格式提取数据
head(ids)

##   probe_id     symbol
## 1  7896746   MTND1P23
## 2  7896754 SEPTIN7P13
## 3  7896759  LINC01128
## 4  7896761     SAMD11
## 5  7896779     KLHL17
## 6  7896798    PLEKHN1

# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有symbol列，也有的GPL平台没有提供注释表格
  b = read.delim("GPL570-55999.txt",
                 check.names = F,
                 comment.char = "#")
  colnames(b)
  ids2 = b[,c("ID","Gene Symbol")]
  colnames(ids2) = c("probe_id","symbol")
  k1 = ids2$symbol!="";table(k1)
  k2 = !str_detect(ids2$symbol,"///");table(k2)
  ids2 = ids2[ k1 & k2,]
  # ids = ids2
}
​
# 方法3 官网下载注释文件并读取
##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") ##保存数据

rm(list = ls())  
load(file = "step1output.Rdata")##导入数据
load(file = "step2output.Rdata")
#输入数据：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))##将数据转置，处理成PCA接受的数据格式。仅PCA适用，所以要赋值给dat。
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

save(pca_plot,file = "pca_plot.Rdata")
​
# 2.top 1000 sd 热图---- 
cg=names(tail(sort(apply(exp,1,sd)),1000))##apply，对exp按行求标准差。
n=exp[cg,]
​
# 直接画热图，对比不鲜明
library(pheatmap)
annotation_col=data.frame(group=Group)
rownames(annotation_col)=colnames(n) ##写一个分组信息的矩阵
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col
)
​
# 按行标准化
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",  ##按行标准化
         breaks = seq(-3,3,length.out = 100)##色带范围改变，随机生成100个数
         ) 

dev.off()

## png 
##   4

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)
​
#为deg数据框添加几列
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
#2.加上探针注释
ids = ids[!duplicated(ids$symbol),]
#其他去重方式在zz.去重方式.R
deg <- inner_join(deg,ids,by="probe_id")
nrow(deg)

## [1] 18857

#3.加change列,标记上下调基因
logFC_t=1
p_t = 0.05
k1 = (deg$P.Value < p_t)&(deg$logFC < -logFC_t)
k2 = (deg$P.Value < p_t)&(deg$logFC > logFC_t)
deg <- mutate(deg,change = ifelse(k1,"down",ifelse(k2,"up","stable")))
table(deg$change)

## 
##   down stable     up 
##    489  17782    586

#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)#人类
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
save(Group,deg,logFC_t,p_t,gse_number,file = "step4output.Rdata")

rm(list = ls()) 
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#1.火山图----
library(dplyr)
library(ggplot2)
dat  = distinct(deg,symbol,.keep_all = T)
​
p <- ggplot(data = dat, 
            aes(x = logFC, 
                y = -log10(P.Value))) +
  geom_point(alpha=0.4, size=3.5, 
             aes(color=change)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",linewidth=0.8) +
  geom_hline(yintercept = -log10(p_t),lty=4,col="black",linewidth=0.8) +
  theme_bw()
p

for_label <- dat%>% 
  filter(symbol %in% c("DUSP6","CD36"))##展示基因名
​
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

#2.差异基因热图----
​
load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,]
rownames(exp) = dat$symbol
if(T){
  #取前10上调和前10下调
  library(dplyr)
  dat2 = dat %>%
    filter(change!="stable") %>% 
    arrange(logFC) 
  cg = c(head(dat2$symbol,10),
         tail(dat2$symbol,10))
}else{
  #全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}
n=exp[cg,]
dim(n)

## [1] 20  6

#差异基因热图
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(-1,1,length.out =100)##色带范围改变，随机生成100个数
) 
heatmap_plot

#拼图
library(patchwork)
volcano_plot + heatmap_plot$gtable

# 3.感兴趣基因的相关性----
library(corrplot)
g = sample(deg$symbol[1:5000],10) # 这里是随机取样，注意换成自己感兴趣的基因
g

##  [1] "BVES"    "RETREG1" "ZBTB2"   "ENPP3"   "EMC8"    "MATN2"   "POLE"    "MRAS"    "DIPK1B" 
## [10] "SMAD5"

M = cor(t(exp[g,])) ###cor计算相关性。转置后基因变成列，按列做相关性
pheatmap(M)

library(paletteer)#配色包
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")
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)

dev.off()

## null device 
##           1

#保存
pdf("deg.pdf",width = 10,height = 10)
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()

## null device 
##           1

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)输入数据
gene_up = deg$ENTREZID[deg$change == 'up'] 
gene_down = deg$ENTREZID[deg$change == 'down'] 
gene_diff = c(gene_up,gene_down)
​
#(2)富集
#以下步骤耗时很长，设置了存在即跳过
f = paste0(gse_number,"_GO.Rdata")
if(!file.exists(f)){
  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参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.
  save(ego,ego_BP,file = f)
}
load(f)
​
#(3)可视化
#条带图
barplot(ego)

barplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

#气泡图
dotplot(ego)

dotplot(ego, split = "ONTOLOGY", font.size = 10, 
        showCategory = 5) + 
  facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y") 

#(3)展示top通路的共同基因，要放大看。
​
#gl 用于设置下图的颜色
gl = deg$logFC
names(gl)=deg$ENTREZID
#Gene-Concept Network,要放大看
cnetplot(ego,
         #layout = "star",
         color.params = list(foldChange = gl),
         showCategory = 3)

# 2.KEGG pathway analysis----
#上调、下调、差异、所有基因
#（1）输入数据
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)
​
#（2）对上调/下调/所有差异基因进行富集分析
f2 = paste0(gse_number,"_KEGG.Rdata")
if(!file.exists(f2)){
  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 = f2)
}
load(f2)
##返回空值别紧张，看看帮助文档~  https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
#(3)看看富集到了吗
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值，在文中说明即可
​
source("kegg_plot_function.R")
g_kegg <- kegg_plot(kk.up,kk.down)

g_kegg

#g_kegg +scale_y_continuous(labels = c(2,0,2,4,6))
​
# 3.能看懂的资料越来越多----
# 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
# 网上的资料和宝藏无穷无尽，学好R语言慢慢发掘~