基因表达分析（中）-富集分析

回顾

首先对基因表达分析（上）做一个简单的回顾

1. 研究基因表达的有如下工具：RNA-Seq，microarray， qRT-PCR等（欢迎补充）

2. RNA-Seq，microarray一般用在探索性阶段，qRT-PCR用于验证

3. RNA-Seq和microarray由于他们的实验方式不同，导致寻找差异表达基因的统计学方法也不同。其中microarray使用寡核苷酸作为探针进行杂交，基因表达量与亮度正相关，而亮度是一个连续型变量，因此大多认为结果是服从正态分布。而RNA-Seq的测序结果是一条条read，是一种离散抽样过程，因此认为是服从泊松分布。

4. ANOVA和简单线性模型都是广义线性模型的特殊情况。ANOVA是研究名义型解释变量和连续型解释变量的关系，简单线性模式是研究连续型解释变量和连续型解释变量的关系。而广义线性模式没特殊要求。

5. 在3,4的背景下，microarray一般用t检验（两个条件），ANOVA分析（多个条件），最常用limma（线性模型）进行检验。RNA-Seq有许多基于count的R包，如DESeq，DESeq2,（基于负二向分布广义线性模型）

6. 以上要求你每个条件都要有3个重复（目前投稿要求），你要是老板穷，一个重复都不给，那你去Google解决方案吧。

7. 用R作差异表达分析大致分为以下几步：1）根据软件包要求导入数据；2）数据预处理，把那些只有0或1计数结果的基因去掉，提高效率。这一步还可以进行探索性数据分析；3）跑程序，得到结果；4）对结果进行可视化，看看基因聚类等结果，这一步不是必须的，但却是展示数据最好的手段了。

如果这些内容都还有印象，那么我们继续上次没有写完的内容的其中一个部分—基因富集分析

为何要基因富集分析

在基因差异表达分析之后，你得到了好多p值特别小（也就是显著性很高）的基因，那么下一步你想做什么？

• 选择一些基因用于验证？

• 对其中基因进行后续研究？

• 在结果中把这些基因都放在后面？

• 尝试着把所有基因相关的文献都都读读看（劝你放弃这个念头）？

• 欢迎补充

这些想法都是非常顺理成章的，但是不要着急。
首先，差异表达找到的基因往往很多，你简单的粗暴去找每一个基因的详细资料，显然不太现实；
其次，如果我们单纯觉得某一个基因和你研究的课题相关，或者说你其实已经找到了一个有可能的基因（或者你只是希望用一些高大上的实验验证一下）那么这个行为是不是有太多主观性，存在一些偏见。
当然，你觉得基因就是你要找的，可是万一它只是碰巧来打酱油的呢，这不就是很尴尬了。
所以为了让审稿人相信你的结果，你就需要做一个基因富集分析哦。

什么是基因富集分析

基因富集分析(gene set enrichment analysis)是在一组基因或蛋白中找到一类过表达的基因或蛋白。一般是高通量实验，如基因芯片，RNA-Seq，蛋白质组学（质谱结果）的后续步骤。

基因富集分析需要我们提供某一类功能基因的集合用于背景，常用的注释数据库如：

• The Gene Ontology Consortium: 描述基因的层级关系

• Kyoto Encyclopedia of Genes and Genomes: 提供了pathway的数据库。

分析方法

在文献Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges（推荐大家看一遍）作者将研究方法归为三种：

其中第三种方法想的很好就是难度很大。而且贴心的把每一种方法有哪些工具都总结出来了，然而图片太大放不出来，麻烦大家去简书（hoptop）或者下载文献了。

Over-Repressentation Analysis(ORA)

ORA是目前商业化最多的方法。为了说明他的基本思想，我要举一个喜闻乐见的例子：读书无用论。

这是我百度找到搜狐财经一篇文章《大数据告诉你真正的有钱人是什么样》的有钱人的学历分布情况，高学历人群（本科以上，因为本科生太多了）所占的比例是9.4%，其他都是一般学历占90.6%。这时候，有些公众号就可以开始不带脑子的说了，读书没什么用呀，有钱人中都是一般学历的呀，以后读书读到大学就行了，甚至也可以不上本科呀（34.2%本科和本科以下）。
你每年回家总能回去看到有人炫耀说，虽然我有钱，可是读书太少了，都不能和你们读书人比的。你总感觉哪里不对劲，但是却又不太方便说出来。

实际上，这就是因为没有考虑到背景。因为高学历本身人数就不多，当然在有钱人里面的人数也就相应不多了。我们要证明有钱人更多是富集高学历这一部分。

类别	有钱人	整体
高学历	10	50
一般学历	90	950
	100	1000

H0: 是否有钱和学历高无关
Ha: 学历高还是有点用的
然后做一个Fisher精确检验，看看p值。

richer.pop <- matrix(data = c(10,90,50,950),nrow=2) fisher.test(richer.pop, alternative = "greater")    Fisher's Exact Test for Count Data data:  richer.pop p-value = 0.03857 alternative hypothesis: true odds ratio is greater than 1 95 percent confidence interval: 1.052584      Inf sample estimates: odds ratio  2.109244

p值小于0.05，看来我读个博士让我以后有钱概率变大了。

现在将我们上面的有钱人改成我们找到的基因，整体改成所有基因。高学历表示属于目标注释基因集，一般学历就是非注释基因组.我们就是要判断我们找到的基因更多是在目标注释集中。所以你需要列出下表，然后再做一个fisher.test()。

类别	gene list	Genome
in anno group	10	50
not in anno group	290	19950
	300	20000

上述的基本思想就是统计学的白球黑球实验：

在一个黑箱里，有确定数量的黑白两种球，你随机抽取（不放回）M个球中，其中两种球的比例分别是多少？

除了用Fisher精确检验，还有其他统计方法：

• Hypergeometric (fisher精确检验用的就是超几何检验)

• Binomial: 二项分布要求是有放回，无放回要求整体足够大大到可以近似。

• Chi-squared chisq.test(counts)

• Z

• Kolmogorov-Smirnov

• Permutation

ORA的方法就是如此的简单，但是有一个问题，就是你如何确定哪些基因是差异表达的，你还是需要设置一个人为的cutoff, 主观能动性成分有点大。

Functional Class Scoring(FCS)

FCS认为，“虽然个体基因表达改变之后会更多在通路中体现，但是一些功能相关基因中较弱但协调的变化也有明显的影响。”

The hypothesis of functional class scoring (FCS) is that although large changes in individual genes can have significant effects on pathways, weaker but coordinated changes in sets of functionally related genes (i.e., pathways) can also have significant effects

FCS分析方法稍微复杂了一点，他要求的输入是一个排序的基因列表和一个基因集合。MIT， Broad Institute 2007年文献就提供了这一方法的软件”GSEA”

有如下特点：

• 计算所有输入基因集合的分数，而不是单个基因

• 不需要设置cutoff

• 找到一组相关的基因

• 提供了更加稳健的统计框架

GSEA是一款图形化的软件，根据他们提供的教程，然后点呀点，就会得到如下结果。下图就是需要好好理解的部分。

中间从蓝色到红色的过渡“带”表示基因从上调到下调排列（排序可以按照fold change,也可以是p-value)。黑色像条形码的竖线表示该位置的基因属于某个指定通路。绿色有波动的曲线表示富集分数，从0开始计算，属于基因通路增加，不属于则减少。最后看下黑色的条形码是不是富集在一端。

那如何做统计检验呢？

The final step in FCS is assessing the statistical significance of the pathway-level statistic. When computing statistical significance, the null hypothesis tested by current pathway analysis approaches can be broadly divided into two categories: i) competitive null hypothesis and ii) self-contained null hypothesis [3], [18], [22], [31]. A self-contained null hypothesis permutes class labels (i.e., phenotypes) for each sample and compares the set of genes in a given pathway with itself, while ignoring the genes that are not in the pathway. On the other hand, a competitive null hypothesis permutes gene labels for each pathway, and compares the set of genes in the pathway with a set of genes that are not in the pathway。

我们要检验的目标是基因富集在一端是因为于目标通路相关的基因都在一端富集。那么空假设就是，你把找到的基因随便摆放也能看到富集现象。用比较专业的话说就是先生成一个零假设的数据分布，然后观察实际数据在这个零假设分布下，是不是在尾端。

对所有零假设下分布情况作图，可以得到大致分布范围。将这两个比较结果就能得到统计学的显著性。

有一个思考题：如果GESA分析结果和Fisher精确检验出现冲突该相信谁？

一些问题

统计检验的能力是有限的，所以还有很多问题存在解决。

• 我们希望找到生物学显著的基因，但是生物学显著和统计显著两者并不是完全相关

• 无论是ORA还是FCS都对背景（也就是这个物种一共有多少基因）有要求，但是随着我们的研究深入，基因数量会改变。有些软件会直接设置一个很大的背景数，从而让p值很显著，然后我们就开心地用他们的结果。

• 有些基因没有注释，也就是注释缺失，处理方法就是扔（欢迎拍砖）。

• 有一些注释项是其他项的子集。

实战：用clusterProfiler做富集分析

clusterProfiler是Y叔良心之作（当然他的良心之作还有很多），目前支持KEGG在线拉数据，支持DAVID，支持Broad iNSTITUTE Molecular Signatures Databases， 支持GSEA。

你问我支不支持，我当然是支持的啦。所以这里放上Y叔的公众号了。

安装：

source("https://bioconductor.org/biocLite.R") biocLite('clusterProfiler')

这里简单举例如何使用，更多内容见, 这个部分你可以直接过掉，因为我只是跟着Y叔的代码敲了一遍而已。

加载差异表达基因，我这里偷懒就随机挑一些基因名（来自于DOSE包，Disease Ontology Semantic and Enrichment analysis ）出来了。

library(org.Hs.eg.db) data(geneList) gene <- names(geneList)[abs(geneList) > 2] head(gene) [1] "4312"  "8318"  "10874" "55143" "55388" "991"

Y叔为了展示他能够处理不同命名方式的ID，用bitr（来自于clusterProfiler)进行生物学ID转换

gene.df <- bitr(gene, fromType = "ENTREZID",                toType = c("ENSEMBL", "SYMBOL"),                OrgDb = org.Hs.eg.db) head(gene.df)  ENTREZID         ENSEMBL SYMBOL 1     4312 ENSG00000196611   MMP1 2     8318 ENSG00000093009  CDC45 3    10874 ENSG00000109255    NMU 4    55143 ENSG00000134690  CDCA8 5    55388 ENSG00000065328  MCM10 6      991 ENSG00000117399  CDC20

clusterProfiler: ORA

然后对不同命名的ID都做ORA富集分析。

ego <- enrichGO(gene          = gene,                universe      = names(geneList),                OrgDb         = org.Hs.eg.db,                ont           = "CC",                pAdjustMethod = "BH",                pvalueCutoff  = 0.01,                qvalueCutoff  = 0.05) ego2 <- enrichGO(gene         = gene.df$ENSEMBL,                 OrgDb         = org.Hs.eg.db,                 keytype       = 'ENSEMBL',                 ont           = "CC",                 pAdjustMethod = "BH",                 pvalueCutoff  = 0.01,                 qvalueCutoff  = 0.05) ego3 <- enrichGO(gene         = gene.df$SYMBOL,                 OrgDb         = org.Hs.eg.db,                 keytype       = 'SYMBOL',                 ont           = "CC",                 pAdjustMethod = "BH",                 pvalueCutoff  = 0.01,                 qvalueCutoff  = 0.05)

结果会有如下
ID： 基因本体论的ID
Description： 描述
GeneRatio： 在GO词条所占的比例
BgRatio ：在背景所占的比例
pvalue： 假设是正确但是被拒绝的概率
p.adjust： 采用BH方法进行多重试验p值校正
qvalue： Q值=被拒绝但却是正确的概率
geneID： 列出在这个GO词条下的我们提供的基因
count： 数量，如果没有约分，就是GeneRatio的分子了。

所以从GeneRatio：24/199， BgRatio：231/11711可以反推出下表：

类别	gene list	Genome
in anno group	24	199
not in anno group	231	11711
	255	11910

最后再画一张美美的点图,根据GeneRatio所做。

dotplot(ego, showCategory=30)

clusterProfiler: GSEC

clusterProfiler支持GSEC，而且画图也很方便。

gsecc <- gseGO(geneList=geneList, ont="CC", OrgDb=org.Hs.eg.db, verbose=F) head(summary(gsecc) gseaplot(gsecc, geneSetID="GO:0000779"

总结

这次重要是根据上次找到差异表达基因做富集分析，主流方法有两种ORA和FCS。相应的工具也有很多，但我推荐使用Y叔的clusterProfiler包。

推荐阅读：

1.  https://zhuanlan.zhihu.com/p/20901089

2. http://www.biotrainee.com/thread-37-1-1.html

苹果打赏通道