Hi Paul,

Thank you very much, your response is helpful and full of insights.
Personally, I think it'd be great if a similar script could be wrapped into
a function and included into a new or existing package, but I'd like to
hear what other people think.

For Reactome, I think a good alternative could be using the ReactomePA
package which does just that, I don't know how I didn't notice it earlier.

Another resource which in my opinion should be considered for future
applications is PathwayCommons [1] which is freely accessible and already
gathers a considerable set of pathways from different sources. It's
basically what Cytoscape uses for pathway retrieval so it should have a
programmer-friendly access interface.

Best,

[1] http://www.pathwaycommons.org



On 11 September 2013 18:00, Martin Morgan <mtmorgan@fhcrc.org> wrote:

> One could also use KEGG pathways from MSigDB (or any other gene set
> collection;
>
>   http://www.broadinstitute.org/**gsea/msigdb/collections.jsp<http://www.broadinstitute.org/gsea/msigdb/collections.jsp>
>
> registration required and provenance unclear [to me, on quick look]) with
>
> > library(GSEABase)
> > sets <- getGmt("~/Downloads/c2.cp.**kegg.v4.0.entrez.gmt")
> > genes.by.pathway <- lapply(sets, geneIds)
> > names(genes.by.pathway) <- names(sets)
>
> and then as below.
>
> Martin
>
>
> On 09/11/2013 09:05 AM, Paul Shannon wrote:
>
>> Hi Enrico,
>>
>> You ask good questions.
>>
>> In answer to your immediate need, pasted in below are two versions of a
>> short KEGG pathway enrichment script.
>> The first emphasizes simplicity and clarity, and uses only the stock R
>> "dhyper" function along with KEGGREST.
>> The second version (courtesy of Martin, who is also responsible for some
>> of the approaches I use in my simpler version)
>> produces a richer data.frame result.
>>
>> Inasmuch as these examples use KEGGREST, you benefit from the current
>> curated data at KEGG.   As you point out, KEGG.db will only become
>> increasingly out of date.  KEGGREST gets gene/pathway memberships
>> dynamically, avoiding the KEGG prohibition on full downloads.
>>
>> We offer these two KEGGREST:dhypher/hyperg examples informally, and in
>> hopes of encouraging discussion about what Bioconductor should offer more
>> formally.
>>
>> It would be illuminating to run the same demo using pathways from
>> Reactome.  Perhaps you could contribute such a demo?
>>
>> Please let us know what you think.
>>
>>   - Paul
>>
>> ----- simpler version
>> library(KEGGREST)
>> library(org.Hs.eg.db)
>>
>>       # created named list, eg:  path:map00010: "Glycolysis /
>> Gluconeogenesis"
>>
>> pathways.list <- head(keggList("pathway", "hsa"), n=10)
>>
>>       # make them into KEGG-style human pathway identifiers
>> pathway.codes <- sub("path:", "", names(pathways.list))
>>
>>     # for demonstration, just use the first ten pathways
>>     # not all pathways exist for human, so TODO: tryCatch the
>>     # keggGet to be robust against those failures
>>
>>     # subsetting by c(TRUE, FALSE) -- which repeats
>>     # as many times as needed, sorts through some
>>     # unexpected packaging of geneIDs in the GENE element
>>     # of each pw[[n]]
>>
>> genes.by.pathway <- sapply(pathway.codes,
>>                             function(pwid){
>>                                 pw <- keggGet(pwid)
>>                                 pw[[1]]$GENE[c(TRUE, FALSE)]
>>                                 })
>>
>> all.geneIDs <- keys(org.Hs.eg.db)
>>
>>     # chose one of these for demonstration.  the first (a whole genome
>> random
>>     # set of 100 genes)  has very little enrichment, the second, a random
>> set
>>     # from the pathways themsevles,  has very good enrichment
>>
>> genes.of.interest <- c("23118", "23119", "23304", "25998", "26001",
>> "51043",
>>                         "55632", "55643", "55743", "55870", "7314",
>>  "56254",
>>                         "7316",  "144193","784",   "8837",  "1111",
>>  "84706",
>>                         "200931","169522","5707",  "5091",  "5901",
>>  "55532",
>>                         "9777")
>>
>>     # the hypergeometric distribution is traditionally explained in terms
>> of
>>     # drawing a sample of balls from an urn containing black and white
>> balls.
>>     # to keep the arguments straight (in my mind at least), I use these
>> terms
>>     # here also
>>
>> pVals.by.pathway <- sapply(names(genes.by.pathway)**,
>>                          function(pathway) {
>>                             pathway.genes <- genes.by.pathway[[pathway]]
>>                             white.balls.drawn <-
>> length(intersect(genes.of.**interest, pathway.genes))
>>                             white.balls.in.urn <- length(pathway.genes)
>>                             total.balls.in.urn <- length(all.geneIDs)
>>                             black.balls.in.urn <- total.balls.in.urn -
>> white.balls.in.urn
>>                             total.balls.drawn.from.urn <-
>> length(genes.of.interest)
>>                             dhyper(white.balls.drawn,
>>                                    white.balls.in.urn,
>>                                    black.balls.in.urn,
>>                                    total.balls.drawn.from.urn)
>>                             })
>>
>> print(pVals.by.pathway)
>>
>> ---- richer version
>>
>> library(KEGGREST)
>> library(org.Hs.eg.db)
>>
>>       # created named list, length 449, eg:
>>       # path:hsa00010: "Glycolysis / Gluconeogenesis"
>>
>> pathways <- keggList("pathway", "hsa")
>>
>>       # make them into KEGG-style human pathway identifiers
>> human.pathways <- sub("path:", "", names(pathways))
>>
>>     # for demonstration, just use the first ten pathways
>>
>> demo.pathway.ids <- head(human.pathways, 10)
>> demo.pathways <- setNames(keggGet(demo.pathway.**ids), demo.pathway.ids)
>>
>> genes.by.pathway <- lapply(demo.pathways, function(demo.pathway) {
>>       demo.pathway$GENE[c(TRUE, FALSE)]
>>        })
>>
>> all.geneIDs <- keys(org.Hs.eg.db)
>>
>>     # chose one of these for demonstration.  the first (a whole genome
>> random
>>     # set of 100 genes)  has very little enrichment, the second, a random
>> set
>>     # from the pathways themselves,  has very good enrichment in some
>> pathways
>>
>> set.seed(123)
>> significant.genes <- sample(all.geneIDs, size=100)
>> #significant.genes <- sample(unique(unlist(genes.by.**pathway)), size=10)
>>
>>     # the hypergeometric distribution is traditionally explained in terms
>> of
>>     # drawing a sample of balls from an urn containing black and white
>> balls.
>>     # to keep the arguments straight (in my mind at least), I use these
>> terms
>>     # here also
>>
>> hyperg <- Category:::.doHyperGInternal
>> hyperg.test <-
>>      function(pathway.genes, significant.genes, all.genes, over=TRUE)
>> {
>>      white.balls.drawn <- length(intersect(significant.**genes,
>> pathway.genes))
>>      white.balls.in.urn <- length(pathway.genes)
>>      total.balls.in.urn <- length(all.genes)
>>      black.balls.in.urn <- total.balls.in.urn - white.balls.in.urn
>>      balls.pulled.from.urn <- length(significant.genes)
>>      hyperg(white.balls.in.urn, black.balls.in.urn,
>>             balls.pulled.from.urn, white.balls.drawn, over)
>> }
>>
>> pVals.by.pathway <-
>>      t(sapply(genes.by.pathway, hyperg.test, significant.genes,
>> all.geneIDs))
>>
>> print(pVals.by.pathway)
>>
>>
>> On Sep 10, 2013, at 9:23 AM, Enrico Ferrero wrote:
>>
>>  Hi Paul,
>>>
>>> Thanks for your suggestions.
>>> I already performed a GO enrichment analysis and I was specifically
>>> looking at the enrichment of pathways. You also suggest GSEA, but,
>>> correct me if I'm wrong, I think I would need some additional value to
>>> rank my list of gene IDs, as the GSEA algorithm requires a ranked list
>>> to start with.
>>> I guess at the moment my best option is to use Category and the
>>> outdated KEGG.db to have an idea of the enriched pathways and then use
>>> some non-Bioconductor tool (any suggestion?).
>>>
>>> I'm not sure this is the right place, but, on a more general note, may
>>> I ask if there are any plans to provide better and more integrated
>>> support for pathway enrichment analysis in Bioconductor?
>>> For example, would it be possible to build something similar to an up
>>> to date KEGG.db using KEGG's REST API? After all, packages such as
>>> gage and ROntoTools get their information from KEGG in that way.
>>> Alternatively, while not as comprehensive as KEGG, Reactome is
>>> completely open source and open access and I would be surprised if
>>> there wasn't a more fruitful way to integrate it with the core
>>> Bioconductor tools an packages.
>>>
>>> I'm probably being naive here, but I'm afraid I don't fully understand
>>> what is the general consensus of the Bioconductor community on current
>>> status and future directions of pathway analysis tools.
>>>
>>> Thank you.
>>> Best,
>>>
>>> On 9 September 2013 16:25, Paul Shannon <paul.thurmond.shannon@gmail.**
>>> com <paul.thurmond.shannon@gmail.com>> wrote:
>>>
>>>> Hi Enrico,
>>>>
>>>> reactome.db is best described, I believe, as simply a bioc rendering of
>>>> Reactome's sql database -- Marc, please correct me if I am wrong.
>>>>
>>>> There is thus a data representation obstacle when using reactome.db:
>>>> molecular relations are described, with pathway/gene mappings not so easy
>>>> to get at.
>>>> In addition, and despite Reactome's many strengths, its coverage is
>>>> incomplete. The canonical wnt pathway, for instance, is (at my last check)
>>>> not included.
>>>>
>>>> If you have a list of geneIDs, exploratory analysis can usefully start
>>>> out with both GO enrichment, KEGG enrichment, and GSEA. Though the
>>>> information in KEGG.db has not been updated in a couple of years, the
>>>> information there is still very useful for exploratory data analysis. Any
>>>> enrichments you discover using these assorted gene/ctaegory associations
>>>> may lead you to a close study of particular functions or pathways, and it
>>>> is this point that you may wish to get the latest and most specific
>>>> information via KEGGREST and Reactome (and, with our next release) the new
>>>> PSICQUIC package (see http://code.google.com/p/**psicquic/<http://code.google.com/p/psicquic/>
>>>> ).
>>>>
>>>> I hope this helps. Let us know if it falls short, or if new questions
>>>> arise.
>>>>
>>>> - Paul
>>>>
>>>>
>>>> On Sep 9, 2013, at 7:53 AM, Enrico Ferrero wrote:
>>>>
>>>>  Dear list,
>>>>>
>>>>> Can anybody suggest how to perform a simple pathway enrichment
>>>>> analysis starting from a list of gene IDs?
>>>>>
>>>>> I know about the gage and ROntoTools packages that use KEGGREST to
>>>>> retrieve an up to date version of the KEGG database, but, as far as I
>>>>> understand, they require a microarray experiment as input (or at least
>>>>> fold changes and pvalues).
>>>>>
>>>>> Since this time around I'm not starting from a microarray experiment
>>>>> but I just have a gene list, I'm looking for a way to perform pathway
>>>>> enrichment analysis using a simple numerical method such as Fisher's /
>>>>> hypergeometric test.
>>>>>
>>>>> I know the Category package still provides a KEGGHyperG class (which
>>>>> would be perfect!), but the results are based on the outdated version
>>>>> of KEGG (via KEGG.db, I guess).
>>>>>
>>>>> Are there any good alternatives available out there? Would it be
>>>>> possible to use reactome.db in conjunction with the Category/GOstats
>>>>> functions for example?
>>>>>
>>>>> Thank you!
>>>>> Best,
>>>>>
>>>>> --
>>>>> Enrico Ferrero
>>>>> Department of Genetics
>>>>> Cambridge Systems Biology Centre
>>>>> University of Cambridge
>>>>>
>>>>> ______________________________**_________________
>>>>> Bioconductor mailing list
>>>>> Bioconductor@r-project.org
>>>>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https://stat.ethz.ch/mailman/listinfo/bioconductor>
>>>>> Search the archives: http://news.gmane.org/gmane.**
>>>>> science.biology.informatics.**conductor<http://news.gmane.org/gmane.science.biology.informatics.conductor>
>>>>>
>>>>
>>>>
>>> --
>>> Enrico Ferrero
>>> Department of Genetics
>>> Cambridge Systems Biology Centre
>>> University of Cambridge
>>>
>>
>> ______________________________**_________________
>> Bioconductor mailing list
>> Bioconductor@r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https://stat.ethz.ch/mailman/listinfo/bioconductor>
>> Search the archives: http://news.gmane.org/gmane.**
>> science.biology.informatics.**conductor<http://news.gmane.org/gmane.science.biology.informatics.conductor>
>>
>>
>
> --
> Computational Biology / Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N.
> PO Box 19024 Seattle, WA 98109
>
> Location: Arnold Building M1 B861
> Phone: (206) 667-2793
>



-- 
Enrico Ferrero
Department of Genetics
Cambridge Systems Biology Centre
University of Cambridge

	[[alternative HTML version deleted]]

