GEO数据挖掘实战全流程1.27万字笔记 by 生信技能树

这个笔记主要是根据生信技能树数据挖掘线上直播课以及B站视频做的，GEO芯片数据分析部分。如果大家也想完整的学习，可以去跟这个马拉松授课，生信入门&数据挖掘线上直播课4月班

一、数据下载与读取

1. 使用R包 GEOquery 下载

推荐用getGEO函数下载，通过GSE号下载后得到的 gset 是一个 ExpressionSet 对象。这个对象被封装在1个 list 中，因为 GSE中可能包含了多个 GPL 平台的数据，不同平台的数据会存放在 list 中，这个 list 的每个对象都是 1 个平台的 ExpressionSet 对象。

这里需要做一下区分，如果通过getGEO读取下载到本地的GSE数据，那么读取进来直接就是ExpressionSet 对象而不是list。

GEOquery包的使用

2. geoChina()

生信技能树健明老师写的 AnnoProbe 包中的 geoChina 函数，提供了一个中国区镜像，可以下载部分 GEO 数据。

library(AnnoProbe)
eSet <- geoChina(gse_number, destdir = '.')

3. NCBI 直接下载

比如 GSE168113_series_matrix.txt。下载后在R中读取txt文件，read.table 函数参数怎么设置提前预览一下txt文件中的内容格式就知道了。

a<- read.table("GSE42872_series_matrix.txt.gz",
               sep = "\t",header = T,
               # row.names=1,   # 把探针名所在的第1列直接读取为行名
               comment.char = "!")

另外是直接下载原始数据。不推荐，因为不同的芯片rawdata 格式不同，需要学习不同的处理方法。比如 GSE168113_RAW.tar

4. 读取芯片数据的方法

三大芯片制造商：affymetrix , illumina , agilent

芯片类型：表达谱芯片（mRNA）、基因分析芯片（SNP,indel）、CNV芯片（拷贝数变异）、甲基化芯片、miRNA芯片

从GEO数据库下载得到表达矩阵一文就够 (qq.com)

用affy包读取 affymetix 的基因表达芯片
用oligo包来读取 affymetix 的基因表达芯片
用R包lumi来处理illumina的bead系列表达芯片

5. 实战(Step 0)

数据集背景了解

数据集: GSE42872 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42872。
数据类型：芯片数据
平台：GPL6244
两个组，实验组 3 ，对照组 3 。宿主为人类

实验设计：用威罗非尼（Vemurafenib）或 DMSO 处理 BRAFV600E A375 黑色素瘤细胞(melanoma)。

Summary: 威罗非尼是一种 BRAF 抑制剂，对黑色素瘤中最常见的 BRAF 突变体（BRAFV600E）具有特异性。威罗非尼通过抑制 MEK / ERK 丝裂原活化蛋白激酶的下游活化来抑制 BRAF 突变人黑色素瘤细胞的增殖。我们使用微阵列检查威罗非尼敏感的BRAFV600E人黑色素瘤细胞系（A375）对威罗非尼的转录反应，以进一步描述 BRAFV600E 驱动人类黑色素瘤细胞增殖和能量代谢的机制。PMID:24469106

6. 实战(Step 1)

数据下载与导入，获取表达矩阵。

rm(list = ls())
library(GEOquery)
gse_number = "GSE42872"
eSet <- getGEO(gse_number, destdir = '.', getGPL = F)
class(eSet)
length(eSet)
eSet = eSet[[1]]
#(1)提取表达矩阵exp
exp <- exprs(eSet)
dim(exp)
boxplot(exp) # 已经log 转化
exp[1:4,1:4]
#(2)提取临床信息
pd <- pData(eSet)
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl_number <- eSet@annotation;gpl_number
save(gse_number,pd,exp,gpl_number,file = "step1output.Rdata")

二、分组信息与探针注释

1. 获取分组信息

pData 函数可以获得临床信息。

pd <- pData(eSet)
# 第一种方法，有现成的可以用来分组的列
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")

2. 获取芯片注释的方法

(1) 捷径 tinyarray

library(tinyarray) # 小洁老师写的R包
find_anno(gpl_number) # 打出找注释的代码
ids <- AnnoProbe::idmap('GPL570') # 曾老师写的R包

(2) bioconductor 包

安装与加载包的时候注意在包名后面加 .db，如 library(org.Hs.eg.db)

GPL与R包对应关系可参考：

用R获取芯片探针与基因的对应关系三部曲-bioconductor | 生信菜鸟团 (bio-info-trainee.com)http://www.bio-info-trainee.com/1399.html
GPL对应的Bioconductor注释包https://blog.csdn.net/weixin_40739969/article/details/103186027

# BiocManager::install("hugene10sttranscriptcluster.db")  # R包安装
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")   # 查看这个包中的对象
ids=toTable(hugene10sttranscriptclusterSYMBOL) # 得到每个探针对应的SYMBOL,这里括号中的根据前面查找的对象可以更改

(3) NCBI 下载文件解析

第一种方法：直接用 GEOquery 下载 soft 文件（soft 文件较大，需要网速）

# 第一种方式：fData 函数提取表格，得到的 fd 表格即含有芯片注释信息
gse42589 <- getGEO('GSE42589', destdir=".",getGPL = T)
eSet = gse42589[[1]] # 下载如果是列表先提取对象
fd <- fData(eSet)
colnames(fd) 
# 第二种方式：Table 函数提取表格，不过是直接从 soft 文件提取的，得到的表与第一种方法一致
gpl6244 <- getGEO('GPL6244', destdir=".")
fd <- Table(gpl6244)
colnames(fd)

第二种方法：可以直接在浏览器上从NCBI里面下载文件来解析，如GEO Accession viewer (nih.gov)https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244

(4) 其它方法

去基因芯片的厂商的官网直接去下载
自主注释

3. 实战(Step 2)

获取分组信息与探针注释

# 载入数据
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)

# 获取分组信息
  Group=ifelse(str_detect(pd$title,"Control"),
               "control",
               "treat")
colnames(exp) <- paste(Group,c(1,2,3),sep = "-");colnames(exp)
# 需要把Group转换成因子，并设置参考水平，指定levels，对照组在前，处理组在后
Group = factor(Group,levels = c("control","treat"))
Group

# 探针注释的获取
library(tinyarray) # 小洁老师写的R包
find_anno(gpl_number) #打出找注释的代码
ids <- AnnoProbe::idmap('GPL6244')
save(exp,Group,ids,gse_number,file = "step2output.Rdata")

三、表达矩阵数据检查

1. 基础检查 tips

log 转化

如何判断下载的表达矩阵是否进行过 log 转化？ => 箱线图的纵坐标如果在 0-20 之间代表是 log 后的表达矩阵

箱线图查看异常样本

不同样本箱线图整体表达水平应该基本一致，因为一般来说差异基因只是少数，不会对整体的数据分布产生太大的影响，下图就存在一个异常样本：

异常样本处理办法

删除异常样本
将样本表达水平拉到同一个水平线：

exp <- limma::normalizeBetweenArrays(exp)

表达矩阵里有负值

进行过 log 转化，正常。
没进行过 log 转化，可能存在错误数据，一般弃用，处理原始数据。
有一半负值，代表进行了 z-score 等标准化，不能用于差异分析，应处理原始数据。

检查某些基因的表达量

如 GAPDH 管家基因

exprSet["GAPDH",]     
boxplot(exprSet[,1])  # 发现GAPDH的表达量很高，因为是管家基因,如果很低反而有问题

2. 绘图

注意, 某个图并没有体现出良好的分组并不是说一定要去除某个离群样本,我们只是用这些图大致检查一下数据质量,最后需不需要去除某个样本还需要结合具体实验设计

(1) 箱线图

箱线图是一种用于展示一组数据分布情况的图表。它可以帮助我们直观地了解数据的中位数、最小值、最大值、四分位数以及异常值等统计指标。具体来说，箱线图由一个矩形和两条线组成。矩形的边界由第一四分位数和第三四分位数组成，而矩形内部的线条代表中位数。箱线图中的两条线分别表示数据的最小值和最大值。如果在最大值和最小值之外存在异常值，这些值将被表示为离群点。

通过观察箱线图，我们可以直观地看出数据的集中趋势和分散程度。如果矩形较小且中位数靠近矩形中央，说明数据分布相对集中；如果矩形较大且中位数偏离矩形中央，说明数据分布相对分散。此外，如果存在离群点，我们需要考虑这些点是否是数据错误或异常值，需要进行进一步的数据分析和处理。

ggplot2 制作箱线图需要对表达矩阵（数据框）进行长宽变换，一列为样本信息作为 x 轴，一列为每个样本每个基因的表达量，作为 y 轴，如果需要添加分组信息，还需要添加一列分组信息。

原始矩阵：行为探针/基因，列为样本。

ggplot2画箱线图输入的数据：

实战(Step 3.1)

对表达矩阵进行长宽变换后用 ggplot2 绘制箱线图。

# 载入数据
rm(list = ls())  
load(file = "step1output.Rdata")
load(file = "step2output.Rdata")

# （1）箱线图
library(reshape2)
library(tibble)
exp <- rownames_to_column(as.data.frame(exp),var = "probe_id") # 行名转换为第1列
exp_L<- melt(exp) # 长宽变换
colnames(exp_L)<- c("probe_id","sample","value");table(exp_L$sample) 
exp_L$group<- rep(Group,each=nrow(exp)) # 添加分组列
head(exp_L)
library(ggplot2)
p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()+
  theme(axis.text.x = element_text(angle=0, vjust=0.6))
print(p)
exp<- column_to_rownames(exp,"probe_id")

(2) 层次聚类

层次聚类笔记

输入：距离矩阵（通过 dist 函数获得）：

> sampleDists
          control-1 control-2 control-3  treat-1  treat-2  treat-3
control-1   0.00000  39.71222  34.36004 81.76124 88.19851 86.35148
control-2  39.71222   0.00000  43.05040 85.04253 90.11748 87.08480
control-3  34.36004  43.05040   0.00000 83.72428 86.22464 89.42490
treat-1    81.76124  85.04253  83.72428  0.00000 40.85336 41.91346
treat-2    88.19851  90.11748  86.22464 40.85336  0.00000 49.70981
treat-3    86.35148  87.08480  89.42490 41.91346 49.70981  0.00000

实战(Step 3.2)

用 dist 函数得到距离矩阵，然后用 factoextra 包的 hcut 函数进行层次聚类，再用 fviz_dend 函数进行可视化。

library(factoextra) 
sampleDists <- dist(t(exp),upper = T,diag=T)  # 得到距离矩阵
# sampleDists1 <- dist(scale(t(exp)),upper = T,diag=T)
res <- hcut(          
  sampleDists,k = 2,   # k表示分为两类
  hc_func = "hclust",     # 功能：聚类
  hc_method = "complete", # 距离定义方式：最大距离
  # hc_method = "average", # 距离定义方式：平均距离
  hc_metric = "euclidean", # 欧式距离
  stand = FALSE,graph = FALSE)                        
fviz_dend(res,rect = TRUE,cex = 0.75, 
          # type="phylogenic",  # 树形图
          type="rectangle", # 矩形图
          horiz = T,
          k_colors = c("#0073C2FF","#CD534CFF"))

(3) 热图

热图是一种可视化方法，用于显示数据中不同部分之间的相对值或强度。热图通常以矩阵形式呈现，其中每个单元格的颜色表示该单元格所代表的值的大小。颜色通常是一个渐变色条，其中深色表示较高的值，浅色表示较低的值。

pheatmap可以绘制具有不同颜色映射方案的热图，并允许用户对行和列进行聚类。在绘制热图之前，pheatmap 通常会对矩阵数据进行预处理，包括对数据进行缩放、标准化、聚类等操作（可通过参数进行设置）。其中，缩放和标准化可用于使数据在不同变量之间具有相同的尺度，以便更好地比较和解释结果。聚类可以将相关变量和行/列聚类在一起，以更好地显示数据中的模式和结构。pheatmap默认聚类、不标准化。pheatmap绘图时通常使用一个渐变色条来表示矩阵数据中的数值大小，其中深色表示高值，浅色表示低值。在 pheatmap 中，默认的层次聚类计算数据点之间距离的算法为欧氏距离euclidean，可更改为其它距离，与 dist() 函数一致；默认的计算簇之间距离的算法为 complete，可更改为其它距离算法，与 hclust() 函数一致。

按行标准化意味着对每一行进行标准化处理，使得每一行的值都满足均值为0，标准差为1的正态分布。按行标准化可以在样本之间比较基因表达量的差异时，将不同样本的数据进行比较，消除了每个样本中表达量大小的影响，使得每个样本的表达量数据更加可比性。

按列标准化意味着对每一列进行标准化处理，使得每一列的值都满足均值为0，标准差为1的正态分布。按列标准化可以在基因之间比较它们在样本中的表达模式时，将不同基因的数据进行比较，消除了基因表达量大小的影响，使得每个基因的表达量数据更加可比性。

另外还可以画样本之间相关性的热图，比如使用 cor 函数计算得到样本相关性，然后画图。

输入：pheatmap 绘制热图直接输入原始矩阵即可

实战(Step 3.3)

筛选方差前1000大的基因，然后用 pheatmap 画基因表达热图。

cg=names(tail(sort(apply(exp,1,sd)),1000)) # 筛选方差前1000大的基因
n=exp[cg,]
library(pheatmap)
group=data.frame(group=Group)
rownames(group)=colnames(n) 
pheatmap(n,
         show_colnames =F,
         show_rownames = F,
         annotation_col=group, # 添加分组信息
         scale = "row", # 按行标准化
         breaks = seq(-3,3,length.out = 100) # 颜色色阶调整
)

(4) PCA

主成分分析

主成分分析 (Principal Component Analysis, PCA) 是一种常用的数据分析技术，用于降维和数据压缩。PCA 通过找到数据的主要成分来实现这一目标，这些成分能够最好地解释数据中的变异性。具体来说，PCA 将原始数据投影到一组新的坐标轴上，使得数据在新坐标系中的方差最大化。这些新的坐标轴称为主成分，是原始数据的线性组合。主成分按照它们解释原始数据的方差的大小排序，因此，可以选择最高排名的主成分，以获得对原始数据的最佳压缩。根据这些主成分对样本进行聚类，代表样本的点在坐标轴上距离越远，说明样本差异越大。

PCA图主要考察：同一分组是否自成一簇（组内相似），中心点之间是否有距离（组间不一致）
PCA 图解释

PCA图是指在PCA过程中绘制的图形，通常用于可视化数据降维的结果。以下是PCA图的一些解释：

散点图：散点图是PCA图中最常见的形式，用于显示数据在新的低维空间中的分布情况。每个数据点都表示为一个点，该点的位置由其在低维空间中的坐标确定。
方差贡献率图：这种图形显示每个主成分解释的总方差的百分比。这可以帮助你确定在低维空间中保留多少主成分才能保留足够的信息。
热图：热图显示原始数据的相关性结构。通常使用颜色来表示相关性的强度，红色表示正相关性，蓝色表示负相关性。
散点矩阵图：散点矩阵图显示原始数据中所有可能的变量之间的散点图。这可以帮助你确定哪些变量对主成分的解释贡献最大。

在散点PCA图中，每个数据点代表一个样本，它的位置由它在主成分分析后的低维空间中的坐标确定。通常，散点PCA图上的不同指标可以表示以下信息：

X轴和Y轴：这些轴代表低维空间中的两个主成分。在散点PCA图上，X轴和Y轴通常代表最具有解释力的两个主成分。
数据点的位置：图上每个点代表一个样本，点与点之间的距离代表样本的差别。每个数据点在散点PCA图上的位置表示它在低维空间中的坐标。通常情况下，距离中心点较远的数据点表示数据中的一个特异点或异常值。而聚集在一起的数据点表示数据中具有相似性质的样本。
颜色：颜色可以用来表示数据点所属的类别或其他特征。例如，在聚类分析中，不同的簇可以用不同的颜色表示。
圆圈：某个簇样本地置信区间画出来地圈。
X轴和Y轴百分比：在PCA中，每个主成分解释的方差贡献率是一个重要的指标，用于评估该主成分对原始数据中方差的解释能力。具体来说，X轴和Y轴百分比表示对应的主成分解释的方差贡献率在总方差中的比例。例如，如果X轴和Y轴百分比分别为30％和20％，则表示X轴和Y轴解释了总方差的30％和20％，而其余主成分解释了剩余的50％。因此，X轴和Y轴的百分比可以帮助我们确定保留多少主成分以达到一定的降维效果。横纵坐标后面的百分比在其它领域可能有要求，比如主成分相加要大于 90% 之类的，不过生信分析中一般不怎么关注这个,因为变量（基因）太多了。

值得注意的是，PCA主成分分析是一种无监督学习方法，因此，X轴和Y轴的标签通常表示主成分的编号，而不是代表任何具体的特征或变量。

输入：转置后表达矩阵（数据框），行为样本，列为基因

实战(Step 3.4)

先得到转置后的表达矩阵，然后进行 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
save(pca_plot,file = "pca_plot.Rdata")

四、差异分析

1. limma 差异分析

limma

Limma 是一个用于分析由微阵列芯片或 RNA-seq 技术产生的基因表达数据的软件包。limma的算法原理基于线性模型和贝叶斯方法。它采用线性模型来描述基因表达量数据中的差异，并使用贝叶斯方法来估计模型参数，如样本间差异和基因间方差。它还利用经验贝叶斯方法来控制假阳性率（false discovery rate，FDR）

差异分析结果解释

Foldchange(FC): 处理组表达水平平均值/对照组平均值
logFC: 处理组 log 后表达量的平均值 - 对照组 log 后表达量的平均值。因为 log(x/y)=log(x)-log(y)
logFC 的合理取值范围：个位数。比如说，logFC = 10 ,代表基因上调表达了 2^10^ 也就是 1024 倍。
logFC >0 , 代表基因表达量上升；logFC <0 , 代表基因表达量下降；
通常所说的上下调基因是指表达量显著上调/下调的基因，显著性可根据 p 值判断
常见logFC取值: 1 2 1.2 1.5 2.2 0.585（没有唯一标准）

实战(Step 4.1)

#载入数据
rm(list = ls()) 
load(file = "step2output.Rdata")
#使用limma进行差异分析
library(limma)
design =model.matrix(~Group) # 构建分组矩阵
fit=lmFit(exp,design) # 线性拟合
fit=eBayes(fit) # 贝叶斯检验
deg=topTable(fit,coef=2,number = Inf) # 提取结果

limma 分析的结果包含 6 列：

logFC：差异倍数（Fold Change）的对数值，是指实验组与对照组基因表达量的比值。正数表示上调表达，负数表示下调表达，0表示没有差异表达。
AveExpr：基因在所有样本中的平均表达水平，通常是对数转换后的表达量平均值。
t：差异检验统计量，用于衡量两组样本基因表达量差异的显著性。t绝对值越大表示差异越显著。
P.Value：差异检验的显著性水平，即基于t统计量计算的P值。P值越小表示差异越显著。
adj.P.Val：校正后的P值，通常使用Benjamini-Hochberg方法进行多重比较校正，控制False Discovery Rate (FDR)。adj.P.Val越小表示差异越显著。
B：Bayes Factor值，可以用于衡量基因的差异表达程度。B值越大表示差异越显著。

2. 补充差异分析表

由于差异分析得到的表只包含探针id号，为了方便后面的分析，我们给差异分析的表添加几列数据，探针 id、基因 symbol、基因差异类型（上调、下调或者不显著）、以及 ENTREZID。其中 ENTREZID 后面富集分析时会用到。

探针 id 列: 表达矩阵的行名即探针 id，只需将其转换成一列数据即可。
基因 symbol 列: 即基因名称。在我们前面 Step2 获得的探针注释信息中, ids 变量中储存了探针 id 与基因 symbol的对应关系。但是需要注意的是多条探针可能会对应1个相同的基因，因为芯片厂商可能会针对某些基因设计多条探针。我们只需要保留其中 1 条探针的差异分析结果即可，所以过滤掉具有重复基因的探针得到的差异表达结果然后再添加基因 symbol列。

去重的方法：随机去重、保留行和或者行平均值最大的探针、取多个探针的平均值。没有唯一标准。

基因差异类型（change）列：我们设置好阈值（logFC 与 p-value）后，将所有基因分为上调up、下调down、或者不显著stable，然后将其作为 1 列数据添加。
ENTREZID 列：使用 clusterProfiler 包的 bitr 函数在相应数据库中（如人类 org.Hs.eg.db）获取基因 symbol 与 ENTREZID 的对应关系。symbol 与 ENTREZID 并不是一一对应的因此，得到的对应关系 symbol数可能会增多一点，也可能会较少一点（比如提示：3.26% of input gene IDs are fail to map...），都是正常的。

对应物种数据库的包名称获取方法：bioconductor链接http://bioconductor.org/packages/release/BiocViews.html


非模式生物没有 OrgDb 需要自己构建：简书链接`https://www.jianshu.com/p/9c9e97167377`

实战(Step 4.2)

# 1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))  

#2.加上探针注释中的symbol
ids = ids[!duplicated(ids$symbol),]  # 去除了ids 中重复的探针
# ids = distinct(ids,symbol,.keep_all = T)  # 与上面语句等价
deg <- inner_join(deg,ids,by="probe_id") # 添加symbol列，同时过滤重复的探针的差异分析结果
nrow(deg)

#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)

#4.加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,  # 根据 deg 的symbol，在相应的数据库中找到对应的ENTREZID
            fromType = "SYMBOL",
            toType = "ENTREZID",
            OrgDb = org.Hs.eg.db) #人类
#其他物种http://bioconductor.org/packages/release/BiocViews.html
deg_plot <- deg  # 因为我不希望后面绘图丢失几个symbol所以这里弄了一个deg_plot用于差异基因绘图 
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL")) # 添加ENTREZID列
save(Group,deg,deg_plot,logFC_t,p_t,gse_number,file = "step4output.Rdata")

最后得到的 deg 表格：

五、差异分析可视化

1. 火山图

火山图（volcano plot）通常用于比较两组样本或实验条件的数据，其中横坐标表示差异度（logFC），纵坐标表示显著性（负对数P值，-log10(P-value)）。每个数据点代表一个基因、蛋白质或其他生物分子，在火山图中呈现为一个点。点的位置表示该基因、蛋白质或分子的显著性和差异度。为什么要用负对数P值？因为 p 值本身很小，直接作为 y 轴画在图上不好看。因此，火山图中横坐标增加 1 ，代表 FoldChange 增加了 2 倍，纵坐标增加 1 ，代表 p 值减小了 10 倍。

一般来说，显著性越高，P值越小，负对数P值越大，数据点越靠近图像顶部；差异度越高，Fold Change越大，数据点越靠近图像的左右两侧。在火山图中，通常会设置一个阈值（如 Fold Change 值或 P 值），以区分显著差异的点和非显著差异的点。一些火山图会进一步使用颜色编码来显示差异的方向，例如使用红色表示高表达，蓝色表示低表达。

输入：差异分析得到的表格（需要 p值、logFC、以及基因差异表达的 change 三列数据）

实战(Step 5.1)

# 载入数据
rm(list = ls()) 
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")

# 画火山图
library(dplyr)
library(ggplot2)
dat <- deg_plot
# dat  = distinct(deg,symbol,.keep_all = T) 如果用deg作图需要对symbol去重
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","ETV5"))

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. 热图(Step 5.2)

根据需求选取一定数量或全部的差异基因，然后用 pheatmap 绘制热图。

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)

#差异基因热图
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

3. 相关性热图(Step 5.3)

使用 cor 函数计算基因之间的相关性，需要注意的是原始的表达矩阵是列为样本，直接用 cor 函数实际上计算的是样本之间的相关性，因此需要转置后才能计算基因之间的相关性。

输入：相关性矩阵，用 pheatmap 或者 corrplot 画图

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

M = cor(t(exp[g,]))  # cor计算列之间的相关性，所以要转置
# pheatmap 绘图
pheatmap(M)
# 用 corrplot 绘图
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)

六、富集分析

1. GO 与 KEGG、富集分析

GO 与 KEGG 数据库简介

GO（Gene Ontology）数据库www.geneontoloy.org是一个关系数据库，包括GO 本体论、基因和基因产物在这些本体论术语中的注释信息。旨在建立一套适用于各种物种的，对基因和蛋白质功能从多个方面进行限定和描述的，并能随着研究不断深入而更新的语义(terms)词汇标准，即基因产物分类标准。注释系统中每一个结点(term)都是基因或蛋白功能的一种命名及描述。

分子功能(Molecular Function, MF)描述在个体分子生物学上的活性，如催化活性或结合活性。
生物学过程(Biological Process,BP)由分子功能有序地组成的，具有多个步骤的一个过程，如细胞周期。
细胞组件(Cellular Component, CC)指基因产物位于何种细胞器或基因产物组中(如糙面内质网，核糖体，蛋白酶体等)，即基因产物在什么地方起作用。

KEGGwww.genome.jp/kegg（Kyoto Encyclopedia of Genes and Genomes）是系统地分析基因功能、链接基因组信息和功能信息的数据库,旨在揭示生命现象的遗传与化学蓝图。

富集分析

通过将差异基因做 GO 富集分析，可以把基因按照不同的功能进行归类，达到对基因进行注释和分类的目的。富集分析可以这么理解：衡量每个通路里的基因在差异基因里是否足够多。

富集的思路：可以选取全部差异基因富集，也可以上下调基因分开富集，或者选择自己感兴趣的基因进行富集。

2. 实战(Step 6.1)

利用 clusterProfiler 包进行 GO 富集分析以及可视化。

输入：输入一组我们感兴趣的差异基因的 ENTREZID

GO 富集分析：enrichGO 函数富集分析，barplot 函数与 dotplot函数画条带图与气泡图

# 载入数据与包
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)GO富集
ego <- enrichGO(gene = gene_diff, # 差异基因的ENTREZID
                OrgDb= org.Hs.eg.db, # 人类对应的数据库
                ont = "ALL", # GO的三个部分都要
                readable = TRUE) # 这个参数可以把结果转换回symbol
#ont参数：One of "BP", "MF", and "CC" subontologies, or "ALL" for all three.

#(3)可视化
#条带图
barplot(ego)
barplot(ego, split = "ONTOLOGY", font.size = 10,  # 按照GO的三个组分分开画图
        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") 
#(4)拼图
a <- barplot(ego)
b <- dotplot(ego) 
library(patchwork)
a+b

G0 富集结果含义

ego@result

ONTOLOGY : GO 的 ONTOLOGY 类别（BP、MF、CC）
ID: 对应于富集分析中找到的功能/通路的唯一标识符。
Description: 功能/通路的描述，通常是文本或注释信息。
GeneRatio: 在功能/通路中富集的差异基因数除以被收录的差异基因总数。比如 20/208 ，20 代表差异基因中有 20 个属于该通路，208 代表差异基因中有 208 个被对应的数据库所收录。
BgRatio: 在功能/通路中的基因总数数除以背景基因集中基因总数。比如 91/8217 ，91 代表该通路中共有 91 个基因，8217 代表对应的数据库中共有8217个基因。
pvalue: 统计分析中得到的原始p值。
p.adjust: 经过多重检验校正（例如，Bonferroni校正或Benjamini-Hochberg校正）后的p值。clusterProfiler包任何可视化都是基于 p.adjust 做的。
qvalue: 经过多重检验校正后的假发现率(FDR)。
geneID: 在相应功能/通路中富集的基因ID列表，以斜杆分隔。
Count: 在相应功能/通路中富集的基因数。

可视化结果

横向条形图

x轴：Count（富集的基因数）
y轴：Description（功能/通路）
气泡颜色：对应通路的显著性水平，按照 p.adjust 分配，从蓝到红显著性逐渐升高

气泡图

x轴：GeneRatio
y轴：Description（功能/通路）
气泡大小：Count（富集的基因数多少）
气泡颜色：对应通路的显著性水平，按照 p.adjust 分配，从蓝到红显著性逐渐升高

其它图形

展示top通路的共同基因。

#gl 用于设置下图的颜色
gl = deg$logFC
names(gl)=deg$ENTREZID
#Gene-Concept Network,要放大看
cnetplot(ego,
         # layout = "star",
         color.params = list(foldChange = gl),
         showCategory = 3)

展示 Go term 之间的关系

有颜色的 term 代表实验中富集到的，灰色的 term 代表挖掘到的与之相关的一些 term

ego_BP <- enrichGO(gene = gene_diff,
                     OrgDb= org.Hs.eg.db,
                     ont = "BP",  # 只选择GO的BP
                     readable = TRUE)
goplot(ego_BP)
# 图的理解可参考 zhuanlan.zhihu.com/p/99789859

3. 实战(Step 6.2)

利用 clusterProfiler 包进行 KEGG 富集分析以及可视化。

输入：输入一组我们感兴趣的差异基因的 ENTREZID

KEGG 富集分析：enrichKEGG 函数富集分析，画图函数与 GO 富集分析相同

#（1）输入数据
gene_up = deg[deg$change == 'up','ENTREZID'] 
gene_down = deg[deg$change == 'down','ENTREZID'] 
gene_diff = c(gene_up,gene_down)

#（2）对上调/下调/所有差异基因进行富集分析
#富集上调基因
kk.up <- enrichKEGG(gene = gene_up,
                    organism = 'hsa') # 指定物种名称，人类为hsa
#富集下调基因
kk.down <- enrichKEGG(gene =  gene_down,
                      organism = 'hsa')
#富集所有基因
kk.diff <- enrichKEGG(gene = gene_diff,
                      organism = 'hsa')

#(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)

结果

按照矫正 p 值小于0.05来看，上调基因没有富集到。

> table(kk.diff@result$p.adjust<0.05)
FALSE  TRUE 
  289    14 
> table(kk.up@result$p.adjust<0.05) 
FALSE 
  282 
> table(kk.down@result$p.adjust<0.05)
FALSE  TRUE 
  211    19

kegg 富集分析结果各指标的含义与 GO 富集结果基本一致，也可以画点图、横向条形图，画图函数与前面GO富集的一致：

KEGG富集时遇到的问题

一开始我以为是没富集到，但是后来有一次重新跑的时候又富集到了，所以如果富集后 kk 直接返回 NULL 的话可能和网络有关系，因为一般来说 kk 都会返回一个列表的，列表中的数据 p 值不显著才算没有富集到。

4. 备注

没有富集到怎么办？

clusterProfiler 是根据 p.adjust 进行富集的，可以试一下用 p 值进行筛选
换一种富集的方式，如 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

七、tinyarray简化分析

1. 双分组数据

前面我们经过六步对芯片数据进行了分析，使用小洁老师的 tinyarray 包可以简化上述流程，一步到位进行芯片数据分析。还是用这个数据集 GSE42872 用 tinyarray 包进行简化分析。我们自己只要完成前面两步得到表达矩阵、分组信息、探针注释信息，就可以直接一步获得差异分析可视化结果了。

# 下载数据
# geo_download()函数直接从国内镜像下载表达矩阵，得到exp，pd,gpl
library(tinyarray)
library(stringr)
gse = "GSE42872"
geo = geo_download(gse)
boxplot(geo$exp) # 发现不需要 log 转化

# 获取分组信息和探针注释
Group=ifelse(str_detect(geo$pd$title,"Control"),
             "control",
             "treat")
Group = factor(Group,levels = c("control","treat"))
find_anno(geo$gpl)
ids <- AnnoProbe::idmap('GPL6244')
head(ids)

#差异分析和它的可视化 get_deg_all函数一步到位
dcp = get_deg_all(geo$exp,Group,ids)
table(dcp$deg$change)
head(dcp$deg)
dcp$plots
library(ggplot2)
ggsave("deg.png",width = 15,height = 5)

#富集分析
genes = dcp$deg$ENTREZID[dcp$deg$change!="stable"]
head(genes)
#有可能因为网络问题报错
g = quick_enrich(genes,destdir = tempdir())
names(g)
g[[1]][1:4,1:4]
library(patchwork)
g[[3]]+g[[4]]
ggsave("enrich.png",width = 12,height = 7)

结果

差异分析可视化结果：热图、PCA图、火山图

富集分析可视化结果，左图为 KEGG 富集气泡图，右图为 GO 富集气泡图。

2. 多分组数据

对于多分组数据，如 1 个对照组、两个以上实验组，或者多个实验组的，可以将原始数据输入用 limma 分析，然后提取两两比较的差异分析结果，再可视化。也可以用 tinyarray 包一步完成分析。以数据集 GSE474 为例：

### 1.获取数据
rm(list = ls())
library(stringr)
library(tinyarray)
gse = "GSE474"
geo = geo_download(gse)
geo$exp = log2(geo$exp+1)

### 2.生成分组向量与探针注释
#注意仍然是对照组放在levels的第一个位置
#三分组一般是两两比较，后面差异分析时看清楚谁比谁就可以了。
Group=str_split(geo$pd$title,"-",simplify = T)[,3]
Group=factor(Group,levels = c("NonObese","Obese","MObese"))
geo$gpl
if(!require(hgu133a.db))BiocManager::install("hgu133a.db")
library(hgu133a.db)
ids <- toTable(hgu133aSYMBOL)
head(ids)

### 3.tinyarray的简化操作
#多分组的数据，get_deg_all仍然可以帮你简化操作，
#目前是三分组就两两差异分析，四个或五个分组的数据是后面几个组与第一个组差异分析，
#暂不支持其他的做法和更多的分组。
dcp = get_deg_all(geo$exp,Group,ids,symmetry = T,logFC_cutoff = 0.585,cluster_cols = F,entriz = F)
dcp$plots
ggplot2::ggsave("deg.png",width = 15,height = 10)

#富集分析的输入数据是差异基因名字，只要你提供一组基因名就能做，
#要分别做三次还是取交集或合集一起做都可以，没有标准答案。

结果

得到三个组的差异基因热图、PCA图、两两比较的差异基因韦恩图、两两比较的火山图

写在后面

GEO数据库里面的表达量芯片数据处理，主要的难点是表达量矩阵获取和探针的基因名字转换，搞定后只需要一定的生物学背景对数据进行合理的分组后就是标准的差异分析，富集分析。主要是参考我八年前的笔记：

因为都是标准的代码，所以每次有学徒和实习生我都会让大家两次十几个数据集，凑成为了一个合辑：《1000个基因芯片表达量矩阵数据处理》：

其实本文已经是一套万能代码，它理论上可以支持GEO数据库的至少5万个表达量芯片数据集，从下载表达量矩阵到后续差异分富集分析一条龙，而且输出大量图表和一个网页报告！但是它并不是傻瓜式的，仍然需要你会R语言，需要生物学背景去修改分组形式，需要人为判断芯片的探针对应基因的关系，其它的图表，比如火山图，热图，GO和KEGG数据库富集图，GSEA图就是自动化的啦。已经是目前我们能想到的最小化干预了。

给大家一个作业本：以GSE16515为例

GEO链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16515
芯片平台：GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
平台链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
样品信息：16个正常样本与36个胰腺导管腺癌（PDAC）样本
文章标题及链接：FKBP51 Affects Cancer Cell Response to Chemotherapy by Negatively Regulating，Akt.Cancer Cell. 2009 Sep 8; 16(3) https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2755578/

需要你完成下面的3个分析图：

首先是分组后查看是否合理：

然后是简单的差异分析：

简单的差异分析

最后是简单的数据库注释：

每一个图表都有背后的统计学原理，也有各自美化的代码，但是都不在我们的万能代码里面哦。

ixxmu / mp_duty

GEO数据挖掘实战全流程1.27万字笔记 #3477

GEO数据挖掘实战全流程1.27万字笔记 by 生信技能树

一、数据下载与读取

1. 使用R包 GEOquery 下载

2. geoChina()

3. NCBI 直接下载

4. 读取芯片数据的方法

5. 实战(Step 0)

6. 实战(Step 1)

二、分组信息与探针注释

1. 获取分组信息

2. 获取芯片注释的方法

(1) 捷径 tinyarray

(2) bioconductor 包

(3) NCBI 下载文件解析

(4) 其它方法

3. 实战(Step 2)

三、表达矩阵数据检查

1. 基础检查 tips

log 转化

箱线图查看异常样本

异常样本处理办法

表达矩阵里有负值

检查某些基因的表达量

2. 绘图

(1) 箱线图

实战(Step 3.1)

(2) 层次聚类

实战(Step 3.2)

(3) 热图

实战(Step 3.3)

(4) PCA

实战(Step 3.4)

四、差异分析

1. limma 差异分析

limma

差异分析结果解释

实战(Step 4.1)

2. 补充差异分析表

实战(Step 4.2)

五、差异分析可视化

1. 火山图

实战(Step 5.1)

2. 热图(Step 5.2)

3. 相关性热图(Step 5.3)

六、富集分析

1. GO 与 KEGG、富集分析

GO 与 KEGG 数据库简介

富集分析

2. 实战(Step 6.1)

G0 富集结果含义

可视化结果

其它图形

3. 实战(Step 6.2)

结果

KEGG富集时遇到的问题

4. 备注

没有富集到怎么办？

画其它图

七、tinyarray简化分析

1. 双分组数据

结果

2. 多分组数据

结果

写在后面

给大家一个作业本：以GSE16515为例