转录组GO富集分析

链接：

https://www.jianshu.com/p/22b7e77b770e

保存的位置：C:/Users/liu/Documents/AppData/.AnnotationHub/annotationhub.sqlite3
下载链接：https://annotationhub.bioconductor.org/metadata/annotationhub.sqlite3

安装与下载：
一：首先安装 AnnotationHub这个包
AnnotationHub是一个包含大量注释信息的数据库，里面有很多物种，以及来源于很多数据库的注释信息。

1. 安装

source("https://bioconductor.org/biocLite.R")
biocLite("AnnotationHub"), AnnotationHub  ##这种方法过时了
BiocManager::install("AnnotationHub")   ##用这种方法安装 ，要是中途下载失败就自己手动下载，然后放到上面保存的位置。 数据地址：https://annotationhub.bioconductor.org

2. 载入这个包

library(AnnotationHub)

3. 建立AnnotationHub对象

ah = AnnotationHub()

到这里安装就算完成了，接下来我们可以按自己的需求来完成。
二： 想要查看AnonotationHub里面包括那些物种：

unique(ah$species）  ##图1，会发现由于包含东西太多，我们可以根据自己想要找的物种，来看这个AnonotationHub里面是否包含我们想要的物种：
ah$species[which(ah$species=="Mus musculus")]   ##图2，筛选我们想要的数据库

图1
图2

可见是包含的，由于我们想要使用clusterProfiler包，这个包只针对含有OrgDb对象的，所以我们继续寻找：
使用R中的qury函数可以看见一些信息：

mu = query(ah,"Mus musculus")

ah

OrgDb属于rdataclass里面的，所以使用

ah[ah$species =='Mus musculus'&ah$rdataclass=='OrgDb']

rdataclass

就找到了，名字为org.Mm.eg.db。
三：我们如何去下载所需信息
我们找到所需要的物种数据库之后，你会得到一个ID，比如说老鼠的OrgDb的注释数据库的ID就是"AH66157"，然后根据这个ID可以进行下载。当然下载方式已经出现过了.

ath <- ah[['AH66157']]  ##下载数据包

ath

小结：
A. AnnotationHub是生物数据的中转站，方面我们搜索目标数据，另一个相似包是biomaRt;

B. 我们通过query,subset等方法（图形界面则是display），逐步从AnnotationHub的metadata筛选到所需数据的ID;

C. 使用[]是查看目标数据的metadata, 使用[[]]用于下载数据;

找到和下载注释数据库只是第一步，学会如何使用这些数据库更加重要。下面附一个bioconductor提供的19个Org包，要是你的物种是其中的,就不用去找了，可以直接用了。

物种库

探索注释数据库：
用class看下它具体是什么类（面向对象编程的概念）

> class(ath)
[1] "OrgDb"
attr(,"package")
[1] "AnnotationDbi"

我们看看AnnotationDbi的5个主要函数

columns(x)： 显示当前对象有哪些数据
keytypes(x)： 有哪些keytypes可以用作select或keys的keytypes参数
keys(x, keytype, ...)：返回当前数据对象的keysselect(x, keys, columns, keytype, ...)：基于keys, columns和keytype以data.frame数据类型返回数据，可以是一对多的关系
mapIds(x, keys, column, keytype, ..., multiVals)： 类似于select，只不过就返回一个列。

columns(ath)

columns

keytypes(ath)

keytypes

keytypes告诉我们可以当做哪些列是keytype类型，那么keys则列出这个keytype下有哪些关键字。

head(keys(ath,keytype = "SYMBOL"))

行

select(ath, keys= "AGO1", columns=c("TAIR","GO"),keytype = "SYMBOL")
'select()' returned 1:many mapping between keys and columns

select则是根据你提供的key值去查找注释数据库，返回你需要的columns信息。

select

可以看到一个Pzp就能返回那么多信息.因为一个基因可以有多个功能（GO注释），当然一个GO注释下也有可以多个基因。

select(ath, keys= "GO:0004866", columns=c("PATH"),keytype = "GO")

select

富集分析就是看不同GO，KEGG注释下，你提供的基因集的分布情况。比如说我随机从l老鼠中抽样200个基因，然后观察这些基因的富集情况。
注：这里用的Y叔的clusterProfiler包

library("clusterProfiler")
tair.sample <- sample(keys(ath,keyType = "ENTREZID"), 100)
test <- enrichGO(gene         = tair.sample,
                 OrgDb         = ath,
                 keyType = "ENTREZID",
                 pAdjustMethod = "none",
                 pvalueCutoff  = 0.1,
                 qvalueCutoff  = 0.2)
summary(test)

summary

mapIds功能和select类似，只不过他返回的是一组向量，而不是数据库。

mapIds(ath, keys = tair.sample, column = c("TAIR"), keytype = "ENTREZID")

到这里就可以画图来展示富集到的GO terms了：
1. 散点图

test <- enrichGO(gene         = tair.sample,
                 OrgDb         = ath,
                 keyType = "ENTREZID",
                 pAdjustMethod = "none",
                 pvalueCutoff  = 0.1,
                 qvalueCutoff  = 0.2)
dotplot(test, showCategory = 10)

image.png

横轴为GeneRatio, 代表该GO term下的差异基因个数占差异基因总数的比例，纵轴为富集到的GO Terms的描述信息， showCategory指定展示的GO Terms的个数，默认展示显著富集的top10个，即p.adjust最小的10个。
图中点的颜色对应p.adjust的值，从小到大，对应蓝色到红色，大小对应该GO terms下的差异基因个数，个数越多，点越大。
2. 柱状图

barplot(test, showCategory = 10)

[图片上传失败...(image-e621bb-1556376670895)]

横轴为该GO term下的差异基因个数，纵轴为富集到的GO Terms的描述信息， showCategory指定展示的GO Terms的个数，默认展示显著富集的top10个，即p.adjust最小的10个。注意的颜色对应p.adjust值，从小到大，对应蓝色到红色。

3.GO有向无环图

调用goplot来实现GO有向无环图的绘制，代码如下

goplot(test)

image.png

4. centpolt图

cnetplot(test, showCategory = 5)

image.png

图中灰色的点代表基因，黄色的点代表富集到的GO terms, 默认画top5富集到的GO terms, GO 节点的大小对应富集到的基因个数。
5. emapplot

对于富集到的GO terms之间的基因重叠关系进行展示，如果两个GO terms系的差异基因存在重叠，说明这两个节点存在overlap关系，在图中用线条连接起来，用法如下

emapplot(test, showCategory = 30)

image.png

每个节点是一个富集到的GO term, 默认画top30个富集到的GO terms, 节点大小对应该GO terms下富集到的差异基因个数，节点的颜色对应p.adjust的值，从小到大，对应蓝色到红色。