大家好,欢迎来到IT知识分享网。
功能注释后如何做富集分析
本文是为了回答知识星球里的一个提问,他为了用clusterProfiler
做富集分析,打算构建一个OrgDb
,也就是物种数据库。
我之前写过用Bioconductor对基因组注释,用Bioconductor/AnnotationHub
对模式植物的基因进行注释。昨天的推送,我讲过新物种的注释基本上都是基于同源相似性搜索数据库完成,最后得到的就是基因名和数据库中注释的对应关系。Orgdb
是Bioconductor
计划中其中一环,通过构建一个物种各个数据库注释条目和基因的对应关系数据库,方便在得到基因后对基因进行注释。
enrichGO
的前三个参数gene
,OrgDb
,keyType
的目的是利用数据库将基因编号转换成GO号。enrichKEGG
的前三个参数gene
, organism
,keyType
的目的也是为了基于物种名和基因编号直接爬取KEGG,将基因编号转换成KO号。
如果你只是为了做GO和KEGG富集分析,有必要构建物种数据库吗?我的答案是没有必要,因为不构建物种数据库也能够用clusterProfiler
做富集分析。
我相信Y叔一定提供了不通过OrgDb
,将转换基因编号为GO/KO编号,然后做富集分析的方法,所以我就去翻了Y叔为clusterProfiler
写的文档。于是我找到这一篇use clusterProfiler as an universal enrichment analysis tool, 这里面提到了一个通用的函数enricher
用于支持新注释物种.
核心参数两个gene
,TERM2GENE
,前者表示的基因编号,后者是GO/KEGG条目和基因编号的对应关系
enricher(gene, pvalueCutoff = 0.05, pAdjustMethod = "BH", universe,
minGSSize = 10, maxGSSize = 500, qvalueCutoff = 0.2, TERM2GENE,
TERM2NAME = NA)
由于我只拿到了我的KEGG注释,GO注释还在运行中,这次就以KEGG富集分析作为例子。
我从KEGG上拿到的注释是下面这种情况,很明显,有些基因没有注释。这些没有注释的基因应该如何注释?Y叔的建议是不要,全部丢掉,原因去参考资料中找。
CAROC969890.1
CAROC969900.1 K12736
CAROC969910.1 K02943
CAROC969920.1 K13356
CAROC969930.1
CAROC969940.1
CAROC969950.1
CAROC969960.1
CAROC969970.1
CAROC969980.1
简单的grep就可以完成这个剔除工作,grep K query.ko > kegg.tsv
,然后将kegg.tsv
导入到我们的R语言中
gene_ko <- read.table("C:/Users/DELL/Desktop/KEGG.tsv", header = FALSE,
sep = "\t")
然后我们随机抽样几个基因作为gene
输入,同时构建TERM2GENE
的输入
term2gene <- data.frame(TERM=gene_ko$V2, GENE=gene_ko$V1)
gene_sample <- sample(gene_ko$V1, 100)
enkegg <- enricher(gene_sample, TERM2GENE = term2gene, pAdjustMethod = "none")
我这里不用多重实验矫正的原因,因为我是随机抽的基因,很有可能是一个富集都找不到。。所以为了后续演示,就把矫正去掉了,真实情况下,你是要的。
参考资料
免责声明:本站所有文章内容,图片,视频等均是来源于用户投稿和互联网及文摘转载整编而成,不代表本站观点,不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益,请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权,请联系删除,如若转载,请注明出处:https://yundeesoft.com/21851.html