Get All The Genes In A Pathway (E.G. Kegg)

Question

[Stephen]

I'm trying out a few methods for enrichment analysis for GWAS data. For this kind of analysis I need a list of pathways (or GO categories, gene families, etc), where I have a gene grouping followed by all the genes in that grouping, for example:

path1   gene1 gene2
path2   gene3 gene4 gene5 gene6
path3   gene7 gene8 gene9

I downloaded the entire hsa-pathway folder from the KEGG FTP site before it closed down, and I thought I had all this information in a file called hsa_gene_map.tab that looks like this:

2793    04062 04724 04725
2796    04912
2797    04912
2798    04080 04912
2799    00531 01100 04142

But it turns out this is the opposite of what I need - the first number is the Entrez gene ID, and the following fields on each line are KEGG pathway IDs, i.e.

gene1   path1 path2
gene2   path3 path4 path5 path6
gene3   path7 path8 path9

For KEGG, how can I get a list of pathways and all the genes in each of those pathways? I don't have a subscription to KEGG here so I can't use the FTP site any more. A list of files in a simple format as in the first example will be the most useful, since I have to read each of these lines into a list in R to do the analysis.

[Steve Moss]

Does it have to be KEGG? I think you should be able to get a list of genes by GO term from EnsEMBL! That is actually something I need to do in future analyses :-)

[Stephen]

I would also like to do this for other gene groupings - e.g. gene ontology, reactome pathways, etc. The procedures described in the first link above use GO by default, but I'd like to start with KEGG first if possible.

Answer 1

[Duff]

In R you can accomplish this as below:

library(org.Hs.eg.db)
kegg <- org.Hs.egPATH2EG
mapped <- mappedkeys(kegg)
kegg2 <- as.list(kegg[mapped])

This will give you a list (in the R sense) where the list names are KEGG Pathway IDs and the list elements are EntrezGene IDs. It does of course depend on the presence of an org.ORG.eg.db package for your species of interest. HTH duff

Answer 2

[Tg]

Since you mention that you need to do enrichment with reactome pathway and several data source, I suggest that you should get data from http://www.broadinstitute.org/gsea/downloads.jsp, it contains pathways-gene data from kegg, reactome, biocarta and some GO as well.

[Stephen]

This is exactly what I wanted. Thanks. Parsed using sed -r 's/\s+/,/' | sed -r 's/\s+/,/' | sed -r 's/\s+/ /g' for input for SNPath http://linchen.fhcrc.org/grass.html

Answer 3

[Will]

Why can't you just figure out the answer from what you have? If you have the list of all genes and which pathways they are in then by definition by have a list of all pathways and which genes are in them.

Simply use a programming language to 'transpose' the data you have. I don't have any R knowledge but this python script can do it very easily:

from collections import defaultdict

with open('path/to/file') as handle:
    gene2pathway = defaultdict(set)
    for line in handle:
        parts = line.strip().split()
        gene2pathway[parts[0]] = set(parts[1:])

pathway2gene = defaultdict(set)
for gene, pathways in gene2pathway.items():
    for path in pathways:
        pathway2gene[path].add(gene)

for gene, pathways in pathway2gene.items():
    print gene, ' '.join(sorted(pathways))

I'm sure someone with more R knowledge then myself could do this just as easily.

Answer 4

[Khader Shameer]

You can use KEGG API to get all genes in a KEGG pathway. You can access KEGG API via Perl, Ruby, Python or Java. See the tutorial section on Pathways here.

You have to use the function get_genes_by_pathway(string:pathway_id)

Details about various functions / methods and examples scripts are available here.

Also check this discussion on gene-pathway association file from KEGG as well.

[Will]

why use the API .. if you already have the hsa_gene_map file then you already have all of the mappings?

[Khader Shameer]

@Will: Simple, if there is an API use it, instead of raw data munging :).

[Khader Shameer]

@Will: Simple, if there is an API, I prefer to use it instead of raw data munging :).

Answer 5

[Mt]

We've just added some features to Biocyc that make it easier to produce these types of lists, for instance here's a group of ecoli pathways and associated genes that I just put together:

This is exportable; it may not be in the format you want but suggestions for extensions are welcome.

Note you can also do enrichment analysis on this group right on the site.

Answer 6

[Guangchuang Yu]

For GO, you can use the un-exported function getGO2ExtID of my package clusterProfiler

For instance:

clusterProfiler:::getGO2ExtID("GO:0044419", organism="human")
$GO:0044419
  [1] "213"    "259"    "274"    "290"    "291"    "292"    "293"    "336" 
  [9] "358"    "373"    "387"    "527"    "664"    "665"    "699"    "708" 
 [17] "735"    "820"    "857"    "904"    "905"    "912"    "919"    "920" 
 [25] "941"    "942"    "948"    "949"    "975"    "995"    "1050"   "1058"
...

For KEGG, you can use bioconductor package KEGG.db

require(KEGG)
Path2ExtID <- mget(PathID, KEGGPATHID2EXTID, ifnotfound=NA)
Path2ExtID <- lapply(Path2ExtID, function(i) unique(i))

[Guangchuang Yu]

out-dated, see A: Collecting Genes with similar GO term

Answer 7

[Chris Evelo]

If you are interested in other resources as well you could use the WikiPathways webservices to get a list of pathways and next a list of genes (names and IDs) per pathway. See: http://www.wikipathways.org/index.php/Help:WikiPathways_Webservice

Answer 8

[Robin Haw]

There is a Reactome gene set available for downloading from our 'Download' page. As mentioned already it is also part of the GSEA gene sets available from the BROAD.

Answer 9

[ff.cc.cc]

I had to do the same job inside a C/C++ program. I downloaded signature files from Broad (as Tg): http://www.broadinstitute.org/gsea/downloads.jsp Then I parsed and analyzed with the following code:

// input:
// geneIds: set of genes to look for
// filename:gsea filename
// cutoff:  min. number of genes to match
// pLimit:  min. desired significance
// output:
// genesetResult
// genesetP

    static inline int overlapGeneSet(const set<string> &geneIds, const string &filename, int cutoff, double pLimit, vector< vector<string> > &genesetResult, vector<double> &genesetP){
        const int BufSize(100000); // oversized input row buffer
        char *buffer = (char *)malloc( BufSize );
        int result(0);
        ifstream gsea(filename.c_str());
        string strVal;    
        char delimiter = '\t';

        int gseaSize;                   // signature size
        string gseaName;                // signature name
        string gseaSource;              // signature desc.
        vector<string> gseaCommonGenes; // number of matching genes

        // foreach row/geneset
        while(!gsea.eof()){
            istringstream strstream;
            gsea.getline(buffer, BufSize);
            strstream.str(buffer);
            gseaCommonGenes.clear();
            gseaSize = -2;
            gseaName = "";
            gseaSource = "";
            // foreach field in geneset
            while(!strstream.eof()){
                strVal = "";
                getline(strstream, strVal, delimiter);
                if(gseaSize==-2)
                    gseaName = strVal;
                if(gseaSize==-1)
                    gseaSource = strVal;
                if(gseaSize >= 0){
                    if(geneIds.find(strVal) != geneIds.end()){
                        gseaCommonGenes.push_back(strVal);
                    }
                }
                gseaSize++;
            }
            if(gseaCommonGenes.size() >= cutoff){
                result++;
                double P = <your enrichment test here>; 
                        // e.g. hypg(NGenes, gseaSize, geneIds.size(), seaCommonGenes.size());
                if(P < pLimit){
                    gseaCommonGenes.push_back(gseaName);
                    genesetResult.push_back(gseaCommonGenes);
                    genesetP.push_back(P);
                }
            }
        }
        free(buffer);
        return result;
    }

You only have to set your preferred enrichment condition at the line where I assign a value to P for each parsed geneset.

Answer 10

[Andrew W]

Perhaps ugly, but here's a Perl one-liner you can run on the command line that will do this type of data transpose:

cat hsa_gene_map.tab | perl -e 'while(<>) { chomp; @d = split(/\s+/, $_); $g = shift @d; map { $map{$_}{$g}++ } @d; } for $p (sort keys %map) { print join("\t", $p, sort keys %{ $map{$p} }), "\n" }' > hsa_gene_map2.tab