最近忙毕业论文,有个植物学好友需要跑水稻的GO分析,我此前没做过,所以就研究了下,记录一下步骤。
朋友提供给了我58个基因。
1、id转换
58个基因为msu的命名格式,需要转换为entrezid,但是没有直接的办法,所以分两步转换。
(1)msu到uniprot的转换
http://structuralbiology.cau.edu.cn/PlantGSEA/
(2)uniprot到entrezid或symbol
https://david.ncifcrf.gov/
结果:共58个基因,第一步剩余54个,第二步剩余28个。
原因分析:朋友提供的基因没被广泛研究,部分并不与库里的基因名匹配。
2、GO分析
(1)数据库下载和选择
#需要的包
library('clusterProfiler')
library(topGO)
library(AnnotationHub)
library(dbplyr)
library(pathview)
#导入数据库
ah <- AnnotationHub()
unique(ah$dataprovider) ##可查看数据注释来源
#AH10561 | hom.Oryza_sativa.inp8.sqlite
#AH80658 | org.Oryza_sativa_(japonica_cultivar-group).eg.sqlite
#AH80659 | org.Oryza_sativa_Japonica_Group.eg.sqlite
#AH80660 | org.Oryza_sativa_subsp._japonica.eg.sqlite
#以上提到的几个数据库,我尝试了后三个,第一个不会用
Ztar_org <- ah[["AH80658"]] ##下载目标物种到org数据
keytypes(Ztar_org)#查看可使用的key
'''
[1] "ACCNUM" "ALIAS" "ENTREZID" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[7] "GID" "GO" "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PMID"
[13] "REFSEQ" "SYMBOL"
'''
ens=read.table("gene.txt")#其中gene.txt使用了david结果里的symbol名,LOC+数字
ens=as.character(ens$V1)#keys要求输入字符向量
#symbol到entrezid的转换
DEG.entrez_id = mapIds(x = Ztar_org, #### 数据库
keys = ens, #####差异基因名字
keytype = "SYMBOL", ####差异基因类型是SYMBOL
column = "ENTREZID")
#GO富集
erich.go.BP = enrichGO(gene = DEG.entrez_id,
OrgDb = Ztar_org,
keyType = "ENTREZID",
ont = "BP",
pvalueCutoff = 0.8,
qvalueCutoff = 0.8)
结果:没有可以map的结果
随便搜寻一个基因
select(Ztar_org, keys= "4344902", columns=c("GID",'REFSEQ','GENENAME','PMID','ALIAS','GO'),keytype = "ENTREZID")
结果
GO为‘NA’