GEO表达芯片数据分析

rm(list = ls())
library(GEOquery)

gse_number = "GSE56649" #自行修改编号！

eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
#下载方法备选：geoChina函数，优点下载速度快，缺点较新数据未收入
# 2018年及以前的芯片数据可以用geoChina快速下载
#library(AnnoProbe)
#eSet <- geoChina(gse_number, destdir = '.') 
class(eSet)
#> [1] "list"
length(eSet)
#> [1] 1
eSet = eSet[[1]] #当一个数据集有多个平台时，分别提取子集，当成多个数据分析

#(1)提取表达矩阵信息exp
exp <- exprs(eSet)
dim(exp)
#> [1] 54675    22
exp[1:4,1:4]
#>           GSM1366348 GSM1366349 GSM1366350 GSM1366351
#> 1007_s_at    279.156    202.866    406.818    232.668
#> 1053_at      487.700    409.018    393.810    397.127
#> 117_at       666.869    388.589    704.633    953.481
#> 121_at       240.646    361.198    305.229    374.757
# 目的：
#1）检查矩阵是否正常：
#如果是空的就会报错，空的和有异常值的矩阵需要处理原始数据。
#如果表达矩阵为空，大多数是转录组数据，不能用这个流程(后面另讲)
#2）判断是否需要取log：
#几百、几千的话就需要取log，一般0-20的说明取过log了，不用再取
#没取过log，有负值（光信号值没负数）→错误数据（弃用或处理原始数据）
#取过log，有负值（数值偏小时）→正常

exp = log2(exp+1)  #自行判断，不需要取log时注释掉

boxplot(exp)

# 目的：
#1）再次判断是否需要取log：
#箱线图如果表现为压得很扁，说明没取log
#2）检查并处理异常样本：
#如果有样本的表达量整体低于其他的一大截，运行如下代码去除异常
#exp = limma::normalizeBetweenArrays(exp)
#boxplot(exp)
#3）表达矩阵的数字一半正一半负：做了标准化，弃用或处理原始数据

#(2)提取临床信息pd
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] "GPL570"
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)
#（1）获取实验分组(Group)-----------------
# 标准流程代码是二分组，多分组数据的分析不适用，后面另讲
# 生成Group向量有三种常规方法，根据数据三选一即可
# 选谁就把第几个逻辑值写成T，另外两个为F
# 如果三种办法都不适用，可以继续往后写else if

if(F){
  # 第一种方法：有现成的可以用来分组的列
  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$source_name_ch1,"control"),
               "control",
               "RA")
}
# 只运行了第三种方法，因为只有第三种方法的if后面是T

Group = factor(Group,levels = c("control","RA")) #修改level！

Group
#>  [1] RA      RA      RA      RA      RA      RA      RA     
#>  [8] RA      RA      RA      RA      RA      RA      control
#> [15] control control control control control control control
#> [22] control
#> Levels: control RA
# Group是一个分类型向量，分类型数据可以用因子factor组织
# 需要把Group转换成因子，并设置参考水平，指定levels
# 注意levels顺序要设置成对照组在前，处理组在后
# 因子正文与levels不对应时，会产生NA

#（2）获取探针注释(ids)-----------------
# 探针注释：探针与基因的对应关系，probe_id和symbol
# 捷径注释的方法(首选，能找到70%-80%，失败后选方法1，以此类推)
library(tinyarray)
find_anno(gpl_number) #打出找注释的代码
#> [1] "`library(hgu133plus2.db);ids <- toTable(hgu133plus2SYMBOL)` and `ids <- AnnoProbe::idmap('GPL570')` are both avaliable"

ids <- AnnoProbe::idmap('GPL570')  #修改芯片平台编号！

head(ids) #调试
#>         probe_id symbol
#> 193731   1053_at   RFC2
#> 193732    117_at  HSPA6
#> 193733    121_at   PAX8
#> 193734 1255_g_at GUCA1A
#> 193735   1316_at   THRA
#> 193736   1320_at PTPN21
#四种探针注释的方法，先看捷径，方法1里找不到就从方法2找，以此类推-----------------
#方法1 Bioconductor R包(最常用，但不是万能的!)
if(F){ 
  #要用方法1时才需要改成T
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")  #看R包有哪些数据，看SYMBOL的出处
ids <- toTable(hgu133plus2SYMBOL) #提取R包的注释结果
head(ids)
}
# 方法2 读取GPL网页的表格文件，按列取子集
##https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
if(F){ 
  #要用方法2时才需要改成T
  #注：表格读取参数、文件列名不统一，活学活用，有的表格里没有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())  
#输入数据：exp和Group
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")

#（1）PCA 图（样本聚类图）----
# Principal Component Analysis
# http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials
# 主成分分析，用于预实验，简单看看组间是否有差别
# 同一分组是否自成一簇（组内重复好），中心点间是否有距离（组间差别大）
# 圆圈代表置信区间，样本数量少时无法计算置信区间故没有圆圈
# 中心点不是样本点，只是代表位置的中心
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

save(pca_plot,file = "pca_plot.Rdata")

#（2）top 1000 sd 热图---- 
cg=names(tail(sort(apply(exp,1,sd)),1000))
n=exp[cg,] #筛选标准差最大的1000个基因（最活跃、离散的基因）

# 直接画热图，对比不鲜明
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
)

# 按行标准化（scale = "row"）
#比较的是一行即一个基因在不同分组间的差别
#放大行内部的差别，抹平行与行之间的差别，图例范围缩小
#使基因在不同分组的表达差异更明显
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         scale = "row",
         breaks = seq(-3,3,length.out = 100) #设置色带范围，避免极端值影响整体配色
         ) 

dev.off()
#> null device 
#>           1

rm(list = ls()) 
load(file = "step2output.Rdata")
#需要表达矩阵exp和分组信息Group，不需要改
library(limma)
design=model.matrix(~Group)           #构建一个模型矩阵
fit=lmFit(exp,design)                 #线性拟合
fit=eBayes(fit)                       #贝叶斯检验
deg=topTable(fit,coef=2,number = Inf) #提取差异分析结果
#可以将上面几句代码运用简化思维组织成一个函数function(exp,Group)

#为deg数据框添加几列
#（1）加probe_id列，把行名变成一列----
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))

#（2）加上探针注释symbol列，随机去重更客观----
# 随机去重,因为存在多个探针对应一个symbol
ids = ids[!duplicated(ids$symbol),]      #方法1
ids = distinct(ids,symbol,.keep_all = T) #方法2
# 随机去重更客观，其他去重方式----
if(F){
  #保留最大值去重
  exp2 = exp[ids$probe_id,]
  identical(ids$probe_id,rownames(exp2))
  ids = ids[order(rowSums(exp2),decreasing = T),]
  ids = ids[!duplicated(ids$symbol),];nrow(ids)
  # 拿这个ids去inner_join
}
if(F){
  exp3 = exp[ids$probe_id,]
  rownames(exp3) = ids$symbol
  exp3[1:4,1:4]
  exp4 = limma::avereps(exp3)
  # 此时拿到的exp4已经是一个基因为行名的表达矩阵
  # 直接差异分析，不再需要inner_join 
}
# 用inner_join加上symbol列----
deg <- inner_join(deg,ids,by="probe_id")
nrow(deg)
#> [1] 20188

#（3）加change列,标记上下调基因----
#logFC_t和p_t可自行调整
logFC_t=1 #差异基因过多，可调大（2.2,2,1.5,1.2,0.585=log2(1.5)）
p_t = 0.05#P值越小，差异越显著，可调小（0.01）
#P.Value可调成adj.P.Val的列名
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 
#>    796  18753    639

#（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)#转换依据的数据库（人类）
# symbol和entrezid并非一一对应，损失或增加一些很正常
# 其他物种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) #按symbol列去重复
#去掉添加ENTREZID后多的行
#防止需要添加到图中的symbol重复出现

#1）作火山图
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

#2）添加火山图的标签：
#① 自己添加感兴趣的基因
#② P.Value的top10
#③ 根据logFC挑也行

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

#（2）差异基因热图----
load(file = 'step2output.Rdata')
# 表达矩阵行名替换
exp = exp[dat$probe_id,]
rownames(exp) = dat$symbol  #转变成以基因为行名的表达矩阵
if(T){
  #取前10上调和前10下调，可以根据logFC（更合适），也可以根据Pvalue
  library(dplyr)
  dat2 = dat %>%
    filter(change!="stable") %>% 
    arrange(logFC)              #将logFC从小到大排序
  cg = c(head(dat2$symbol,10),  #logFC最小的10个基因
         tail(dat2$symbol,10))  #logFC最大的10个基因
}else{ #没有运行
  #全部差异基因
  cg = dat$symbol[dat$change !="stable"]
  length(cg)
}
n=exp[cg,]
dim(n)
#> [1] 20 22
#绘制差异基因热图
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

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

#$gtable提取热图的主体部分，才能完成拼图的操作

#（3）感兴趣基因的相关性热图----
# 注意：这里的流程只适用于多个基因的相关性
# 两个基因的相关性需要用相关性点图表示
library(corrplot)
g = sample(deg$symbol[1:500],10) 
# 这里是随机取样，注意换成自己感兴趣的基因
g
#>  [1] "YWHAB"   "P3H2"    "COIL"    "ARL6IP5" "MYDGF"   "R3HDM4" 
#>  [7] "ALYREF"  "RIF1"    "PBX2"    "TDG"
M = cor(t(exp[g,]))  #相关系数矩阵
#计算的是列与列的相关性，需要转置才能计算基因相关性
pheatmap(M)        

library(paletteer) #配色R包
my_color = rev(paletteer_d("RColorBrewer::RdYlBu"))
#选出一组颜色，通过逆转配色，使红色代表正相关，蓝色代表负相关
my_color = colorRampPalette(my_color)(10) #切10种颜色
corrplot(M, type="upper",
         method="pie",
         order="hclust", 
         col=my_color,
         tl.col="black", 
         tl.srt=45)

library(cowplot)         #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)
#width,height可以调整导出图片的比例
plot_grid(pca_plot,cor_plot,
          volcano_plot,heatmap_plot$gtable)
dev.off()
#> null device 
#>           1

# 通过判断某个基因的箱线图和logFC是否匹配
check_gene = deg$symbol[1]
boxplot(exp[check_gene,]~Group) #处理组<对照组→logFC<0

deg$logFC[1]                    #logFC<0成立
#> [1] -1.633518

# 富集分析：衡量每个通路里的基因在差异基因里是否足够多
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) #直接转换为symbol
  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)可视化
#1）条带图
barplot(ego)

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

#分面，根据自己的选择决定是否运行

#2）气泡图
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)

#4）展示GO term 间的关系
goplot(ego_BP) #只能画一个子类

# https://zhuanlan.zhihu.com/p/99789859
# 有颜色的点是富集到的通路
# 灰色的点是根据有颜色的点挖掘到的通路

# 2.KEGG pathway analysis----
#上调、下调、差异、所有基因
#（1）输入数据（与GO完全一致）
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)

#(3)看看富集到了吗？https://mp.weixin.qq.com/s/NglawJgVgrMJ0QfD-YRBQg
table(kk.diff@result$p.adjust<0.05)
#> 
#> FALSE 
#>   317
table(kk.up@result$p.adjust<0.05)
#> 
#> FALSE 
#>   292
table(kk.down@result$p.adjust<0.05)
#> 
#> FALSE 
#>   285
# 富集不到的话，很正常，但有补救方法：
# ①p.adjust改为P.Value
# ②换一种富集的方式GSEA
# ③调整差异基因的阈值

#(4)双向富集图
# 富集分析所有图表默认都是用p.adjust,富集不到可以退而求其次用p值，在文中说明即可
# 调用kegg_plot函数
library(ggthemes)
library(ggplot2)
kegg_plot <- function(kk.up,kk.down,p_t = 0.05,p.adjust = F){
  down_kegg <- kk.down@result %>%
    filter(pvalue<p_t) %>% #筛选行
    mutate(group=-1) #新增列
  
  up_kegg <- kk.up@result %>%
    filter(pvalue<p_t) %>%
    mutate(group=1)
  dat=rbind(up_kegg,down_kegg)
  if(p.adjust) dat$pvalue = dat$p.adjust
  colnames(dat)
  dat$pvalue = -log10(dat$pvalue)
  dat$pvalue=dat$pvalue*dat$group 
  
  dat=dat[order(dat$pvalue,decreasing = F),]
  
  gk_plot <- ggplot(dat,aes(reorder(Description, pvalue), y=pvalue)) +
    geom_bar(aes(fill=factor((pvalue>0)+1)),stat="identity", width=0.7, position=position_dodge(0.7)) +
    coord_flip() +
    scale_fill_manual(values=c("#0072B2", "#B20072"), guide="none") +
    labs(x="", y="" ) +
    theme_pander()  +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          #axis.ticks.x = element_blank(),
          axis.line.x = element_line(linewidth = 0.3, colour = "black"),#x轴连线
          axis.ticks.length.x = unit(-0.20, "cm"),#修改x轴刻度的高度，负号表示向上
          axis.text.x = element_text(margin = margin(t = 0.3, unit = "cm")),##线与数字不要重叠
          axis.ticks.x = element_line(colour = "black",linewidth = 0.3) ,#修改x轴刻度的线                         
          axis.ticks.y = element_blank(),
          axis.text.y  = element_text(hjust=0),
          panel.background = element_rect(fill=NULL, colour = 'white')
    )
}
# 作图：
g_kegg <- kegg_plot(kk.up,kk.down,p_t = 0.01,p.adjust = F)
g_kegg

# 注意：
#①p.adjust可改为T，按照p.adjust作图
#②可以在括号里自己用tab加参数
#③根据不同的横坐标改变，因为理论上横坐标不可能为负数
#g_kegg +scale_y_continuous(labels = c(2,0,2,4,6))