得到了差异基因并进行了一顿操作可视化后，我们可以开始富集分析了，Don't say so much，要富集当然首推Y叔的成功之作-- Clusterprofiler, 因为我的数据物种比对的牛的基因组，也是属于模式物种，用该包去做富集是更为方便~~ ，当然这个包也不仅仅限于模式物种，开发者当然会考虑的比较全面，提供了几个函数去做非模式物种或无参的富集，后面我也会说到~~ OK，开始富集吧！！
R包(OrgDb注释包)的准备
-
我们先去看下牛这个物种的缩写是什么，全部的物种及其简写（三个字母)。可以看到简写为bta
-
然后我们我们去Bioconductor上去搜这个物种的Org类型的注释包，选择最新的包安装就行了
图片.png
-
搜不到的时候怎么办？别急，AnnotationHub包，一个包含大量注释信息的数据库，里面有很多物种及来源于很多数据库的注释信息，我们通过这个包去找到你物种的Org注释信息，制作为标准注释库，就可和模式生物一样使用了
#BiocManager::install("AnnotationHub")
packageVersion('AnnotationHub') #查看版本
library('AnnotationHub')
ah <- AnnotationHub() #这步时间较长，耐心等待，也有时候会加载失败。
#AnnotationHub的数据来源有哪些？
unique(ah$dataprovider)
> #AnnotationHub的数据来源有哪些？
> unique(ah$dataprovider)
[1] "UCSC" "Ensembl" "RefNet"
[4] "Inparanoid8" "NHLBI" "ChEA"
[7] "Pazar" "NIH Pathway Interaction Database" "Haemcode"
.......
#AnnotationHub目前支持哪些物种？我想找的物种在这里面么
unique(ah$species)
[1] "Homo sapiens" "Vicugna pacos"
[3] "Dasypus novemcinctus" "Otolemur garnettii"
[5] "Papio hamadryas" "Papio anubis"
[7] "Felis catus" "Pan troglodytes"
.....
#查看目标物种
ah$species[which(ah$species == "Bos taurus")]
#查看该物种信息
mu = query(ah,"Bos taurus")
mu
#可以看到信息OrgDb属于rdataclass里面的，所以使用 以下方式精确搜索
ah[ah$species == "Bos taurus" & ah$rdataclass =='OrgDb']
AnnotationHub with 1 record
# snapshotDate(): 2019-10-29
# names(): AH75735
# $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
# $species: Bos taurus
# $rdataclass: OrgDb
# $rdatadateadded: 2019-10-29
# $title: org.Bt.eg.db.sqlite
# $description: NCBI gene ID based annotations about Bos taurus
# $taxonomyid: 9913
# $genome: NCBI genomes
# $sourcetype: NCBI/ensembl
# $sourceurl: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, ftp://ftp.ensembl.org/pub/current_fasta
# $sourcesize: NA
# $tags: c("NCBI", "Gene", "Annotation")
# retrieve record with 'object[["AH75735"]]'
关注结果里的names 和 title 俩行 可以看到AH75375对应牛的编号，且存在OrgDb对象
#制作为标准注释库，就可和模式生物一样使用了,使用[[]]用于下载数据 结果就是一个OrgDb类的对象了
Bta.OrgDb <- ah[["AH75375"]]
OK，根据上面最后得到这个Bta.OrgDb这个对象我们可以用之做富集了
- 还是搜不到，但还是想用Org类型的注释包，那就麻烦点看下这个帖子吧，自己构建物种包OrgDb，然后用clusterProfiler富集分析，原理是基于eggnog-mapper结果中的基因ID和Go的对应信息文件制作的，非模式富集原理类似只不过用包的形式把对应文件封装起来，其实个人认为没必要单单为了富集去做个一个Org的包，因为提供了函数可以直接基于gene和Go|Kegg对应文件做富集。
富集分析
转换ID，准备数据文件
我这里就直接在Bioconductor上去下载牛的Orgdb包去做富集了，富集分析可以基于ENSEMBL的ID号，但是建议都转换为ENTREZ ID。
#两列信息就够了
diff_gene_Deseq2 <- diff_gene %>% dplyr::select(GeneID,log2FoldChange)
#或者拿整个差异基因一起进行富集
All_gene_id <- bitr(diff_gene_Deseq2$GeneID, fromType = "ENSEMBL", toType = c("ENTREZID", "SYMBOL"), OrgDb = org.Bt.eg.db)
#保留转换成功的ENTREZID"
diff_gene_Deseq2 <- diff_gene_Deseq2[diff_gene_Deseq2$GeneID%in%All_gene_id$ENSEMBL,]
#转换id
tran_id <- bitr(up_id, fromType = "ENSEMBL", toType = c("ENTREZID", "SYMBOL"), OrgDb = org.Bt.eg.db)
另一种转换转换id方法， toTable函数会提取该包中的ensemble和geneid对应关系
EG2Ensembl=toTable(org.Bt.egENSEMBL) #提取该包的ensemble和geneid对应关系
colnames(diff_gene_Deseq2)[1] <- 'ensembl_id'
results1=merge(diff_gene_Deseq2,EG2Ensembl,by='ensembl_id',all.x=T)
id <- na.omit(results1)
GO功能富集
All <- enrichGO(OrgDb="org.Bt.eg.db",
gene = diff_gene_Deseq2$GeneID,
ont = "all",
keyType = 'ENSEMBL',
pvalueCutoff = 1,
qvalueCutoff = 1,
pAdjustMethod = "BH",
readable= TRUE)
names(attributes(All)) #查看对象包含元素
result <- All@result #其中的 result，即为所有基因的富集结果数据框形式
y = All[All$p.adjust < 0.5, asis = T]
#提供asis参数可以直接过滤enrich对象，这样方便我们不用先提取数据框再去做筛选了
我们也可以直接对于这个对象去使用barplot，dotplot函数作图了
可视化
内置了许多可视化的函数
##常用到的dotplot
dotplot(All, split="ONTOLOGY",showCategory=20,color = 'pvalue')+
facet_grid(ONTOLOGY~.,scale="free")+
scale_size(range=c(2, 10)) +
scale_color_continuous(low='gold', high='green')
dotplot(All, x = ~GeneRatio/BgRatio)
dotplot(All, x = ~ -log(qvalue)) ## 注意我们可以指定x轴的内容
#其他一些可视化
barplot(All,showCategory=20,title="EnrichmentGO")
cnetplot(All, node_label="category") #网络图展示富集功能和基因的包含关系
emapplot(All) #网络图展示各富集功能之间共有基因关系
heatplot(All) #热图展示富集功能和基因的包含关系
goplot(ego) #igraph布局的DAG 有向无环图
plotGOgraph(All) #矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。
建议大家要多去Y叔公众号去学习该包的专题，阅读帖子，了解更细致的作图细节，比如上面dotplot函数默认映射的是按照p.adjust值的大小，加上color = 'pvalue'参数之后就可以按照p值大小去映射，在你的矫正p值都普遍较高时候可以用到，还有我们可以通过scale_size/color等去调点大小，图颜色，就不需要再去把富集结果的表弄出来再用ggplot去画 一步到位省时省力，当然其他的centplot，emapplot函数里也有一些好用的参数(比如node_label参数控制展示的内容)，大家请自行了解，不再多讲~~
图片.png
图片.png
Goplot包富集
除了内置函数可视化，这里再推荐另外一个包Goplot对结果文件进行可视化，图形总体感觉还是比较美观的~
该包文档地址：https://wencke.github.io/
这个包内置了一个EC数据集 方便我们构造数据格式，需要注意的是该包对列名，数据格式要求严格，我们要跟着模仿修改，示例文件中是ENSEMBLID，我这里用基因Symbol ID试一下，这里值简单演示，更多类型多请详细看文档~~
data(EC)
head(EC$david)
head(EC$genelist)
library(GOplot)
## 用自己的数据演示
gene_list <- id[,c(3,4)] #一列Symbol ID 一列 Log(foldchange)
colnames(gene_list) <- c('ID','logFC') #修改列名和内置数据一样
enrich_result <- All@result #提取该基因列表富集的结果
enrich_result <- dplyr::select(enrich_result,ONTOLOGY,ID,Description,geneID,p.adjust)
names(enrich_result) <- c('Category','ID','Term','Genes','adj_pval') #修改列名和内置数据一样
enrich_resulte$Genes <- gsub('/',',',enrich_resulte$Genes)
rownames(enrich_resulte) <- seq(1,nrow(enrich_resulte))
##
circ <- circle_dat(enrich_resulte,gene_list)
GOBubble(circ, title = 'Bubble plot', colour = c('orange', 'darkred', 'gold'), display = 'multiple', labels = 3)
# Colour the background according to the category
GOBubble(circ, title = 'Bubble plot with background colour', display = 'multiple', bg.col = T, labels = 3)
图片.png
KEGG功能富集
#要尽量使用Entrezid
KEGG <- enrichKEGG(gene= tran_id$ENTREZID,
organism = 'bta',
keyType = "kegg",
pAdjustMethod = "BH",
#minGSSize = 10,
#maxGSSize = 500,
qvalueCutoff = 1,
pvalueCutoff = 1)
KEGG1 <- as.data.frame(KEGG)
#barplot
barplot(KEGG, font.size=12, showCategory=20, title="Enrichment Top20")
#标题，字体，泡泡大小
dotplot(KEGG, font.size=8, showCategory=10, title="Enrichment KEGG Top10") + scale_size(rang=c(5,20))
#pathway关系网络图（通过差异基因关联）
emapplot(kk, showCategory = 30)
#pathway与差异基因关系网络图
cnetplot(kk, showCategory = 5)
#pathway映射
browseKEGG(kk, "hsa04934") #在pathway通路图上标记富集到的基因，会链接到KEGG官网
图片.png
这里可视化用到的函数和Go富集分析的基本一致，另外提一点，KEGG除了pathway这个数据库，另一个Module数据库也是比较常用的，大家也可以做下富集看有没有什么新的发现
## KEGG Moudule富集分析
xx <- enrichMKEGG(gene = id$gene_id,
organism='bta',
minGSSize=1)
GSEA富集分析
一般拿全谱数据跑，步骤简单讲就是可以取DEseq2结果中的ID列和log(FC)列数据取出来，然后转换为ENTREZID，从大到小排序，最后用 gseGO和gseKEGG做基因集富集
library(dplyr)
gene_list <- dplyr::select(resdata,Gene,log2FoldChange)
names(gene_list) <- c('ensembl_id','Log2FC')#排序数据， 可以根据log2foldchange, 也可以是pvalues
#转换id
tran_id <- bitr(gene_list$ensembl_id, fromType = "ENSEMBL", toType = c("ENTREZID", "SYMBOL"), OrgDb = org.Bt.eg.db)
#查看没有转换成功的id 也可以选择foldchange
gene_list <- merge(gene_list[gene_list$ensembl_id%in%tran_id$ENSEMBL,],tran_id,by.x = 'ensembl_id',by.y = 'ENSEMBL')%>%
dplyr::select(.,c('ENTREZID','Log2FC')) %>% arrange(desc(Log2FC))
gene_list <- gene_list[!duplicated(gene_list$ensembl_id),] #第一列必须没有重复
genelist <- gene_list$Log2FC
names(genelist) <- gene_list$ENTREZID
head(genelist)
Go富集
ego <- gseGO(
geneList = genelist,
OrgDb = org.Bt.eg.db,
keyType = "ENTREZID",
ont = "ALL",
nPerm = 1000, #置换检验的置换次数
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 1,
verbose = FALSE)
gsea_GO <- ego@result
#查看 enrichResult 类对象中的元素
names(attributes(ego))
#同上文提到的，clusterProfiler 包里的一些默认作图方法，例如
dotplot(ego,color = 'pvalue') #富集气泡图
emapplot(ego) #网络图展示各富集功能之间共有基因关系
gseaplot(ego, 'GO:0002376') #基因集富集得分图
cnetplot(ego, showCategory = 5)#GO term与差异基因关系网络图
图片.png
Kegg富集
kk <- gseKEGG(
geneList = genelist,
keyType = 'kegg',
organism = 'bta',
nPerm = 1000,
minGSSize = 10,
maxGSSize = 500,
pvalueCutoff = 0.05,
pAdjustMethod = "BH"
)
dotplot(kk,color = 'pvalue') #富集气泡图
#pathway映射
browseKEGG(kk, "bta04064") #在pathway通路图上标记富集到的基因，会链接到KEGG官网
关于KEGG的可视化，如果我们想去看具体的通路上基因的表达变化，除了用上面browseKEGG函数去链接官网查看，我还想到一个类似功能且十分强大的R包pathview
library(pathview)
pathview(gene.data = genelist, pathway.id = 'bta04110',species="bta",
limit=list(gene=max(abs(genelist)), cpd=1),
kegg.native = T
)
图片.png
然后在目录里会直接生成俩张图，一张是该目标物种在这条pathway里存在的基因，一张是差异基因的表达变化图，颜色的改变代表foldchange的变化，方便我们直观的观察基因的上下调~
OK!! 常规的富集分析或者基因集富集(GSEA)流程基本如上所示， 但是各个环节参数调用或更个性化的作图还需自己多练习摸索，另外无参物种或者说非模式物种的基因富集 ，我们留到下次再说吧~~
参考链接：1.https://www.jianshu.com/p/bb4281e6604e?utm_campaign=haruki
2.https://www.cnblogs.com/jessepeng/p/12159139.html
网友评论