需要做一个非模式生物的Go富集，在目前能搜索到的教程里面，大多是通过AnnotationHub检索并抓取OrgDb。但是本次发现这个方法获得的数据库有点问题，
具体过程是下载的org.db和扒下来的某物种全基因组交集太小，导致最后go富集的时候p值太大
按说go富集本来p值就是很容易小的，结果做出的结果0.01都费劲
于是本次发现可以通过uniprot数据库如今强大的注释信息来跳过复杂的过程，直接搞定go富集
这种方法不止可以应用于某物种的独立研究，可以任意设计自己的背景基因信息
首先回顾一下org.db依赖的go富集方法
1.Org.db依赖的GO富集
以鸭子(Anas,Taxa_id:8835)为例
library(clusterProfiler)# 3.2版本以上的R需要biomanager安装
library(biomaRt)
library(AnnotationHub)
hub <- AnnotationHub(localHub=FALSE)
#AnnotationHub在最开始使用的时候会建立本地的缓存
anas_result<-query(hub,'Anas')
unique(anas_result$rdataclass)
>output:
[1] "GRanges" "Inparanoid8Db" "TwoBitFile" "ChainFile" "EnsDb" "SQLiteFile"
[7] "OrgDb"
#这里显示了annotationHub所储存到不同类别的数据库
anas_result<-anas_result[anas_result$rdataclass=='OrgDb']#提取OrgDb类型
anas_result$description
>output:
[1] "NCBI gene ID based annotations about Ananas comosus"
[2] "NCBI gene ID based annotations about Ananas comosus_var._comosus"
[3] "NCBI gene ID based annotations about Ananas lucidus"
[4] "NCBI gene ID based annotations about Anas boschas"
[5] "NCBI gene ID based annotations about Anas domesticus"
[6] "NCBI gene ID based annotations about Anas platyrhynchos_f._domestica"
[7] "NCBI gene ID based annotations about Anas platyrhynchos"
[8] "NCBI gene ID based annotations about Anas olor"
[9] "NCBI gene ID based annotations about Anas atrata"
[10] "NCBI gene ID based annotations about Anas fuligula"
[11] "NCBI gene ID based annotations about Drosophila ananassae"
[12] "NCBI gene ID based annotations about Drosophila annanassae"
#这里就是显示了所收录的orgdb，其中6到10都是鸭子，具体可以调取其中的内容来看分类号，产生的时间
anas_result[6]
>output:
'''
AnnotationHub with 1 record
# snapshotDate(): 2022-04-25
# names(): AH101743
# $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
# $species: Anas platyrhynchos_f._domestica
# $rdataclass: OrgDb
# $rdatadateadded: 2022-04-21
# $title: org.Anas_platyrhynchos_f._domestica.eg.sqlite
# $description: NCBI gene ID based annotations about Anas platyrhynchos_f._domestica
# $taxonomyid: 8839
# $genome: NCBI genomes
# $sourcetype: NCBI/UniProt
# $sourceurl: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, ftp://ftp.uniprot.org/pub/databases/uniprot/current_rele...
# $sourcesize: NA
# $tags: c("NCBI", "Gene", "Annotation")
# retrieve record with 'object[["AH101743"]]'
'''
到了这一步，我们算是找到了可以用的数据库，之后用这个来做go富集就OK
anas_org<-anas_result[["AH101742"]]
随后开始进行富集分析，gene_symbol是已经导入进去的基因名
gene_symbol = mapIds(x = anas_org,keys = gene_symbol,keytype = "ENTREZID",column = 'PMID')
erich.go.BP = enrichGO(gene = GO,
OrgDb = anas_org,
keyType ="GO",
ont = "BP",
pvalueCutoff =0.05,
qvalueCutoff = 0.05)
erich.go.CC = enrichGO(gene = GO,
OrgDb = anas_org,
keyType ="GO",
ont = "CC",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)
erich.go.MF = enrichGO(gene = GO,
OrgDb = anas_org,
keyType ="GO",
ont = "MF",
pvalueCutoff =0.05 ,
qvalueCutoff = 0.05)这两步的意思是：
1. 首先通过orgdb进行更名，统一gene_symbol的名称。orgdb的本质就是一个巨大的命名表，里面是一个基因在各种数据库中一共24个各个类型的名称。
> columns(anas_org)
[1] "ACCNUM" "ALIAS" "CHR" "ENSEMBL" "ENTREZID" "EVIDENCE" "EVIDENCEALL"
[8] "GENENAME" "GID" "GO" "GOALL" "ONTOLOGY" "ONTOLOGYALL" "PMID"
[15] "REFSEQ" "SYMBOL" 2.org.db同时提供了索引背景基因的作用，因为GO富集作为超几何分布检验，无非是个背景基因和前景基因数量进行的一个计算。
问题就是出现在这里，这个非模式生物的org.db仅仅在我所导入的全物种基因中识别了几百个，仅仅是一个交集，这就很尴尬了。在我检查了自己之后，认定这个org.db数据库有点问题。（当然各位大佬要是发现我的步骤哪里不对的话请帮忙指出啊TAT）
2.通过uniprot自定义背景基因
当然，首先推荐一下刘小泽老师的博客，讲自己构建org.db的过程
https://www.jianshu.com/p/9c9e97167377/
刘小泽老师的步骤要下载并比对大量的基因来确定ID，但是我的数据量太大网络也太卡，一个简单的NR数据库一周都搞不定（就是下不完整怎么破）。
这里我们可以绕过org.db，方法是通过uniprot数据库的注释
1.首先，通过分类学找到自己需要的物种
2.点击Browse all浏览所有的基因信息
3.点击Customize columns来调整显示的信息
4.然后选择上GO注释
5.这下就可以看到所有基因的GO ID了
6.之后我们download下所有的信息（以表格的方式可以得到GO注释，别的好像不行）
7.打开R，开始编辑
接下来就是个命名的游戏了
整体思路：
通过enricher这个函数来解决这个问题
enricher函数需要输入三个数据：差异基因，名称转换表，描述信息
1.差异基因：一条基因id
2.名称转换表：两列数据框，第一列是基因id，第二列是对应的goid，只要保证大于差异基因列表就行
3.描述表：三列数据框：第一列是go id，第二列是基因描述，第三列是go分类（MF，CC，BP），是全部背景基因的集合。
注：需要实现根据go分类将名称转换表和描述表进行分类，因为enricher相当于一个做超几何分布的函数，他的运行逻辑就是：读入gene_id——转化goid——描述表中有多少基因富集到某个功能？——差异基因中有多少基因富集到某个功能？——超几何分布
根据这个逻辑，他是不会对goid背后的go分类进行解析的，所以需要事先通过go分类将uni2go和description来分成BP组，CC组，MF组
具体代码如下：
library(tidyr)
library(clusterProfiler)
library(stringr)
###########1.准备gene_symbol###########
gene_symbol <- read.csv('C:/Users/gah20/Desktop/xxx/uniprot-txid8835_CDS_100979.csv',header = F)
gene_symbol<-gene_symbol$V17
###########2.准备uni2go###########
uni2go<-read.csv('C:/Users/gah20/Desktop/xxx/my_data.csv',header = T)
View(head(uni2go,10))
colnames(uni2go)<-c('gene_id','ID')
# length(unique(uni2go$gene_id))
#>[1] 100979
uni2go <-uni2go %>% as_tibble() %>%
separate_rows(ID, sep = ";")
remove_blank<-function(x){str_replace_all(x," ", "")}#这里是发现有些GO前面是带有空格，删除空格
uni2go$ID<-lapply(uni2go$ID,remove_blank)
###########3.准备description文件########
dp<-read.csv('C:/Users/gah20/Desktop/xxx/description.csv',header = T)
dp<-unlist(strsplit(dp$锘縂ene.Ontology..GO.,split = ';'))#这个下下来，GOID那行的名称是乱码，就将乱就乱了，总之目的是把go_id1;go_id2;go_id3...这样的格式拆开
dp<-as.data.frame(dp)
dp$GO<-str_extract(dp$dp,"(?<=\\[).*?(?=\\])")
dp$description<-str_extract(dp$dp,".*?(?=\\[)")
dp<-dp[,-1]
###########4.根据BP、MF、CC分割description与MF############
go_class<-read.csv('C:/Users/gah20/Desktop/xxx/description.csv',header = T)
View(head(go_class,10))
#uniprot数据库上下载的数据，是将MF，BP，CC分成了三列，然后把Go_id中的内容分配到其中，所以要用这种方法注释
MF<-go_class$Gene.Ontology..molecular.function.
MF<-MF[which(MF!="")]
MF <-MF %>% as.matrix() %>% as_tibble() %>%
separate_rows(V1, sep = ";")
BP<-go_class$Gene.Ontology..biological.process.
BP<-BP[which(BP!="")]
BP <-BP %>% as.matrix() %>% as_tibble() %>%
separate_rows(V1, sep = ";")
CC<-go_class$Gene.Ontology..cellular.component.
CC<-CC[which(CC!="")]
CC <-CC %>% as.matrix() %>% as_tibble() %>%
separate_rows(V1, sep = ";")
MF<-str_extract(MF$V1,"(?<=\\[)GO.*?(?=])")%>% as.data.frame()
BP<-str_extract(BP$V1,"(?<=\\[)GO.*?(?=])")%>% as.data.frame()
CC<-str_extract(CC$V1,"(?<=\\[)GO.*?(?=])")%>% as.data.frame()
colnames(MF)<-c("GO")
colnames(BP)<-c("GO")
colnames(CC)<-c("GO")
MF2<-subset(dp,GO%in%MF$GO)
BP2<-subset(dp,GO%in%BP$GO)
CC2<-subset(dp,GO%in%CC$GO)
uni2go_MF<-subset(uni2go,ID%in%MF2$GO)
uni2go_BP<-subset(uni2go,ID%in%BP2$GO)
uni2go_CC<-subset(uni2go,ID%in%CC2$GO)
###########4.运行############
Go_rich_MF <- enricher(gene = gene_symbol,
TERM2GENE = uni2go_MF[c('ID','gene_id')],
TERM2NAME = MF2[c('GO','description')],
pvalueCutoff = 0.05,
pAdjustMethod = 'BH',
qvalueCutoff = 0.2,
maxGSSize = 200)
Go_rich_BP <- enricher(gene = gene_symbol,
TERM2GENE = uni2go_BP[c('ID','gene_id')],
TERM2NAME = BP2[c('GO','description')],
pvalueCutoff = 0.05,
pAdjustMethod = 'BH',
qvalueCutoff = 0.2,
maxGSSize = 200)
Go_rich_CC <- enricher(gene = gene_symbol,
TERM2GENE = uni2go_CC[c('ID','gene_id')],
TERM2NAME = CC2[c('GO','description')],
pvalueCutoff = 0.05,
pAdjustMethod = 'BH',
qvalueCutoff = 0.2,
maxGSSize = 200)相关的注释直接写进代码里面了
结语：本人也是处在一个摸索的阶段，各位路过的大佬有什么指点，请不吝赐教~
3333
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
