RNA-seq学习：No.5富集分析--ORF过表达
1、基础知识
（1）基本概念
富集分析(enrichment analysis)简单来说就是将成百上千个基因、蛋白或者其他分子分到不同的类中，以减少分析的复杂度。比如之前差异分析得到几百个显著差异基因，如果一个一个单独研究未免太复杂，若按照一定的准则将差异基因归类即可较为快速，方便的了解某一类基因的变化情况。
（2）分类标准
分类的标准即人们根据目前研究建立的基因注释库，目前常用的有两个：GO与KEGG；
- GO简单来说GO term共有三种类型 ①细胞组成（cellular component，CC）；②生物过程（biological process，BP）；③分子功能（Molecular Function，MF）。每个go term都由相应的GO annotation来说明该term的详细信息，例如人的go注释文件为org.Hs.eg.db（该类型文件均为.db结尾）。通过GO富集分析可以了解差异基因富集在哪些生物学功能、途径或细胞定位。
- KEGG，Kyoto Encyclopedia of Genes and Genomes，京都基因和基因组百科全书，是系统分析基因功能，联系基因组信息和功能信息的知识库,其中包含有大量的通路（PATHWAY）图。人的KEGG注释文件为“hsa”。KEGG分析的最终结果就是把判断某些基因是都富集到某一通路上。
（3）过表达分析 over-representation analysis
举一例：首先获得一组感兴趣的基因(一般是差异表达基因，前景基因)，假设有10个；其中有4个都归类到某一GO term中或者落在某一通路中（pathway）；而在整个基因组中（假设为100个，背景基因）有30个都对应该GO term中或者落在该通路中（pathway）中。基于此来研究4/10与30/100间是否有统计学差异，即观察的计数值是否显著高于随机，即待测功能集在基因列表中是否显著富集。
（4）分析平台
目前有蛮多不错的网站在线富集分析软件，当然也可以通过R语言的R包实现。这里以众多人推荐的clusterProfiler包为例进行学习
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
library(clusterProfiler)
差异基因数据使用前期基于airway包分析得到的结果(筛选条件为
padj<0.1&abs(log2FoldChange)>2)
mydata=read.table("results.csv",header=TRUE,
sep=",",stringsAsFactors=FALSE)
genelist=mydata$X
genelist #共有316个差异基因
genelist
注意此时的基因名为
ENSEMBL格式，特征是以ENSG00000字段开头的；其它常见的格式还有ENTREZID，为纯数字序列；SYMBOL为字母为主的字符串。
2、go分析
BiocManager::install("org.Hs.eg.db")
library("org.Hs.eg.db")
# 安装加载人类go注释包
go.all <- enrichGO(genelist, OrgDb = org.Hs.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05,
qvalueCutoff = 0.2,keyType = 'ENSEMBL')
# 需要稍等一会
-
enrichGO函数可以支持多种基因名格式，使用keyType =指定下即可。
dim(go.all)
head(go.all)
go.all.df=as.data.frame(go.all)
由返回结果可以看出共归类到97个 go term中，其中BP类较多；关于表格中的部分变量意义--
-
Description：Gene Ontology功能的描述信息; -
GeneRatio：差异基因中与该Term相关的基因数与整个差异基因总数的比值； -
BgRation：所有背景基因中与该Term相关的基因数与所有基因的比值； - 三个value值：一般情况下， pvalue < 0.05 该功能为富集项；p.adjust 矫正后的pvalue；qvalue为对p值进行统计学检验的q值；
-
Count：差异基因中与该Term相关的基因数。
go.all.df
genelist里的差异基因明明有316个，表格中显示好像只有221个？？奇怪
基于上述分析，将结果可视化
（1）柱形图
barplot(go.all,showCategory=20,drop=T)
- 筛选20个p.adjust值最小的GO term；
-
统计量为Count值；颜色程度根据p.adjust；纵轴标签为Description
barplot
（2）点图
dotplot(go.aal,showCategory=20)
- 筛选20个分类到某一term中基因数最多的GO term
-
注意下GeneRatio是与Count值成正比的，可查看之前的定义
dotplot
（3）有向无环图DAG
由于这里的数据只能是富集一个GO通路（BP、CC或MF）的数据，因此重新针对某一类go，再分析一下。
go.BP <- enrichGO(genelist,
OrgDb = org.Hs.eg.db,
ont='BP',
pAdjustMethod = 'BH',
pvalueCutoff = 0.05,
qvalueCutoff = 0.2,
keyType = 'ENSEMBL')
plotGOgraph(go.BP)
# 之前失败，提示说一下两个相关包未找到
# BiocManager::install("topGO")
# BiocManager::install("Rgraphviz")
-
得到的结果图比较大，建议保存成pdf文件查看
plotGOgraph(go.BP)
#觉得默认文字有点小，调大一下~
opar <- par(no.readonly=TRUE)
par.axis(cex=4)
plotGOgraph(go.BP)
par(opar)
局部放大图
如上图DAG图结果--
- 颜色越深，代表该GO term越显著（p值越小），函数默认将最显著的10个term设置成方形；
- 图形内标注信息分别是GOterm号、Description、P.adjust以及差异基因注释到该term数与背景基因注释到该term数的比值；
- 超接近根结点的GO term的概括性越强，越往下，分支的GO term表示的结果更细。
3、kegg分析
由于enrichKEGG()需要输入的基因名格式为ENTREZID，所以需要转换一下，这里使用clusterProfiler包的bitr()函数
gene=bitr(genelist,fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
gene=gene$ENTREZID
基因数有316变成了267，没匹配到？
bitr()
kegg <- enrichKEGG(gene,
organism = 'hsa',
keyType = 'kegg',
pvalueCutoff = 0.05,
pAdjustMethod = 'BH',
minGSSize = 10,
maxGSSize = 500,
qvalueCutoff = 0.2
)
head(kegg)
dim(kegg)
kegg.df=as.data.frame(kegg)
-
结果好像只有富集到两条通路上，由于第一次做，不知道对于我的数据来说结果是否正常。
kegg.df
列的意义可参考go分析得到的表格
结果可视化
（1）柱状图、点图
# 画法同前
barplot(kegg)
dotplot(kegg)
barplot
dotplot
（2）通路图--针对每一个富集到的通路图画的
browseKEGG(kegg,'hsa04213')
-
红色边框的为上调的，绿色边框的为下调的。
部分截图
-
以上就是基于ORP方法，利用两种注释库（go、kegg）进行的分析。其实过表达分析有不少缺点的，比如
（1）仅使用了基因数目信息，而没有利用基因表达水平或表达差异值；而基因筛选条件基于人为选定差异水平（比如log2FoldChange）；
（2）忽略差异不显著，但比较关键的基因；
（3）将基因同等对待，ORA法假设每个基因都是独立的，忽视了基因在通路内部生物学意义的不同(如调控和被调控基因的不同)及基因间复杂的相互作用；
（4）ORA假设通路与通路间是独立的，但这个前提假设是错误的。 -
由于上述的不足，GSEA方法也常为人们的选择，将会在下一次中继续总结。
