转录组分析—再谈GSEA

转录组分析—再谈GSEA

之前一直对GSEA的分析朦朦胧胧，这里再重新梳理下相关的知识点

1 相关概念

Gene Set Enrichment Analysis (GSEA) 是一种生物信息学方法，用于确定基因集合（gene sets）在基因表达数据中的显著性变化。它广泛应用于基因表达数据的功能解释，帮助研究者理解在特定实验条件下哪些生物学通路或功能类别是活跃的。以下是GSEA的相关知识点：

1. 基本概念和背景
2. 基因集合（Gene Set）：GSEA分析的基础是<span style="color: red;">预先定义的基因集合</span>，这些集合可以是与某一特定生物学过程、功能或通路相关的基因。例如，KEGG通路、GO（Gene Ontology）生物过程、反应富集等。
3. 富集分析（Enrichment Analysis）：目的是确定这些基因集合在某一条件下是否富集，即基因集合内的基因在实验条件与对照条件之间是否表现出显著差异。
4. GSEA的原理

GSEA的基本步骤包括：

1. 排序基因列表：首先，根据某种度量（如差异表达的统计量、相关系数等）对所有基因进行排序。
2. 计算富集分数（Enrichment Score, ES）：对每一个基因集合，GSEA计算其富集分数，这个分数反映了集合中基因在排序列表中的分布情况。通常使用“步行统计量”（walk statistics），沿着排序列表移动，碰到集合中的基因时增加分数，否则减少分数。
3. 随机化和统计检验：为了评估富集分数的显著性，GSEA通过在原始数据上多次随机化（置换）生成富集分数的空分布。然后将观察到的富集分数与随机生成的富集分数进行比较，计算出p值和FDR（False Discovery Rate）。
4. 显著性评估：确定基因集合在基因排序中的显著性，并识别在特定条件下显著上调或下调的通路或功能。
5. GSEA的特点和优点
6. 无需预先筛选基因：与传统的富集分析不同，GSEA不需要预先筛选出显著差异的基因。这减少了因阈值选择而可能导致的信息损失。（<span style="color: red;">为什么流程是用的差异基因呢</span>>）
7. 捕捉微弱但一致的信号：GSEA可以识别单个基因表达变化不显著但作为一个集合具有一致变化趋势的基因集合。
8. 基于所有基因的排序：通过使用排序列表的所有基因，GSEA可以利用全部数据的信息，而不仅仅是差异显著的基因。
9. GSEA的应用和局限性

应用

• 生物学通路和功能注释：GSEA广泛用于解释基因表达数据，尤其是在理解特定条件下的生物学通路的活跃情况。
• 药物作用机制研究：通过比较药物处理和对照样品的基因表达数据，可以识别药物作用的潜在机制。
• 疾病机制研究：用于比较健康和疾病样本，揭示疾病相关的基因通路。

局限性

• 依赖于基因集合的定义：GSEA的结果强烈依赖于预定义的基因集合，因此高质量和适当的基因集合是关键。
• 敏感性较低：由于GSEA利用的是排序列表的全局信息，对个别基因的微小变化敏感性较低。
• 假阳性和假阴性：尽管使用了FDR校正，但假阳性和假阴性结果仍可能存在。
• GSEA的实现和工具

常用的GSEA工具包括：

• Broad Institute的GSEA软件：经典的GSEA工具，提供了丰富的预定义基因集合和直观的可视化。
• R包：clusterProfiler：R语言中的GSEA实现，提供了KEGG、GO等多种数据库的支持，具有较强的可定制性。
• MSigDB：提供了广泛的基因集合数据库，可与GSEA工具结合使用。
• GSEA的结果解释
• 富集分数（ES）：表示基因集合在排序列表中的富集程度。正值表示基因集合的基因倾向于在排序列表的前端（上调），负值则表示在排序列表的后端（下调）。
• 正态化富集分数（NES）：标准化的ES，用于跨基因集合比较。
• p值和FDR：用于评估富集的显著性。FDR校正用于多重检验。

通过GSEA，研究者可以更全面地理解复杂生物系统中的基因调控网络和通路活动情况，进而揭示潜在的生物学机制。

2 全部代码

全部代码

只能用一个包的差异分析结果来做GSEA? 这里直接导入了limma包差异分析后的结果

rm(list = ls()) 
load( file =  'DEG_limma_voom.Rdata' )
colnames(DEG_limma_voom)

nrDEG=DEG_limma_voom[,c("logFC", "adj.P.Val"  )]
head(nrDEG) 
colnames(nrDEG)=c('logFC','P.Value')

library(org.Hs.eg.db)
library(clusterProfiler)
gene <- bitr(rownames(nrDEG), fromType = "SYMBOL",
             toType =  "ENTREZID",
             OrgDb = org.Hs.eg.db)
head(gene)
head(nrDEG)
nrDEG = nrDEG[rownames(nrDEG) %in% gene$SYMBOL,]
nrDEG$ENTREZID = gene$ENTREZID[match(rownames(nrDEG) , gene$SYMBOL)]
  # https://www.ncbi.nlm.nih.gov/gene/?term=SPARCL1

geneList=nrDEG$logFC
names(geneList)=nrDEG$ENTREZID
geneList=sort(geneList,decreasing = T)
head(geneList)

library(clusterProfiler)
# 一本书，很多数据库，很多可视化
#～～～这里需要替换物种～～～
#mmu
#mmu
str(geneList)
kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',#按需替换
                  nPerm        = 1000,
                  minGSSize    = 10,
                  pvalueCutoff = 0.99,
                  verbose      = FALSE)

tmp=kk_gse@result
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',
                     keyType='ENTREZID')#按需替换
kk@result$Description=gsub(' - Mus musculus \\(house mouse\\)',
                           '',kk@result$Description )
head(kk@result$Description)
tmp=kk@result 


pro='test'
write.csv(kk@result,file = file.path(paste0(pro,'_kegg.gsea.csv')))
save(kk,file = file.path( 'gsea_kk.Rdata'))

# 假如这里找不到显著下调的通路，可以选择调整阈值，或者其它。
# down_kegg<-kk_gse[kk_gse$pvalue<0.01 & kk_gse$enrichmentScore < -0.6,];down_kegg$group=-1
# up_kegg<-kk_gse[kk_gse$pvalue<0.01 & kk_gse$enrichmentScore > 0.3,];up_kegg$group=1
# down_k <- down_kegg[head(order(down_kegg$pvalue,decreasing = F)),]
# up_k <- up_kegg[head(order(up_kegg$pvalue,decreasing = F)),]

#～～～取前6个上调通路和6个下调通路～～～
up_k <- kk[head(order(kk$enrichmentScore,decreasing = T)),];up_k$group=1
down_k <- kk[tail(order(kk$enrichmentScore,decreasing = T)),];down_k$group=-1

dat=rbind(up_k,down_k)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]

p7 <- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low="#34bfb5",high="#ff6633",guide = FALSE) + 
  scale_x_discrete(name ="Pathway names") +
  scale_y_continuous(name ="log10P-value") +
  coord_flip() + 
  theme_ggstatsplot()+
  theme(plot.title = element_text(size = 15,hjust = 0.5),  
        axis.text = element_text(size = 12,face = 'bold'),
        panel.grid = element_blank())+
  ggtitle("Pathway Enrichment") 
p7
ggsave(p7,filename = 'kegg_top_gsea.png' ,width = 8,height = 4)

#ggsave(p7,filename = 'step5_kegg_gsea.pdf',width = 8,height = 4)

# geneSetID对应表格的id，排序后取前6个和后六个
p8up <- enrichplot::gseaplot2( kk, geneSetID = head( rownames(up_k)) )  
pdf(file=file.path("step5_kegg_up_gseaplot.pdf"),width = 7,height = 6)
print(p8up)
dev.off()
# ggsave(p8up[[1]],filename = 'step5_kegg_up_gseaplot.png',width = 7,height = 4)#直接p8up 保存是不行的

p8down <-   enrichplot::gseaplot2(kk, geneSetID = head( rownames(down_k)))
pdf(file=file.path("step5_kegg_top_down_gseaplot.pdf"),width = 7,height = 6)
print(p8down)
dev.off()

3 解析及可视化

nrDEG$ENTREZID = gene$ENTREZIDmatch(rownames(nrDEG), gene$SYMBOL)

match(rownames(nrDEG), gene$SYMBOL)：在gene$SYMBOL中查找nrDEG中的每个基因符号的位置，返回一个整数向量，该向量中的每个元素表示nrDEG中的基因符号在gene$SYMBOL中的位置。

gene$ENTREZID...：使用上一步得到的整数向量从gene$ENTREZID向量中提取对应位置的Entrez基因ID。

geneList = nrDEG$logFC names(geneList) = nrDEG$ENTREZID geneList = sort(geneList, decreasing = TRUE) head(geneList)

<span style="color: red;">这是一个数值向量，包含基因的排序信息（例如，log2FoldChange或其他评分）。向量的名称通常是Entrez基因ID</span>

kk_gse <- gseKEGG(geneList = geneList, organism = 'hsa',#按需替换 nPerm = 1000, minGSSize = 10, pvalueCutoff = 0.99, verbose = FALSE)

使用clusterProfiler包中的gseKEGG函数进行基因集富集分析，具体来说是针对KEGG通路的基因集富集分析（GSEA）。以下是各个参数的详细解释：

1. geneList:
   • 。gseKEGG函数使用这个列表来计算基因集的富集分数。
2. organism = 'hsa':
   • 这是物种代码，用于指定所分析的生物物种。'hsa'表示人类（Homo sapiens）。可以根据研究对象的不同，替换为其他物种的代码，例如'mmu'表示小鼠（Mus musculus）。
3. nPerm = 1000:
   • 表示重采样（置换）次数，用于计算p值。GSEA算法通过多次随机置换基因列表来估计富集分数的显著性。这里设置为1000次。
4. minGSSize = 10:
   • 最小基因集大小。只有包含至少10个基因的基因集才会被考虑在分析中。这有助于避免分析中包含过小的基因集。
5. pvalueCutoff = 0.99:
   • p值截断阈值。只有p值小于或等于0.99的基因集才会被保留在结果中。通常会选择更小的值（如0.05）来筛选显著的结果，但这里设置为0.99可能是为了保留更多的基因集用于后续筛选。
6. verbose = FALSE:
   • 是否显示运行过程的详细信息。设置为FALSE表示不显示详细信息。

tmp = kk_gse@result kk = DOSE::setReadable(kk_gse, OrgDb = 'org.Hs.eg.db', keyType = 'ENTREZID') # 按需替换 kk@result$Description = gsub(' - Mus musculus (house mouse)', '', kk@result$Description) head(kk@result$Description) tmp = kk@result

kk_gse@result是一个数据框，包含基因集富集分析的详细结果。

使用DOSE包中的setReadable函数，将结果中的Entrez基因ID转换为更加易读的基因符号.

OrgDb = 'org.Hs.eg.db指定了人类的基因注释数据库，keyType = 'ENTREZID表示输入的基因ID类型是Entrez基因ID。转换后的结果存储在kk对象中。

使用gsub函数，从kk@result$Description中删除包含的“ - Mus musculus (house mouse)”字符串，通常是在描述中包含物种信息时使用。gsub函数的作用是查找并替换字符串，这里将匹配到的字符串替换为空字符串''。

tmp=kk@result

<span style="color: red;">整理好的GSEA矩阵</span>

up_k <- kk[head(order(kk$enrichmentScore,decreasing = T)),];up_k$group=1
down_k <- kk[tail(order(kk$enrichmentScore,decreasing = T)),];down_k$group=-1

根据基因集富集分析结果中的富集分数（enrichment score，简称ES）将基因集分为上调和下调两组。具体来说：

1. up_k：包含富集分数最高的几个基因集（上调）。
2. down_k：包含富集分数最低的几个基因集（下调）。

kk：包含GSEA结果的对象，其中包括每个基因集的富集分数等信息。

order(kk$enrichmentScore, decreasing = TRUE)：按富集分数（enrichmentScore）降序排列基因集，返回排序后的索引。

head(...)：返回排序后前几项（默认前6项）的索引，即富集分数最高的基因集的索引。

kk...：根据索引提取富集分数最高的基因集信息，存储在up_k对象中。

up_k$group = 1：为up_k添加一列group，值为1，用于标记这些基因集为上调，同理，-1表下调。

dat=rbind(up_k,down_k)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]

将up_k和down_k数据框按行绑定，合并成一个新的数据框dat。这一步将上调和下调的基因集结果合并到一起。

将dat中的p值（pvalue）转换为其负对数值（-log10），

dat$pvalue = dat$pvalue * dat$group：根据dat中的group列，将变换后的p值乘以1或-1。group列的值为1表示上调，-1表示下调。这样，上调的基因集的p值保持为正，下调的基因集的p值变为负。这种处理方式便于后续分析中区分上调和下调的基因集。

dat = datorder(dat$pvalue, decreasing = FALSE),：按dat中的p值列对数据进行排序。从最小的p值（负数）到最大的p值（正数）。排序后，下调的基因集（负数）会排在前面，上调的基因集（正数）会排在后面。

kegg_top_gsea.png

step5_kegg_top_down_gseaplot.png

step5_kegg_up_gseaplot.png

PS:

GSEA不是应该拿全部的基因来做分析吗？为什么这里使用limma包差异分析后得到的基因来做差异分析呢?

GSEA（Gene Set Enrichment Analysis）通常是基于所有基因的排序结果进行分析，而不是仅仅使用差异表达基因。然而，在实际应用中，有时会出现使用差异表达分析结果进行后续分析的情况。

完整基因集分析：传统的GSEA是基于全基因表达数据的排序来评估基因集的富集情况。这种方法不要求预先筛选出差异表达基因，而是通过对基因表达数据的排序，计算每个基因集的富集得分。

特定基因集分析：有时，研究者可能更关心特定的基因集（如DEGs）的功能或通路富集情况。在这种情况下，使用差异表达分析后的基因（如nrDEG）来进行富集分析可以集中探讨这些显著变化的基因是否在特定的生物学通路或功能类别中有富集倾向。

补充：之前的理解有误，这个流程并非用差异基因来做的，而是用全部的基因，导入的DEG等是为全部的基因，而非筛选过的基因，脑子瓦特了。