文章目录

描述
学习目标
一、安装
- 数据集
- 读取数据文件
- 加载R包
二、基因组注释
- 数据库
- - 通用数据库
  - 注释用数据库
- 基因组构建
- 访问数据库的工具
- AnnotationHub
- AnnotationDbi
- - org.Hs.eg.db
  - - 为什么有如此多重复？
    - 为什么会出现基因有gene symbol但没有与之相关的ID？
  - EnsDb.Hsapiens.v75
三、功能分析
- Over-representation analysis
- - 超几何检验
- Gene Ontology project
- - GO Ontologies
  - GO术语层级
- clusterProfiler
- - 下载注释文件
  - 创建基因列表
  - 使用ClusterProfiler
  - 可视化clusterProfiler结果
- 其他功能分析方法
- Functional class scoring 工具
- - 使用clusterProfiler和Pathview进行GSEA
- Pathway topology 工具
- 其他工具
- - 共表达聚类
- 功能分析资源
推荐阅读

描述

功能分析方法能帮助我们一窥所得基因背后所蕴含的生物学意义。这些基因列表 (gene list) 可能来自差异表达分析、GWAS分析、蛋白组学分析。不论是何种来源，功能分析可以探索特定的通路或生物学过程是否富集在一列基因中。
本次学习中，将使用过表征分析 (over-representation analysis, ORA) 和功能分类评分 (functional class scoring, FCS) 法来鉴定出与基因列表相关的潜在通路。我们将会使用clusterProfiler R包确定基因列表中是否存在任何 基因本体（gene ontology, GO） 过程的富集，并根据结果生成图。我们还将简要介绍使用 clusterProfiler 通过 基因集富集分析 (gene set enrichment analysis, GSEA) 执行 FCS，然后使用 Pathview R 包进行可视化。

学习目标

使用GO术语来探索富集的过程
解释富集检验的结果
描述不同类别的功能分析工具 (over-representation analysis, functional class scoring, and pathway topology methods)
利用功能分析工具产生关于富集过程/途径的假设

一、安装

首先下载该课程提供的R project
下好后可看到文件名为Functional_analysis.zip，解压然后移动文件夹到合适的位置，打开文件夹后可见以下：

双击.Rproj file，将在RStudio中打开“Functional_analysis” project，如下：

数据集

来自文章 Kenny PJ et al, Cell Rep 2014
研究目的：研究脆性X综合征中多种基因之间的相互作用。脆性X综合征患者FMRP蛋白异常产生，导致认知受损和自闭症样特点。

FMRP“最常见于大脑中，对正常的认知发育和女性生殖功能至关重要。该基因的突变可导致脆性 X 综合征、智力低下、卵巢早衰、自闭症、帕金森病、发育迟缓和其他认知缺陷。” - 来自wikipedia
MOV10 是一种假定的 RNA 解旋酶，在 microRNA 通路的背景下也与 FMRP 相关。

本篇文章的科学假设：FMRP和MOV10具有关联性，可调节部分RNA的翻译。

对本课程，使用从DESeq2 R包分析得到的差异表达数据（过表达Mov10蛋白后相比对照组的差异表达基因）。数据在下载好的Mov10oe_DE_results.csv文件里。

读取数据文件

打开已有的空R script文件Functional_analysis.R，读取差异表达结果文件，将对象命名为res_tableOE （OE为过表达之意）。

## Load libraries
library(tidyverse)## Read in differential expression results
res_tableOE <- read.csv("data/Mov10oe_DE_results.csv", row.names = 1)## Create a tibble
res_tableOE_tb <- res_tableOE %>%rownames_to_column(var="gene") %>% as_tibble()

第一步read.csv()读取后View()一下，row.names = 1参数即将第一列视为行名（观测）；列变量包括"baseMean", “log2FoldChange”, “lfcSE”, “stat”, “pvalue”, “padj”

第二步，将行名单独作为一列，变量命名为"gene"

加载R包

需先下载，再进行加载

# Install CRAN packages
install.packages(c("BiocManager", "devtools", "tidyverse"))# Install Bioconductor packages
BiocManager::install(c("clusterProfiler", "DOSE", "org.Hs.eg.db", "pathview", "AnnotationDbi", "EnsDb.Hsapiens.v75"))# Optional for the lesson:
BiocManager::install(c("gProfileR", "treemap", "SPIA", "stephenturner/annotables"))

## Load libraries one at a time
library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)
library(pathview)
library(tidyverse)
library(AnnotationDbi)
library(EnsDb.Hsapiens.v75)# Optional for the lesson
library(gProfileR)
library(treemap)
library(SPIA)
library(annotables)

二、基因组注释

本节目标：

使用基因组注释数据完成对基因列表的功能分析
理解在不同数据库中储存的信息类型
探索检索基因组注释信息不同R包的优缺点

二代测序结果需要关联基因、转录本、蛋白等与功能或调控信息。为了对基因列表进行功能分析，我们经常需要获得与我们希望使用的工具兼容的基因标识符，这是很关键的。在这里，我们将讨论获取基因注释信息的方法以及每种方法的一些优缺点。

数据库

我们从存储信息的必要数据库中检索有关过程、通路等（涉及基因）的信息。选择的数据库将取决于尝试获取的信息类型。经常查询的数据库包括：

通用数据库

提供有关基因组特征、特征坐标、同源性、变异信息、表型、蛋白质结构域/家族信息、相关生物过程/途径、相关 microRNA 等的综合信息：

Ensembl (use Ensembl gene IDs)
NCBI (use Entrez gene IDs)
UCSC
EMBL-EBI

注释用数据库

Gene Ontology (GO): database of gene ontology biological processes, cellular components and molecular functions - based on Ensembl or Entrez gene IDs or official gene symbols
KEGG: database of biological pathways - based on Entrez gene IDs
MSigDB: database of gene sets
Reactome: database of biological pathways
Human Phenotype Ontology: database of genes associated with human disease
CORUM: database of protein complexes for human, mouse, rat
…

除此之外，还有其他很多未列出的数据库。

基因组构建

在开始搜索这些数据库之前，应该知道使用哪个基因组构建 (genome build) 来生成基因列表，并确保在功能分析期间使用相同的build进行注释。当获得新的基因组build时，基因组特征（基因、转录本、外显子等）的名称和/或坐标位置可能会发生变化。因此，关于基因组特征（基因、转录本、外显子等）的注释是特定于基因组duild的，我们需要确保我们的注释是从适当的资源中获得的。
例如，如果我们使用人类基因组的 GRCh38 构建来量化用于差异表达分析的基因表达，那么我们应该使用相同的genome build GRCh38 构建来在基因 ID 之间进行转换并识别每个基因的注释。

访问数据库的工具

在 R 中，有许多流行的包用于基因/转录水平的注释。这些软件包提供的工具可以利用提供的基因列表，并使用上面列出的一个或多个数据库检索每个基因的信息。

注释工具：用于访问/查询来自特定数据库的注释

Tool	Description	Pros	Cons
org.Xx.eg.db	Query gene feature information for the organism of interest	gene ID conversion, biotype and coordinate information	only latest genome build available
EnsDb.Xx.vxx	Transcript and gene-level information directly fetched from Ensembl API (similar to TxDb, but with filtering ability and versioned by Ensembl release)	easy functions to extract features, direct filtering	Not the most up-to-date annotations, more difficult to use than some packages
TxDb.Xx.UCSC.hgxx.knownGene	UCSC database for transcript and gene-level information or can create own TxDb from an SQLite database file using the GenomicFeatures package	feature information, easy functions to extract features	only available current and recent genome builds - can create your own, less up-to-date with annotations than Ensembl
annotables	Gene-level feature information immediately available for the human and model organisms	super quick and easy gene ID conversion, biotype and coordinate information	static resource, not updated regularly
biomaRt	An R package version of the Ensembl BioMart online tool	all Ensembl database information available, all organisms on Ensembl, wealth of information	-

接口工具：用于访问/查询来自多个不同注释源的注释

AnnotationHub: queries large collection of whole genome resources, including ENSEMBL, UCSC, ENCODE, Broad Institute, KEGG, NIH Pathway Interaction Database, etc.
AnnotationDbi: queries the OrgDb, TxDb, Go.db, EnsDb, and BioMart annotations.

AnnotationHub

AnnotationHub 包为存储在 AnnotationHub Web 服务中的资源提供客户端接口。它允许直接在 RStudio 窗口中查询我们上面提到的大多数大型数据库。虽然我们不会在本课中讨论它，但如果你有兴趣自己探索它，我们会在此处链接材料。

AnnotationDbi

AnnotationDbi 是一个 R 包，它提供了一个接口，用于使用 SQLite 数据存储连接和查询各种注释数据库。 AnnotationDbi 包，作为Bioconductor 包（即 OrgDb、TxDb、EnsDb、Go.db 和 BioMart），提供了一个接口来查询注释资源。在本课中，我们将通过两个不同数据库的示例。为了从此处未涵盖的任何其他数据库中提取数据，可仔细阅读此帮助文档。

org.Hs.eg.db

有大量特定于生物体的 orgDb 包，例如用于人类的 org.Hs.eg.db 和用于鼠的 org.Mm.eg.db，并且可以在此处找到生物数据库的列表。这些数据库最适合转换基因 ID 或获取当前基因组构建 (current genome builds) 的 GO 信息，但不适用于较旧的基因组构建。这些包提供了与包的发布日期相对应的当前版本，并且每 6 个月更新一次。如果某个包不适用于你所感兴趣的生物体，你可以使用 AnnotationHub 创建自己的包。

# Load libraries
library(org.Hs.eg.db)
library(AnnotationDbi)# Check object metadata
org.Hs.eg.db

我们可以通过输入数据库的名称来查看数据库的metadata，包括物种、不同来源信息的最后更新以及来源 url。请注意，该数据库中的 KEGG 数据最后一次更新是在 2011 年，因此可能不是该信息的最佳站点。

OrgDb object:
| DBSCHEMAVERSION: 2.1
| Db type: OrgDb
| Supporting package: AnnotationDbi
| DBSCHEMA: HUMAN_DB
| ORGANISM: Homo sapiens
| SPECIES: Human
| EGSOURCEDATE: 2021-Sep13
| EGSOURCENAME: Entrez Gene
| EGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
| CENTRALID: EG
| TAXID: 9606
| GOSOURCENAME: Gene Ontology
| GOSOURCEURL: http://current.geneontology.org/ontology/go-basic.obo
| GOSOURCEDATE: 2021-09-01
| GOEGSOURCEDATE: 2021-Sep13
| GOEGSOURCENAME: Entrez Gene
| GOEGSOURCEURL: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA
| KEGGSOURCENAME: KEGG GENOME
| KEGGSOURCEURL: ftp://ftp.genome.jp/pub/kegg/genomes
| KEGGSOURCEDATE: 2011-Mar15
| GPSOURCENAME: UCSC Genome Bioinformatics (Homo sapiens)
| GPSOURCEURL:
| GPSOURCEDATE: 2021-Jul20
| ENSOURCEDATE: 2021-Apr13
| ENSOURCENAME: Ensembl
| ENSOURCEURL: ftp://ftp.ensembl.org/pub/current_fasta
| UPSOURCENAME: Uniprot
| UPSOURCEURL: http://www.UniProt.org/
| UPSOURCEDATE: Wed Sep 15 18:21:59 2021

我们可以使用 AnnotationDbi 轻松地从该数据库中提取信息，方法是：columns、keys、keytypes 和 select。例如，我们将使用我们的 org.Hs.eg.db 数据库来获取信息，但要知道相同的方法适用于 TxDb、Go.db、EnsDb 和 BioMart 注释。

# Return the Ensembl IDs for a set of genes
annotations_orgDb <- AnnotationDbi::select(org.Hs.eg.db, # databasekeys = res_tableOE_tb$gene,  # data to use for retrievalcolumns = c("ENSEMBL", "ENTREZID","GENENAME"), # information to retreive for given datakeytype = "SYMBOL") # type of data given in 'keys' argument

这很容易将我们想要的信息返回给我们，但请注意返回的 warning：‘select()’ returned 1:many mapping between keys and columns. 我们可以检查我们的数据，看看有多少重复的基因条目返回。

which(duplicated(annotations_orgDb$SYMBOL))  %>%length()

为什么有如此多重复？

让我们看一个例子。如果你在 annotations_orgDb 中搜索基因符号 AATF，你将看到一个重复基因的示例。可以看到，除了 Ensembl ID 之外，两者的所有信息字段都是相同的。一个基因映射到 ENSG00000275700，另一个映射到 ENSG00000276072。不同之处在于第一个是对首次基因组组装的映射，而第二个映射是提供替代等位基因的新补丁序列的结果。此外，还有可能导致重复的修复补丁和伪基因。（此处原文：The difference is that the first one is a mapping to the primary assembly, whereas the second is a mapping that is a result of a novel patch sequence providing alternate alleles. There are also fix patches and pseduogenes that can contribute to the duplicates.）
我们没有时间单独检查每个基因，因此对于这些重复，我们将只保留第一个映射。

# Determine the indices for the non-duplicated genes
non_duplicates_idx <- which(duplicated(annotations_orgDb$SYMBOL) == FALSE)# Return only the non-duplicated genes using indices
annotations_orgDb <- annotations_orgDb[non_duplicates_idx, ]

请注意，每列中（ENSEMBL列）还有一些基因具有 NA 值。这些是未注释的基因。

# Check number of NAs returned
is.na(annotations_orgDb$ENSEMBL) %>%which() %>%length()

为什么会出现基因有gene symbol但没有与之相关的ID？

这是因为我们没有使用与分析期间用作参考的构建相匹配的基因组构建(genome build)！我们的数据集是基于人类基因组的 GRCh37/hg19 构建创建的，而 orgDb 仅使用最新的构建 GRCh38。可能某些基因已更改名称，因此在此版本的数据库中不存在。 在整个分析过程中与基因组构建保持一致是非常重要的。 为了获得正确的注释集，我们将使用 EnsDb 数据库。

EnsDb.Hsapiens.v75

要为特定构建生成 Ensembl 注释，也可以使用 AnnotationDbi 轻松查询 EnsDb 数据库。你需要决定要查询的 Ensembl 版本并安装该软件包。此处列出了所有 Ensembl 版本。我们知道我们的数据是针对 GRCh37 的，对应版本 75，所以我们可以安装这个版本的 EnsDb 数据库。
加载库后，我们可以像之前一样查看数据库：

# Load the library
library(EnsDb.Hsapiens.v75)# Check object metadata
EnsDb.Hsapiens.v75# Explore the fields that can be used as keys
keytypes(EnsDb.Hsapiens.v75)

使用 AnnotationDbi 查询数据库，我们可以使用相同的 select 函数为我们的基因列表返回所有gene IDs。请注意，列名略有不同，以匹配数据库中存储的内容。

# Return the Ensembl IDs for a set of genes
annotations_edb <- AnnotationDbi::select(EnsDb.Hsapiens.v75,keys = res_tableOE_tb$gene,columns = c("GENEID", "ENTREZID","GENEBIOTYPE"),keytype = "SYMBOL")

现在，如果我们检查重复项，我们会看到仍然存在许多一对多的映射。

# Check number of duplicates returned
which(duplicated(annotations_edb$SYMBOL))  %>%length()

我们可以删除这些条目，只保留其中一个重复项：

# Determine the indices for the non-duplicated genes
non_duplicates_idx <- which(duplicated(annotations_edb$SYMBOL) == FALSE)# Return only the non-duplicated genes using indices
annotations_edb <- annotations_edb[non_duplicates_idx, ]

可以看到此数据框中的基因少于使用 OrgDb 返回的基因。那我们还有没有注释的gene symbols 吗？

# Check number of NAs returned
is.na(annotations_edb$GENEID) %>%which() %>%length()

此时不再有那些 NA 条目。通过使用正确的注释构建，我们确保基因正确映射到标识符。
许多注释包包含的信息比功能分析所需的信息多得多，我们提取的信息主要用于对不同工具的基因 ID 转换。但是，当对这些包越来越熟悉时鼓励你对它们进行更多探索。

三、功能分析

学习目标：

使用 Gene Ontology 术语确定功能如何归因于基因
了解功能富集工具如何产生有统计意义的富集功能或相互作用的理论
使用over-representation analysis, functional class scoring 和 pathway topology方法讨论功能分析
探索功能分析工具

RNA-seq差异表达分析的输出结果是一个具有显著性的差异表达基因(DEGs)列表。为了获得对DEGs更好的生物学解读，可以进行以下分析：

明确是否能够富集到已知的生物学功能、相互作用或通路
通过基于相似趋势将基因分组在一起来识别基因参与的新通路或调控网络
通过在外部交互数据的背景下可视化所有基因显著上调或下调来使用基因表达的全局变化

通常对于任何差异表达分析，使用免费提供的基于网络和 R 的工具来解释生成的基因列表。虽然功能分析工具涵盖多种技术，但它们可以大致分为三种主要类型：over-representation analysis, functional class scoring, and pathway topology [1]。

功能分析的目标是提供生物学见解，因此有必要在实验假设的背景下分析结果：FMRP 和 MOV10 相互作用并调节部分RNA 的翻译。因此，根据作者的假设，我们可以预期富集到与翻译、剪接和 mRNA 调节相关的过程/通路，后续需要通过实验验证。
请注意，下面描述的所有工具都是验证实验结果和做出假设的好工具。这些工具可能分析出你感兴趣的有关基因/通路；但是，不应该由此对实验过程中涉及的通路直接做出结论。需要对任何由分析表明的通路进行实验验证。

Over-representation analysis

有大量的功能富集工具通过查询包含基因功能和相互作用信息的数据库来执行某种类型的“over-representation”（过表征）分析。
这些数据库通常根据共同的功能，或参与某一途径，或存在于特定的细胞位置，或其他分类，如功能途径等，将基因归为一组（基因集）。从本质上讲，已知的基因被分到已被一致命名的类别中（控制词汇），这些类别是基于基因如何被注释为功能的。这些类别是独立于任何生物体的，然而每个生物体都有不同的分类可用。
为了确定是否有任何类别被过度表征，你可以根据背景集中与同一类别相关的基因比例（相应生物体的基因分类）来确定在你的基因列表中与特定类别相关的基因的观察到比例的概率。

判断某物是否真的过表征的统计检验是超几何检验。

超几何检验

使用上面第一个功能类别的例子，超几何分布是一个概率分布，它描述了 25 个基因 (k) 与“功能类别 1”相关联的概率，对于我们基因列表 (n=1000) 中的所有基因，来自整个基因组中所有基因的种群（N = 13,000），其中包含与“功能类别 1”相关的 35 个基因 (K)。
k成功的概率的计算公式如下：

该检验将得出每个测试类别的校正后P值（经过多重检验校正）。

Gene Ontology project

最广泛使用的分类法之一是由基因本体项目建立的基因本体（GO）。
“基因本体项目是一项协作努力，旨在解决各数据库对基因产物一致描述的需求”。基因本体论联盟维护着GO术语，这些GO术语被纳入许多流行的动物、植物和微生物基因组库的基因注释中。
研究生物功能或相互作用富集的工具通常使用基因本体 (GO) 分类，即 GO 术语来确定在给定基因列表中是否有任何显著变化。因此，为了最好地使用和解释这些功能分析工具的结果，对 GO 术语本身及其组织有一个很好的理解是有帮助的。

GO Ontologies

为了描述基因和基因产物的作用，GO 术语以与物种无关的方式组织成三个独立的控制词汇表（本体），GO的组织是一个无环树：

生物学过程（Biological process）：是指涉及基因或基因产物的生物学作用，可以包括“转录”、“信号转导”和“细胞凋亡”。生物过程通常涉及起始材料或输入的化学或物理变化。
分子功能（Molecular function）：代表基因产物的生化活性，这些活性可以包括“配体”、“GTP酶”和“转运蛋白”。
细胞组分（Cellular component）：是指基因产物在细胞中的位置。细胞成分可以包括“细胞核”、“溶酶体”和“细胞膜”。

每个 GO 术语都有一个术语名称（例如 DNA repair）和一个唯一的术语登录号（GO:0005125），并且单个基因产物可以与许多 GO 术语相关联，因为单个基因产物“可能在多个过程中起作用，包含执行多种分子功能的结构域，并参与与细胞中其他蛋白质、细胞器或位置的多种相互作用”。

GO术语层级

一些基因产物经过充分研究，有大量关于其生物学过程和功能的数据。然而，其他基因产物关于它们在细胞中的作用的可用数据非常少。
例如，蛋白质“p53”将包含有关其在细胞中的作用的大量信息，而另一种蛋白质可能仅被称为“membrane-bound protein”，没有其他可用信息。
GO 本体被开发用于描述和查询具有不同可用信息级别的生物学知识。为此，GO 本体是松散的层次结构，从一般的“父”术语到更具体的“子”术语。 GO 本体是“松散”分层的，因为“子”术语可以有多个“父”术语。
一些信息较少的基因可能只与一般的“父”术语相关联或根本没有术语，而其他信息量大的基因与许多术语相关联。

clusterProfiler

我们将使用 clusterProfiler 对与我们的重要基因列表相关的 GO 术语进行过表征分析。该工具将差异表达基因列表和背景基因列表作为输入，并使用超几何检验进行统计富集分析。基本参数允许用户选择适当的有机体和 GO 本体（BP、CC、MF）进行检验。
加载包：

# Load libraries
library(DOSE)
library(pathview)
library(clusterProfiler)

下载注释文件

通过右键单击此处并选择“将链接另存为…”将注释下载到data文件夹中。
读取文件：

annotations_ahb <- read.csv("data/annotations_ahb.csv")

为了运行 clusterProfiler GO 过表征分析，我们将把我们的基因名称更改为 Ensembl IDs，因为该工具使用 Ensembl IDs 会更容易一些。此外，对于功能分析中的不同步骤，我们需要 Ensembl 和 Entrez IDs。我们将使用新的基因注释与我们的差异表达结果合并。

## Merge the AnnotationHub dataframe with the results
res_ids <- inner_join(res_tableOE_tb, annotations_ahb, by=c("gene"="gene_name"))

创建基因列表

为了进行过表征分析，我们需要一个背景基因列表和一个显著基因列表。对于我们的背景数据集，我们将使用所有用于测试差异表达的基因（结果表中的所有基因）。对于显著的基因列表，我们将使用 p.adj < 0.05 的基因（如果有许多 DEGs，也可以设置倍数变化阈值）。

## Create background dataset for hypergeometric testing using all genes tested for significance in the results
allOE_genes <- as.character(res_ids$gene_id)## Extract significant results
sigOE <- dplyr::filter(res_ids, padj < 0.05)sigOE_genes <- as.character(sigOE$gene_id)

使用ClusterProfiler

下面开始进行GO富集分析并保存结果：

## Run GO enrichment analysis
ego <- enrichGO(gene = sigOE_genes, universe = allOE_genes,keyType = "ENSEMBL",OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05, readable = TRUE)

注意：可以在此处找到具有可用于 OrgDb 参数的注释数据库的不同生物体。
此外，keyType 参数可能在不同版本的 clusterProfiler 中写为 keytype。
最后，ont 参数可以接受“BP”（生物过程）、“MF”（分子功能）和“CC”（细胞组分）子本体，或所有这三个子本体的“ALL”。

## Output results from GO analysis to a table
cluster_summary <- data.frame(ego)write.csv(cluster_summary, "results/clusterProfiler_Mov10oe.csv")

注意：除了保存ego对象的结果摘要之外，保存对象本身也是很有用的。 save() 函数使可以将其保存为 .rda 文件，例如保存save(ego, file="results/ego.rda")。 save() 的互补函数是函数 load()，例如ego <- load(file="results/ego.rda")。
这是一组有用的功能，因为它使我们能够在特定阶段保留分析并在需要时重新加载它们。有关这些功能的更多信息，请参见此处和此处。

可视化clusterProfiler结果

clusterProfiler 有多种选项可用于查看过表征的 GO 术语。我们将探索点图 (dotplot)、富集图 (enrichment plot) 和类别网络图 (category netplot)。
点图显示了与前 50 个术语（大小）相关的基因数量以及这些术语的 p.adj（颜色）。此图按基因比率显示前 50 个 GO 术语（与 GO 术语相关的基因数 / sig 基因的总数），而不是 p.adj。

## Dotplot
dotplot(ego, showCategory=50)

要保存图形，请单击 RStudio Plots 选项卡中的 Export 按钮并 Save as PDF...。在弹出窗口中，更改：
Orientation: to Landscape
PDF size to 8 x 14，为文本标签提供适当大小的图形

接下来，类别网络图显示了与前五个最重要的 GO 术语相关的基因与与这些术语相关的显著基因的倍数变化（颜色）之间的关系。 GO 术语的大小反映了术语的 p 值，越重要的越大。这个图可以帮助我们提出假设，确定可能对几个受影响最大过程很重要的基因。

## To color genes by log2 fold changes, we need to extract the log2 fold changes from our results table creating a named vector
OE_foldchanges <- sigOE$log2FoldChangenames(OE_foldchanges) <- sigOE$gene## Cnetplot details the genes associated with one or more terms - by default gives the top 5 significant terms (by padj)
cnetplot(ego, categorySize="pvalue", showCategory = 5, foldChange=OE_foldchanges, vertex.label.font=6)## If some of the high fold changes are getting drowned out due to a large range, you could set a maximum fold change value
OE_foldchanges <- ifelse(OE_foldchanges > 2, 2, OE_foldchanges)
OE_foldchanges <- ifelse(OE_foldchanges < -2, -2, OE_foldchanges)cnetplot(ego, categorySize="pvalue", showCategory = 5, foldChange=OE_foldchanges, vertex.label.font=6)

NOTE: You might get an error saying Error in loadNamespace(name) : there is no package called 'ggnewscale', in that case install it using `install.packages(“ggnewscale”) and try plotting the Cnetplot again.

同样，要保存图形，请单击 RStudio Plots 选项卡中的 Export 按钮并 Save as PDF...。将 PDF size更改为 12 x 14 以提供适合文本标签大小的图形。

如果你对不在前五名中的重要生物学过程感兴趣，可以对ego数据集取子集来专门展示这些过程：

## Subsetting the ego results without overwriting original `ego` variable
ego2 <- egoego2@result <- ego@result[c(1,3,4,8,9),]## Plotting terms of interest
cnetplot(ego2, categorySize="pvalue", foldChange=OE_foldchanges, showCategory = 5, vertex.label.font=6)

其他功能分析方法

过表征分析只是一种功能分析方法，可用于梳理你所感兴趣的状况来说很重要的生物过程。其他类型的分析可能同样重要或者有信息价值，包括功能分类评分（functional class scoring）和通路拓扑（pathway topology）方法。

Functional class scoring 工具

FCS工具，比如GSEA最常使用差异表达结果中所有基因的基因水平统计或 log2 倍数变化，然后查看特定生物途径的基因集是否在较大的正或负倍数变化中富集。

FCS 方法的假设是，虽然单个基因的大幅度变化会对通路产生显着影响（并将通过 ORA 方法检测到），但功能相关基因集（即通路）中较弱但协同的变化也会产生显着影响。因此，与其设置任意阈值来识别“重要基因”，不如在分析中考虑所有基因。来自数据集的基因水平统计数据被聚合以生成单个通路水平的统计数据，并报告每个通路的统计显着性。如果差异表达分析只得到一小部分重要的 DE Gs，这种类型的分析可能特别有用。

使用clusterProfiler和Pathview进行GSEA

使用从每个基因的差异表达分析中获得的 log2 倍变化，可以通过 clusterProfiler 和 Pathview 工具进行基因集富集分析和通路分析。
对于使用 clusterProfiler 的基因集或通路分析，测试的是基因集上的协同差异表达，而不是单个基因的变化。 “基因集是预先定义的一群基因，它们在功能上是相关的。常用的基因集包括来自 KEGG 通路、GO术语、MSigDB、Reactome 或共享一些其他功能注释的基因群等。这些基因集的一致性变化常表明机制的改变”。
GSEA帮助生物学家在两种不同的生物学状态 (biological states) 中，判断某一组有特定意义的基因集合的表达模式更接近于其中哪一种。[2]

为了对 KEGG 基因集进行 GSEA 分析，clusterProfiler 要求用我们结果数据集中所有基因的 Entrez IDs 来识别基因。我们还需要在分析之前删除 NA 值和重复项（由于基因 ID 转换）：

## Remove any NA values (reduces the data by quite a bit)
res_entrez <- dplyr::filter(res_ids, entrezid != "NA")## Remove any Entrez duplicates
res_entrez <- res_entrez[which(duplicated(res_entrez$entrezid) == F), ]

最后，提取和命名fold changes：

## Extract the foldchanges
foldchanges <- res_entrez$log2FoldChange## Name each fold change with the corresponding Entrez ID
names(foldchanges) <- res_entrez$entrezid

接下来，我们需要按降序排列fold changes。为此，我们将使用 sort() 函数，该函数将向量作为输入。这与 Tidyverse 的arrange() 形成对比，后者需要数据框作为输入。

## Sort fold changes in decreasing order
foldchanges <- sort(foldchanges, decreasing = TRUE)head(foldchanges)

使用 KEGG 基因集执行 GSEA：

## GSEA using gene sets from KEGG pathways
gseaKEGG <- gseKEGG(geneList = foldchanges, # ordered named vector of fold changes (Entrez IDs are the associated names)organism = "hsa", # supported organisms listed belownPerm = 1000, # default number permutations # 注：新版函数会报错，故无需设置此参数minGSSize = 20, # minimum gene set size (# genes in set) - change to test more sets or recover sets with fewer # genespvalueCutoff = 0.05, # padj cutoff valueverbose = FALSE)## Extract the GSEA results
gseaKEGG_results <- gseaKEGG@result

KEGG pathway的organisms信息列在此处。

有多少通路被富集出来？ 查看：

## Write GSEA results to file
View(gseaKEGG_results)write.csv(gseaKEGG_results, "results/gseaOE_kegg.csv", quote=F)

此处同样可通过save()保存一个.rda文件
注意：不同人会得到不同的 GSEA 结果，因为执行的排列 (permutations) 使用随机重排序。如果我们希望每次运行函数时都用相同的排列（即每次运行函数时都希望得到相同的结果），那么我们可以在运行之前使用 set.seed(123456) 函数，即随机种子。 set.seed() 的输入可以是任何数字，但如果想要相同的结果，则需要使用与输入相同的数字。

探索 GSEA 图富集排名列表中的一种pathway：

## Plot the GSEA plot for a single enriched pathway, `hsa03040`
gseaplot(gseaKEGG, geneSetID = 'hsa03040')

使用 Pathview R package将来自 clusterProfiler 的 KEGG pathway数据集整合到通路图像中：

detach("package:dplyr", unload=TRUE) # first unload dplyr to avoid conflicts## Output images for a single significant KEGG pathway
pathview(gene.data = foldchanges,pathway.id = "hsa03040",species = "hsa",limit = list(gene = 2, # value gives the max/min limit for foldchangescpd = 1))

注意：打印所有重要pathways的 Pathview 图像可以很容易地执行如下：

## Output images for all significant KEGG pathways
get_kegg_plots <- function(x) {pathview(gene.data = foldchanges, pathway.id = gseaKEGG_results$ID[x], species = "hsa", limit = list(gene = 2, cpd = 1))
}purrr::map(1:length(gseaKEGG_results$ID), get_kegg_plots)

除了探索 KEGG 基因集的富集，我们还可以使用GSEA来探索 BP 基因本体术语的富集：

# GSEA using gene sets associated with BP Gene Ontology terms
gseaGO <- gseGO(geneList = foldchanges, OrgDb = org.Hs.eg.db, ont = 'BP', nPerm = 1000,  # 注：同样省略minGSSize = 20, pvalueCutoff = 0.05,verbose = FALSE) gseaGO_results <- gseaGO@resultgseaplot(gseaGO, geneSetID = 'GO:0007423')

GSEA usage:

Uer input: All the genes, ranked by DE, no cutoff
Set of genes with specific annotation involved in coordinated up or down regulation
– Significant as a whole
KS test, multiple hypotheses testing correction
Need to define the set before looking at the data
FDR?

除了对特定gene set的分析，反过来GSEA也可以用于发现两组样本从表达或其它度量水平分别与哪些特定生物学意义的基因集有显著关联，或者发现哪些基因集的表达模式或其他模式更接近于表型A、哪些更接近于表型B。这些特定的基因集合可以从GO、KEGG、Reactome、hallmark或MSigDB等基因集中获取，其中MSigDB数据库整合了上述所有基因集。研究者也可自定义gene set (即新发现的基因集或其它感兴趣的基因的集合)。[2]
clusterProfiler 中还有其他可用于 GSEA 分析的基因集（Disease Ontology, Reactome pathways等）。此外，可以使用特殊的 clusterProfiler 函数提供你自己的基因集 GMT 文件，例如用于 MSigDB 的 GMT，MSigDB能帮助你找到可能与你实验相关的已发表研究，如下所示：

BiocManager::install("GSEABase")
library(GSEABase)# Load in GMT file of gene sets (we downloaded from the Broad Institute for MSigDB)c2 <- read.gmt("/data/c2.cp.v6.0.entrez.gmt.txt")msig <- GSEA(foldchanges, TERM2GENE=c2, verbose=FALSE)msig_df <- data.frame(msig)

GSEA分析似乎与GO分析类似但又有所不同。GO分析更加依赖差异基因，实则是对一部分基因的分析 (忽略差异不显著的基因)，而GSEA是从全体基因的表达矩阵中找出具有协同差异 (concordant differences)的基因集，故能兼顾差异较小的基因。GSEA分析可以将那些GO/KEGG富集分信息中容易遗漏掉的差异表达不显著却有着重要生物学意义的基因包含在内。因此二者的应用场景略有区别。另外GO富集是定性的分析，GSEA考虑到了表达或其它度量水平的值的影响。另外，对于时间序列数据或样品有定量属性时，GSEA的优势会更明显，不需要每个分组分别进行富集，直接对整体进行处理。[2]

Pathway topology 工具

最后一种主要的功能分析技术是通路拓扑分析。通路拓扑分析通常考虑基因相互作用信息以及差异表达分析中的倍数变化和调整后的 p 值，以识别失调的通路。根据工具的不同，通路拓扑工具探索基因如何相互作用（例如激活、抑制、磷酸化、泛素化等）以确定通路水平的统计数据。基于通路拓扑的方法利用基因产物（我们的 DEGs）和其他基因产物之间相互作用的数量和类型来推断基因功能或通路关联。
例如，SPIA (Signaling Pathway Impact Analysis) 工具可用于整合差异表达基因列表、它们的倍数变化和通路拓扑结构，以识别受影响的通路。

其他工具

共表达聚类

共表达聚类通常用于通过基于相似的表达趋势将基因分组在一起来识别参与新途径或网络的基因。当基因参与通路和/或通路本身未知时，这些工具可用于识别通路中的基因。这些工具对具有相似表达模式的基因进行聚类，以创建共同表达基因的“模块”，这些基因通常反映功能相似的基因群。然后可以不同条件或在时间过程实验中比较这些“模块”，以识别生物学相关的途径或网络信息。
可以使用热图可视化共表达聚类，当然应该被视为仅具有suggestive；基因的严格分类需要更好的方法。
工具执行聚类的方式是获取整个表达矩阵并计算成对的共表达式值。然后生成一个网络，我们从中探索拓扑结构以推断基因共调节。例如， WGCNA R 包是一种更复杂的共表达聚类方法。

功能分析资源

g:Profiler - http://biit.cs.ut.ee/gprofiler/index.cgi
DAVID - http://david.abcc.ncifcrf.gov/tools.jsp
clusterProfiler - http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html
GeneMANIA - http://www.genemania.org/
GenePattern - http://www.broadinstitute.org/cancer/software/genepattern/ (need to register)
WebGestalt - http://bioinfo.vanderbilt.edu/webgestalt/ (need to register)
AmiGO - http://amigo.geneontology.org/amigo
ReviGO (visualizing GO analysis, input is GO terms) - http://revigo.irb.hr/
WGCNA - https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/
GSEA - http://software.broadinstitute.org/gsea/index.jsp
SPIA - https://www.bioconductor.org/packages/release/bioc/html/SPIA.html
GAGE/Pathview - http://www.bioconductor.org/packages/release/bioc/html/gage.html

基因功能分析——哈佛大学