update pkg_resources

src/sspa/__init__.py

@@ -1,5 +1,6 @@
-from pkg_resources import get_distribution
-__version__ = get_distribution('sspa').version
+from importlib.metadata import version
+
+__version__ = version("sspa")
 
 from .process_pathways import process_reactome, process_kegg, process_gmt, process_pathbank
 from .sspa_cluster import sspa_ssClustPA
@@ -11,4 +12,4 @@
 from .sspa_gsea import sspa_gsea
 from .sspa_ssGSEA import sspa_ssGSEA
 from .download_pathways import download_KEGG, download_reactome
-from .identifier_conversion import identifier_conversion, map_identifiers
+from .identifier_conversion import identifier_conversion, map_identifiers

src/sspa/download_pathways.py

@@ -10,6 +10,30 @@
 import requests
 import io
 
+def extract_kegg_field(lines, field):
+    import re
+    """
+    Extract one KEGG flat-file field, including continuation lines.
+    Example field: 'MODULE', 'GENE', 'COMPOUND', 'KO_PATHWAY'
+    """
+    values = []
+    in_field = False
+    # KEGG field names usually occupy the first 12 characters
+    field_pattern = re.compile(r"^[A-Z_]+")
+    for line in lines:
+        key = line[:12].strip()
+        value = line[12:].strip()
+        if key == field:
+            in_field = True
+            values.append(value)
+        elif in_field and key == "":
+            # continuation line
+            values.append(value)
+        elif in_field and key != "":
+            # next field starts
+            break
+    return values
+
 def download_KEGG(organism, filepath=None, omics_type='metabolomics'):
     '''
     Function for KEGG pathway download
@@ -33,38 +57,34 @@ def download_KEGG(organism, filepath=None, omics_type='metabolomics'):
     for path in pathways:
         path = path.split("\t")
         name = path[1]
-        pathid = re.search(r"(.*)", path[0]).group(1)
+        pathid = re.search(r"(.*)", path[0]).group(1).strip(f"{organism}:")
         pathway_dict[pathid] = name
 
     # get compounds for each pathway
-    base_url = 'http://rest.kegg.jp/get/'
+    base_url = 'https://rest.kegg.jp/get/'
 
     pathway_ids = [*pathway_dict]
     pathway_names = list(pathway_dict.values())
 
     # get release details
-    release_data = requests.get('http://rest.kegg.jp/info/kegg')
+    release_data = requests.get('https://rest.kegg.jp/info/kegg')
     version_no = release_data.text.split()[9][0:3]
 
     if omics_type == 'metabolomics':
         pathway_compound_mapping = dict()
 
         for index,i in enumerate(tqdm(pathway_ids)):
             complist = []
-            current_url = base_url + "pathway:" +i
+            current_url = base_url + organism + i
             # parse the pathway description page
             page = requests.get(current_url)
             lines = page.text.split("\n")
 
-            try:
-                cpds_start = [lines.index(i) for i in lines if i.startswith("COMPOUND")][0]
-                reference_start = [lines.index(i) for i in lines if i.startswith("REFERENCE") or i.startswith("REL_PATHWAY")][0]
-                cpds_lines = lines[cpds_start:reference_start]
-                first_cpd = cpds_lines.pop(0).split()[1]
-                complist.append(first_cpd)
-                complist = complist + [i.split()[0] for i in cpds_lines]
-                pathway_compound_mapping[i] = list(set(complist))
-            except IndexError:
+            cpds_lines = extract_kegg_field(lines, "COMPOUND")
+            complist = [i.split()[0] for i in cpds_lines]
+            pathway_compound_mapping[i] = list(set(complist))
+            if complist==[]:
+                print(current_url + " has no compounds, skipping...")
                 pathway_compound_mapping[i] = []
 
         # remove empty pathway entries
@@ -90,26 +110,17 @@ def download_KEGG(organism, filepath=None, omics_type='metabolomics'):
         for index,i in enumerate(tqdm(pathway_ids)):
             complist = []
             genelist = []
-            current_url = base_url + "pathway:" +i
+            current_url = base_url + organism + i
             # parse the pathway description page
             page = requests.get(current_url)
             lines = page.text.split("\n")
-
-            try:
-                genes_start = [lines.index(i) for i in lines if i.startswith("GENE")][0]
-                cpds_start = [lines.index(i) for i in lines if i.startswith("COMPOUND")][0]
-                reference_start = [lines.index(i) for i in lines if i.startswith("REFERENCE") or i.startswith("REL_PATHWAY")][0]
-                genes_lines = lines[genes_start:cpds_start]
-                cpds_lines = lines[cpds_start:reference_start]
-
-                first_cpd = cpds_lines.pop(0).split()[1]
-                complist.append(first_cpd)
-                complist = complist + [i.split()[0] for i in cpds_lines]
-                first_gene = genes_lines.pop(0).split()[1]
-                genelist.append(first_gene)
-                genelist = genelist + [i.split()[0] for i in genes_lines]
-                pathway_mapping[i] = list(set(complist)) + list(set(genelist))
-            except IndexError:
+            gene_lines = extract_kegg_field(lines, "GENE")
+            cpds_lines = extract_kegg_field(lines, "COMPOUND")
+            gene_list = [i.split()[0] for i in gene_lines]
+            compound_list = [i.split()[0] for i in cpds_lines]
+            pathway_mapping[i] = list(set(compound_list)) + list(set(gene_list))
+            if compound_list==[] and gene_list==[]:
+                print(current_url + " has no compounds or genes, skipping...")
                 pathway_mapping[i] = []
 
         # remove empty pathway entries
@@ -441,4 +452,4 @@ def download_pathbank(organism, filepath=None, omicstype='metabolomics'):
         print("Complete!")
 
         return multiomics_pathways_gmt
-
+

src/sspa/process_pathways.py

@@ -1,8 +1,8 @@
 import pandas as pd
-import pkg_resources
+from importlib import resources as resources
 import sspa.download_pathways 
 
-def process_reactome(organism, infile=None, download_latest=False, filepath=None, omics_type='metabolomics', identifiers=None):
+def process_reactome(organism, infile=None, download_latest=False, filepath=None, omics_type='metabolomics', identifiers=None, sep="\t"):
     '''
     Function to load Reactome pathways 
     Args:
@@ -12,7 +12,7 @@ def process_reactome(organism, infile=None, download_latest=False, filepath=None
         filepath (str): filepath to save pathway file to, default is None - save to variable
         omics_type(str): If using download_latest, specify type of omics pathways to download. Options are 'metabolomics', 'proteomics', 'transcriptomics', or 'multiomics'
         identifiers (list): list of identifiers to download for multi-omics pathways, default is None (download all). Options are 'chebi', 'uniprot', 'gene_symbol'
-
+        sep (str): separator for the input file, default is '\t'
     Returns: 
         GMT-like pd.DataFrame containing Reactome pathways
     '''
@@ -27,68 +27,77 @@ def process_reactome(organism, infile=None, download_latest=False, filepath=None
         if omics_type != 'metabolomics':
             raise ValueError('Proteomics/multi-omics pathways only accessible when download_latest=True')
         if infile == None or infile == "R78":
-            stream = pkg_resources.resource_stream(__name__, 'pathway_databases/ChEBI2Reactome_All_Levels_R78.txt')
-            f = pd.read_csv(stream, sep="\t", header=None, encoding='latin-1')
+            stream = resources.files(__name__).joinpath('pathway_databases/ChEBI2Reactome_All_Levels_R78.txt').open('rb')
+            f = pd.read_csv(stream, sep=sep, header=None, encoding='latin-1')
+            f.columns = ['CHEBI', 'pathway_ID', 'link', 'pathway_name', 'evidence_code', 'species']
+            f_filt = f[f.species == organism]
+            name_dict = dict(zip(f_filt['pathway_ID'], f_filt['pathway_name']))
+
+            groups = f_filt.groupby(['pathway_ID'])['CHEBI'].apply(list).to_dict()
+            groups = {k: list(set(v)) for k, v in groups.items()}
+            df = pd.DataFrame.from_dict(groups, orient='index', dtype="object")
         else:
-            f = pd.read_csv(infile, sep="\t", header=None)
-
-        f.columns = ['CHEBI', 'pathway_ID', 'link', 'pathway_name', 'evidence_code', 'species']
-        f_filt = f[f.species == organism]
-        name_dict = dict(zip(f_filt['pathway_ID'], f_filt['pathway_name']))
+            df = pd.read_csv(infile, index_col=0,sep=sep, header=None)
 
-        groups = f_filt.groupby(['pathway_ID'])['CHEBI'].apply(list).to_dict()
-        groups = {k: list(set(v)) for k, v in groups.items()}
-        df = pd.DataFrame.from_dict(groups, orient='index', dtype="object")
+
         pathways_df = df.dropna(axis=0, how='all', subset=df.columns.tolist()[1:])
         pathways_df = df.dropna(axis=1, how='all')
 
-        pathways_df["Pathway_name"] = pathways_df.index.map(name_dict)
-        pathways_df.insert(0, 'Pathway_name', pathways_df.pop('Pathway_name'))
+        # Remove duplicated compounds
+        mask = pathways_df.apply(pd.Series.duplicated, 1) & pathways_df.astype(bool)
+        pathways_df = pathways_df.mask(mask, None)
+
+        pathways_df = pathways_df.rename(columns={1: "Pathway_name"})
 
         return pathways_df
 
-def process_kegg(organism, infile=None, download_latest=False, filepath=None, omics_type='metabolomics'):
+def process_kegg(organism, infile=None, download_latest=False, filepath=None, omics_type='metabolomics',sep="\t"):
     '''
     Function to load KEGG pathways 
     Args:
         organism (str): KEGG organism code
         infile (str): default None, provide a KEGG pathway file to process into the GMT-style dataframe 
         download_latest (Bool): Downloads the latest version of KEGG metabolic pathways
         filepath (str): filepath to save pathway file to, default is None - save to variable
+        sep (str): separator for the input file, default is '\t'
     Returns: 
         GMT-like pd.DataFrame containing KEGG pathways
     '''
     if download_latest:
         pathways_df = sspa.download_pathways.download_KEGG(organism, filepath, omics_type)
+
+        pathways_df = pathways_df.rename(columns={1: "Pathway_name"})
         return pathways_df
 
     else:
         if omics_type != 'metabolomics':
             raise ValueError('Proteomics/multi-omics pathways only accessible when download_latest=True')
         if infile == None or infile == "R98":
-            stream = pkg_resources.resource_stream(__name__, 'pathway_databases/KEGG_human_pathways_compounds_R98.csv')
+            stream = resources.files(__name__).joinpath('pathway_databases/KEGG_human_pathways_compounds_R98.csv').open('rb')
             pathways_df = pd.read_csv(stream, index_col=0, encoding='latin-1')
         else:
-            pathways_df = pd.read_csv(infile, index_col=0)
+            pathways_df = pd.read_csv(infile, index_col=0, sep=sep, header=None)
 
         pathways_df = pathways_df.dropna(axis=0, how='all', subset=pathways_df.columns.tolist()[1:])
         pathways_df = pathways_df.dropna(axis=1, how='all')
 
         # Remove duplicated compounds
         mask = pathways_df.apply(pd.Series.duplicated, 1) & pathways_df.astype(bool)
         pathways_df = pathways_df.mask(mask, None)
+        pathways_df = pathways_df.rename(columns={1: "Pathway_name"})
 
         return pathways_df
 
 
-def process_pathbank(organism, infile=None, download_latest=False, filepath=None, omics_type='metabolomics'):
+def process_pathbank(organism, infile=None, download_latest=False, filepath=None, omics_type='metabolomics', sep="\t"):
     '''
     Function to load PathBank pathways 
     Args:
         infile (str): default None, provide a PathBank pathway file to process into the GMT-style dataframe 
         download_latest (Bool): Downloads the latest version of PathBank metabolic pathways
         filepath (str): filepath to save pathway file to, default is None - save to variable
         omics_type(str): If using download_latest, specify type of omics pathways to download. Options are 'metabolomics', 'proteomics', or 'multiomics'
+        sep (str): separator for the input file, default is '\t'
     Returns: 
         GMT-like pd.DataFrame containing PathBank pathways
     '''
@@ -98,21 +107,22 @@ def process_pathbank(organism, infile=None, download_latest=False, filepath=None
 
     else:
         if infile:
-            pathways_df = pd.read_csv(infile, index_col=0)
+            pathways_df = pd.read_csv(infile, sep=sep,index_col=0)
         else: 
             print('Set download_latest=True to download latest version of PathBank pathways or provide a saved PathBank pathway .csv/.gmt file to load')
 
 
-def process_gmt(infile):
+def process_gmt(infile, sep=","):
     '''
     Function to load pathways from a custom GMT-like file
     Args:
         infile (str): default None, provide a GMT pathway file to process into the GMT-style dataframe, file ending can be .csv or .gmt
+        sep (str): separator for the input file, default is ','
     Returns: 
         GMT-like pd.DataFrame containing pathways
     '''
     if infile[-4:] == ".csv":
-        pathways_df = pd.read_csv(infile, index_col=0, dtype='object')
+        pathways_df = pd.read_csv(infile, index_col=0, dtype='object', sep=sep)
     elif infile[-4:] == ".gmt":
         input_gmt = []
         with open(infile, "r") as f:

src/sspa/sspa_gsea.py

@@ -14,13 +14,19 @@ def sspa_gsea(mat, metadata, pathway_df, ranking_metric='signal_to_noise', min_e
             Other options are 't_test' and see GSEApy package https://github.com/zqfang/GSEApy/blob/2b5419e14615b6fd19a575ff065256dc7099bbec/gseapy/gsea.py#L135 for more options. 
         min_entity (int, optional): minimum number of molecules mapping to pathways for GSEA to be performed. Defaults to 2.
     """
-
+    if "Gene" not in mat.columns:
+        #check if  gene names is set as a column, if not, transpose and reset index
+        mat = mat.T.reset_index()
+        mat = mat.rename(columns={"index": "Gene"})
+    else:
+        mat = mat.copy()
+
     pathway_names = pathway_df["Pathway_name"].to_dict()
     pathways = utils.pathwaydf_to_dict(pathway_df)
-    compounds_present = mat.columns.tolist()
+    compounds_present = mat["Gene"].tolist()
     pathways = {k: v for k, v in pathways.items() if len([i for i in compounds_present if i in v]) >= min_entity}
 
-    gsea_res = gseapy.gsea(data=mat.T, 
+    gsea_res = gseapy.gsea(data=mat, 
                  gene_sets=pathways, 
                  cls=metadata,
                  min_size=min_entity,

src/sspa/utils.py

@@ -1,6 +1,6 @@
 import pandas as pd
 import numpy as np
-import pkg_resources
+from importlib import resources as resources
 import scipy.stats as stats
 import statsmodels.api as sm
 
@@ -22,11 +22,13 @@ def load_example_data(omicstype="metabolomics", processed=True):
 
     if omicstype == "metabolomics":
         if processed:
-            stream = pkg_resources.resource_stream(__name__, 'example_data/Su_covid_metabolomics_processed.csv')
+            stream = resources.files(__name__).joinpath('example_data/Su_covid_metabolomics_processed.csv').open('rb')
+            #stream = pkg_resources.resource_stream(__name__, 'example_data/Su_covid_metabolomics_processed.csv')
             f = pd.read_csv(stream, index_col=0, encoding='latin-1')
             return f
         else:
-            stream = pkg_resources.resource_stream(__name__, 'example_data/Su_metab_data_raw.csv')
+            stream = resources.files(__name__).joinpath('example_data/Su_covid_metabolomics_processed.csv').open('rb')
+            #stream = pkg_resources.resource_stream(__name__, 'example_data/Su_metab_data_raw.csv')
             f = pd.read_csv(stream, index_col=0, encoding='latin-1')
             return f
 
@@ -81,7 +83,7 @@ def pathwaydf_to_dict(df):
     pathway_dict = {}
 
     for pathway in pathways_df.index:
-        pathway_compounds = list(set(pathways_df.loc[pathway, :].tolist()))
+        pathway_compounds = list(set(pathways_df.loc[pathway].dropna().tolist()))
         pathway_compounds = [str(i) for i in pathway_compounds if str(i) not in ["None", np.nan, 'nan']]
 
         if len(pathway_compounds) > 1: