GO enrichGO using input from g:orth and using biomaRt

[Guerande29]

Hello YuLab-SMU

I am trying to obtain GO for genes from Brassica oleracea, but somethings is wrong.

I was obtained orthologous for Brassica using the tool g:orth from gprofiler.
With this input I am trying to do GO using clusterprofiler, but something is wrong. I think there is lacking some step to fix anything in my input (genes from Brassica) or maybe another step that I unknown.

My input for GO look like this:
pal_BOL green_BOL
Bo8g117760 Bo5g002710
Bo8g114550 Bo8g117670
Bo5g010000 Bo8g117820
Bo8g058660 Bo5g002450
Bo8g107870 Bo8g117890
Bo5g020520 Bo5g002380
Bo5g020860 Bo8g117940
Bo8g106040 Bo8g118140
Bo5g021040 Bo8g118180
Bo8g106000 Bo5g001330
Bo8g066810 Bo8g118290
Bo8g067030 Bo5g003250
Bo5g027600 Bo8g001000
Bo5g034170 Bo5g003340
Bo5g034280 Bo5g003370
Bo7g059440 Bo8g001100
Bo01055s030 Bo8g117320
Bo24568s010 Bo2g158800
Bo5g008910 Bo5g003550

PDTA: gprofiler inf: the "ortholog name" is the preferred gene name in relevant Ensembl Biomart for that gene. The "ortholog ENSG" is the ensembl gene id. It's poorly named due to historical reasons, but we haven't changed it to keep it as compatible with already created pipelines as possible.

This is my code

library(lattice)
library(BiocFileCache)
library(biomaRt)
library(dplyr)
library(ggplot2)
library(RSQLite)
library(devtools)
library(dbplyr)

#listMarts()

head(biomaRt::listMarts(host = "https://plants.ensembl.org"))

https://plants.ensembl.org

I need use "plants_mart"

#create an object that stores the "plants_mart" information
db <- biomaRt::useMart(biomart = "plants_mart", dataset = "boleracea_eg_gene", host = "https://plants.ensembl.org")

#get the list of data sets available in the mart:
datasets <- biomaRt::listDatasets(db)

dataset #### description #### version

boleracea_eg_gene #### Brassica oleracea genes (BOL) #### BOL

attributes <- biomaRt::listAttributes(db)
head(attributes)

keytype: ensembl_gene_id (Gene stable ID)

setwd()

#***** gen ontology *****#
library(AnnotationDbi)
library(clusterProfiler)
library(forcats)
library(enrichplot)
library(pathview)
library(data.table)
library(ggplot2)
library(GOsummaries)
library(DOSE)
library(tidyverse)

uploading universe file

here my universe https://bit.ly/40t8gMG

universe = read.delim("universe_Bol.txt", header = T)
universe<-as.character(universe[,1])
universe <- sort(universe, decreasing = TRUE)
head(universe)

uploading File with genes to analyze

gene <- read.delim("GO.txt", header = T)
head(gene)

ORA

ORA <- compareCluster(geneClusters=gene,
enrichGO,
OrgDb = db,
keyType = "ensembl_gene_id",
ont = "BP",
pAdjustMethod = "BH",
universe = universe,
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = FALSE)

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘mapIds’ for signature ‘"Mart"’

It seems that there is no match between the input and the "db" (from the ensembl plant base).
I also think the problem may be in the "keyType".
I'm not sure how to solve the problem
Could you please help me?

Thank you very much in advance

[guidohooiveld]

Please note that because of the lack of formatting your post is difficult to read. Please put your code in code sections; after selecting all code in your post use this button <> to reformat (and check the preview). Moreover, your universe is not available for download, and it is also unclear what exactly your input is....

Yet, based on the last part of your post, in which you perform an ORA using the function enrichGO, i believe your problem can be solved by using the universal enrichment function enricher together with the arguments TERM2GENE and TERM2NAME. This allows you to use your input data, because the more specific function enrichGO (and enrichKEGG) work only with NCBI-based org.xx.eg.db annotation packages.

This post of mine may be helpful: #588 (comment)

[Guerande29]

Thank you @guidohooiveld for your reply.
My apologies for entering the code incorrectly. You are right, the universal enricher solved the problem.
I used the post 588# and it works very well.
I will leave my corrected code here for whoever needs it.
Thank you so much :)

library(lattice)`
library(BiocFileCache)
library(biomaRt)
library(dplyr)
library(ggplot2)
library(RSQLite)
library(devtools)
library(dbplyr)

## *** PART 1: check available plant marts
listMarts( host="https://plants.ensembl.org" )

## connect to the mart database
EPgenes = useEnsembl(biomart="plants_mart",  host="https://plants.ensembl.org")

## find names of available plant data sets
dsets = listDatasets(EPgenes)
head(dsets)

## in this case, find the brassica oleracea one
dsets[grep("Brassica oleracea", dsets$description),]

## take a note of the dataset name 'boleracea_eg_gene'
## *** PART 2: check available filters and attributes
EPgenes <- biomaRt::useMart(biomart = "plants_mart",  
                            dataset = "boleracea_eg_gene", 
                            host = "https://plants.ensembl.org")

head( listFilters(EPgenes) )
head( listAttributes(EPgenes) )

go.all = getBM(attributes=c("ensembl_gene_id", "go_id", "name_1006", "namespace_1003"), mart=EPgenes) 
head(go.all)
dim(go.all)
#[1] 194991      4

## To get a feeling about the GO annotations:
## check how many of the 194991 entries have a GO annotation
sum(go.all$go_id != "")
#[1] 174260

## check how many unique genes are represented in the 194991 entries?
sum(!duplicated( go.all$ensembl_gene_id) )
#[1] 60586

## remove the genes that don't have a GO annotation.
go.all <- go.all[go.all$go_id != "", ]
dim(go.all)
#[1] 174260      4

## check how many unique genes are represented in the 174260 entries?
sum(!duplicated( go.all$ensembl_gene_id) )
#[1] 39855

### **** PART 4: perform GO ORA analysis GO overrepresentation analysis (ORA) using my input
# I should use the universal enrichment function enricher together with the arguments TERM2GENE and TERM2NAME
# because the more specific function enrichGO (and enrichKEGG) work only with NCBI-based org.xx.eg.db annotation packages.
# note the use of arguments TERM2GENE and TERM2NAME. Their column
# order is important!

#setwd
library(clusterProfiler)    
library(forcats)            
library(enrichplot)         
library(pathview)           
library(data.table)         
library(ggplot2)            
library(GOsummaries)        
library(DOSE)     
library(tidyverse)

# for universe you will use this line:
# universe=go.all$ensembl_gene_id,  
# OR or make your own universe, according to your interest

#### uploading universe file
universe = read.delim("universe_Bol.txt", header = T)
universe<-as.character(universe[,1])
universe <- sort(universe, decreasing = TRUE)
head(universe)

#### uploading File with genes to analyze
gene <- read.delim("GO.txt", header = T)
head(gene)

ORA <- compareCluster(
  geneClusters=gene,
  enricher,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  universe= universe,
  minGSSize = 10,
  maxGSSize = 500,
  qvalueCutoff = 0.05,
  TERM2GENE = go.all[go.all$namespace_1003 == "biological_process",  c("go_id","ensembl_gene_id")],
  TERM2NAME = go.all[go.all$namespace_1003 == "biological_process",  c("go_id","name_1006")] )

## check
as.data.frame(ORA)[1:15,]