涛哥文章系列（25）：识别和比较基因列表的功能图谱

近年来，基因芯片、RNA-Seq和质谱等高通量实验技术可以在系统水平上检测细胞分子。这类分析产生了大量的数据，需要给予生物学上的解释。一种常用的方法是通过在基因维度上进行聚类，根据基因的相似性对不同的基因进行分组。

为了寻找基因之间共享的功能，一种常见的方法是结合生物学知识，如基因本体学(GO)和京都基因和基因组百科全书(KEGG)，以确定一组基因的主要生物学主题。

在聚类分析之后，研究人员不仅要确定特定基因簇是否有共同的主题，而且要比较基因簇之间的生物主题。选择感兴趣的聚类，然后对每个选定的聚类进行富集分析的手动步骤是缓慢而乏味的。为了弥补这一差距，我们设计了clusterProfiler，用于比较和可视化基因簇之间的功能图谱。

使用clusterProfiler时请引用以下文章。

GYu,LGWang,YHan,::,16(5),284-287.

目前，clusterProfiler支持大约有20个物种，如拟南芥、按蚊、牛、犬、鸡、黑猩猩、大肠杆菌K12株、大肠杆菌Sakai株、苍蝇、人、疟原虫、小鼠、猪、大鼠、恒河猴、蠕虫、非洲爪蟾、酵母菌、斑马鱼，这些物种都得到了GO和KEGG分析的支持。分析也支持天蓝色链霉菌和刚地弓形虫。

在clusterProfiler中，groupGO是为在特定级别上基于GO分布进行基因分类而设计的。

require(DOSE)data(geneList)gene-names(geneList)[abs(geneList)2]head(gene)ggo-groupGO(gene=gene,organism="human",ont="BP",level=3,readable=TRUE)head(summary(ggo))

5.1超几何模型

富集分析是识别生物主题的一种广泛使用的方法。在这里，我们使用超几何模型来评估与疾病相关的选定基因的数量是否比预期的要多。为了确定是否有任何术语以高于偶然预期的频率注释指定的基因列表，clusterProfiler使用超几何分布计算p值：

在该等式中，N是背景分布中的基因总数，M是(直接或间接)注释到感兴趣节点的该分布内的基因的数量，n是感兴趣的基因列表的大小，k是该列表中注释到该节点的基因的数量。默认情况下，背景分布是所有有注释的基因。P值按照多重比较来调整，并计算q值用于FDR控制。

5.2基因集富集分析

分析基因表达谱的一种常见方法是识别被认为有趣的差异表达基因。我们之前展示的富集分析是基于这些差异表达的基因。这种方法将发现差异很大的基因，但它不会检测到差异很小的情况，而是在一组相关基因中以协同的方式得到证明。基因集富集分析(GSEA)直接解决了这一限制。所有基因都可以在GSEA中使用；GSEA聚合了一个基因集中的每个基因的统计数据，因此可以检测到预定义的集合中的所有基因以一种小但协调的方式发生变化的情况。因为很可能许多相关的表型差异表现为一组基因的微小但一致的变化。

·富集分数的计算。

·估计ES的显著性水平。

ES的p值是用排列检验计算的。具体地说，我们对基因列表L的基因标签进行置换，并为置换后的数据重新计算基因集合的ES，从而生成ES的零分布。然后，相对于该零分布计算观察到的ES的p值。

·多重假设检验的调整。

当评估整个GO或KEGG基因集时，clusterProfiler调整估计的显著性水平，以考虑多个假设检验，并计算用于FDR控制的q值。

5.3GO富集分析

ego-enrichGO(gene=gene,universe=names(geneList),'',ont="CC",pvalueCutoff=0.01,readable=TRUE)head(summary(ego))

5.4KEGG通路富集分析

kk-enrichKEGG(gene=gene,organism="human",pvalueCutoff=0.01)head(summary(kk))

5.5DO富集分析

疾病本体论(DO)富集分析是以DOSE实施的，请参考说明文件。enrichDO函数对于识别感兴趣基因的疾病关联性非常有用，函数gseAnalyzer函数是为DO的基因集富集分析而设计的。

5.6Reactome途径富集分析

5.7函数调用

groupGO,enrichGO,enrichKEGG,enrichDO和enrichPathway的函数调用是一致的。基因的输入参数为entrezgene向量(人、鼠)或ORF(酵母菌)的IDs，organism必须是支持的物种(如上所述)。

对于基因集富集分析，gseGO、gseKEGG、gseAnalyzer和gsePathway函数需要额外的参数nPerm来指定排列数。

对于GO分析，ont必须被分别指定为“BP”、“MF”和“CC”中的一个，分别代表生物过程、分子功能和细胞成分。在groupGO中，该级别指定用于基因投影的GO级别。

在富集分析中，pvalueCutoff是根据它们的p值和调整后的p值来限制结果。计算Q值以控制假发现率(FDR)。

readable是指示输入基因ID是否映射到基因符号的逻辑参数。

5.8可视化

groupGO、enrichGO和enrichKEGG的输出可以通过条形图、富集图和类别-基因-网络图来可视化。在条形图或饼图中可视化富集结果是非常常见的。我们认为饼图具有误导性，仅提供条形图。

5.8.1条形图

barplot(ggo,drop=TRUE,showCategory=12)

barplot(ego,showCategory=8)

5.8.2富集图

富集图目前用enrichplot包中的emapplot函数可视化，支持超几何检验和基因集富集化分析的结果。

library(enrichplot)emapplot(ego)

5.8.3cnetplot

为了考虑一个基因可能属于多个注释类别的潜在生物学复杂性，并提供数字变化的信息(如果有的话)，我们开发了cnetplot函数来提取复杂的关联。

cnetplot(ego,categorySize="pvalue",foldChange=geneList)

5.8.4gseaplot

基因集富集分析的运行分数及其与表型的关联性可通过gseaplot可视化。

gseaplot(kk2,geneSetID="hsa04145")

5.8.5来自pathview包的pathview

clusterProfiler用户还可以使用来自pathview包的pathview函数来可视化KEGG路径。下面的例子说明了如何可视化“hsa04110”途径，我们在前面的分析中富集了这一通路。

require(pathview)hsa04110-pathview(=geneList,="hsa04110",species="hsa",limit=list(gene=max(abs(geneList)),cpd=1))

欲了解更多信息，请参阅pathview的说明文档。

clusterProfiler还开发了生物主题比较，它提供了一个函数compareCluster，用于自动计算每个基因簇富集的功能类别。

data(gcSample)xx-compareCluster(gcSample,fun="enrichKEGG",organism="hsa",pvalueCutoff=0.05)dotplot(xx)

默认情况下，仅绘制每个群集的前5个(最重要的)类别。用户可以更改参数showCategory以指定要绘制每个集群的多少个类别，如果showCategory设置为NULL，则将绘制所有结果。

默认情况下，点的大小基于其对应的行百分比，用户可以将参数设置为”count”，以基于基因计数进行比较。参数by还可以设置为“rowPercentage”来标准化点大小，因为某些类别可能包含大量基因，并且使那些小类别的点的大小太小而无法比较。默认参数by设置为“geneRatio”，它对应于输出的“geneRatio”列。为了提供完整的信息，我们还提供了当“by”设置为“rowPercentage”时每个类别中已识别的基因数量(圆括号中的数字)和当“by”设置为“geneRatio”时每个簇标签中的基因簇数(圆括号中的数字)，如图3所示。如果点大小基于“count”，则不会显示行数。

p值表示哪些类别更有可能具有生物学意义。绘图中的点根据其相应的p值进行颜色编码。从红色到蓝色的颜色渐变对应于p值递增的顺序。也就是说，红色表示低p值(高度富集)，蓝色表示高p值(低富集)。用参数pvalueCutoff给出的阈值过滤掉p值和调整后的p值，由qvalue估计FDR。用户可以参考[2]中的例子；我们分析了200名患者的公开可用的乳腺肿瘤组织表达数据集(GSE11121，GeneExpressionOmnibus)[6]。我们从差异表达的基因中鉴定出8个基因簇，并利用compareCluster对这些基因簇进行了丰富的生物学过程的比较。

另一个例子如[7]所示，我们使用[8]中描述的方法计算病毒miRNAs之间的功能相似性，并使用compareCluster比较不同病毒调控的重要KEGG通路。

比较功能被设计成一个通用软件包，用于比较任何种类的本体关联的基因簇，不仅该软件包提供的groupGO、enrichGO和enrichKEGG，而且还包括其他生物和生物医学本体，例如，DOSE包中的enrichDO函数和ReactomePA包中的enrichPathway函数与compareCluster在疾病和反应组途径角度的生物学主题比较中工作得很好。更多细节可以在DOSE和ReactomePA的说明文件中找到。