GO enrichment analysis on Solanum lycopersicum proteomics dataset (UniProt IDs)
1
0
Entering edit mode
Alberto • 0
@Alberto-24625
Last seen 11 days ago
Netherlands

Hello! I am analyzing a proteomic dataset of Solanum lycopersicum that is organized as following:

Treated vs Control - in a time series: TP1; TP2; TP3...

On this dataset I need to perform a GeneOntology enrichment analysis in R, and I was planning to do it with clusterProfiler.

I want to perform both ORA, (enrichment analysis on over represented tomato proteins in treated vs control on each time point) and GESA (Gene set enrichment analysis, that in this case would be the whole proteome set for each time point). Now I have two problems that are kind of related to each other.

1) Solanum lycopersicum is not among the org.db model organism supported by clusterprofiler. I tried to import the AnnotationHub Solanum as following but it doesn't seem to work with clusterprofiler.

hub <- AnnotationHub()
query(hub, c("Solanum"))
org.Sl.eg.db <- hub[["AH80808"]]

2) Since I am starting with a proteomics dataset my query IDs are Uniprot-ID and these seem not supported in the org.db of solanum lycopersicum I am trying to make from AnnotationHub

AnnotationDbi::keytypes(org.Sl.eg.db)
 [1] "ACCNUM"      "ALIAS"       "ENTREZID"    "EVIDENCE"    "EVIDENCEALL" "GENENAME"    "GID"         "GO"          "GOALL"       "ONTOLOGY"    "ONTOLOGYALL" "PMID"  "REFSEQ"      "SYMBOL"      "UNIGENE"

I was trying to play around and/or find the way to use clusterProfiler from non model organism, and see if i can start from this kind of table format

uniprotID   GO  ONTOLOGY

Please any suggestion/tip is very welcome.

clusterProfiler AnnotationHub GO uniprot • 845 views
ADD COMMENT
0
Entering edit mode

Minor comment: the AnnotationHub ID for org.Sl.eg.db is "AH85737" (and not "AH80808")! But this will unfortunately not change your 'problem'...

> query(hub, c("Solanum"))
AnnotationHub with 8 records
# snapshotDate(): 2020-10-27
# $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, Inparanoid8
# $species: Solanum tuberosum, Solanum lycopersicum, Solanum pennellii, Sola...
# $rdataclass: OrgDb, Inparanoid8Db
# additional mcols(): taxonomyid, genome, description,
#   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
#   rdatapath, sourceurl, sourcetype 
# retrieve records with, e.g., 'object[["AH10593"]]' 

            title                                             
  AH10593 | hom.Solanum_lycopersicum.inp8.sqlite              
  AH10606 | hom.Solanum_tuberosum.inp8.sqlite                 
  AH85602 | org.Solanum_pennellii.eg.sqlite                   
  AH85603 | org.Solanum_pennelli.eg.sqlite                    
  AH85660 | org.Solanum_tuberosum.eg.sqlite                   
  AH85736 | org.Solanum_esculentum.eg.sqlite                  
  AH85737 | org.Solanum_lycopersicum.eg.sqlite                
  AH85738 | org.Solanum_lycopersicum_var._humboldtii.eg.sqlite
> 
ADD REPLY
2
Entering edit mode
Guido Hooiveld ★ 3.1k
@guido-hooiveld-2020
Last seen 2 days ago
Wageningen University, Wageningen, the …

Some time ago I had a similar challenge (i.e. having Uniprot ids as main identifier, and wanting to perform GO analysis with clusterProfiler), and in the end I got it working. The key thing is to know about the existence of the Gene Ontology Annotation (GOA) group at the EBI. These are the people that provide the GO annotation to the UniProt database. Very conveniently, they also provide annotation sets for full (inferred) proteomes for download, including tomato!

To find the file for your organism of interest, you should check here (website; https://www.ebi.ac.uk/GOA/proteomes) or here (FTP server; http://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/ ); file for tomato is called 284179.S_lycopersicum_genome_cv_Heinz_1706.goa (direct link (http://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/284179.S_lycopersicum_genome_cv_Heinz_1706.goa)). Note: ftp links don't work anymore in Chome or Firefox!?

The next challenge could be the import of such GOA file, but luckily you can make use of some infrastructure already existing in Bioconductor, for example the package mgsa here.

With respect to clusterProfiler (here), the section on 'Supported organisms' states that when there is no suitable OrgDb object is available, functions (enricher() and gseGO()) can also be run when GO annotation is available in data frame format (with 2 columns). Also notice the statement regarding the difference between direct and indirect GO annotation, and in case of the former the need to use the function buildGOmap() to annotate ancestor GO nodes. Since the function readGAF() also does this (see ?readGAF), there was no need to run buildGOmap() in the code below.

Below find some code that I used (adapted for tomato). Realize that this is one way of doing it, but certainly not the only one. Also check whether you may need to change (default) values of the arguments for each function.

## Prepare GO annotation
library(mgsa)

# read GO annotation file
# Note that I manually downloaded the *.goa file, and is in working directory.
GAF <- readGAF(filename="284179.S_lycopersicum_genome_cv_Heinz_1706.goa")

# extract relevant info.
# unfortunately could not find accessor functions for all required
# info, thus sometimes had to utilize object slots directly (with @)
mapping.index <-  GAF@itemName2ItemIndex
ID.annotations <- itemAnnotations(GAF)

GO.sets <- GAF@sets
GO.annotation <- setAnnotations(GAF)


# create a 2-column data frame with GOID and ID index
# after little further processing, this will be used as input for clusterProfiler
GO.df <- data.frame("GOID" = rep(names(GO.sets), sapply(GO.sets, length)),
        "ID.index" = unlist(GO.sets),  row.names = NULL)


# do some processing for objects GO.annotation and GO.df
# in both remove category 'all',
# and to GO.df also add column with Uniprot ids

# GO.annotation
GO.annotation <- GO.annotation[GO.annotation[,"term"] != "all", ]
GO.annotation[,"GOID"] <- rownames(GO.annotation)

# GO.df
GO.df <- GO.df[GO.df[,"GOID"] != "all", ]
GO.df[,"UNIPROTKB"] <- names(mapping.index [GO.df[,"ID.index"] ])


## objects are now ready for use with clusterProfiler!


## Functional analysis with clusterProfiler.
library(clusterProfiler)

# For function analysis some test input is needed. Let's generate that first.
# For ORA analysis, generate the 'background/universe' list (= all
# proteins identified in the tomato proteome), plus some test input
# (='foreground' list), which usually is identified after performing
# statistics on a data set.

# all proteins
proteins.all <- rownames(ID.annotations)
# random selection of 1750 proteins
proteins.random <- base::sample(proteins.all, 1750)


# Perform GO ORA analysis
# Note the use of arguments TERM2GENE and TERM2NAME. Their column
# order is important!
# Also: to show proof-of-principle no cutoff on significance
# is applied! i.e Cutoff=1

res.GO.ora <- enricher(
    gene=proteins.random,
    pvalueCutoff = 1,
    pAdjustMethod = "none",
    universe=proteins.all,
    minGSSize = 10,
    maxGSSize = 500,
    qvalueCutoff = 1,
    TERM2GENE = GO.df[ ,c("GOID","UNIPROTKB")],
    TERM2NAME = GO.annotation[ ,c("GOID", "term")]
)

# check
as.data.frame(res.GO.ora)[1:15,]

# generate plots
library(enrichplot)

barplot(res.GO.ora, showCategory=10)
dotplot(res.GO.ora, showCategory=10)
cnetplot(res.GO.ora, showCategory=10)
heatplot(res.GO.ora, showCategory=10)
emapplot( pairwise_termsim(res.GO.ora), node_scale=1.5,layout="kk" )


# Also perform KEGG-based analyses.
# note: only 10841 (~52%) uniprot ids have been annotated by KEGG!
uni2kegg <- bitr_kegg( proteins.all, fromType="uniprot",
    toType='kegg', organism='sly')

# KEGG ORA.
kk <- enrichKEGG(
    gene =  proteins.random,
    organism = "sly",
    keyType = "uniprot",
    pvalueCutoff = 1,
    pAdjustMethod = "none",
    universe = proteins.all,
    minGSSize = 10,
    maxGSSize = 500,
    qvalueCutoff = 1
)

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

# object kk can be used to generate same plots as above!


# GSEA with KEGG pathways.
# for GSEA a per-gene metric is required that is used for ranking.
# assume dataset consists of 17.5k uniprot ids.
# make up fold changes (= ranking metric) ranging from -15 to 15.

protList <- runif(n = 17500, min = -15, max = 15)
names(protList) <- base::sample( proteins.all, 17500)
# for GSEA, sort on FC
protList <- sort(protList, decreasing = TRUE)

# KEGG GSEA
kk2 <- gseKEGG(geneList = protList,
    organism     = "sly",
    keyType      = "uniprot",
    eps          = 0.0,
    minGSSize    = 10,
    maxGSSize    = 500,
    pAdjustMethod = "none",
    pvalueCutoff = 0.9,
    verbose      = FALSE)

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

# again, object kk2 can be used to generate same plots as above!

# Lastly, compare multiple clusters using KEGG.
# First generate list with multiple clusters
clusterData <- list(
    "Cluster 1" = base::sample(proteins.all, 19500),
    "Cluster 2" = base::sample(proteins.all, 18900),
    "Cluster 3" = base::sample(proteins.all, 22000),
    "Cluster 4" = base::sample(proteins.all, 20500) )

# run compareCluster
xx <- compareCluster(clusterData, fun="enrichKEGG",
                     organism="sly",  keyType = "uniprot",
                     pAdjustMethod = "none", pvalueCutoff=1.0)

# plot
emapplot( pairwise_termsim(xx) )
ADD COMMENT
0
Entering edit mode

@guido I am having trouble downloading the goa annotation file, I have tried wget ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/4144447 but it not downloading. I just get the html file. Could you help with this

ADD REPLY
0
Entering edit mode

Could you also help me with SARS-COV-2 goa ftp link? I tried to get it from https://www.ebi.ac.uk/GOA/index but still I can locate the exact file

ADD REPLY
0
Entering edit mode

That should be: http://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/4378222.W_seafood_market_pneumonia_virus_1.goa.

Found through the GOA proteome website (https://www.ebi.ac.uk/GOA/proteomes [note: takes few second to load!]), and then searched for the taxonomy id of SARS-CoV-2, which is 2697049.

ADD REPLY
0
Entering edit mode

I managed to download the goa file now. I have the following question, I have differentially expressed proteins following DE analysis. I cannot see where to pass this information because ID.annotations information is from the goa database.

ADD REPLY
0
Entering edit mode

Have a look at the code after the text ## Functional analysis with clusterProfiler. You basically will need to define two inputs; which I called proteins.all and proteins.random.

proteins.all are all proteins that have a GO annotation (=background), and proteins.random are the proteins you selected/find interesting (e.g. after performing a statistical test) (=foreground). FYI: I selected these 'foreground' proteins randomly, hence its naming.

ADD REPLY

Login before adding your answer.

Traffic: 376 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6