R包(OrgDb注释包)的准备
- 我们先去看下牛这个物种的缩写是什么，全部的物种及其简写（三个字母)。可以看到简写为bta
- 然后我们我们去Bioconductor上去搜这个物种的Org类型的注释包，选择最新的包安装就行了
- 搜不到的时候怎么办？别急，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参数控制展示的内容)，大家请自行了解，不再多讲~~
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)
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官网
这里可视化用到的函数和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与差异基因关系网络图
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 )
然后在目录里会直接生成俩张图，一张是该目标物种在这条pathway里存在的基因，一张是差异基因的表达变化图，颜色的改变代表foldchange的变化，方便我们直观的观察基因的上下调~
OK!! 常规的富集分析或者基因集富集(GSEA)流程基本如上所示， 但是各个环节参数调用或更个性化的作图还需自己多练习摸索，另外无参物种或者说非模式物种的基因富集 ，我们留到下次再说吧~~
参考链接：1.https://www.jianshu.com/p/bb4281e6604e?utm_campaign=haruki
