生信分析生信绘图科研学习进步

GSEA富集分析+R分析+多类美化工具

2022-09-23  本文已影响0人  粥粥zz

一.GSEA富集分析定义

GSEA(Gene Set Enrichment Analysis),是一种基于基因集的富集分析方法。具体介绍查阅Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles
GSEA的适用场景是:在两种不同的生物学状态下,可以理解为处理组与对照组,判断某一组基因集其表达模式更接近于那种过程或者通路,从而推断这些基因对这个生物学过程的重要贡献。在分析结果上GSEA也是更加可靠,可信度较高的。

二.GSEA富集分析优势

  基因表达分析的传统策略侧重于识别在两种感兴趣状态之间表现出差异的单个基因。尽管有用,但它们无法检测分布在整个基因网络中并且在单个基因水平上差异小的生物过程,例如代谢途径、转录程序和压力反应。但基因富集分析(GSEA)不需要指定明确的差异基因阈值,把基因按照在两组样本中的差异表达程度进行排序,然后采用统计学方法检验预先设定的基因集合是否在排序表的顶端或低段富集

  与单基因方法相比,GSEA 具有许多优势。
    首先,它通过识别路径和过程简化了对大规模实验的解释。研究人员可以专注于基因组,而不是专注于高分基因(可能注释不佳且可能无法重现),这往往更具重现性和可解释性。
    其次,当基因组的成员表现出很强的互相关性时,GSEA 可以提高信噪比,并使检测单个基因的适度变化成为可能。
    第三,前沿分析可以帮助定义基因子集以阐明结果。
    第四,常规的超几何富集分析富集到某个通路时会出现这种情况,这条通路既有上调的差异基因又有下调的差异基因,那么这条通路总体是被抑制还是激活是不清楚的。GSEA是基于基因集的富集分析方法,对基因表达量数据进行分析时,选择一个或多个功能基因集进行分析(以某个kegg通路为例,可以认为一个基因集),基因表达量的数据与表型或者不同样本的关联度进行排序,然后判断每个基因集内的基因是否位于基于表型相关度排序后的基因列表的上部或者下部,来判断词基因集内的基因协同变化对表型或不同处理样本的影响。

三.数学原理

1.输入文件

    GSEA的输入是一个基因表达量矩阵,其中的样本分成了A和B两组,找到两组之间差异表达的基因,然后根据foldchange进行排序,用来表示基因在两组间表达量的变化趋势。排序之后的基因列表其顶部可以看做是上调的差异基因,其底部是下调的差异基因。GSEA分析的是一个基因集下的所有基因是否在这个排序列表的顶部或者底部富集,如果在顶部富集,我们可以说,从总体上看,该基因集是上调趋势,反之,如果在底部富集,则是下调趋势。


image.png

2.三个关键步骤

GSEA主要包括三个步骤:计算富集得分(Enrichment Score);估计富集得分的显著性水平;多重假设检验。

第 1 步:计算富集分数。我们计算了一个浓缩分数(ES),它反映了一个集合S在整个排序列表L的极端(顶部或底部)处被过度代表的程度。分数是通过遍历列表L来计算的,当我们遇到 S 中的基因时增加运行总和统计量,当遇到不在 S 中的基因时减少它。增量*的大小取决于基因与表型的相关性. ES即为该过程中的x的最大取值,对应加权的Kolmogorov-Smirnov样的统计量。

第 2 步:估计 ES 的显着性水平。我们通过使用基于经验表型的置换测试程序来估计ES的统计显着性(标称P值),该程序保留了基因表达数据的复杂相关结构。
    具体来说,构建零分布:对每个样本重新分配表型标签、重新排序所有基因、重新计算基因集S的ES值;以上过程重复1000次,该1000个ES值构成零分布(null distribution)。
    计算P值:绘制直方图,p值即为ESnull distribution中≥或≤ESobserved的比例。该p值为经验名义p值。
    结果解读:小于α值(如0.05),则拒绝零假设,认为基因集S在排序列表L的top端或bottom端富集;若≥α值,则接受零假设,认为兴趣基因集S内基因在排序列表L中随机分布。

第 3 步:调整多重假设检验。当评估整个基因集数据库时,我们会调整估计的显着性水平以考虑多重假设检验。
   我们首先基于基因集S大小,对ES进行标准化处理,从而产生归一化富集分数 ( NES )。然后我们计算每个NES对应的错误发现率 (FDR) 来控制误报的比例。
   FDR值代表,给定NES值对应的基因集为假阳性(富集)的估计概率;值越大,则假阳性的概率越大。它是通过比较NES*的观测分布和空分布的尾部来计算的。

四.结果解读

1.富集通路的分析

富集分析的结果结果,要综合考虑p值及ES值两个指标


image.png

S1:基因集S1主要分布在排序列表的top端,ES分值较高,p值显著;
S2:基因集S2在排序列表中随机分布,ES值低,p值不显著;
S3:基因集S3非随机分布,但也并不在top or bottom呈现集中分布模式,ES值较高,但p值不显著。

2.分析结果

image.png

五.R语言分析

常用的方法有两种,一种是clusterProfiler包中的GSEA函数,另一种则是fgsea包中的fgsea方法

1.制作输入文件

#加载包
library(clusterProfiler)
library(gggsea)
library(ggplot2)
library(GSVA)
library(enrichplot)
library('GSEABase')
library(fgsea)
library(org.Hs.eg.db)
library(tidyverse)
library(dplyr)

 
#读取基因列表(包含基因名称以及 log2 转化后的 Fold Change 信息)
#根据logfc降序排列基因
load('deg.RData')
alldiff <- deg[order(deg$logFC,decreasing = T),]
genelist <- alldiff$logFC
names(genelist) <- rownames(alldiff)

#读取基因集(本示例使用的 MSigDB 数据库中的 hallmark)
#使用clusterProfiler中的read.gmt函数读取下载的gmt基因集文件
hallmark <- read.gmt("h.all.v2022.1.Hs.symbols.gmt")

2.clusterProfiler分析

gsea.re1<- GSEA(genelist,  #待富集的基因列表
    TERM2GENE = hallmark,  #基因集
    pvalueCutoff = 1,  #指定 p 值阈值(可指定 1 以输出全部)
    pAdjustMethod = 'BH')  #指定 p 值校正方法
     
#GSEA(geneList, exponent = 1, nPerm = 1000, minGSSize = 10,
#       maxGSSize = 500, pvalueCutoff = 0.05, pAdjustMethod = "BH", TERM2GENE,
#       TERM2NAME = NA, verbose = TRUE, seed = FALSE, by = "fgsea")
##exponent: weight of each step
##nPerm置换检验的次数,默认为1000
##minGSSize富集到某个条目的最小包含基因数,如果基因数小于该值则这个条目将被过滤掉,默认为10
##maxGSSize富集到某个条目的最大包含基因数,如果基因数大于该值则这个条目将被过滤掉,默认为500
##verbose是否输出提示信息
##seed是否使结果具有可重复性
##by选择使用的统计学方法,默认为fgsea

#输出
write.table(gsea.re1, 'gsea_human_hallmark.txt', sep = '\t', row.names = FALSE, quote = FALSE)

#提取显著富集的基因集
g1<-as.data.frame(gsea.re1)
g1<-subset(g1,p.adjust<0.05)
g1<-g1[order(g1$NES,decreasing = T),]

3.fgsea分析

## 这里去掉了基因集前缀
hallmark$ont <- gsub('HALLMARK_','',hallmark$ont)
#[[代表提取元素,比如数据框中data[1]提取第一列,数据类型还是data.frame,data[[1]]提取第一列,数据类型是character
#下面表示取第二列,并表示为character
hallmark.list <- hallmark %>% split(.$ont) %>% lapply( "[[", 2)

gsea.re2 <- fgsea(pathways = hallmark.list,#基因集列表
                  stats = genelist,#排序后的基因level,这里是logFC
                  nperm=1000,#置换检验的次数
                  minSize=1,#富集到某个条目的最小包含基因数,如果基因数小于该值则这个条目将被过滤
                  maxSize=10000,#富集到某个条目的最大包含基因数,如果基因数大于该值则这个条目将被过滤
#                 nproc = 0#如果不等于零,则将 BPPARAM 设置为使用的nproc (默认值 = 0)。
#                 gseaParam = 1,#GSEA 权重参数(0 为未加权,建议值为 1)
#                 BPPARAM  = NULL,#Bplapply中使用的并行化参数。可用于指定要运行的群集。如果没有显式初始化或通过设置“nproc”默认值“,bpparam()”被使用
                    )

#提取显著富集的基因集
colnames(gsea.re2)
g2 <- gsea.re2[gsea.re2$padj<0.05,]
g2 <- g2[order(g2$NES,decreasing = T),]
#输出
save(gsea.re1,g1,gsea.re2,g2,file = 'gsea.RData')

pathway – name of the pathway as in 'names(pathway)';
pval – an enrichment p-value;
padj – a BH-adjusted p-value;
ES – enrichment score, same as in Broad GSEA implementation;
NES – enrichment score normalized to mean enrichment of random samples of the same size;
nMoreExtreme' – a number of times a random gene set had a more extreme enrichment score value;
size – size of the pathway after removing genes not present in 'names(stats)'.
leadingEdge – vector with indexes of leading edge genes that drive the enrichment

4.可视化展示

可以绘制富集气泡图、得分图、条形图等
GSEA集的结果的可视化有一些R包和函数,比如gseaplot、gseaplot2、gggsea、GseaVis

load('gsea.RData')
num<-g1[,c(1,11)]
#添加富集的核心基因个数
num<-num%>% separate_rows(core_enrichment, sep = "/")%>%group_by(ID)%>%count()
num<-num[match(g1$ID,num$ID),]
sum(num$ID==g1$ID)
g1$Count<-num$n
data<-g1%>%mutate(GeneRatio = Count/setSize)
data$sign<-ifelse(data$NES>0,"activated","suppressed")
data$Description<- gsub('HALLMARK_','',data$Description)

library(RColorBrewer)
library(wesanderson)
ggplot(data)+geom_point(aes(x=GeneRatio ,y=Description,colour=p.adjust,size=Count))+facet_grid(~sign,scales = "free") +scale_colour_gradientn(colors=wes_palette("Zissou1",80,type="continuous"))
image.png
library(ggsci)
col_gsea1<-pal_simpsons()(16)

num1=1
gseaplot2(gsea.re1,geneSetID = rownames(g1)[1:num1],
              title = "",#标题
              color = col_gsea1[1:num1],#颜色
              base_size = 14,#基础大小
              rel_heights = c(1, 0.2, 0.4),#小图相对高度
              subplots = 1:3,#展示小图
              pvalue_table = FALSE,#p值表格
              ES_geom = "line"#line or dot
          )

#展示多个基因集,设置多个名称即可,比如这里展示4个,如果有特定基因集,定义完成后传入geneSetID参数即可
num2=4
gseaplot2(gsea.re1,geneSetID = rownames(g1)[1:num2],
          title = "",#标题
          color = col_gsea1[1:num2],#颜色
          base_size = 14,#基础大小
          rel_heights = c(1, 0.2, 0.4),#小图相对高度
          subplots = 1:3,#展示小图
          pvalue_table = FALSE,#p值表格
          ES_geom = "line"#line or dot
)
image.png
image.png
plotGseaTable(hallmark.list[g2$pathway],
              genelist , 
              gsea.re2,gseaParam = 0.5,
              colwidths = c(0.5,0.2,0.1,0.1,0.1)
)
image.png
names(hallmark.list)
se_hall<-c(head(g2$pathway,2),tail(g2$pathway,2))

#3 所选中的基因集
sig1<-g2[g2$pathway%in%se_hall,]

hallmark.se<-hallmark.list[sig1$pathway]
df.new <- gseaCurve(genelist , hallmark.se, sig1)

#主要由三部分组成
# ggplot() + geom_gseaLine(df.new) + theme_gsea() #曲线
# ggplot() + geom_gseaTicks (df.new) + theme_gsea()#竖线
# ggplot() + geom_gseaGradient(df.new) + theme_gsea()#色块

pal_line<-pal_lancet()(9)#各曲线颜色

ggplot() + 
  geom_gsea(df.new,
            prettyGSEA=T,
            tickcolor='grey30',#竖线颜色
            ticksize=0.1,#线条粗细
            #colour='grey',#色块边框
            linecolor=pal_line,
            linesize=2,lty=1,
            #多个图形时使用,设定每行列图形个数
            nrow = 2,
            ncol = 2
  ) + 
  theme_bw()+
  theme(strip.text = element_text(size = 7,face = 'italic'),#增大分面标题字体
        strip.background = element_rect(fill = 'white'))+
  xlab(bquote(italic('Rank')))+ylab(bquote(italic('Enrichment Score')))+
  theme(axis.text.x = element_text(size=8,angle = -30,face = 'plain',hjust = 0.5),
        axis.text.y = element_text(size=8,angle = 0,face = 'plain',vjust = 0.5))
image.png
library(GseaVis)
num1=1
#指定名称:图纸
p1<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1])
# 只保留曲线
p2<-gseaNb(object =  gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       subPlot = 1)
#曲线保留和热图retain curve and heatmap
p3<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       subPlot = 2)
# 名称太长截断换行
p4<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       subPlot = 2,
       termWidth = 30)
#标记里的一些名称:标记里
# add gene in specific pathway
mygene <- c("BIRC5","HMGB3","UBE2T","CCNB2")

# plot
p5<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       subPlot = 2,
       addGene = mygene)

#基因改变标签颜色和箭头类型:
p6<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       subPlot = 2,
       addGene = mygene,
       arrowType = 'open',
       geneCol = 'black')

#剩下所有图形:
p7<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       subPlot = 3,
       addGene = mygene,
       rmSegment = TRUE)

# clsaasic with pvalue
p8<-gseaNb(object = gsea.re1,
       geneSetID = rownames(g1)[1:num1],
       addGene = mygene,
       addPval = T,
       pvalX = 0.75,pvalY = 0.8,
       pCol = 'black',
       pHjust = 0)


image.png

修改主题,自定义颜色
添加统计指标、水平线
改变排列方式
https://mp.weixin.qq.com/s/TzczjJ1KfU3omux9yXH9oA
https://mp.weixin.qq.com/s/z5UrUMAXMxbCE_c4SB4JgA

上一篇下一篇

猜你喜欢

热点阅读