gsea富集分析结果怎么看_基因集富集分析(Gene Set Enrichment Analysis, GSEA)

89513a10f97ca938e216c136a1e54fe3.png

下面这篇文章中简单介绍了基因矩阵转置文件(gmt)

基因矩阵转置文件格式(* .gmt)​mp.weixin.qq.com
5cb89d24a0539e89e7388308fdbb7b6b.png

而下面这篇文章展示了如何使用R来读取gmt格式的文件

R读取gmt文件​mp.weixin.qq.com
88476ac313abe5f091d8e1636b0b2570.png

今天我们来看看如何做GSEA(Gene Set Enrichment Analysis,基因集富集分析)以及GSEA的结果如何解读。

首先我们需要了解一下GSEA跟传统的基因富集分析有什么区别,有什么优势。我相信大家在做传统的基因功能富集分析的时候肯定遇到这样的情况,一条富集到的通路中,既有上调的差异表达基因,也有下调的差异表达基因,那么这条通路总体是被抑制还是被激活呢?那么这条通路中的基因表达水平在实验组相比于对照组究竟是上升了呢,还是下降了呢?

在传统的富集分析时,我们其实根本不关心这些差异表达的基因究竟是上调还是下调。这是因为传统的富集分析根本不考虑基因表达量的变化趋势,其算法的核心只关注这些差异表达基因的分布是否跟随机抽样得到的分布一致,即使在后续可视化时,我们在通路图上用不同颜色标记了上调和下调的基因,但是由于没有采用有效的统计学方法去分析这条通路中所有差异基因的总体变化趋势,这使得传统的富集分析结果无法回答上述的问题。

即使有些文章里面根据差异表达基因的上下调将差异表达基因分成两组分别进行基因富集分析,这样得到的结果也会有失偏颇,并不能反应差异表达基因的整体情况。有时甚至会出现自相矛盾的情况,上调的基因和下调的基因富集到相同的一条通路中,这时就很难解释结果了。

GSEA(Gene Set Enrichment Analysis),该方法发表于2005年的Gene set enrichment analysis: a knowledge-based approach forinterpreting genome-wide expression profiles,是一种基于基因集的富集分析方法,在对基因表达数据分析时,首先确定分析的目的,即选择MSigDB中的一个或多个功能基因集进行分析(基因矩阵转置文件格式(* .gmt)中已经介绍过),然后基于基因表达数据与表型的关联度(也可以理解为表达量的变化)的大小进行排序。然后判断每个基因集内的基因是否富集于表型相关度排序后基因列表的上部或下部,从而判断此基因集内基因的协同变化对表型变化的影响。以上其实就是GSEA的分析原理。下面我们就借助一张图来帮助大家更好的理解GSEA的分析原理。

5c35b5256f94d7f0a1489810803e17c5.png

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

以上就是GSEA的分析原理,那么进行GSEA的结果怎样解读呢?

GSEA分析结果最常见的是下图:

e8c300f1a35fa09d709c5503ca0a2bbf.png

1、图最上面部分展示的是富集分数(ES,enrichment score)值计算过程,从左至右每到一个基因,计算出一个ES值,连成线。在最左侧或最右侧有一个特别明显的峰值就是基因集表型上的ES值。图中间部分每一条线代表基因集中的一个基因,及其在基因列表中的排序位置。

2、最下面部分展示的是基因与表型关联的矩阵,红色为与第一个表型(class A)正相关,在class A中表达高,蓝色与第二个表型(class B)正相关,在class B中表达高。

3、Leading-edge subset 对富集得分贡献最大的基因成员。若富集得分为正值,则是峰左侧的基因;若富集得分为负值,则是峰右侧的基因。

4、FDR GSEA默认提供所有的分析结果,并且设定FDR<0.25为可信的富集,最可能获得有功能研究价值的结果。但如果样品数目少,而且选择了gene_set作为Permumation type则需要使用更为严格的标准,比如FDR<0.05。

下面我们来看看如何使用R语言来进行GSEA分析,这里跟大家分享两种方法,一个是fgsea包,另一个是clusterProfiler包。

一、fgsea包,fgsea是Fast Gene Set Enrichment Analysis的缩写

1. 首先需要准备好rank文件,就是排好序的基因列表文件。一般做完差异表达分析都能得到这样一个文件。第一列是entrez gene id号,第二列可以是t值,也可以是foldchange。

> 

2. 同时还需要准备gmt文件,gmt文件的下载我们在基因矩阵转置文件格式(* .gmt)一文中已经讲解过,读入到R之后如下

> 

3. 看每一套基因集中的基因是否显著富集在差异表达基因列表的上部或者下部

#进行富集分析
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize=500, nperm=1000)

4. 显示显著富集在上部和下部的各10条通路的NES,pval和padj

由于ES是根据分析的数据集中的gene是否在一个功能gene set中出现来计算的,但各个功能gene set中包含的gene数目不同,且不同功能gene set与data之间的相关性也不同,因此,比较data set在不同功能gene set中的富集程度要对ES进行标准化处理,也就是NES
NES=某一功能gene set的ES/数据集所有随机组合得到的ES平均值

eff495029b3dce9e36419c15af1f8994.png

5. 柱形图展示显著富集在上部和下部的各20条通路的富集分值

8199fc81156f217a61e679ed84a91d1b.png

6. 显示特定通路上的富集结果

显著富集在cell cycle通路上部

5e3771e494e3592e485ddf893206533e.png

显著富集在chromation organization通路下部

60d003e477d5fb9d235b1135e4749f39.png

二、clusterProfiler

步骤跟fgsea包差不多

气泡图展示显著富集的前20条通路

9f5600ca93c2bbe7dc6d9146e3791a92.png

绘制cell cycle通路,显著富集在上部

2071c6d2a8569fc48ac1d70adbd03c6d.png

具体的R脚本参考下面这篇文章

基因集富集分析(Gene Set Enrichment Analysis, GSEA)​mp.weixin.qq.com
07accdff29fa1f743f1856accfeb188d.png
  • 3
    点赞
  • 35
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
以下是一个简单的基因GSEA分析代码示例,供参考: ```python import pandas as pd import numpy as np from scipy.stats import norm # 读取基因表达数据 gene_expression_data = pd.read_csv("gene_expression_data.csv") # 读取基因合数据 gene_set_data = pd.read_csv("gene_set_data.csv") # 计算每个基因的表达值的均值 gene_expression_mean = gene_expression_data.mean(axis=1) # 基于基因表达均值计算基因的排名 gene_ranks = gene_expression_mean.rank(method='min') # 初始化GSEA参数 n_genes = len(gene_ranks) n_perms = 1000 gene_set_size_min = 15 gene_set_size_max = 500 alpha = 0.05 # 计算基因合的ES值 def calculate_es(gene_set): gene_ranks_in_set = gene_ranks[gene_ranks.index.isin(gene_set)] gene_set_size = len(gene_ranks_in_set) if gene_set_size < gene_set_size_min: return np.nan else: es = 0 running_sum = 0 for gene_rank in gene_ranks: if gene_rank in gene_ranks_in_set: running_sum += 1 else: running_sum -= 1 es = max(es, running_sum) return es / gene_set_size # 计算基因合的NES值 def calculate_nes(gene_set, gene_set_es, gene_set_size): gene_set_mean = gene_expression_mean[gene_expression_mean.index.isin(gene_set)] gene_set_sd = np.std(gene_set_mean, ddof=1) if gene_set_sd == 0: return np.nan else: nes = gene_set_es / gene_set_sd return nes # 计算随机排列基因合的ES值 def calculate_permutation_es(gene_set): permuted_ranks = gene_ranks.sample(frac=1).reset_index(drop=True) permuted_es = calculate_es(gene_set, permuted_ranks) return permuted_es # 计算基因合的p值和FDR校正的p值 def calculate_p_value(gene_set_es, permutation_es): p_value = (np.sum(permutation_es >= gene_set_es) + 1) / (len(permutation_es) + 1) fdr_p_value = p_value * n_perms return p_value, fdr_p_value # 进行GSEA分析 results = [] for index, row in gene_set_data.iterrows(): gene_set = row['gene_set'].split(',') gene_set_es = calculate_es(gene_set) if np.isnan(gene_set_es): continue gene_set_size = len(gene_set) if gene_set_size < gene_set_size_min or gene_set_size > gene_set_size_max: continue permutation_es = [calculate_permutation_es(gene_set) for i in range(n_perms)] p_value, fdr_p_value = calculate_p_value(gene_set_es, permutation_es) if fdr_p_value <= alpha: nes = calculate_nes(gene_set, gene_set_es, gene_set_size) results.append({ 'gene_set': row['gene_set'], 'es': gene_set_es, 'nes': nes, 'p_value': p_value, 'fdr_p_value': fdr_p_value, 'gene_set_size': gene_set_size }) # 输出结果 results_df = pd.DataFrame(results) results_df = results_df.sort_values(by=['fdr_p_value'], ascending=True) print(results_df) ``` 这个代码是一个简单的基因GSEA分析代码示例,具体实现过程如下: 1. 读取基因表达数据和基因合数据,计算每个基因的表达值的均值。 2. 基于基因表达均值计算每个基因的排名。 3. 初始化GSEA参数,包括基因合的大小范围、置信水平等。 4. 定义计算基因合的ES值、NES值、随机排列基因合的ES值、p值和FDR校正的p值的函数。 5. 对每个基因合进行GSEA分析,计算其ES值、NES值、p值和FDR校正的p值,筛选出FDR校正的p值小于等于置信水平的基因合,并将结果保存在结果列表中。 6. 输出结果,按照FDR校正的p值升序排列。 需要注意的是,这个代码只是一个简单的示例,实际应用中可能需要根据具体情况进行修改和优化。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值