搜索
查看: 3679|回复: 1

[R] 非模式基因GO富集分析:以玉米为例+使用OrgDb

[复制链接]

3

主题

4

帖子

38

积分

新手上路

Rank: 1

积分
38
发表于 2017-8-31 20:44:38 | 显示全部楼层 |阅读模式
模式生物做什么都简单,非模式生物则很多缺少注释,没有注释你就没法做,只能是借助于各种软件比如blastgo,自己跑电子注释。但今天要讲的不是这种情况,很多物种还是有注释的,只是你有时候不知道该去那里下载,或者你有数据,却不知道该怎么用!很多的软件都是针对模式生物的,或者针对某一些类型的非模式生物,能够支持多种非模式生物,能够支持用户自己的注释文件的软件相对来讲,就非常少有了,然而clusterProfiler就是这类少有的软件之一。
获得OrgDb
今天要讲的是通过OrgDb来做GO分析,这是clusterProfiler的enrichGO函数所支持的背景注释,Bioconductor自带20个OrgDb可供使用,多半是模式生物,难道我们要做的物种不在这20个里面就不行了吗?显然不是的,clusterProfiler能支持的物种我自己都数不过来。
我们可以通过AnnotationHub在线检索并抓取OrgDb,比如这里以玉米为例:
> require(AnnotationHub)> hub <- AnnotationHub()> query(hub, "zea")AnnotationHub with 2 records# snapshotDate(): 2017-04-25 # $dataprovider: Inparanoid8, ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/# $species: Gibberella zeae, Zea mays# $rdataclass: Inparanoid8Db, OrgDb# additional mcols(): taxonomyid, genome, description,#   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,#   rdatapath, sourceurl, sourcetype # retrieve records with, e.g., 'object[["AH10514"]]'             title                            AH10514 | hom.Gibberella_zeae.inp8.sqlite  AH55736 | org.Zea_mays.eg.sqlite
通过检索,org.Zea_mays.eg.sqlite就是我们所要的OrgDb,可以通过相应的accession number, AH55736抓取文件,并存入了maize对象中,它包含了51097个基因的注释:
> maize <- hub[['AH55736']]> length(keys(maize))[1] 51097
这个OrgDb,包含有以下一些注释信息:
> columns(maize) [1] "ACCNUM"      "ALIAS"       "CHR"         "ENTREZID"    "EVIDENCE"    [6] "EVIDENCEALL" "GENENAME"    "GID"         "GO"          "GOALL"      [11] "ONTOLOGY"    "ONTOLOGYALL" "PMID"        "REFSEQ"      "SYMBOL"     [16] "UNIGENE"转换ID
我们可以使用bitr来转换ID,甚至于直接检索GO注释:
> require(clusterProfiler)> bitr(keys(maize)[1], 'ENTREZID', c("REFSEQ", "GO", "ONTOLOGY"), maize)   ENTREZID         REFSEQ         GO ONTOLOGY1    541612 XP_008648268.1 GO:0009507       CC2    541612 XP_008648268.1 GO:0051537       MF3    541612 XP_008648268.1 GO:0009055       MF4    541612 XP_008648268.1 GO:0046872       MF5    541612 XP_008648268.1 GO:0022900       BP6    541612 NP_001104837.2 GO:0009507       CC7    541612 NP_001104837.2 GO:0051537       MF8    541612 NP_001104837.2 GO:0009055       MF9    541612 NP_001104837.2 GO:0046872       MF10   541612 NP_001104837.2 GO:0022900       BP11   541612 XM_008650046.2 GO:0009507       CC12   541612 XM_008650046.2 GO:0051537       MF13   541612 XM_008650046.2 GO:0009055       MF14   541612 XM_008650046.2 GO:0046872       MF15   541612 XM_008650046.2 GO:0022900       BP16   541612 NM_001111367.2 GO:0009507       CC17   541612 NM_001111367.2 GO:0051537       MF18   541612 NM_001111367.2 GO:0009055       MF19   541612 NM_001111367.2 GO:0046872       MF20   541612 NM_001111367.2 GO:0022900       BPGO富集分析> sample_genes <- keys(maize)[1:100]> head(sample_genes)[1] "541612" "541613" "541614" "541615" "541617" "541618"
这里我只是简单地使用ID列表中前100个ENTREZ基因ID,也可以使用其它的ID,通过借助于bitr进行转换,或者通过给enrichGO指定ID类型(keyType参数)。
有了OrgDb,使用起来,就跟文档中使用人类基因做为例子一样,用法一致,并且也可以通过clusterProfiler所提供的各种可视化函数对结果进行展示:
> require(clusterProfiler)> res = enrichGO(sample_genes, OrgDb=maize, pvalueCutoff=1, qvalueCutoff=1)> res## over-representation test##...@organism      Zea mays #...@ontology      MF #...@keytype      ENTREZID #...@gene      chr [1:100] "541612" "541613" "541614" "541615" "541617" "541618" "541619" ...#...pvalues adjusted by 'BH' with cutoff <1 #...114 enriched terms found'data.frame':    114 obs. of  9 variables: $ ID         : chr  "GO:0004871" "GO:0000155" "GO:0004673" "GO:0016775" ... $ Description: chr  "signal transducer activity" "phosphorelay sensor kinase activity" "protein histidine kinase activity" "phosphotransferase activity, nitrogenous group as acceptor" ... $ GeneRatio  : chr  "9/80" "5/80" "5/80" "5/80" ... $ BgRatio    : chr  "81/14230" "22/14230" "23/14230" "23/14230" ... $ pvalue     : num  6.65e-10 1.21e-07 1.54e-07 1.54e-07 1.90e-07 ... $ p.adjust   : num  7.58e-08 3.60e-06 3.60e-06 3.60e-06 3.60e-06 ... $ qvalue     : num  6.37e-08 3.03e-06 3.03e-06 3.03e-06 3.03e-06 ... $ geneID     : chr  "541618/541625/541627/541634/541636/541638/541641/541642/541663" "541627/541634/541641/541642/541663" "541627/541634/541641/541642/541663" "541627/541634/541641/541642/541663" ... $ Count      : int  9 5 5 5 4 4 4 5 5 5 ...#...Citation  Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.  clusterProfiler: an R package for comparing biological themes among  gene clusters. OMICS: A Journal of Integrative Biology  2012, 16(5):284-287
如果你有表达谱数据,你也可以使用gseGO进行GSEA分析,这里我懒得上网找数据来演示了,用法反正跟文档里的一样,只不过换成了你自己的数据和相应物种的OrgDb对象而已。
如果没有OrgDb怎么办?
必须也是可以做的,这个将在以后讲解!
相关阅读
[size=0em]​

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:GCE计算基因组重复率
下一篇:求大神指导安装EMBOSS
回复

使用道具 举报

0

主题

4

帖子

41

积分

新手上路

Rank: 1

积分
41
发表于 2018-8-30 21:44:57 | 显示全部楼层
居然是Y叔!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 02:09 , Processed in 0.040572 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.