基因功能富集方法和基因注释数据库介绍
1. 两种基因功能富集方法介绍
名词解释：
1. 基因功能富集分析：包括过表达分析ORA、功能分类打分FCS、基于同路拓扑结构PT和基于网络拓扑结构NT和基因富集分析GSEA。常用的是ORA和GSEA。
2. ORA（Over Representation Analysis）：也就是我们熟知的传统的基因富集方法，是一种超几何分布检验方法（Fisher's Exact Test），也就是常见的2x2方法。ORA是一种广泛使用的方法，用于确定已知的生物学功能或过程是否在实验衍生的基因列表中被过度表达（丰富），例如，差异表达基因列表（DEGs）。
3. GSEA（Gene Set Enrichment Analysis）：ORA方法只考虑显著差异基因，这种策略会因过高的阈值而忽略变化较小的基因，GSEA直接解决了这一局限性。所有的基因均可用于GSEA，GSEA聚合了一个基因集内每个基因统计数据，因此能够以一种小而协调的方式检测预定义的基因集中所有基因的变化情况。
4. 基因功能注释数据库：收集了各种物种基因的功能注释的数据库，包括GO、KEGG、Reactome和MSigDB等常见功能注释数据库
1.1 传统的功能富集方法：
N为所有基因中具有pathway/GO term注释的基因数目
n为N中差异表达基因的数目
M为所有基因中注释到某特定pathway/GO term条目的基因数目
m为注释到某特定pathway/GO term条目的差异表达基因数目
超几何分布（hypergeometric）是统计学上的一种离散概率分布。它描述了由有限物体中抽出n个物体，成功抽出制定种类的物件的个数（不归还）。超几何分布和Fisher's Exact Test是一摸一样的原理，只是两个不同的称谓。
1.2 GSEA
传统方法输入的是显著富集的差异基因list，需要经常尝试不同的cutoff值，并将上调和下调基因分开富集。
而GSEA输入的是所有的差异基因和对应的fold change值，基于fc去做统计学检验，计算p值。
1.2.1 GSEA介绍
GSEA是一种基于基因集（也可以理解为某个term）的富集分析方法。
基于基因表达数据与表型的关联度（也可以理解为表达量的变化，比如log2FC）的大小进行排序。
判断没个基因集内的基因是否富集于表型相关度排序后的基因列表的上部或下部，从而判断此基因集内基因的协调变化对表型变化的影响。
1.2.2GSEA原理和步骤
- 基因排序
如上图左边的热图所示，GSEA分析的第一步就是利用所有基因的表达数据，然后计算每个基因在两个分组（或者表型）ClassA、ClassB中的差异度（GSEA提供了6种算法，默认方法是signal2 noise，GSEA官网有提供公式），然后按照在两个表型中的差异度从大到小排序，形成一个排好序的基因列表。
- 分析基因集是否富集
这里的基因集，是事先根据功能或者其他的一些原理把很多的基因分类成不同的基因集合，具体一个基因集可以是某一个通路或者go term中的所有基因，也可以是一个miRNA靶标对应的多个基因等等。GSEA提供了多个分类基因集，在分析数据时，只需要选择不同基因集就可以，当然也可以自己制作基因集。我们可以对每一个小的基因集（GeneSet ）里面的基因对应一下上一步排序表里面的位置，例如上图中的GeneSet1 （一个箭头代表一个基因），看基因集里的成员在基因列表里面的分布情况是否均匀，例如GeneSet1就在基因列表中均匀分布，GeneSet2里面的成员主要分布在基因列表的顶部，GeneSet3里面的成员主要分布在基因列表的底部。如果基因集中的成员在基因列表中均匀分布，说明这个基因集不在这两个表型中富集。如果基因集中的成员在基因列表的顶端例如图中的GeneSet2，说明这个基因集在第一个表型ClassA中富集。如果基因集中的成员在基因列表的底部例如图中GeneSet3，说明这个基因集在第二个表型ClassB中富集。
- 计算基因集的ES值
GSEA分析的第三步就是使用加权法计算基因集的ES值（enrichment score），对位于中部（与性状相关性低）的部分采用较小的权值，所以越集中在两端，与表型的相关性越高。ES曲线最大值为富集分数（Enrichment Score）。
- Permutation test
评估富集得分ES的显著性，通过基于表型而不改变基因之间关系的排列检验Permutation test计算观察到的富集得分ES出现的可能性，计算p值和FDR值。
Permutation test原理：
1）假设全部基因集L有10000个基因，包含某通路基因50个（目标基因集S），GSEA计算得到目标基因集富集分数ES=3。
2）采用随机抽样法，从10000个基因中抽取50个基因，构建虚拟目标基因集A，计算A的ES值。随机抽样10000次，从而得到虚拟ES值10000个。
3）根据真实ES值在10000个虚拟ES值中的位置，如果位于两端极值位置（例如，前top1%或后top1%），则说明S所计算出的ES值不是随机获得，而是显著偏大（说明基因集X基因的排位靠前）或显著偏小（说明基因集X基因的排位靠后）。真实的ES值对应位置 ，例如属于前第50名，p=50/10000
1.2.3 GSEA计算中几个关键概念：
| 概念 | 含义 |
|---|---|
| 计算富集得分 (ES, enrichment score) | ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是，从基因集L的第一个基因开始，计算一个累计统计值。当遇到一个落在s里面的基因，则增加统计值。遇到一个不在s里面的基因，则降低统计值。每一步统计值增加或减少的幅度与基因的表达变化程度（更严格的是与基因和表型的关联度）是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集，负值ES表示基因集在列表的底部富集。 |
| 评估富集得分(ES)的显著性 | 通过基于表型而不改变基因之间关系的排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少，也可基于基因集做排列检验 (permutation test)，计算p-value。 |
| 多重假设检验矫正 | 首先对每个基因子集s计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score (NES)。随后针对NES计算假阳性率。（计算NES也有另外一种方法，是计算出的ES除以排列检验得到的所有ES的平均值） |
| Leading-edge subset | 对富集得分贡献最大的基因成员。 |
| false discovery rate | 错误发现率 |
1.3 GSEA和ORA方法的比较
GO富集分析是先筛选差异基因，再判断差异基因在哪些注释的通路存在富集，这就涉及到阈值的设定，存在一定主观性并且只能用于表达变化较大的基因，即我们定义的显著差异基因。
GSEA不局限于差异基因，而是从基因集的富集角度出发，理论上更容易囊括细微变化对生物通路的影响。另外，对于时间序列数据或样品有定量属性时，GSEA的优势会更明显，不需要每个分组分别进行富集，直接对整体进行处理（可以类比WGCNA分析）。
| 方法 | GSEA | ORA |
|---|---|---|
| 原理 | 统一通路所有的基因是否在某些指标中靠前或靠后 | 差一基因集中某条通路所占比例是否与背景间存在显著区别 |
| 特点 | 没有固定阈值，不受限制；原理略复杂 | 筛选的阈值必须合适；原理简单，适合大规模流程化分析 |
| 显著性检验方法 | Premutation test（p值会有随意性，可通过增加抽样次数减少差异） | 超几何分布检验（Fisher's Exact Test） |
| ⚠️使用的基因 | 全部基因列表 | 差异基因集 |
| 分析是否依赖差异基因集的选取 | 不依赖 | 依赖 |
2. 多种基因注释数据库介绍
2.1 GO数据库
网址：http://geneontology.org/
GO（Gene Ontology）数据库由基因本体论联合会建立，该数据库将全世界所有与基因相关的研究结果进行分类汇总。对不同数据库中关于基因和基因产物的生物学术语进行标准化，对基因和蛋白功能进行统一的限定和描述。基因本体论定义了用来描述基因功能的概念/类，以及这些概念之间的关系。GO术语组织在一个有向无环图中，其中术语之间的边表示父子关系。
它将功能分为三个方面：BP（Biological Process，生物学过程）、CC（Cellular Component，细胞元件）和MF（Molecular Function，分子功能）。在这三个大分支下面又分很多小层级（level），level级别数字越大，功能描述越细致。通过GO注释，可以大致了解 某个物种的全部基因产物的功能分类情况。
2.2 KEGG数据库
网址：http://www.genome.jp/kegg
KEGG（Kyoto Encyclopedia of Genes and Genomes）是由日本京都大学和东京大学联合开发的数据库，是基因组测序和其他高通量实验技术生成的大规模分子数据集的整合和解读的参考知识库，可以用来查询代谢途径、酶（或编码酶的基因）、产物等，也可以通过BLAST比对查询未知序列的代谢途径信息。
KEGG是一组人工绘制的代表分子相互作用和反应网络的通路图。这些途径涵盖了广泛的生化过程，可分为7大类：新陈代谢、遗传和环境信息处理、细胞过程、机体系统、人类疾病和药物开发。
2.3 Reactome数据库
网址：https://reactome.org/
Reactome数据库是一个免费开源的通路数据库，提供直观的生物信息学工具，用于可视化，解释和分析途径相关知识，以支持基础研究，基因组分析，建模，系统生物学研究等。Reactome利用PSIQUIC Web服务来覆盖Reactome功能交互网络和外部交互数据库（如IntAct，ChEMBL，BioGRID和Reflndex）的分子交互数据。Pathway注释由生物学专家和Reactome编辑人员合作编写，并交叉引用许多生物信息学数据库。
该数据库目前覆盖19个物种的通路研究，包括经典的代谢通路、信号转导、基因转录调控、细胞凋亡与疾病。数据库引用了100多个不同的生物信息学资源库，包括NCBI、Ensembl、UniProt、UCSC基因组浏览器、ChEBI小分子数据库和PubMed文献数据库等。
2.4 MSigDB数据库
MSigDB数据库网址：https://www.gsea-msigdb.org/gsea/index.jsp
2015年的文献
MSigDB数据库是由Broad Institute研究所的科学家提出GSEA富集方法的同时提供的基因集数据库。
该数据库从位置、功能、代谢途径和靶标结合等多种角度出发，构建出许多的基因集合。目前包括H和C1-C8这九个系列的基因及，可供下载以及R软件包（msigdbr）载入，以用于富集分析。
msigdbr官网使用说明：https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html
| 基因集 | 说明 |
|---|---|
| H | 癌症特征基因集合，共50 gene sets，最常用。 |
| C1 | 染色体位置基因集合，共299 gene sets，用的很少。 |
| C2 | 包含了已知数据库，文献和专家支持的基因集信息，包含5529 gene sets。 |
| C3 | 调控靶基因集合，包括miRNA-target以及转录因子-target调控模式，3735 gene sets。 |
| C4 | 计算机软件预测出来的基因集合，主要是和癌症相关的基因，858 gene sets。 |
| C5 | Gene Ontology对应的基因集合，10192 gene sets。 |
| C6 | 致癌基因集合，大部分来源于NCBI GEO发表芯片数据，189 gene sets。 |
| C7 | 免疫相关基因集，4872 gene sets。 |
| C8 | single cell identitiy gene sets, 302 gene sets |
msigdbr的使用：
# 1. 查看包含的物种（包含了11个物种）
library(msigdbr)
msigdbr_species()
# 2. 查看特定物种的基因集，比如小鼠
all_gene_sets = msigdbr(species = "Mus musculus")
head(all_gene_sets)
# 3. 检索某个物种特定的基因集，如小鼠的hallmark基因集
h_gene_sets = msigdbr(species = "Mus musculus", category = "H")
head(h_gene_sets)
# 4. 检索某个物种特定的基因集中的某个亚集，如检索鼠类C2 (curated) CGP (chemical and genetic perturbations)基因集
cgp_gene_sets = msigdbr(species = "Mus musculus", category = "C2", subcategory = "CGP")
head(cgp_gene_sets)
# 5. 查看所有可以用的基因集
msigdbr_collections()
# 6. msigdbr()的输出，常规的R数据框操作都可以进行
m_df = msigdbr(species = "Mus musculus") %>% dplyr::filter(gs_cat == "H")
head(m_df)
msigdbr通路富集分析
# clusterProfiler (for genes as Entrez Gene IDs).
m_t2g = m_df %>% dplyr::select(gs_name, entrez_gene) %>% as.data.frame()
enricher(gene = gene_ids_vector, TERM2GENE = m_t2g, ...)
# clusterProfiler (for genes as gene symbols).
m_t2g = m_df %>% dplyr::select(gs_name, gene_symbol) %>% as.data.frame()
enricher(gene = gene_symbols_vector, TERM2GENE = m_t2g, ...)
-# fgsea
m_list = m_df %>% split(x = .$gene_symbol, f = .$gs_name)
fgsea(pathways = m_list, ...)
注意：
msigdbr使用的是MSigDB v7.1 (released March 2020)。
可以从MSigDB 直接下载基因集，GSEABase包getGmt()函数导入数据；但是从MSigDB 的数据主要是人源的，其他物种少一些。
如果通过改变基因大小写来进心人源和鼠的基因转化，会有一小部分基因失败。
不使用msigdbr包，biomaRt 可以进行物种间同源基因转换；但是会面对一对多的情况。
Ge Lab Gene Set Files 也提供多个物种的GMT文件供下载； WEHI也提供人类和鼠类的MSigDB 基因集，但是数据只是 Entrez ID以及每个基因集是单独的文件。 MSigDF 是基于WEHI 数据的。
2.5 数据库总结
| 数据库 | 适用范围 | 特点 |
|---|---|---|
| GO | 几乎所有物种 | 分为3个方面：BP(生物学过程)、CC(细胞元件)和MF(分子功能)。 |
| KEGG | 几乎所有物种 | 分为7大类：新陈代饿、遗传和环境系信息处理、细胞过程、机体系统、人类疾病和药物开发。 |
| Reactome | 19个物种 | 包括经典的代谢通路、信号转导、基因转录调控、细胞凋亡和疾病 |
| MSigDB | 11个物种 | 包括H和C1-C8九个部分(Collection)，可以直接使用，也可以用其中的一部分进行分析。 |
