# 转录组入门（8）：富集分析

PS: 下面的内容我直接从我之前的文章里摘录过来了。

## 为何要基因富集分析

• 选择一些基因用于验证？
• 对其中基因进行后续研究？
• 在结果中把这些基因都放在后面？
• 尝试着把所有基因相关的文献都都读读看（劝你放弃这个念头）？
• 欢迎补充

## 什么是基因富集分析

• The Gene Ontology Consortium: 描述基因的层级关系
• Kyoto Encyclopedia of Genes and Genomes: 提供了pathway的数据库。

## 分析方法

### Over-Repressentation Analysis(ORA)

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

100900

H0: 是否有钱和学历高无关
Ha: 学历高还是有点用的

richer.pop <- matrix(data = c(10,90,50,850),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，看来我读个博士让我以后有钱概率变大了。

in anno group1050
not in anno group29019950
30020000

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

the install screen for GSEA

• 计算所有输入基因集合的分数，而不是单个基因
• 不需要设置cutoff
• 找到一组相关的基因
• 提供了更加稳健的统计框架

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

GESA

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。

## 一些问题

• 我们希望找到生物学显著的基因，但是生物学显著和统计显著两者并不是完全相关
• 无论是ORA还是FCS都对背景（也就是这个物种一共有多少基因）有要求，但是随着我们的研究深入，基因数量会改变。有些软件会直接设置一个很大的背景数，从而让p值很显著，然后我们就开心地用他们的结果。
• 有些基因没有注释，也就是注释缺失，处理方法就是扔（欢迎拍砖）。
• 有一些注释项是其他项的子集。

## 富集分析

• Over-Representation Analysis
• Gene Set Enrichment Analysis
• Biological theme comparison

### 前提准备：

deseq2.sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

1. 安装包下载注释数据（rog.HS.eg.db)
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite（"AnnotationHub")
library(AnnotationHub)
ah <- AnnotationHub()
org.hs <- ah[['AH53766']]

### GO富集和GESA

GO分析参考Y叔写的GO analysis using clusterProfiler文档：

enrichGO(gene, OrgDb, keytype = "ENTREZID", ont = "MF",
pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,
minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)

• OrgDb： 物种注释数据库，一般是org开头
• keytpe: gene的命名格式
• ont： 是BP（Biological Process）, CC（Cellular Component）, MF（Molecular Function）。一个基因的功能可以从生物学过程，所属细胞部分，和分子功能三个角度定义

ego <- enrichGO(
gene = row.names(deseq2.sig),
OrgDb = org.hs,
keytype = "ENSEMBL",
ont = "MF"
)

dotplot(ego,font.size=5)
plotGOgraph(ego)

GO图

gseGO(geneList, ont = "BP", OrgDb, keyType = "ENTREZID", exponent = 1,
nPerm = 1000, minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05,
pAdjustMethod = "BH", verbose = TRUE, seed = FALSE, by = "fgsea")

geneList： 排序数据， 可以根据log2foldchange, 也可以是pvalues

nPerm： 重抽取次数

minGSSize： 每个基因集的最小数目

maxGSSize： 用于测试的基因注释最大数目

genelist <- sig.deseq2\$log2FoldChange
names(genelist) <- rownames(sig.deseq2)
genelist <- sort(genelist, decreasing = TRUE)
gsemf <- gseGO(genelist,
OrgDb = org.hs,
keyType = "ENSEMBL",
ont="MF"
)


gseaplot(gsemf, geneSetID="GO:0004871")


GO:0004871

### KEGG富集分析

KEGG富集分析那家强，必须是Y叔的clusterProfiler。因为它能爬取最新的KEGG在线版数据库，而不是用不再更新的KEGG.db。

enrichKEGG(gene, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05,
pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500,
qvalueCutoff = 0.2, use_internal_data = FALSE)
library(clusterProfiler)
gene_list <- mapIds(org.hs, keys = row.names(deseq2.sig),
column = "ENTREZID", keytype = "ENSEMBL" )

kk <- enrichKEGG(gene_list, organism="hsa",
keyType = "ncbi-geneid",
qvalueCutoff=0.1)
# 气泡图等图限于篇幅不画了，流量有限
dotplot(kk)

## 总结

[1] Guangchuang Yu., Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

[2] Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges

• 版权声明 本文源自 hoptop 整理整理发表
• 除非特殊声明，本站文章均为原创，转载请务必保留本文链接

• hyp

转录组入门（8）：富集分析 学习