什么是富集分析?手把手教你理解富集分析

1. 为什么要做富集分析

生物医学研究领域,做完差异分析,获得大量差异基因只是研究的起点。面对成百上千的差异基因,我们不可能一个一个查找这些差异基因的生物学功能,如何快速定位它们的生物学功能和潜在意义?富集分析便是打开这一奥秘之门的关键钥匙,它能帮助我们从基因的 “汪洋大海” 中,打捞起有价值的生物学线索。

2. 什么是富集分析

富集分析是生物医学研究中,对差异基因进行功能解读的关键手段。它通过将差异分析得到的基因,与已知的生物学功能和通路信息进行比对,找出这些基因在哪些生物学过程、细胞组分或分子功能中显著聚集 ,从而揭示基因与生命活动、疾病发生发展的关联。

这一分析基于众多专业数据库,其中 GO 数据库从分子功能、细胞组分、生物过程三个层面描述基因功能;KEGG 数据库聚焦代谢通路和信号转导通路;此外还有 Reactome 等提供更丰富视角。

在算法上主要基于超几何分布原理,通过统计计算判断差异基因在特定通路或功能类别中的富集程度是否显著高于随机分布,为生物学结论提供可靠证据 。

3. 富集分析:核心数据库解析

科学家们耗费大量精力构建了丰富的基因功能注释数据库,为富集分析提供了坚实的基础。其中,GO(Gene Ontology)和KEGG(Kyoto Encyclopedia of Genes and Genomes) 是应用最为广泛的两大数据库。

GO 数据库:它是一个结构化、受控的词汇表,从三个不同层面描述基因功能。

    • 分子功能(MF,Molecular Function关注基因产物在分子水平的活性,比如某个基因是否具有酶催化活性;
    • 细胞组分(CC,Cellular Component界定基因产物在细胞中的位置,像在线粒体、细胞膜还是细胞核发挥作用;
    • 生物过程(BP,Biological Process则描绘基因参与的一系列生物学事件,例如细胞增殖、凋亡等过程。

这三个层面相互补充,能全面探究基因功能的改变,就像从多个角度观察一个复杂的机器零件如何工作。

KEGG 数据库:侧重于代谢通路和信号转导通路,是细胞内代谢和调控网络的 “地图集”。它可以用图形化的方式展示基因在代谢通路中的上下游关系,比如糖代谢通路中,不同基因如何协同作用将葡萄糖转化为能量。通过 KEGG 富集分析,我们可以清晰看到差异基因在细胞代谢通路图上汇聚于哪些关键路段,从而找到潜在的药物靶点或疾病机制。

除了这两大 “巨头”,还有Reactome数据库,GO Slim数据库等为富集分析提供了更多元化的选择。

4. 富集分析原理算法:超几何分布

富集分析背后的核心算法是超几何分布。

假设我们将全基因组比作一个装满彩球的大箱子,每个彩球代表一个基因,不同颜色的彩球对应不同的功能通路。

我们从箱子里随机抓取一把彩球(这些彩球就相当于差异分析得到的差异基因),然后统计抓取的彩球中,某一种颜色(某一功能通路)彩球的比例。

超几何分布算法就是:计算这个比例与箱子里该颜色彩球的原始比例相比,是否显著更高。如果显著更高,就意味着我们抓取的这把彩球(差异基因)在这个颜色对应的功能通路(功能或通路)中出现了 “富集”,即这些差异基因在这个通路中过度表达,很可能在该通路中发挥着关键作用。

下面是超几何分布分公式:

注:

  • N为注释数据库中的总基因数;
  • M为数据库中属于某个注释条目的基因数;
  • n为我们得到的需要进行富集分析的基因的总数目;
  • m为n中属于M的数目。

5. 富集分析结果表格解读

以下就是GO富集分析的核心结果表格。KEGG的结果也类似这样。

这个表格包含了富集分析结果的相关数据,其中每一列代表不同的指标。

下面逐个解释各列的含义:

  • Class: 该列表示所分析的基因分类的类型。例如GO中,分为了BP, CC, MF
  • TermID: 这是每个富集条目对应的唯一标识符。
  • Term: 这是与"TermID"对应的富集条目的名称或描述。
  • bg_pro: 表示背景基因集中基因总数。(背景基因可以用鉴定到的基因当背景,也可以用该物种全部的有注释的基因当背景)
  • bg_term: 背景基因集中实际与该富集条目相关的基因数目。
  • fg_pro: 表示前景基因集中基因总数。(前景基因个数,一般是差异基因的个数,当然也可以只放上调基因,或只放下调基因做富集分析)
  • fg_term: 前景基因集中实际与该富集条目相关的基因数目。
  • Ratio: 该列表示前景基因集中富集条目相关基因数量与背景基因集中富集条目相关基因数量的比值。也就是fg_term/bgtrem
  • PValue: P值是用于评估该富集条目是否具有统计显著性的值。通过超几何分布算法计算出P值。约定俗成,p小于0.05,认为显著富集。
  • Enrichment: 富集度,计算方法为-log10(PValue),用于评估富集的强度。
  • FDR:校正后的P值,用于控制多重假设检验中的假阳性率。
  • GeneSymbol:富集到该条目下具体是哪些基因

6. 富集分析结果图片的解读

6.1. 富集气泡图

气泡图共有四个维度来描述数据,x轴,y轴,点大小,点颜色。

富集气泡图就是把富集结果的ratio,Term名称,富集条目个数,富集度或PValue,分别映射到气泡图的四个维度上。

图形解释:

1. 横坐标:Ratio。表示前景基因集中富集条目相关基因数量与背景基因集中富集条目相关基因数量的比值。

2. 纵坐标: Term。显示富集条目的名称或描述,每一行代表一个富集条目。通常条目按照富集程度排序,方便观察哪些富集项更为显著。

3. 气泡:通过气泡的阿颜色,位置、大小,来传达不同富集条目的显著性、富集程度以及相关基因的数量。气泡的大小代表富集到该条目的基因数量。气泡越大,表示与该富集条目相关的基因数量越多。气泡的颜色反映了富集度或P值的大小。

6.2. GO注释三联图

GO(Gene Ontology)注释三联图,也叫 GO 分类柱状图 ,用于展示GO注释结果的常见可视化形式,由BP, CC, MF 三个维度的柱状图组成,用于呈现不同 GO 分类下基因的基因富集情况,帮助研究者快速了解基因在功能、定位和参与过程等方面的富集信息 。

图形解释:

1. 横坐标:是每个GO term,前景基因集中实际与该富集条目相关的基因数目。

2. 纵坐标:富集条目。不同颜色条代表GO的三个子类(BP、CC、MF),每个子类展示前n个富集度最高的 GO 词条。

6.3. 富集弦图

富集弦图:分别展示出每个Term富集到哪些基因以及这些基因的上下调情况。

图形解释:

1. 基因:图形左侧为差异基因,颜色变化表示logFC从小到大的变化,红色为上调,蓝色为下调。

2. 富集条目:图形右侧为富集结果的Term,

3. 连线:连接基因和富集条目,表示这些基因富集到相应的富集条目中。

富集弦图一般显示3-10个通路较为合适,太多通路,会导致太过密集,反而看不清。

6.4. 富集圈图

富集圈图(Enrichment Circle Plot)是一种展示基因富集分析结果的可视化工具。下图由四个圈组成,展示了从外到内的富集分析的多维度信息。

图形解释:

第一圈:富集条目与分类这一圈显示的是基因富集的不同条目与对应的分类。每个类别由不同的颜色表示,以便清晰区分。

第二圈:每个条目背景基因的数量及P值这一圈显示每个条目中背景基因的数目,同时标注该分类的统计显著性值。基因数目较多的条形图较长,而P值越小,富集度越高的部分颜色越红,表明该分类的富集显著性越强。这一圈的作用是展示该分类在背景基因中的频次及其统计显著性。

第三圈:上下调基因比例条形图这一圈展示了每个分类中上调和下调基因的比例。条形图上,深紫色代表上调基因,浅紫色代表下调基因。条形的长度表示该分类中上下调基因的比例。如果输入的绘图数据没有区分上下调基因(,第三圈将显示前景基因的总数目。

第四圈:Ratio值这一圈展示了每个条目的Ratio值,即该条目中前景基因的数量与背景基因的比例。

6.5. 桑基气泡图

桑基气泡图将桑基图的流动路径和气泡图的数量表现结合起来,能够同时展示数据之间的流动关系和每个类别的相对大小。

相比于传统的气泡图,还能展示富集条目下具体富集了哪些基因。

该图结合了桑基图和气泡图的特点。

左侧为桑基图,展示每个pathway中包含的基因信息,流动路径表示基因在不同pathway之间的分布情况。

右侧为常规的气泡图,气泡的大小表示每个pathway所属的基因数量,气泡的颜色则反映了该pathway的P值,通过颜色深浅区分统计显著性。这样,图表直观展示了基因在不同pathway中的分布及其统计显著性。

7. 分析绘图的数据准备

2种格式,差异蛋白列表或者鉴定蛋白列表。

示例数据可以在BioLadder v2.0-生物信息在线分析云平台找到并下载。

7.1. 鉴定蛋白列表

只需要一列ID,可以使Accession,Symbol,uniprot Entry Name等常见ID。会以鉴定蛋白为前景,背景为该物种uniprot数据库。

7.2. 差异蛋白格式

需要有3列,第一列是ID,第二列是差异分析完后的log2FC,第三列是上下调情况,1表示上调,-1表示下调。最后富集分析时,会把上调,下调,和所有的差异蛋白分别充当前景,背景是所有的差异蛋白。

8. 做富集分析的方法

8.1. BioLadder v2.0 云平台在线分析绘图

BioLadder v2.0 是一个在线绘图平台,我们只需要上传数据,系统自动生成图表,不需要写代码,适合0代码基础的同学。网址:BioLadder v2.0-生物信息在线分析云平台

8.2. R语言做富集分析

喜欢自己写代码的同学,也附上R语言代码可供参考。

代码源自https://zhuanlan.zhihu.com/p/488948980

8.2.1. 安装clusterProfiler包及所需文件
# 安装clusterProfiler包。
# 可能会报错,可以下载到本地安装,可参考https://zhuanlan.zhihu.com/p/436671645
BiocManager::install("clusterProfiler")

# 安装所需的物种数据库注释文件
BiocManager::install("org.Hs.eg.db") # 人的注释数据库。其他物种的可在https://bioconductor.org/packages/3.6/data/annotation/找到
8.2.2. ID Mapping

将ID转换为需要的格式

# 通过bitr函数,做Symbol,Gene ID,uniprot ID 等之间的转换
# 支持如下类型之间的相互转换:ACCNUM, ALIAS, ENSEMBL, ENSEMBLPROT, ENSEMBLTRANS, ENTREZID, ENZYME, EVIDENCE, EVIDENCEALL, GENENAME, GO, GOALL, IPI, MAP, OMIM, ONTOLOGY, ONTOLOGYALL, PATH, PFAM, PMID, PROSITE, REFSEQ, SYMBOL, UCSCKG, UNIGENE, UNIPROT
library(clusterProfiler)
symbolList <- c("GPX3",  "GLRX",   "LBP",   "CRYAB", "DEFB1", "HCLS1",   "SOD2",   "HSPA2", 
                 "ORM1",  "IGFBP1", "PTHLH", "GPC3",  "IGFBP3","TOB1",    "MITF",   "NDRG1", 
                 "NR1H4", "FGFR3",  "PVR",   "IL6",   "PTPRM", "ERBB2",   "NID2",   "LAMB1", 
                 "COMP",  "PLS3",   "MCAM",  "SPP1",  "LAMC1", "COL4A2",  "COL4A1", "MYOC",  
                 "ANXA4", "TFPI2",  "CST6",  "SLPI",  "TIMP2", "CPM",     "GGT1",   "NNMT",
                 "MAL",   "EEF1A2", "HGD",   "TCN2",  "CDA",   "PCCA",    "CRYM",   "PDXK",  
                 "STC1",  "WARS",  "HMOX1", "FXYD2", "RBP4",   "SLC6A12", "KDELR3", "ITM2B")
eg <- bitr(symbolList, 
           fromType="SYMBOL",    # 定义输入类型
           toType=c("ENTREZID","ENSEMBL","UNIPROT"),  # 输出何种类型
           OrgDb="org.Hs.eg.db") # (人)物种数据库

8.2.3. GO富集
genelist <- eg$ENTREZID
genelist = genelist[-duplicated(genelist)] # 去重

go <- enrichGO(genelist,
               OrgDb = org.Hs.eg.db,   # 选择对应物种的数据库
               ont='ALL',              # 选择GO富集的种类"BP", "MF", "CC" OR "ALL"
               pAdjustMethod = 'BH',   # FDR算法
               pvalueCutoff = 1,       # 过滤p值小于多少的条目,1为不过滤
               qvalueCutoff = 1,       # 过滤q值小于多少的条目,1为不过滤
               keyType = 'ENTREZID',   # 定义输入类型为NCBI的Gene ID
               readable = T            # 是否将结果表格中的富集条目改为Gene Symbol,增加可读性
)
head(go)
barplot(go,showCategory=20,drop=T) # 绘制富集条形图
dotplot(go,showCategory=20)        # 绘制富集气泡图
8.2.4. KEGG富集
# 由于kegg官网的限制,需要用下面的命令将下载方式更改为自动,也许新版本的包,解决了这个问题。2022年11月16日更新
R.utils::setOption( "clusterProfiler.download.method",'auto' )
kegg <- enrichKEGG(genelist, 
                   organism = 'hsa',     # 选择对应物种的数据库,通过如下网址,获取对应物种的三字母缩写http://www.genome.jp/kegg/catalog/org_list.html
                   keyType = 'kegg',     # 定义输入类型
                   pvalueCutoff = 1,     # 过滤p值小于多少的条目,1为不过滤
                   qvalueCutoff = 1,     # 过滤q值小于多少的条目,1为不过滤
                   pAdjustMethod = 'BH', # FDR算法
                   use_internal_data = F # 是否使用KEGG.db本地数据库
)
head(kegg)
barplot(kegg,showCategory=20,drop=T)     # 绘制富集条形图
dotplot(kegg,showCategory=20)            # 绘制富集气泡图

更多图形的绘制

富集气泡图:bubble – R2Omics

富集弦图:enrichChord – R2Omics

富集圈图:R语言绘制富集圈图 – R2Omics

桑基气泡图:R语言绘制桑基气泡图 – R2Omics

### 关于单细胞转录组测序数据分析的详细程 #### 数据预处理阶段 在进行任何深入的数据分析之前,原始读取数据需要经过严格的质控和过滤过程。这一过程中会去除低质量序列以及可能存在的污染片段[^1]。 #### 细胞捕获与表达矩阵构建 通过特定算法识别并分离来自不同细胞的独特分子标签(UMI),从而建立每个样本中各个细胞对应的基因表达量矩阵。此步骤对于后续所有类型的生物学解释至关重要[^2]。 #### 质量控制(QC) 为了确保下游分析的有效性和准确性,在获得初步的表达谱之后还需要执行严格的质量评估。通常包括但不限于计算每条记录下的总reads数、检测到的特征数量等指标来判断哪些单元格应该被保留用于进一步研究[^3]。 #### 归一化(Normalization) 由于实验条件差异和技术噪声的影响,直接比较不同样品间绝对计数值并不合适;因此需采用适当的方法调整这些偏差以便更好地反映真实的生物变化情况。常见的做法有log转换法或是基于大小因子缩放等方式。 #### 可视化降维(Dimensionality Reduction & Visualization) 高维度空间中的点难以直观理解其分布特性,故而引入诸如PCA(主成分分析)、t-SNE 或 UMAP 技术降低特征向量数目至二维或三维平面内展示出来,便于观察群体结构及潜在亚群的存在形式。 #### 集群(Clustering) 利用无监督聚类算法如K-means、层次聚类(Hierarchical Clustering)或者Louvain社区发现模型对相似度较高的对象分组归类,进而揭示出隐藏在内的异质性模式及其背后所代表的功能分化状态。 #### 差异性表达分析(Differential Expression Analysis) 针对已知类别之间开展统计测试找出显著上调/下调的目标列表,并借助GO富集检验等功能预测手段挖掘关联通路机制,最终达到解析疾病发生发展规律的目的。 #### 可变剪接事件探测(Alternative Splicing Detection) 专门设计用来量化个体水平上的RNA加工变异程度工具——BRIE能够实现全转录本范围内的拼接比例测定工作,有助于探索复杂调控网络下蛋白质多样性产生的源头所在。 ```python import scanpy as sc adata = sc.read_10x_h5('filtered_gene_bc_matrices.h5') sc.pp.filter_cells(adata, min_genes=200) sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000) sc.tl.pca(adata, svd_solver='arpack') sc.pl.pca_variance_ratio(adata, log=True) sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40) sc.tl.umap(adata) sc.tl.leiden(adata) sc.pl.umap(adata, color=['leiden']) ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值