GEO数据挖掘-基于芯片
原创
GEO数据挖掘-基于芯片
1 00_pre_install.R
1.1 代码
options("repos"="https://mirrors.ustc.edu.cn/CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
cran_packages <- c('tidyr',
'tibble',
'dplyr',
'stringr',
'ggplot2',
'ggpubr',
'factoextra',
'FactoMineR',
'devtools',
'cowplot',
'patchwork',
'basetheme',
'paletteer',
'AnnoProbe',
'ggthemes',
'VennDiagram',
'survminer',
"tinyarray")
Biocductor_packages <- c('GEOquery',
'GO.db',
'hgu133plus2.db',
'ggnewscale',
"limma",
"impute",
"GSEABase",
"GSVA",
"clusterProfiler",
"org.Hs.eg.db",
"preprocessCore",
"enrichplot")
for (pkg in cran_packages){
if (! require(pkg,character.only=T,quietly = T) ) {
install.packages(pkg,ask = F,update = F)
require(pkg,character.only=T)
}
}
for (pkg in Biocductor_packages){
if (! require(pkg,character.only=T,quietly = 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)
}
#没有任何提示就是成功了，如果有warning xx包不存在，用library检查一下。
#library报错，就单独安装。在require()函数中，如果直接传递包的名称作为参数，不需要加引号；如果包的名称以字符串形式存储在变量中，则需要使用character.only = TRUE来指定这个变量是一个字符串
1.2 解析
1.2.1 require(pkg,character.only=T,quietly = T)
直接传递包名称（不加引号）
require(ggplot2) # 加载ggplot2包包名称存储在字符串变量中（需要加引号并使用character.only = TRUE）
package_name <- "ggplot2"
require(package_name, character.only = TRUE) # 加载ggplot2包为什么不加引号
当你直接传递包的名称时，R会把它视为一个标识符，而不是一个字符串。例如，require(ggplot2)等同于告诉R直接加载名为ggplot2的包。
为什么需要character.only = TRUE
当包名称存储在一个变量中时，比如package_name <- "ggplot2"，变量package_name包含的是一个字符串。因此，你需要告诉require()函数这是一个字符串，并且需要解释成包的名称。通过设置character.only = TRUE，require()函数会正确地将字符串变量解释为包的名称。
require()函数中的quiet参数用于控制加载包时的消息输出：
quiet = FALSE（默认值）：输出加载包的消息。quiet = TRUE：抑制加载包的消息，保持输出简洁。
2 01_start_GEO.R
2.1 代码
rm(list = ls())
#打破下载时间的限制,改前60秒，改后10w秒
options(timeout = 100000)
options(scipen = 20)#不要以科学计数法表示
#传统下载方式
library(GEOquery)
eSet = getGEO("GSE7305", destdir = '.', getGPL = F)
#网速太慢，下不下来怎么办
#1.从网页上下载/发链接让别人帮忙下，放在工作目录里
#2.试试geoChina,只能下载2019年前的表达芯片数据
#library(AnnoProbe)
#eSet = geoChina("GSE7305") #选择性代替第8行
#研究一下这个eSet
class(eSet)
length(eSet)
eSet = eSet[[1]]
class(eSet)
#(1)提取表达矩阵exp
exp <- exprs(eSet)
#⭐第一个要检查的地方👇，表达矩阵行列数，正常是几万行，列数=样本数，
#如果0行说明不是表达芯片或者是遇到特殊情况，不能用此流程分析
dim(exp)
#⭐二个要检查的地方👇
range(exp)#看数据范围决定是否需要log，是否有负值，异常值，如有负值，结合箱线图进一步判断
#⭐可能要修改的地方👇
exp = log2(exp+1) #需要log才log，不需要log要注释掉这一句
#⭐第三个要检查的地方👇
boxplot(exp,las = 2) #看是否有异常样本
#(2)提取临床信息
pd <- pData(eSet)
#⭐多分组中提取两分组的代码示例，二分组不需要
if(F){
#因为现在这个例子不是多分组，所以编造一列做示例。
pd$fake = paste0(rep(c("a","b","c","d"),each = 5),1:5)
k1 = str_detect(pd$fake,"b");table(k1)
k2 = str_detect(pd$fake,"c");table(k2)
pd = pd[k1|k2,]
}
#(3)让exp列名与pd的行名顺序完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) {
s = intersect(rownames(pd),colnames(exp))
exp = exp[,s]
pd = pd[s,]
}
#(4)提取芯片平台编号，后面要根据它来找探针注释
gpl_number <- eSet@annotation;gpl_number
save(pd,exp,gpl_number,file = "step1output.Rdata")
# 原始数据处理的代码，按需学习
# https://mp.weixin.qq.com/s/0g8XkhXM3PndtPd-BUiVgw2.2 解析
2.2.1 rm(list=ls())
在单个R脚本中，开头都应该是rm(list = ls()),然后需要什么数据再重新加载，养成良好的代码习惯，可以确保结果的可重复性，减少问题的产生。
2.2.2 options(timeout = 100000)
在R语言中，options(timeout = 100000) 是用来设置全局选项的，特别是调整网络连接和下载操作的超时时间。timeout选项控制的是当R进行网络操作（如下载文件或访问网络资源）时等待响应的最长时间（以秒为单位）。
默认情况下，R的timeout值可能设置得较低（如60秒），这意味着如果网络操作在该时间内未完成，R会抛出一个超时错误。通过设置一个较大的timeout值，可以避免网络操作因超时而失败。
2.2.3 options(scipen = 20)
scipen 是 scientific penalty 的缩写。它是一个数值，用于影响R在打印数值时选择是否使用科学计数法的倾向。scipen 的值越大，R越倾向于使用普通的定点数表示法而不是科学计数法。反之，scipen 的值越小（或为负值），R越倾向于使用科学计数法表示数值。
2.2.4 getGEO("GSE7305", destdir = '.', getGPL = F)
getGEO()函数是Bioconductor包GEOquery中的一个函数，用于从Gene Expression Omnibus (GEO)数据库下载GEO数据集。
GSE7305：
- 这是GEO数据集的访问编号（GEO Series accession number），指定了你要下载的数据集。在这个例子中，你下载的是编号为
GSE7305的数据集。
destdir = '.'：
- 这个参数指定下载文件的保存目录。
.表示当前工作目录。你可以将其更改为任何你希望保存文件的目录路径。
getGPL = FALSE：
- 这个参数决定是否下载平台注释文件（GEO Platform file）。如果设置为
FALSE（如示例中），平台注释文件将不会被下载。如果设置为TRUE，则会下载这些文件。平台注释文件包含关于实验所用平台的信息，如芯片上的探针序列等。
2.2.5 eSet = eSet[1] ；class(eSet)；
由于getGEO()返回的eSet是一个包含一个或多个ExpressionSet对象的列表，所以你需要提取列表中的第一个元素，即eSet[[1]]，并将其赋值回eSet。查看class(eSet)可以确认它现在是一个ExpressionSet对象。
补充知识：<font color=red size=5>ExpressionSet</font>
ExpressionSet对象是Bioconductor框架中的一个核心类，用于存储高通量基因表达数据及其相关的元数据。它主要用于微阵列和RNA-Seq数据分析。ExpressionSet对象整合了表达矩阵、样本信息和特征信息，提供了一个一致的数据结构，使得后续的数据分析和可视化更加方便和一致。
主要组成部分
一个典型的ExpressionSet对象包含以下几个主要组成部分：
- 表达矩阵（Expression Matrix）：
- 存储基因表达数据的矩阵。行通常表示基因（探针、特征），列表示样本。矩阵中的每个元素表示某个基因在某个样本中的表达量。
- 可以通过
exprs()函数提取。
exp <- exprs(eSet)
- 样本元数据（Sample Metadata）：
- 描述样本的元数据（例如，样本的分组信息、处理条件等），存储在
phenoData中。 - 可以通过
pData()函数提取。
sample_info <- pData(eSet)
- 特征元数据（Feature Metadata）：
- 描述基因或探针的元数据（例如，基因的注释信息、探针的序列等），存储在
featureData中。 - 可以通过
fData()函数提取。
feature_info <- fData(eSet)
- 实验描述（Experiment Description）：
- 描述实验的相关信息，如实验名称、实验设计等，存储在
experimentData中。 - 可以通过
experimentData()函数提取。
experiment_info <- experimentData(eSet)
2.2.6 p = identical(rownames(pd),colnames(exp));p
p = identical(rownames(pd),colnames(exp));p
identical() 是 R 语言中的一个函数，用于比较两个对象是否完全相同。与 == 不同，identical() 比较对象的内容和属性，确保两个对象在所有方面都完全相同。即identical() 用于比较表达矩阵（exp）的列名和临床信息数据框（pd）的行名，以确保它们完全一致。
如果p为false,执行
if(!p) {
s = intersect(rownames(pd),colnames(exp))
exp = exp[,s]
pd = pd[s,]
}- 取行名和列名的交集：
s = intersect(rownames(pd), colnames(exp))- 这行代码取临床信息数据框
pd的行名和表达矩阵exp的列名的交集。交集s包含了同时出现在pd和exp中的样本名称。
- 根据交集重新排序表达矩阵和临床信息数据框：
exp = exp[, s]- 重新排列表达矩阵
exp的列，使其顺序与交集s中的样本顺序一致。 pd = pd[s, ]- 重新排列临床信息数据框
pd的行，使其顺序与交集s中的样本顺序一致。
这样做的目的是确保在后续分析中，每个样本的表达数据和临床信息能够正确对应。
2.2.7 gpl_number <- eSet@annotation;gpl_number
eSet[[1]]@annotation这种语法用于访问S4类对象的槽（slot）。在R语言中，ExpressionSet对象是S4类对象，S4类对象的槽通过@操作符来访问。下面是详细的解释。
S4类和槽（Slot）:S4类是R中一种更严格和复杂的类定义方式，适用于需要更严格数据结构的情况。S4类对象包含一个或多个槽，每个槽存储特定类型的数据。
ExpressionSet对象:Bioconductor框架中用于存储高通量基因表达数据及其元数据的S4类对象。这个对象有多个槽，例如assayData、phenoData、featureData、experimentData和annotation`。
annotation槽:存储芯片平台编号（例如GPL编号），用于指定该数据集使用的微阵列或测序平台。
3 02_group_ids.R
3.1 代码
# Group(实验分组)和ids(探针注释)
rm(list = ls())
load(file = "step1output.Rdata")
# 1.Group----
library(stringr)
# 标准流程代码是二分组，多分组数据的分析后面另讲
#⭐要修改的地方：分组信息,必须学会ifelse和str_detect
k = str_detect(pd$title,"Normal");table(k) #不在title就在pd的其他列
Group = ifelse(k,"Normal","Disease")
# 需要把Group转换成因子，并设置参考水平，指定levels
#⭐要修改的地方，对照组在前，处理组在后
Group = factor(Group,levels = c("Normal","Disease"))
Group
#⭐检查自己得到的分组是否正确
data.frame(pd$title,Group)
#2.探针注释的获取-----------------
#四种方法，方法1里找不到就从方法2找，以此类推。
#方法1 BioconductorR包(最常用)
#⭐要操作的地方
library(tinyarray)
gpl_number #首先看看编号是多少
#View(pkg_all)
#然后在pkg_all里搜索gpl编号，找到对应的R包前缀(第二列)，没搜到就是没有R包,再看方法2。
#也可以用代码直接得到对应的R包前缀：
pkg_all[pkg_all$gpl==gpl_number,2]
#⭐要操作的地方
#用上面找到的前缀替换下面所有的hgu133plus2，共5处
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db",ask = F,update = F)
library(hgu133plus2.db)
ls("package:hgu133plus2.db") #列出R包里都有啥
ids <- toTable(hgu133plus2SYMBOL) #把R包里的注释表格变成数据框
# 方法2 下载并读取GPL网页的表格文件，按列取子集
#⭐要操作的地方
library(tinyarray)
get_gpl_txt(gpl_number) #获取表格文件的下载链接
# 接下来是复制网址去浏览器下载、放在工作目录下、读取、提取探针id和基因symbol(没有现成的需要拆分和转换)，不同文件代码不统一，等看同学们的例子。
# 注意:最终的数据ids只能有两列，第一列列名是probe_id,第二列列名是symbol,且都是字符型，否则后面代码要报错咯。
# 方法3 官网下载注释文件并读取
# 方法4 自主注释，了解一下
#https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
save(exp,Group,ids,file = "step2output.Rdata")
#比较复杂的探针注释参考资料
#资料1：拆分取列https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5
#资料2：多种id的转换https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/pn0s1mmsaxocfynb?singleDoc# 《又一个有点难的探针注释（多种id的转换）》
#资料3：其他id转换 https://www.jianshu.com/p/f4e799f06b524 03_pca_heatmap.R
4.1 代码
rm(list = ls())
#⭐不用改代码，会看结果就行
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))
library(FactoMineR)
library(factoextra)
dat.pca <- PCA(dat, graph = FALSE)
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"
)
#如果少于4个点就不会画出圈，是正常的。因为圈是置信区间，样本太少无法计算，不是必须的。
# 2.top 1000 sd 热图----
g = names(tail(sort(apply(exp,1,sd)),1000)) #day7-apply的思考题
n = exp[g,]
library(pheatmap)
annotation_col = data.frame(row.names = colnames(n),
Group = Group)
pheatmap(n,
show_colnames =F,
show_rownames = F,
annotation_col=annotation_col,
scale = "row", #按行标准化，只保留行内差别，不保留行间差别，会把数据范围缩放到大概-5~5之间
breaks = seq(-3,3,length.out = 100) #设置色带分布范围为-3~3之间，超出此范围的数字显示极限颜色
)
# 关于scale的进一步学习：zz.scale.R4.2 解析
4.2.1 dat = as.data.frame(t(exp))
将表达矩阵 exp 转置后转换为数据框。
在基因表达数据分析中，表达矩阵 exp 通常是一个二维矩阵，其中：
- 行代表基因。
- 列代表样本。
为了进行主成分分析（PCA）等分析，需要将矩阵转置，以便样本成为行，基因成为列。
4.2.2 fviz_pca_ind(...)
以下是 fviz_pca_ind 函数的详细解释和代码示例，它用于绘制主成分分析（PCA）图，并按组别进行颜色区分。
fviz_pca_ind 函数参数说明
dat.pca：PCA分析的结果对象。geom.ind：表示样本点的几何形状，这里设置为 "point" 表示仅显示点。col.ind：指定样本点的颜色，这里根据Group进行颜色区分。palette：指定颜色调色板，这里使用了蓝色和黄色。addEllipses：是否添加浓度椭圆，这里设置为TRUE。legend.title：图例的标题，这里设置为 "Groups"。
4.2.3 g = names(tail(sort(apply(exp,1,sd)),1000))
apply(exp, 1, sd)：对表达矩阵 exp 的每一行（即每个基因）计算标准差。
sort()：将这些标准差按升序排序。
tail(..., 1000)：取出排序后的最后1000个值，即标准差最大的1000个基因（基因探针编号）。
names()：获取这些基因的名称（基因探针编号）。
4.2.4 pheatmap(...)
pheatmap(...)：使用 pheatmap 包绘制热图。
show_colnames = FALSE：不显示列名。
show_rownames = FALSE：不显示行名。
annotation_col = annotation_col：添加列注释，即样本的分组信息。
scale = "row"：按行标准化，使每行数据的均值为0，标准差为1。
breaks = seq(-3, 3, length.out = 100)：设置颜色梯度的分布范围为 -3 到 3。
<font color=red>为什么选择标准差最大的1000个基因并绘制热图？</font>
- 识别差异：标准差最大的基因通常是表达变化最大的基因，这些基因更有可能在不同的样本或组别之间显示出显著的差异。
- 降低数据维度：基因表达数据通常非常高维，选择1000个基因可以降低数据的维度，使得可视化和分析更为可行和清晰。
5 04_DEG_R
5.1 代码
rm(list = ls())
# 加载之前保存的数据，包括 exp（表达矩阵）、Group（样本分组）和 ids（探针注释）。
load(file = "step2output.Rdata")
#差异分析
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 = distinct(ids,symbol,.keep_all = T)
#其他去重方式在zz.去重方式.R
deg = inner_join(deg,ids,by="probe_id")
#⭐检查
nrow(deg) #如果行数为0就是你找的探针注释是错的。
#3.加change列,标记上下调基因
#⭐阈值，可按需修改
logFC_t = 1
p_t = 0.05
#⭐思考，如何使用padj而非p值
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)
#火山图
library(ggplot2)
ggplot(data = deg, 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()
# 差异基因热图----
# 表达矩阵行名替换为基因名
exp = exp[deg$probe_id,]
rownames(exp) = deg$symbol
diff_gene = deg$symbol[deg$change !="stable"]
n = exp[diff_gene,]
library(pheatmap)
annotation_col = data.frame(group = Group)
rownames(annotation_col) = colnames(n)
pheatmap(n,show_colnames =F,
show_rownames = F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col,
breaks = seq(-3,3,length.out = 100)
)
#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)#⭐Hs是人类,注意物种
#一部分基因没匹配上是正常的。只要物种正确，<30%的失败都没事。
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
#⭐检查行数，减少太多说明不正常
nrow(deg)
deg = inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
#多了一些行少了一些行都正常，SYMBOL与ENTREZID不是一对一的。
nrow(deg)
save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")5.2 解析
5.2.1 design = model.matrix(~Group)
创建设计矩阵，用于差异表达分析。在差异基因表达分析中，设计矩阵是一个非常重要的步骤。设计矩阵描述了实验设计和样本分组信息，为后续的线性模型拟合提供基础。
注：
- 因子变量 Group
Group 是一个因子变量，表示实验分组。它有两个水平："Normal" 和 "Disease"。前10个样本属于 "Disease" 组，后10个样本属于 "Normal" 组。
- 创建设计矩阵
model.matrix(~Group) 创建了一个包含分组信息的设计矩阵。对于20个样本，设计矩阵可能如下：
(Intercept) 列表示截距项，每个样本的值都为1。
GroupDisease 列表示 "Disease" 组样本。如果样本属于 "Disease" 组，值为1；如果属于 "Normal" 组，值为0。
5.2.2 deg = topTable(fit,coef = 2,number = Inf)
fit = lmFit(exp, design)：使用线性模型拟合表达数据。
fit = eBayes(fit)：使用贝叶斯方法计算统计量。
deg = topTable(fit, coef = 2, number = Inf)：提取所有基因的差异表达结果，coef = 2 表示第二个因子的系数（通常是对照组和处理组之间的比较）。
注：
topTable 函数用于从线性模型拟合的结果中提取基因表达的统计信息。
fit：这是前面用 lmFit 和 eBayes 函数得到的线性模型拟合结果。它包含了表达矩阵和设计矩阵的信息，以及通过贝叶斯方法计算的统计量。
topTable：这是 limma 包中的一个函数，用于提取差异表达分析的结果。
coef = 2：指定要提取的系数。在设计矩阵 design 中，每个因子（即实验组）都有一个对应的系数。coef = 2 表示我们要提取的是设计矩阵中第二个因子的系数（在这种情况下，通常是对照组与处理组的比较）。
number = Inf：指定要提取的基因数量。Inf 表示提取所有基因的结果。如果你只想提取前 n 个基因，可以将 Inf 替换为具体的数字，比如 100 表示提取前100个基因。
5.2.3 deg = mutate(deg,probe_id = rownames(deg))
使用 dplyr 包中的 mutate 函数为数据框 deg 添加一列 probe_id，该列的值为数据框 deg 的行名。
5.2.4 ids = distinct(ids,symbol,.keep_all = T)
使用 dplyr 包中的 distinct 函数，从数据框 ids 中移除重复的行，并保留每个 symbol 列唯一的行，同时保留所有其他列。
ids：要处理的数据框。
symbol：指定根据哪一列进行去重（这里是 symbol 列）。
.keep_all = TRUE：表示在去重时，保留所有列的数据。
注：如果不写 .keep_all = TRUE，distinct 函数只会返回去重列，并且不会保留去重列之外的其他列。具体来说，在默认情况下，distinct 函数只返回去重后的 symbol 列，不会保留 probe_id 等其他列的数据。
5.2.5 差异基因热图
- 过滤和重命名表达矩阵
exp = exp[deg$probe_id,]：将 exp 矩阵的行过滤为 deg 数据框中 probe_id 列对应的行。这一步确保表达矩阵 exp 只包含差异表达基因分析结果中的探针。
rownames(exp) = deg$symbol：将表达矩阵 exp 的行名设置为 deg 数据框中的 symbol 列。这一步将表达矩阵中的探针 ID 替换为对应的基因符号，使得矩阵更加易读。
- 提取差异基因
diff_gene = deg$symbol[deg$change != "stable"]：从 deg 数据框中提取非稳定状态（即有差异表达）的基因符号。
- 提取差异基因的表达数据
n = exp[diff_gene,]：从表达矩阵 exp 中提取差异基因的表达数据。
- 准备注释数据框
annotation_col = data.frame(group = Group)：创建一个包含分组信息的注释数据框。
rownames(annotation_col) = colnames(n)：将注释数据框 annotation_col 的行名设置为表达数据矩阵 n 的列名，确保注释信息与样本数据对齐。
- 绘制热图
library(pheatmap)：加载 pheatmap 包，用于绘制热图。
pheatmap(n, show_colnames = F, show_rownames = F, scale = "row", annotation_col = annotation_col, breaks = seq(-3, 3, length.out = 100))：绘制热图，参数解释如下：
show_colnames = F：不显示列名。
show_rownames = F：不显示行名。
scale = "row"`：按行标准化数据，使得每个基因的表达值在同一范围内进行比较。
annotation_col = annotation_col：使用注释数据框 annotation_col 添加列注释，标注样本的分组信息。breaks = seq(-3, 3, length.out = 100)：设置色带的分布范围为 -3 到 3，超出此范围的值显示为极限颜色。
6 05_anno.R
6.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)输入数据
gene_diff = deg$ENTREZID[deg$change != "stable"]
#(2)富集
#⭐下面的三句都要注意物种
ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
#其他物种https://www.genome.jp/kegg/catalog/org_list.html
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")
#如果ekk是空的，这句就会报错，因为没富集到任何通路。
ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,
ont = "ALL",readable = TRUE)
#setReadable和readable = TRUE都是把富集结果表格里的基因名称转为symbol
class(ekk)
#(3)可视化
#barplot可以换成dotplot
barplot(ego, split = "ONTOLOGY") +
facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")
barplot(ekk)
#如果ekk中没有padj<0.05的通路，就会报错，因为只画padj<0.05,没有参数
# 更多资料---
# GSEA：https://www.yuque.com/docs/share/a67a180f-dd2b-4f6f-96c2-68a4b86fe862?#
# Y叔的书：http://yulab-smu.top/clusterProfiler-book/index.html
# GOplot：https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew
# 网上的资料和宝藏无穷无尽，学好R语言慢慢发掘~6.2 解析
6.2.1 ekk <- enrichKEGG(gene = gene_diff,organism = 'hsa')
使用 clusterProfiler 包中的 enrichKEGG 函数对差异基因进行KEGG通路富集分析。
gene = gene_diff：提供要进行富集分析的基因列表。
organism = 'hsa'：指定分析的人类基因（hsa代表Homo sapiens）。
ekk <- setReadable(ekk,OrgDb = org.Hs.eg.db,keyType = "ENTREZID")使用 setReadable 函数将KEGG富集结果中的基因ID转换为更容易理解的基因符号。
6.2.2 ego <- enrichGO(gene = gene_diff,OrgDb= org.Hs.eg.db,ont = "ALL",readable = TRUE)
使用 clusterProfiler 包中的 enrichGO 函数对差异基因进行GO富集分析。
gene = gene_diff：提供要进行富集分析的基因列表。
OrgDb = org.Hs.eg.db：指定用于注释的基因数据库。
ont = "ALL"：指定进行所有GO分类（生物过程BP、分子功能MF、细胞组分CC）的富集分析。
readable = TRUE：将富集结果中的基因ID转换为基因符号。
6.2.3 barplot(ego,split...)
使用 barplot 函数绘制GO富集结果的柱状图。
split = "ONTOLOGY"：按GO分类（BP、MF、CC）进行分割。
facet_grid(ONTOLOGY ~ ., space = "free_y", scales = "free_y")：使用 ggplot2 包中的 facet_grid 函数将不同GO分类的结果分开显示。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
