diff --git a/src/sspa/__init__.py b/src/sspa/__init__.py index edbb0ad..8cc2b7a 100644 --- a/src/sspa/__init__.py +++ b/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 \ No newline at end of file +from .identifier_conversion import identifier_conversion, map_identifiers diff --git a/src/sspa/download_pathways.py b/src/sspa/download_pathways.py index 0aab603..b7faa17 100644 --- a/src/sspa/download_pathways.py +++ b/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,17 +57,17 @@ 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': @@ -51,20 +75,16 @@ def download_KEGG(organism, filepath=None, omics_type='metabolomics'): 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 - \ No newline at end of file + diff --git a/src/sspa/process_pathways.py b/src/sspa/process_pathways.py index 24e0858..eb6ad76 100644 --- a/src/sspa/process_pathways.py +++ b/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,27 +27,31 @@ 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: @@ -55,21 +59,24 @@ def process_kegg(organism, infile=None, download_latest=False, filepath=None, om 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') @@ -77,11 +84,12 @@ def process_kegg(organism, infile=None, download_latest=False, filepath=None, om # 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: @@ -89,6 +97,7 @@ def process_pathbank(organism, infile=None, download_latest=False, filepath=None 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: diff --git a/src/sspa/sspa_gsea.py b/src/sspa/sspa_gsea.py index 7c67f0b..53a83d1 100644 --- a/src/sspa/sspa_gsea.py +++ b/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, diff --git a/src/sspa/utils.py b/src/sspa/utils.py index f85ae6e..4347616 100644 --- a/src/sspa/utils.py +++ b/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: