KEGG富集分析

KEGG（Kyoto Encyclopedia of Genes and Genomes）是一个全球公认的生物信息学数据库，用于解释高级功能和生物系统功能的分子水平信息，特别关注基因组、基因功能和生物化学通路等方面。KEGG数据库提供了对生物学、生物化学和生物信息学领域的多种信息资源，包括基因和蛋白质的功能注释、代谢通路、信号通路、药物目标等信息。

library(ggplot2) # 绘图使用
library(clusterProfiler) #KEGG富集分析使用
library(org.Hs.eg.db) #转换基因ID使用
library(stats) #数据处理使用
library(data.table) #数据读取使用
library(dplyr) #数据处理使用
##差异基因筛选，这里选择abs(logFC) > 2,FDR < 0.05的值作为差异基因。也可以根据自己的情况进行设定，abs(logFC) < -2为下调基因富集。
DEG_data <- LUNG_Match_DEG %>% dplyr::filter(abs(logFC) > 2 & FDR < 0.05)
#Gene名转化为GeneID
gene.df <- bitr(DEG_data$gene_name, fromType = "SYMBOL", #fromType是指你的数据ID类型是属于哪一类的
toType = c("ENTREZID", "SYMBOL"), #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = org.Hs.eg.db) #Orgdb是指对应的注释包是哪个
#合并
colnames(gene.df)[1] <- "gene_name"
DEG_data1 <- left_join(gene.df,DEG_data)

差异基因的选择也可以单独选择上调或下调基因分别进行富集分析。KEGG富集分析时，需要将GeneID从SYMBOL转换成ENTREZID，由于数据是人类数据因此选择org.Hs.eg.db数据集。如果是研究其它物种，需改成对应物种的数据集。具体物种信息可以在http://bioconductor.org/packages/release/BiocViews.html#___OrgDb中查看。

## KEGG富集分析
kegg <- enrichKEGG(DEG_data1$ENTREZID, organism = 'human', pvalueCutoff = 0.05, 
                   pAdjustMethod = 'BH', minGSSize = 3, maxGSSize = 500, qvalueCutoff = 0.05, 
                   use_internal_data = FALSE)

kegg_result <- data.frame(kegg)

organism = 'human'：这是指定要进行KEGG富集分析的生物学物种。
pvalueCutoff = 0.05：用于过滤富集分析结果的原始p-value的阈值。
pAdjustMethod = 'BH'：指定用于校正多重假设检验的方法。在这里，使用的是Benjamini-Hochberg方法，通常用于控制假阳性率（False Discovery Rate，FDR）。
minGSSize = 3：富集分析中考虑的最小基因集合大小。
maxGSSize = 500：富集分析中考虑的最大基因集合大小。
qvalueCutoff = 0.05：用于过滤富集分析结果的校正后q-value的阈值。
use_internal_data = FALSE：是否使用内部数据集进行富集分析。如果设置为TRUE，将使用内部数据集，否则将从KEGG网站下载数据。
富集分析结果共有9列，这里解释一下各列代表的含义。
ID：KEGG通路或功能类别的唯一标识符。
Description：KEGG通路或功能类别的描述，即通路或功能的名称或注释。
GeneRatio：输入的差异基因中，映射到此KEGG通路的基因数目/映射到所有KEGG通路的基因总数目。
BgRatio：在整个基因数据库（或背景基因集合）中映射到KEGG通路的基因数目，以及该通路中的总基因数目。
pvalue、p.adjust和qvalue：p值、校正后p值和q值；
geneID：映射到该通路的基因的标识符列表，通常是ENTREZID。
Count：在富集分析中映射到该通路或功能类别的基因数目，即GeneRatio列中的第一个值。

KEGG结果绘制

自带绘图函数

## 条形图
barplot(kegg, drop = TRUE, showCategory = 15,color = "p.adjust",title = "KEGG Pathway")
## 气泡图
dotplot(kegg, showCategory=15)

根据自己需要画图

# 根据p-value值排序，选择显著的通路或功能（这里设定阈值为0.05）
significant_pathways <- subset(kegg_result, p.adjust < 0.05)
# 绘制基于ggplot2的条形图
ggplot(significant_pathways, aes(x = reorder(Description, -log10(p.adjust)), y = -log10(p.adjust))) +
geom_bar(stat = "identity", fill = "skyblue") +
labs(title = "KEGG Pathway Enrichment Analysis",
x = "Pathway",
y = "-log10(P-value)") +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
coord_flip()

气泡图

# 设置颜色的阶梯和对应的颜色
significant_pathways$pvalue_group <- cut(significant_pathways$p.adjust, breaks = c(0, 0.001, 0.01, 0.05), labels = c("p < 0.001", "0.001 <= p < 0.01", "0.01 <= p < 0.05"))

# 绘制基于ggplot2的气泡图
ggplot(significant_pathways, aes(x = reorder(Description, -log10(p.adjust)), y = -log10(p.adjust))) +
  geom_point(aes(size = Count, color = pvalue_group), alpha = 0.7) +
  scale_size_continuous(range = c(1, 10)) +
  labs(title = "KEGG Pathway Enrichment Analysis",x = "Pathway",y = "-log10(P-value)",size = "Count",color = "P-value") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  coord_flip()+
  scale_color_manual(values = c("p < 0.001" = "#DC0000B2", "0.001 <= p < 0.01" = "#F39B7FB2", "0.01 <= p < 0.05" = "#4DBBD5B2"))

GO分析

GO数据库（Gene Ontology database）是一个由国际合作组织维护的生物信息学资源库，用于对基因和蛋白质的功能进行注释和分类，以便更好地理解生物学过程和分析实验数据。
GO数据库的主要内容包括三个主要领域（也称为三个GO本体）：
生物学过程（Biological Process）：用于描述基因和蛋白质在生物学过程中的活动，例如细胞分裂、信号传导、代谢过程等。
分子功能（Molecular Function）：用于描述基因和蛋白质在分子水平上的功能，例如酶活性、DNA结合、受体活性等。
细胞组分（Cellular Component）：用于描述基因和蛋白质在细胞中所处的位置或组分，例如细胞核、细胞膜、线粒体等。

#GO富集分析
GO_all <- enrichGO(gene = DEG_data1$ENTREZID,  #基因列表(转换的ID)
               keyType = "ENTREZID",  #指定的基因ID类型，默认为ENTREZID
               OrgDb=org.Hs.eg.db,  #物种对应的org包
               ont = "ALL",   #CC细胞组件，MF分子功能，BF生物学过程，ALL以上三个
               pvalueCutoff = 0.01,  #p值阈值
               pAdjustMethod = "fdr",  #多重假设检验校正方式
               minGSSize = 10,   #注释的最小基因集，默认为10
               maxGSSize = 500,  #注释的最大基因集，默认为500
               qvalueCutoff = 0.01,  #q值阈值
               readable = TRUE)  #基因ID转换为基因名

GO_result <- data.frame(GO_all)

ONTOLOGY：表示BP、CC、MF`中的哪一类
ID：GO通路或功能类别的唯一标识符。
Description：GO通路的描述，即通路或功能的名称或注释。
GeneRatio：输入的差异基因中，映射到此GO通路的基因数目/映射到某一过程(BP/CC/MF之一)所有GO通路的基因总数目。
BgRatio：在整个基因数据库（或背景基因集合）中映射到GO通路的基因数目，以及该通路中的总基因数目。
pvalue、p.adjust和qvalue：p值、校正后p值和q值；
geneID：映射到该通路的基因的标识符列表，这里转换成基因名SYMBOL了。
Count：在富集分析中映射到该通路或功能类别的基因数目，即GeneRatio列中的第一个值。

## 自带的绘图
## 条形图
barplot(GO_all,showCategory = 20)
dotplot(GO_all, showCategory=15)

纵状柱状图

# GO三种类别，每种选择显著性最高的12个展示出来
go_enrichment_pathway <- GO_result %>% group_by(ONTOLOGY) %>% top_n(n = 12, wt = -p.adjust)

###纵向柱状图
ggplot(go_enrichment_pathway, aes(x=reorder(Description, -Count), y=Count, fill=ONTOLOGY)) +
  geom_bar(stat="identity", position=position_dodge(width=0.8), width=0.6) +
  theme_minimal() +
  labs(x="GO Term", y="Gene_Number", title="Top 10 Enriched GO Terms") +
  facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y')+
  coord_flip() +  #让柱状图变为纵向
  scale_fill_manual(values=c("CC"="skyblue","BP"="pink","MF"="lightgreen"))+
  theme_bw()

横装柱状图

ggplot(go_enrichment_pathway, 
       aes(x=reorder(Description, -Count),y=Count, fill=ONTOLOGY)) +  #x、y轴定义；根据ONTOLOGY填充颜色
  geom_bar(stat="identity", width=0.8) +  #柱状图宽度
  scale_fill_manual(values=c("CC"="skyblue","BP"="pink","MF"="lightgreen")) + #柱状图填充颜色
  facet_grid(.~ONTOLOGY, scale = 'free_x', space = 'free_x')+
  labs(x="GO Term", y="Gene_Number", title="Top 10 Enriched GO Terms")+
  theme_bw() + 
  theme(axis.text.x=element_text(family="sans",face = "bold", color="gray50",angle = 70,vjust = 1, hjust = 1 )) #对字体样式、颜色、还有横坐标角度（）

气泡图

# 绘制气泡图
ggplot(go_enrichment_pathway, aes(x=reorder(Description, Count), y=Count)) +
  geom_point(aes(size=Count,color=-log10(p.adjust))) +
  scale_size_continuous(range=c(1, 10)) +
  facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y')+
  coord_flip() +  #让柱状图变为纵向
  theme_minimal() +
  scale_color_gradient(low = "pink",high ="red")+
  labs(color=expression(-log10(p.adjust),size="Count"), 
       x="Gene Ratio",y="Gene_Number",title="GO Enrichment")+
  theme_bw()

GSEA富集分析

在常规的转录组分析中，往往会得到大量的差异表达基因，但如何将这些差异基因与生物学功能结合在一起成为一个挑战。传统的KEGG和GO分析可以提取差异基因序列，与预设的通路进行比较，并得到通路富集的结果，但是它们通常只提供了基因集的整体富集信息，没有考虑基因的调控方向（上调或下调）。因此，我们无法了解差异基因如何影响整个通路或生物学功能。
GSEA不依赖预定义的基因列表（一般指差异基因），而是使用整个基因表达数据集。它通过对基因排序和富集得分的计算，能够全面考察基因集的富集情况，并区分富集在基因表达上调或下调的基因集。这样，我们可以更好地理解差异基因在整个通路或生物学功能中的集体影响，而不仅仅是关注单个差异基因的表达情况。
做GSEA富集分析有两个输入，一个是基因列表，另一个是基因集。
基因列表就是用做完差异基因之后得到的基因列表（跟KEGG/GO不同，不需要根据条件筛选出差异基因）。基因集可以使用数据库中提供的基因集，当然也自己制作感兴趣的基因集。
基因集的下载
可以使用GSEA官网MSigDB提供的基因集，这里选择的是hallmark gene set，自己可以根据情况选择，MSigDB链接：https://www.gsea-msigdb.org/gsea/msigdb/index.jsp

#加载基因集（自行下载）
geneSet <- read.gmt("h.all.v2023.2.Hs.symbols.gmt") 
# 1、直接使用logFC作为排序依据
geneList <- LUNG_Match_DEG$logFC #获取GeneList
names(geneList) <- LUNG_Match_DEG$gene_name #使用转换好的ID，对GeneList命名
geneList <- sort(geneList, decreasing = T) #从高到低排序
GSEA_enrichment <- GSEA(geneList, #排序后的gene
TERM2GENE = geneSet, #基因集
pvalueCutoff = 0.05, #P值阈值
minGSSize = 10, #最小基因数量
maxGSSize = 500, #最大基因数量
eps = 0, #P值边界
pAdjustMethod = "BH") #校正P值的计算方法
result <- data.frame(GSEA_enrichment)
dim(GSEA_enrichment@result)

每一列代表的含义：
ID: 这一列代表基因集的标识符或名称。
Description: 这一列表示基因集的描述或注释，通常是对基因集功能、通路或生物学过程的解释。
setSize: 这一列表示基因集中包含的基因数量。
enrichmentScore: 这一列表示富集得分（Enrichment Score），它是一个反映基因集在基因表达数据中富集程度的统计量。
NES: 这一列表示标准化富集得分（Normalized Enrichment Score），它是将富集得分标准化后的值，使得不同基因集的富集得分可比较。
pvalue: 这一列表示富集得分的显著性水平（p-value），用于衡量基因集在基因表达数据中的显著性富集。
p.adjust: 这一列表示多重比较校正后的p-value，通常使用FDR或其他方法进行校正。
qvalue: 这一列表示估计的FDR（False Discovery Rate），用于控制多重假设检验引起的假阳性。
rank: 这一列表示基因集在排序后的基因列表中的排名。
leading_edge: 这一列指示哪些基因在计算富集得分时对富集结果产生了主要贡献。
tags=60%: 这表示60%的基因集中的基因在富集分析中对结果产生了影响。
list=10%: 这表示在富集分析中使用的整体基因列表（gene list）中，有10%的基因在该基因集中。
signal=55%: 这表示在整个基因集中，有55%的基因在样本中显示出富集信号，即在表达数据中呈现出差异表达的特征。(由tags和list计算得来)。
core_enrichment: 这一列指示哪些基因是核心富集基因，对于形成富集得分起关键作用。

GSEA结果绘图

## 展示最显著的15个通路
dotplot(GSEA_enrichment,showCategory=15,color="p.adjust") 
## 将通路分为激活和抑制两个部分
dotplot(GSEA_enrichment,split = ".sign")+facet_grid(~.sign)+
  theme(plot.title = element_text(size = 10,color="black",hjust = 0.5),
        axis.title = element_text(size = 10,color ="black"), 
        axis.text = element_text(size= 10,color = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1 ),
        legend.position = "top",
        legend.text = element_text(size= 10),
        legend.title= element_text(size= 10))

绘制特定通路富集情况

BiocManager::install("enrichplot")
library(enrichplot)
#看一下MYC通路
gseaplot2(GSEA_enrichment,"HALLMARK_MYC_TARGETS_V1",color="red",pvalue_table = T)

同时绘制多个通路

gseaplot2(GSEA_enrichment,c("HALLMARK_MYC_TARGETS_V1","HALLMARK_G2M_CHECKPOINT"),color=c("red","blue"),pvalue_table = T)