Parsing databases for creating a comprehensive metabolite-gene interaction network
This tutorial demonstrates how to parse and integrate data from the Human Metabolome Database (HMDB), Reactome, and KEGG to create a comprehensive metabolite-gene interaction network (in AnnData format). The resulting network can be used for various bioinformatics analyses, including enrichment analysis and subsequent graph-based analyses.
Installation
From PyPI
pip install cospa
From Source
pip install git+https://github.com/compbioclub/COSPA.git@main
Importing Required Libraries
import os
import pandas as pd
import numpy as np
import cospa.parsing_database as pdb
Human Metabolome Database (HMDB)
Data Download
HMDB XML files can be downloaded from https://hmdb.ca/downloads
e.g., All Metabolites: https://hmdb.ca/system/downloads/current/hmdb_metabolites.zip
Parsing HMDB XML Files
Once we have downloaded and extracted the HMDB XML files, we can parse them to extract compound mapping information and metabolite-gene interactions by using the pdb.parse_hmdb_xml() function.
# filename = "hmdb_serum_metabolites.xml"
filename = "./hmdb/hmdb_all_metabolites.xml"
df = pdb.parse_hmdb_xml(filename)
The df DataFrame contains information about metabolites, including their IDs (HMDB, KEGG, CHEBI), names, chemical formulas, and associated genes.
df
We can save the parsed data to a CSV file for future reference (CHEBI-KEGG mapping).
df.to_csv("hmdb_all_metabolites.csv", index=False)
Extracting Metabolite-Gene Interactions from HMDB
We can extract metabolite-gene interactions gene_compound_df (pairwise), gene-gene interactions gene_relation_df (pairwise, second-order neighbors connected via metabolites), and compound-compound interactions compound_relation_df (pairwise, second-order neighbors connected via genes) from the parsed HMDB data using the pdb.extract_hmdb_relations() function.
gene_compound_df, gene_relation_df, compound_relation_df = pdb.extract_hmdb_relations(df)
gene_compound_df
gene_relation_df
compound_relation_df
Mapping Uniprot gene id to gene symbol
We can map the Uniprot gene IDs to gene symbols using the pdb.apply_mapping_and_cleanup() function. This function takes a DataFrame, a mapping file path, and other parameters to perform the mapping and cleanup.
# all_gene_compound_df, all_gene_relation_df
gene_compound_df = pdb.apply_mapping_and_cleanup(
df=gene_compound_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="gene",
target_col="symbol",
origin_ref_col="uniprot_ids",
target_ref_col="symbol",
origin_col_dtype=str,
origin_ref_dtype=str,
)
gene_relation_df = pdb.apply_mapping_and_cleanup(
df=gene_relation_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="gene1",
target_col="symbol1",
origin_ref_col="uniprot_ids",
target_ref_col="symbol",
)
gene_relation_df = pdb.apply_mapping_and_cleanup(
df=gene_relation_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="gene2",
target_col="symbol2",
origin_ref_col="uniprot_ids",
target_ref_col="symbol",
)
gene_compound_df
gene_relation_df
gene_relation_df
Creating three corresponding AnnData objects
After obtaining the cleaned DataFrames for gene-compound, gene-gene, and compound-compound interactions, we can convert them into AnnData objects using the pdb.create_adata() function. These AnnData objects can be used for various downstream analyses.
adata_gc, adata_gg, adata_cc = pdb.create_adata(gene_compound_df, gene_relation_df, compound_relation_df)
adata_gc, adata_gg, adata_cc
KEGG pathway database
Data Download
KEGG pathway data can be fetched using the KEGG REST API. We can download pathway XML files for a specific species (e.g., human) using the pdb.download_all_pathway_xml() function.
This function takes the species code (e.g., "hsa" for human) and the directory path where the XML files will be saved. If no species code is provided, it defaults to downloading KEGG Orthology (KO) pathways.
The XML files will be saved in the specified directory ($kgml_path/$species_code).
species_code = "hsa" # Change to desired species code or None for all pathways
kgml_path = "./kegg"
pdb.download_all_kegg_pathway_xml(species_code, kgml_path)
Parsing one KEGG pathway XML file
Once we have downloaded the KGML XML files, we can parse a specific pathway XML file to extract metabolite-gene interactions using the pdb.parse_kegg_relations() function. This function takes the path to a KGML XML file and returns a DataFrame containing gene-compound interactions.
species_code = "hsa" # Change to desired species code or None for all pathways
pathway_id = "hsa01100" # Example pathway ID for Glycolysis / Gluconeogenesis
kgml_path = f"./kegg/{species_code}/{pathway_id}.xml"
gene_compound_df = pdb.parse_kegg_relations(kgml_path)
The resulting DataFrame contains KEGG gene IDs and compound IDs (KEGG CIDs).
gene_compound_df
Parsing all KEGG pathway XML files for a species
We can parse all KGML XML files for a specific species using the pdb.parse_all_kegg_pathways() function.
This function takes the species code and the directory path where the XML files are stored. It returns a DataFrame containing all gene-compound interactions for the specified species.
species_code = "hsa" # Change to desired species code or None for all pathways
all_gene_compound_df = pdb.parse_all_kegg_pathways(species_code, "./kegg")
all_gene_compound_df
Mapping KEGG gene id & compound id to gene symbol & CHEBI id
For the consistency of all data formats, we need to map the KEGG gene IDs to gene symbols and KEGG compound IDs to CHEBI IDs using the pdb.apply_mapping_and_cleanup() function.
all_gene_compound_df = pdb.apply_mapping_and_cleanup(
df=all_gene_compound_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="kegg_gene_id",
target_col="symbol",
origin_ref_col="entrez_id",
target_ref_col="symbol",
cleanup_subset_cols=None,
origin_col_dtype=str,
origin_ref_dtype=str,
)
all_gene_compound_df = pdb.apply_mapping_and_cleanup(
df=all_gene_compound_df,
mapping_file_path="./hmdb/hmdb_all_metabolites.csv",
mapping_sep=",",
origin_col="kegg_compound_id",
target_col="compound",
origin_ref_col="kegg_id",
target_ref_col="chebi_id",
cleanup_subset_cols=None,
origin_col_dtype=str,
origin_ref_dtype=str,
)
all_gene_compound_df
Constructing gene/compound relations
We can construct gene-gene and compound-compound interaction DataFrames from the gene-compound interaction DataFrame using the pdb.construct_kegg_relations() function. This function takes the gene-compound DataFrame and returns two DataFrames: one for gene-gene interactions gene_relation_df (pairwise, second-order neighbors connected via metabolites), and another for compound-compound interactions compound_relation_df (pairwise, second-order neighbors connected via genes)
all_gene_relation_df, all_compound_relation_df = pdb.construct_kegg_relations(all_gene_compound_df, gene_col="symbol", compound_col="compound", output_gene_col="symbol", output_compound_col="compound")
all_gene_relation_df
all_compound_relation_df
Creating three corresponding AnnData objects
Same as HMDB, after obtaining the cleaned DataFrames for gene-compound, gene-gene, and compound-compound interactions, we can convert them into AnnData objects using the pdb.create_adata() function. These AnnData objects can be used for various downstream analyses.
adata_gc, adata_gg, adata_cc = pdb.create_adata(all_gene_compound_df, all_gene_relation_df, all_compound_relation_df)
adata_gc, adata_gg, adata_cc
Reactome pathway database
Data Download
The Reactome SBML files can be fetched from: https://reactome.org/download-data
e.g., for human pathways: https://download.reactome.org/94/homo_sapiens.3.1.sbml.tgz, for all pathways: https://download.reactome.org/94/all_species.3.1.sbml.tgz
Parsing Reactome Pathway Files
Once we have downloaded and extracted the Reactome pathway files, we can parse them to extract pathway information and gene-metabolite interactions by using the pdb.parse_reactome_sbml() function.
df_compartments, df_species, df_reactions = pdb.parse_reactome_sbml('./reactome/reactome_hsa/R-HSA-15869.sbml')
df_compartments
df_species
df_reactions
Extracting Metabolite-Gene Interactions from Reactome
From df_species and df_reactions these two DataFrames, we can extract metabolite-gene interactions gene_compound_df (pairwise), gene-gene interactions gene_relation_df (pairwise, second-order neighbors connected via metabolites), and compound-compound interactions compound_relation_df (pairwise, second-order neighbors connected via genes) from the parsed Reactome data using the pdb.extract_reactome_relations() function.
gene_compound_df, gene_relation_df, compound_relation_df = pdb.extract_reactome_relations(df_species, df_reactions)
Parsing all Reactome Pathways
We can parse all KGML XML files for a specific species using the pdb.parse_all_kegg_pathways() function.
This function takes the species code and the directory path where the XML files are stored. It returns a DataFrame containing all gene-compound interactions for the specified species.
all_gene_compound_df, all_gene_relation_df, all_compound_relation_df = pdb.process_all_reactome_pathways('hsa', pathway_dir="./reactome/reactome_hsa")
Mapping Reactome uniprot gene id to gene symbol
For the consistency of all data formats, we need to map the uniprot gene IDs to gene symbols using the pdb.apply_mapping_and_cleanup() function.
all_gene_compound_df = pdb.apply_mapping_and_cleanup(
df=all_gene_compound_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="gene",
target_col="symbol",
origin_ref_col="uniprot_ids",
target_ref_col="symbol",
origin_col_dtype=str,
origin_ref_dtype=str,
)
all_gene_relation_df = pdb.apply_mapping_and_cleanup(
df=all_gene_relation_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="gene1",
target_col="symbol1",
origin_ref_col="uniprot_ids",
target_ref_col="symbol",
origin_col_dtype=str,
origin_ref_dtype=str,
)
all_gene_relation_df = pdb.apply_mapping_and_cleanup(
df=all_gene_relation_df,
mapping_file_path="./hgnc_complete_set.txt",
mapping_sep="\t",
origin_col="gene2",
target_col="symbol2",
origin_ref_col="uniprot_ids",
target_ref_col="symbol",
origin_col_dtype=str,
origin_ref_dtype=str,
)
Creating three corresponding AnnData objects
Same as HMDB & KEGG, after obtaining the cleaned DataFrames for gene-compound, gene-gene, and compound-compound interactions, we can convert them into AnnData objects using the pdb.create_adata() function. These AnnData objects can be used for various downstream analyses.
adata_gc, adata_gg, adata_cc = pdb.create_adata(all_gene_compound_df, all_gene_relation_df, all_compound_relation_df)
adata_gc, adata_gg, adata_cc