diff --git a/README.md b/README.md index c80926a..5cc5ad0 100644 --- a/README.md +++ b/README.md @@ -29,6 +29,8 @@ This repository provides accessory tools associated with papers and activities o - **ThermoMPNN_trainingdata_overlap**: Contains a python script to find matching proteins in a MAVISp dataset and the training datasets of ThermoMPNN (see paper: https://doi.org/10.1073/pnas.2314853121) with the purpose of defining a benchmarking dataset for the ThermoMPNN tool integrated into MAVISp. The directory also includes scripts for data cleanup steps that do **not** have to be rerun before using the comparison script. - **enzymes_annotation**: This Python script retrieves enzyme annotations starting from a UniProt accession (UniProt AC) by querying the M-CSA (https://www.ebi.ac.uk/thornton-srv/m-csa/) database and, when needed, the UniProt REST API. Its main purpose is to obtain enzyme classification (EC number), catalytic site residues, and information about whether the annotation is manually curated or predicted by homology. Depending on database availability and selected flags, the script can return manually curated catalytic residues from M-CSA, predicted catalytic residues based on homologous alignments, or, if no catalytic site information is available, only the EC classification retrieved from UniProt. The script uses REST API requests to collect data, processes results into structured tables using pandas, and outputs the results as CSV files. - **mavisp_isoform_support**: A collection of tools and utilities that provide the necessary support to run **MAVISp** on **non-canonical isoforms**. The module helps extend standard MAVISp analyses beyond canonical UniProt entries by enabling identifier mapping, isoform-aware annotation, and preparation of the required input data and workflows for isoform-specific variant interpretation. +- **reactome_to_mavisp.py**: This pipeline runs an automated Reactome analysis for one or more UniProt accessions. For each target protein, the workflow retrieves human Reactome pathways, identifies reactions that contain the target, parses BioPAX-level reaction/protein/complex annotations, and writes cleaned CSV outputs in MAVISp supported format for downstream analysis. + ## Citations diff --git a/tools/reactome_to_mavisp/.DS_Store b/tools/reactome_to_mavisp/.DS_Store new file mode 100644 index 0000000..313e3bf Binary files /dev/null and b/tools/reactome_to_mavisp/.DS_Store differ diff --git a/tools/reactome_to_mavisp/README.md b/tools/reactome_to_mavisp/README.md new file mode 100644 index 0000000..44598dc --- /dev/null +++ b/tools/reactome_to_mavisp/README.md @@ -0,0 +1,477 @@ +# Reactome UniProt Reaction Workflow + +This pipeline runs an automated Reactome analysis for one or more UniProt accessions. For each target protein, the workflow retrieves human Reactome pathways, identifies reactions that contain the target, parses BioPAX-level reaction/protein/complex annotations, and writes cleaned CSV outputs in MAVISp supported format for downstream analysis. + +The code is organized as a small Python package. The command-line entry point is `reactome_to_mavisp.py`, while the core logic is split across the `reactome_pipeline/` modules. + +--- + +## Requirements + +### Python + +```text +Python >= 3.8 +``` + +### Python packages + +Required packages: + +```text +reactome2py +pybiopax +pandas +numpy +networkx +requests +``` + +The workflow also uses standard-library modules such as `argparse`, `os`, `shutil`, `re`, `time`, `copy`, `contextlib`, `io`, `collections`, `pathlib`, and `typing`. + +Example installation: + +```bash +pip install reactome2py pybiopax pandas numpy networkx requests +``` + +--- + +## Description + +### Main files + +| File | Role | +|---|---| +| `reactome_to_mavisp.py` | Main command-line entry point. Parses arguments, checks whether a UniProt accession is mapped in Reactome, runs the workflow, and writes `entries_not_in_reactome.csv` when needed. | +| `reactome_pipeline/workflow.py` | Contains the `ReactomeScript` class and orchestrates the full Reactome workflow for one UniProt accession. | +| `reactome_pipeline/reactome_analysis.py` | Contains Reactome/BioPAX helper functions, safe retry logic for Reactome calls, BioPAX parsing, disease-link extraction, and target-reaction checks. | +| `reactome_pipeline/data_processing.py` | Flattens nested reaction/pathway/protein annotations into a cleaned `pandas.DataFrame`, reorders columns, removes duplicates, and propagates complex information. | +| `reactome_pipeline/graph_utils.py` | Builds and processes directed reaction graphs used for optional pathway-ordering analysis. | +| `reactome_pipeline/uniprot_utils.py` | Contains UniProt helper functions for gene-to-accession conversion, accession-to-protein-name retrieval, and parsing protein names. | +| `reactome_post_process.py` | Optional post-processing script that merges individual `result.csv` files into summary CSV tables. | +| `reactome_pipeline/__init__.py` | Marks `reactome_pipeline/` as a Python package. It can remain empty. | + +The scripts are organized as follow: + +```text +project/ +├── reactome_to_mavisp.py +├── reactome_post_process.py +├── reactome_pipeline/ +│ ├── __init__.py +│ ├── workflow.py +│ ├── reactome_analysis.py +│ ├── data_processing.py +│ ├── graph_utils.py +│ └── uniprot_utils.py +└── README.md +``` + +Python cache files such as `__pycache__/` or `*.pyc` are generated automatically and should not be edited or tracked manually. + +For each UniProt accession, the workflow performs the following steps. + +### 1. Initial Reactome mapping check + +Before running the full workflow, `reac_classified.py` first checks whether the input UniProt accession is mapped to at least one human Reactome pathway. + +This check is performed using: + +```python +rc.content.mapping( + id=uniprot_ac, + resource="UniProt", + species="9606", + by="pathways" +) +``` + +This initial step is used only to decide whether the accession should be processed further. + +If Reactome returns no mapped human pathways, the accession is not analysed and is written to: + +```text +entries_not_in_reactome.csv +``` + +with the status: + +```text +not_found_in_reactome +``` + +This prevents the workflow from spending time on accessions for which Reactome has no pathway-level annotation. + +--- + +### 2. Retrieve Reactome pathways + +For accessions that pass the initial mapping check, the workflow retrieves all human Reactome pathways associated with the UniProt accession using: + +```python +rc.content.mapping( + id=uniprot_ac, + resource="UniProt", + species="9606", + by="pathways" +) +``` + +Each retrieved pathway is later expanded into its full Reactome hierarchy using ancestor information. + +This allows the final output to report not only the specific low-level pathway containing a reaction, but also the broader pathway context in which that reaction occurs. + +--- + +### 3. Optional pathway ordering + +Unless `--skip_pathway_order` is used, the workflow attempts to infer the order of reactions within each retrieved pathway. + +For each pathway, the workflow downloads the corresponding BioPAX model and extracts next/previous reaction relationships. These relationships are used to build a directed graph where: + +- nodes represent reactions; +- edges represent next/previous relationships between reactions; +- inferred reaction paths are written to `ordered_paths.csv`. + +The pathway-ordering step is optional because it can be slow for proteins associated with many Reactome pathways. This is due to the need for multiple BioPAX downloads and graph operations. + +When `--skip_pathway_order` is used, the workflow skips this step and still produces the main `result.csv` output. + +--- + +### 4. Resolve target information + +Before filtering reactions, the workflow collects target-level information for the input UniProt accession. This is done in: + +```python +ReactomeScript.resolve_target_information() +``` + +The purpose of this step is to prepare the identifiers that will later be used to decide whether a candidate Reactome reaction actually contains the input protein. + +The workflow collects: + +- the readable name of the target protein, which will be written in `result.csv`; +- optional Reactome stable identifiers for the target protein; +- alternative Reactome forms of the same protein; +- Reactome complexes associated with the UniProt accession; +- the identifiers needed for target-reaction detection. + +The final `target_name` is selected using the following priority: + +1. protein name retrieved from UniProt; +2. Reactome display/name retrieved from `search_fireworks`, if that optional call works; +3. the original UniProt accession as fallback. + +The `search_fireworks` call is used only as an optional extra source of Reactome-specific protein identifiers. These identifiers can help detect reactions containing specific Reactome protein forms of the target. + +However, `search_fireworks` is not required for the workflow to continue. If Reactome returns a server error or no valid entries, the workflow continues using: + +- Reactome complexes associated with the UniProt accession; +- direct UniProt accession matching inside BioPAX reaction models. + +The target information dictionary contains: + +| Field | Meaning | +|---|---| +| `target_name` | Final target name written in `result.csv`. Preferentially retrieved from UniProt, otherwise from Reactome, otherwise set to the input UniProt accession. | +| `target_protein_stId` | Reactome stable IDs corresponding to the target protein, when available from `search_fireworks`. Used as an optional way to recognise reactions containing the target. | +| `target_protein_stId_other` | Alternative Reactome forms of the target protein, retrieved from Reactome when target stable IDs are available. | +| `all_target_protein_stId` | Combined list of direct and alternative Reactome target IDs. Used later during reaction filtering. | +| `target_protein_complexes_names` | Reactome complexes associated with the UniProt accession. Used later to identify reactions where the target appears as part of a complex. | + +--- + +### 5. Collect candidate reactions + +After retrieving the Reactome pathways, the workflow collects the reactions that may potentially involve the input protein. + +For each lowest-level pathway, the workflow retrieves the contained events and keeps only true Reactome reactions. + +Pathway containers are excluded to avoid analysing higher-level pathway objects as if they were individual biochemical reactions. + +At this stage, the reactions are still considered **candidate reactions**. This means that they belong to Reactome pathways associated with the input UniProt accession, but they have not yet been confirmed to directly contain the target protein. + +--- + +### 6. Identify target reactions + +The script takes: + +- the candidate reactions collected from Reactome pathways; +- the target information prepared by `resolve_target_information()`. + +For each candidate reaction, the workflow checks whether the reaction actually contains the input protein. + +A reaction is considered a **target reaction** if at least one of the following checks is true: + +- the reaction contains a Reactome complex associated with the target UniProt accession; +- the reaction contains a Reactome protein stable ID corresponding to the target protein or one of its alternative forms; +- the BioPAX model of the reaction directly contains the input UniProt accession. + +The third check is the most robust one because it searches BioPAX protein and entity-reference cross-references directly for the input UniProt accession. + +Only reactions that pass this filtering step are kept for detailed annotation extraction. + +--- + +### 7. Extract BioPAX annotations + +For each reaction confirmed to contain the target protein, the workflow downloads or reuses the corresponding BioPAX model and extracts detailed reaction-level annotations. + +The extracted information includes: + +- protein display name; +- UniProt accession; +- cellular location; +- sequence intervals; +- sequence sites; +- modification type; +- complex membership; +- stoichiometry; +- parent protein family or physical entity; +- pathway hierarchy; +- reaction name and Reactome stable ID; +- biochemical left/right participants; +- conversion direction; +- regulatory controllers; +- disease links. + +This step produces nested annotation dictionaries describing the target-containing reactions. + +--- + +### 8. Build and write output tables + +The nested annotation dictionaries are flattened into a `pandas.DataFrame`. + +The final table is then: + +- cleaned; +- deduplicated; +- column-ordered; +- filtered to remove protein-family rows; +- optionally reordered using pathway-ordering files; +- written to `result.csv`. + +Protein-family rows are removed so that the final output focuses on individual protein entries rather than broad family-level Reactome entities. + +The final `result.csv` therefore contains only Reactome reactions that passed the target-reaction filtering step and for which BioPAX-level annotations could be extracted. + + +### Optional post-processing + +After running the main workflow for multiple UniProt accessions, the optional script `reactome_post_process.py` can merge individual `result.csv` files. It concatenates the available result tables, harmonizes the target columns when needed, removes duplicate rows, and writes three post-processed output files: + + +**merged_reaction.csv** contains a compact reaction-level summary across all analysed UniProt accessions. It keeps the target accession, target name, pathway hierarchy, disease annotation, lowest-level pathway, reaction name, reaction ID, reaction participants, and reaction direction. This file is useful for quickly comparing which reactions are associated with each protein across the full Reactome output. + +**merged_highest_pathways.csv** contains a simplified pathway-level summary. It reports the highest-level Reactome pathways associated with each target protein, together with the corresponding Reactome pathway ID, UniProt accession, and target name. This file is useful for obtaining a non-redundant overview of the major biological areas covered by the analysed proteins. + +**disease_single_sequence_site.csv** contains a filtered subset of the concatenated results. It keeps only rows belonging to the highest-level Reactome pathway Disease and where the sequence-site annotation corresponds to a single numeric residue position. This output is useful for downstream inspection of disease-associated reactions involving specific modified or annotated residue sites + +--- + +## Input + +Input for reactome_to_mavisp.py script: + +| Argument | Description | +|---|---| +| `-u`, `--uniprot_ac` | UniProt accession to analyze. Default: `Q8N726`. Ignored when `--uniprot_file` is supplied. | +| `-uf`, `--uniprot_file` | Text file containing one UniProt accession per line. Blank lines and lines starting with `#` are ignored. | +| `-o`, `--output_dir` | Main output directory. Default: `reactome_outputs`. A subfolder is created for each accession. | +| `-s`, `--skip_pathway_order` | Skip pathway-order inference. This speeds up the analysis and avoids writing `pathways_order/` files. | + +Here an example of file with a list of uniprot ac + + +```text +P04637 +Q8N726 +Q9Y2X3 +``` +Input for reactome_post_process.py script: + +Arguments: + +| Argument | Description | +|---|---| +| `-i`, `--input_dir` | Main Reactome output directory containing UniProt-specific folders. | +| `-o`, `--output_dir` | Output directory for merged tables. Default: same as `--input_dir`. | +| `--result_filename` | Name of the result file inside each UniProt folder. Default: `result.csv`. | + +--- + +## Output + +By default, outputs are written under: + +```text +reactome_outputs/ +``` + +For a single accession such as `P04637`, the output structure is: + +```text +reactome_outputs/ +├── entries_not_in_reactome.csv # Only created when at least one accession fails +└── P04637/ + ├── result.csv + ├── skipped_reactions.csv # Only created when reactions are skipped + └── pathways_order/ # Only created when pathway ordering is enabled + └── / + ├── graph_edges.csv + ├── graph_nodes.csv + └── ordered_paths.csv +``` + + +### `result.csv` + +`result.csv` is the main output of the workflow. Each row corresponds to a protein entry in a Reactome reaction involving the target protein. + +Common columns include: + +| Column | Description | +|---|---| +| `target_uniprot_ac` | Input UniProt accession. | +| `target_name` | Final target name. Preferentially from UniProt, otherwise from Reactome, otherwise the UniProt accession. | +| `highest_pathway` | Highest-level Reactome pathway. | +| `highest_pathway_id` | Reactome stable ID of the highest-level pathway. | +| `pathway_1`, `pathway_2`, ... | Intermediate pathway hierarchy levels. | +| `pathway_1_id`, `pathway_2_id`, ... | Reactome stable IDs of intermediate pathway levels. | +| `lowest_pathway` | Lowest-level pathway containing the reaction. | +| `lowest_pathway_id` | Reactome stable ID of the lowest-level pathway. | +| `reaction_name` | Reactome reaction name. | +| `reaction_id` | Reactome stable reaction ID. | +| `protein` | Protein entry parsed from BioPAX. | +| `uniprot_ac` | UniProt accession associated with the parsed protein entry. | +| `cellular_location` | Cellular location of the protein entry. | +| `SequenceInterval` | Sequence interval annotation, when available. | +| `SequenceSite` | Sequence site annotation, when available. | +| `Modification_type` | Protein modification annotation, when available. | +| `is_a_protein_family` | Boolean flag indicating whether the BioPAX entry represents a protein family. Protein-family rows are removed from the final output. | +| `complex_of` | Complexes in which the protein participates. | +| `stoichiometry` | Stoichiometric coefficient of the protein within the corresponding complex. | +| `member_physical_entity_of` | Parent protein family or physical entity, when available. | +| `reaction_Left` | Left-side participants of the biochemical reaction. | +| `reaction_Right` | Right-side participants of the biochemical reaction. | +| `reaction_Conversion_Direction` | Biochemical reaction directionality. | +| `Controller_of_reaction_*` | Controllers, activators, or inhibitors associated with the reaction. | +| `disease_name` | Disease annotation/link when available. | +| `ordered` | Boolean flag indicating whether the reaction could be ordered using pathway-ordering files. | + +Additional columns can appear depending on the BioPAX content returned by Reactome. + + +### `skipped_reactions.csv` + +This file is written inside a UniProt-specific output folder when one or more reactions could not be processed. + +Typical reasons include: + +| Reason | Meaning | +|---|---| +| `query_id_returned_none_after_retries` | Reactome did not return usable metadata for a reaction after retries. | +| `biopax_unavailable_after_retries` | The BioPAX model for that reaction could not be downloaded after retries. | + +Common columns: + +| Column | Description | +|---|---| +| `uniprot_ac` | UniProt accession being analyzed. | +| `reaction_id` | Reactome stable reaction ID. | +| `pathway_id` | Reactome stable ID of the pathway associated with the reaction. | +| `pathway_name` | Pathway name associated with the reaction. | +| `reason` | Reason why the reaction was skipped. | + + +### `entries_not_in_reactome.csv` + +This file is written at the global output-directory level when at least one input accession does not produce a valid output. + +Possible statuses: + +| Status | Meaning | +|---|---| +| `not_found_in_reactome` | The UniProt accession had no mapped human Reactome pathways. | +| `no_valid_reactome_output` | Reactome contained mapped pathways for the accession, but no valid reactions remained after filtering. | + +Common columns: + +| Column | Description | +|---|---| +| `uniprot_ac` | UniProt accession. | +| `status` | Failure/status category. | +| `reason` | Explanation of why no final output was produced. | + + +### `pathways_order/` + +This directory is created only when pathway ordering is enabled. + +For each pathway, the workflow writes: + +| File | Description | +|---|---| +| `graph_edges.csv` | Directed edges between reaction nodes. | +| `graph_nodes.csv` | Reaction nodes with `is_start` and `is_end` flags. | +| `ordered_paths.csv` | Ordered reaction paths inferred from the directed graph. | + +When `--skip_pathway_order` is used, this directory is not created and reactions in `result.csv` are marked as unordered. + +### `post process analysis` + +The post-processing script writes: + +| File | Description | +|---|---| +| `merged_reaction.csv` | Deduplicated reaction-level summary across all analyzed UniProt accessions. | +| `merged_highest_pathways.csv` | Deduplicated table of highest-level pathways per target. | +| `disease_single_sequence_site.csv` | Subset of disease-pathway rows where the sequence-site annotation is a single numeric residue position. | + +--- + +## Run + +### Run one UniProt accession + +```bash +python reactome_to_mavisp.py -u P04637 +``` + +### Run one UniProt accession and skip pathway ordering + +```bash +python reactome_to_mavisp.py -u P04637 -s +``` +Skipping pathway ordering is faster because it avoids building pathway-level reaction graphs. + +### Run with input list + +```bash +python reactome_to_mavisp.py -uf uniprot_list.txt -o reactome_outputs +``` +### Run post process analysis + +```bash +python reactome_post_process.py -i reactome_outputs -o reactome_outputs/summary +``` + +--- + +## example + +```bash +# Single protein, faster run without pathway ordering +python reactome_to_mavisp.py -u P04637 -s + +# Multiple proteins +python reactome_to_mavisp.py -uf uniprot_list.txt -o reactome_outputs -s + +# Merge results +python reactome_to_mavisp.py -i reactome_post_process.py -o reactome_outputs/summary +``` + diff --git a/tools/reactome_to_mavisp/example/readme.txt b/tools/reactome_to_mavisp/example/readme.txt new file mode 100644 index 0000000..0a7b042 --- /dev/null +++ b/tools/reactome_to_mavisp/example/readme.txt @@ -0,0 +1 @@ +bash run.sh diff --git a/tools/reactome_to_mavisp/example/run.sh b/tools/reactome_to_mavisp/example/run.sh new file mode 100644 index 0000000..ea2fc4d --- /dev/null +++ b/tools/reactome_to_mavisp/example/run.sh @@ -0,0 +1,2 @@ +source /data/user/marnaudi/reactome/reactome2py/bin/activate +python ../reactome_to_mavisp.py -uf uniprot_list.txt -o reactome_outputs -s diff --git a/tools/reactome_to_mavisp/example/uniprot_list.txt b/tools/reactome_to_mavisp/example/uniprot_list.txt new file mode 100644 index 0000000..a40fe71 --- /dev/null +++ b/tools/reactome_to_mavisp/example/uniprot_list.txt @@ -0,0 +1,2 @@ +Q8N726 +P04637 diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__init__.py b/tools/reactome_to_mavisp/reactome_pipeline/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/__init__.cpython-38.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/__init__.cpython-38.pyc new file mode 100644 index 0000000..7d83f43 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/__init__.cpython-38.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/__init__.cpython-39.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/__init__.cpython-39.pyc new file mode 100644 index 0000000..9f66854 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/__init__.cpython-39.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/data_processing.cpython-38.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/data_processing.cpython-38.pyc new file mode 100644 index 0000000..c314668 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/data_processing.cpython-38.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/data_processing.cpython-39.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/data_processing.cpython-39.pyc new file mode 100644 index 0000000..a535a39 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/data_processing.cpython-39.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/graph_utils.cpython-38.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/graph_utils.cpython-38.pyc new file mode 100644 index 0000000..18d615d Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/graph_utils.cpython-38.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/graph_utils.cpython-39.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/graph_utils.cpython-39.pyc new file mode 100644 index 0000000..c48ff45 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/graph_utils.cpython-39.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/reactome_analysis.cpython-38.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/reactome_analysis.cpython-38.pyc new file mode 100644 index 0000000..5eacaf4 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/reactome_analysis.cpython-38.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/reactome_analysis.cpython-39.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/reactome_analysis.cpython-39.pyc new file mode 100644 index 0000000..02dc1ed Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/reactome_analysis.cpython-39.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/uniprot_utils.cpython-38.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/uniprot_utils.cpython-38.pyc new file mode 100644 index 0000000..56a1b76 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/uniprot_utils.cpython-38.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/uniprot_utils.cpython-39.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/uniprot_utils.cpython-39.pyc new file mode 100644 index 0000000..f321b5b Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/uniprot_utils.cpython-39.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/workflow.cpython-38.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/workflow.cpython-38.pyc new file mode 100644 index 0000000..cbf1c57 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/workflow.cpython-38.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/workflow.cpython-39.pyc b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/workflow.cpython-39.pyc new file mode 100644 index 0000000..60a56b3 Binary files /dev/null and b/tools/reactome_to_mavisp/reactome_pipeline/__pycache__/workflow.cpython-39.pyc differ diff --git a/tools/reactome_to_mavisp/reactome_pipeline/data_processing.py b/tools/reactome_to_mavisp/reactome_pipeline/data_processing.py new file mode 100644 index 0000000..065abce --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_pipeline/data_processing.py @@ -0,0 +1,341 @@ +import re +from typing import Any, Dict, List, Tuple + +import numpy as np +import pandas as pd + +try: + import reactome2py as rc + from reactome2py import analysis, content, utils +except ImportError: + rc = None # type: ignore + +from reactome_pipeline.reactome_analysis import ReactomeAnalysisFunctions + +class DataProcessingFunctions: + """Functions for cleaning and manipulating DataFrames.""" + + @staticmethod + def replace_empty_with_nan(value: Any) -> Any: + if isinstance(value, (set, list)) and len(value) == 0: + return np.nan + return value + + @staticmethod + def process_data(data: List[Dict[str, Any]]) -> pd.DataFrame: + """Flatten nested pathway/reaction annotations into a tabular DataFrame. + + Each row represents one protein entry within one Reactome reaction and includes + pathway hierarchy, reaction metadata, protein features, complex membership, + stoichiometry, regulatory information, and disease annotations when available. + """ + rows: List[Dict[str, Any]] = [] + for entry in data: + pathways_keys = [int(k) for k in entry.keys() if k != 'reaction'] + pathways_keys.sort() + for reaction in entry['reaction']: + proteins = reaction.get('proteins', []) + for protein in proteins: + row: Dict[str, Any] = {} + row["protein"] = protein['display_name'] + row["uniprot_ac"] = protein['uniprot_ac'] + row["cellular_location"] = protein['cellular_location'] + prot_stId = protein['stId'] + + # Protein sequence features and residue modifications + + sequence_interval: List[str] = [] + sequence_site: List[str] = [] + modification_type: List[str] = [] + if 'feature' in protein: + features = protein['feature'] + for feature in features: + + if 'modification_type' in feature: + mod_type_match = re.search(r'"(.*?)"', str(feature['modification_type'])) + if mod_type_match: + modification_type.append(mod_type_match.group(1)) + elif mod_type_match is None: + modification_type.append("could be altered by a mutation event") + else: + raise ValueError("Unexpected modification_type annotation") + position = re.findall(r'\((\d+)\)', str(feature['feature_location'])) + position = list(map(str, position)) + + if "SequenceInterval" in str(feature['feature_location']): + sequence_interval.append("-".join(position)) + if "SequenceSite" in str(feature['feature_location']) and not "SequenceInterval" in str(feature['feature_location']): + sequence_site.append("-".join(position)) + if sequence_interval: + row["SequenceInterval"] = "_".join(sequence_interval) + if sequence_site: + row["SequenceSite"] = "_".join(sequence_site) + if modification_type: + row["Modification_type"] = "_".join(modification_type) + + # Protein family annotation + if not protein['member_physical_entity_of'] and protein['member_physical_entity']: + row['is_a_protein_family'] = True + else: + row['is_a_protein_family'] = False + + # Complex membership and stoichiometry + + # Get the complexes in which the protein participates + complex_of = protein['component_of'] + # Get all the complexes of the reaction + reaction_complexes = reaction.get('complex', {}) + # Annotate all the complexes in which the protein participates + complexes_per_protein: List[str] = [] + # Annotate the stoichiometry of the protein in the corresponding complex + stoichiometry: List[str] = [] + protein_complex_ids = { + prot_complex["stId"] + for prot_complex in complex_of + if prot_complex.get("stId") + } + + prot_reaction_complexes = [ + rcx for rcx in reaction_complexes + if protein_complex_ids.intersection(set(rcx.get("complexes_ids", []))) + ] + for prc in prot_reaction_complexes: + for entity in prc["complex_components"]: + entity_id = entity["stId"] + + if prot_stId == entity_id: + complexes_per_protein.append(str(prc["complex_name"])) + stoichiometry.append(str(entity["stoc_coefficient"])) + + else: + if entity_id not in ReactomeAnalysisFunctions.complex_entities_cache: + ReactomeAnalysisFunctions.complex_entities_cache[entity_id] = ( + ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.entities_complex, + entity_id, + default=[] + ) if rc else [] + ) + + components_in_complex = ReactomeAnalysisFunctions.complex_entities_cache[entity_id] + + for component in components_in_complex: + if component and prot_stId == component.get("stId"): + if str(prc["complex_name"]) not in complexes_per_protein: + complexes_per_protein.append(str(prc["complex_name"])) + stoichiometry.append(str(entity["stoc_coefficient"])) + if complex_of: + row["complex_of"] = "_".join(complexes_per_protein) + row["stoichiometry"] = "_".join(stoichiometry) + + # Parent protein family, when available + + match = re.search(r"\((.*?)\)", str(protein['member_physical_entity_of'])) + if match: + result = match.group(1) + else: + result = None + + row['member_physical_entity_of'] = result + + # Pathway hierarchy annotation + + row['highest_pathway'] = entry[pathways_keys[0]].get('name') + row['highest_pathway_id'] = entry[pathways_keys[0]].get('id') + row['lowest_pathway'] = entry[pathways_keys[-1]].get('name') + row['lowest_pathway_id'] = entry[pathways_keys[-1]].get('id') + for key, index in zip(pathways_keys[1:-1], range(1, len(pathways_keys) - 1)): + row[f'pathway_{index}'] = entry[key].get('name') + row[f'pathway_{index}_id'] = entry[key].get('id') + + # Reaction identifiers + row['reaction_id'] = reaction['stID'] + row['reaction_name'] = reaction["Display_Name"] + + # Biochemical reaction metadata + biochemical = reaction.get('biochemical', {}) + if len(biochemical) > 1: + raise ValueError("Unexpected multiple biochemical entries") + for bio in biochemical: + for key, value in bio.items(): + if isinstance(value, list): + str_value: List[str] = [] + for val in value: + str_value.append(str(val)) + if str_value: + row[f'reaction_{key}'] = " ".join(str_value) + else: + row[f'reaction_{key}'] = [] + else: + row[f'reaction_{key}'] = value + + # Reaction regulation, such as activation or inhibition + + control_information = reaction.get('control_information', []) + + # Store the regulatory role of each controller, such as activation or inhibition. + + for reaction_control in control_information: + controllers: List[str] = [] + + for controller in reaction_control['controller']: + match = re.search(r"\((.*?)\)", str(controller)) + if match: + controllers.append(match.group(1)) + + control_type = reaction_control['control_type'] + row[f'Controller_of_reaction_{control_type}'] = "_".join(controllers) + + # Catalytic metadata, when available + catalytic = reaction.get('catalytic', {}) + for key, value in catalytic.items(): + if isinstance(value, list): + str_value: List[str] = [] + for val in value: + str_value.append(str(val)) + if str_value: + row[f'catalytic_{key}'] = " ".join(str_value) + else: + row[f'catalytic_{key}'] = [] + else: + row[f'catalytic_{key}'] = value + + # Disease annotation, when available + row["disease_name"] = reaction["disease"] + rows.append(row) + df = pd.DataFrame(rows) + df = df.applymap(DataProcessingFunctions.replace_empty_with_nan) + df.drop_duplicates(inplace=True) + return df + + @staticmethod + def get_column_type(col: str) -> str: + """Classify a DataFrame column into a semantic category.""" + + # Pathways: Columns that start with "highest_pathway", "pathway_", or "lowest_pathway" + if "pathway" in col: + if "highest_pathway" in col: + return "highest_pathway" + elif "lowest_pathway" in col: + return "lowest_pathway" + elif col.startswith("pathway_"): + return "pathway" + if col.startswith("disease_name"): + return "disease_name" + + # Reaction: Columns that start with "reaction" + if col.startswith("reaction"): + if "reaction_name" in col: + return "reaction_name" + elif "reaction_id" in col: + return "reaction_id" + elif "reaction_Left" in col: + return "reaction_left" + elif "reaction_Right" in col: + return "reaction_right" + elif "reaction_Conversion_Direction" in col: + return "reaction_conversion_direction" + elif "reaction_controller" in col: + return "reaction_controller" + return "reaction" + # Catalytic: Columns that start with "catalytic" + + if col.startswith("catalytic"): + return "catalytic" + + # Complex or other physical entities: Columns that start with "Complex", "SmallMolecule", or "Protein" + if col.startswith("Complex") or col.startswith("SmallMolecule") or col.startswith("Protein"): + return "Physical_entity" + return "Physical_entity" # Default category for uncategorized columns + + @staticmethod + def column_sort_key(col: str) -> Tuple[int, str]: + """Return a key for sorting columns.""" + column_type = DataProcessingFunctions.get_column_type(col) + # Assign priority for each column type + type_priority = { + "highest_pathway": 0, + "disease_name": 1, + "pathway": 2, + "lowest_pathway": 3, + "reaction_name": 4, + "reaction_id": 5, + "reaction_right": 7, + "reaction_left": 6, + "reaction_conversion_direction": 8, + "reaction_controller": 9, + "reaction": 10, + "catalytic": 11, + "complex": 12, + "Physical_entity": 13, + } + # First sort by category type, then alphabetically within each category + return (type_priority.get(column_type, 99), col) + + @staticmethod + def reorder_dataframe_from_ordered_paths( + df: pd.DataFrame, + lowest_pathway_id: str, + ordered_paths_df: pd.DataFrame + ) -> pd.DataFrame: + + df_filtered = df[df["lowest_pathway_id"] == lowest_pathway_id].copy() + + if df_filtered.empty or ordered_paths_df.empty: + return df_filtered + + reaction_set = set(df_filtered["reaction_name"].dropna().tolist()) + + best_path_id = None + best_overlap = set() + + for path_id, path_df in ordered_paths_df.groupby("path_id"): + path_reactions = set(path_df["reaction"].dropna().tolist()) + overlap = reaction_set.intersection(path_reactions) + + if len(overlap) > len(best_overlap): + best_overlap = overlap + best_path_id = path_id + + if best_path_id is None or not best_overlap: + df_filtered["ordered"] = False + return df_filtered + + best_path_df = ordered_paths_df[ + ordered_paths_df["path_id"] == best_path_id + ].sort_values("step") + + ordered_reactions = [ + r for r in best_path_df["reaction"].tolist() + if r in reaction_set + ] + + missing_reactions = [ + r for r in reaction_set + if r not in ordered_reactions + ] + + final_order = ordered_reactions + missing_reactions + + df_filtered["ordered"] = df_filtered["reaction_name"].isin(ordered_reactions) + + df_filtered["reaction_name"] = pd.Categorical( + df_filtered["reaction_name"], + categories=final_order, + ordered=True + ) + + df_filtered = df_filtered.sort_values("reaction_name") + + return df_filtered + + @staticmethod + def fill_complex_info(group: pd.DataFrame) -> pd.DataFrame: + """Fill missing complex and stoichiometry values in grouped DataFrame.""" + for index, row in group.iterrows(): + if pd.isna(row["complex_of"]): + matching_row = group[group["protein"] == row["member_physical_entity_of"]] + if not matching_row.empty: + group.at[index, "complex_of"] = matching_row.iloc[0]["complex_of"] + group.at[index, "stoichiometry"] = matching_row.iloc[0]["stoichiometry"] + return group \ No newline at end of file diff --git a/tools/reactome_to_mavisp/reactome_pipeline/graph_utils.py b/tools/reactome_to_mavisp/reactome_pipeline/graph_utils.py new file mode 100644 index 0000000..d937c0f --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_pipeline/graph_utils.py @@ -0,0 +1,95 @@ + +from typing import Any, Dict, Iterable, List, Set, Tuple + +import networkx as nx + +class PathwayGraphFunctions: + """Utility functions for pathway graph construction and analysis.""" + + @staticmethod + def make_pathway_graph( + reactions_lists: Iterable[Tuple[str, str, str]], + pathway_id: str + ) -> Tuple[nx.DiGraph, List[str], List[str], List[Dict[str, Any]], List[Dict[str, Any]]]: + + G = nx.DiGraph() + starting_nodes: List[str] = [] + ending_nodes: List[str] = [] + edge_rows: List[Dict[str, Any]] = [] + + for current, next_reaction, previous in reactions_lists: + + # Skip invalid nodes to avoid adding empty nodes to the graph. + if not current: + print(f"[WARNING] Skipping invalid current: {current}, next={next_reaction}, previous={previous}") + continue + + if next_reaction: + G.add_edge(current, next_reaction) + edge_rows.append({ + "pathway_id": pathway_id, + "source": current, + "target": next_reaction + }) + + if previous: + G.add_edge(previous, current) + edge_rows.append({ + "pathway_id": pathway_id, + "source": previous, + "target": current + }) + + if not previous: + starting_nodes.append(current) + + if not next_reaction: + ending_nodes.append(current) + + # Remove duplicated edges while preserving row structure. + edge_rows = list({(e["pathway_id"], e["source"], e["target"]): e for e in edge_rows}.values()) + + node_rows = [] + for node in G.nodes: + node_rows.append({ + "pathway_id": pathway_id, + "node": node, + "is_start": node in starting_nodes, + "is_end": node in ending_nodes + }) + + return G, starting_nodes, ending_nodes, edge_rows, node_rows + + @staticmethod + def remove_duplicates_order(lists: Iterable[List[Any]]) -> List[List[Any]]: + """Remove duplicate lists while preserving order.""" + seen: Set[Tuple[Any, ...]] = set() + result: List[List[Any]] = [] + for sub in lists: + t = tuple(sub) + if t not in seen: + seen.add(t) + result.append(list(sub)) + return result + + @staticmethod + def is_subsequence(sub: List[Any], main: List[Any]) -> bool: + """Return True if ``sub`` is a subsequence of ``main``.""" + it = iter(main) + return all(x in it for x in sub) + + @staticmethod + def find_non_subsequences(lists: Iterable[List[Any]]) -> List[List[Any]]: + """Return those lists that are not subsequences of any other list.""" + lists = list(lists) + non_subsequences: List[List[Any]] = [] + for i, sub in enumerate(lists): + is_sub = False + for j, main in enumerate(lists): + if i != j and PathwayGraphFunctions.is_subsequence(sub, main): + is_sub = True + break + if not is_sub: + non_subsequences.append(sub) + return non_subsequences + diff --git a/tools/reactome_to_mavisp/reactome_pipeline/reactome_analysis.py b/tools/reactome_to_mavisp/reactome_pipeline/reactome_analysis.py new file mode 100644 index 0000000..17427af --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_pipeline/reactome_analysis.py @@ -0,0 +1,659 @@ +import re +import copy +import contextlib +import io +import time +from collections import defaultdict +from typing import Any, Dict, List, Tuple + +import requests + +try: + import reactome2py as rc + from reactome2py import analysis, content, utils +except ImportError: + rc = None # type: ignore + +try: + import pybiopax +except ImportError: + pybiopax = None # type: ignore + +from reactome_pipeline.uniprot_utils import UniprotFunctions + + +class ReactomeAnalysisFunctions: + """Functions for extracting information from BioPAX objects.""" + + biopax_cache: Dict[str, Any] = {} + reactome_call_cache: Dict[Any, Any] = {} + complex_entities_cache: Dict[str, Any] = {} + disease_link_cache: Dict[str, Any] = {} + reaction_target_cache: Dict[Tuple[str, str], bool] = {} + reaction_information_cache: Dict[Tuple[str, Tuple[str, ...]], List[Dict[str, Any]]] = {} + + @staticmethod + def safe_model_from_reactome_cached(st_id: str, retries: int = 5, sleep: int = 3) -> Any: + import time + import requests + + if st_id in ReactomeAnalysisFunctions.biopax_cache: + return ReactomeAnalysisFunctions.biopax_cache[st_id] + + reactome_id = st_id.split("-")[2] + + for attempt in range(1, retries + 1): + try: + model = pybiopax.api.model_from_reactome(reactome_id) + ReactomeAnalysisFunctions.biopax_cache[st_id] = model + return model + + except requests.exceptions.RequestException as e: + print( + f"[WARNING] BioPAX download failed for {st_id} " + f"({attempt}/{retries}): {e}" + ) + + if attempt < retries: + time.sleep(sleep * attempt) + + print(f"[ERROR] Skipping {st_id}: BioPAX model unavailable") + ReactomeAnalysisFunctions.biopax_cache[st_id] = None + return None + + @staticmethod + def safe_reactome_call(func, *args, retries: int = 8, sleep: int = 5, default=None, **kwargs): + import time + import requests + import io + import contextlib + import copy + + func_key = f"{getattr(func, '__module__', '')}.{getattr(func, '__name__', str(func))}" + + cache_key = ( + func_key, + repr(args), + repr(sorted(kwargs.items())) + ) + + if cache_key in ReactomeAnalysisFunctions.reactome_call_cache: + return copy.deepcopy(ReactomeAnalysisFunctions.reactome_call_cache[cache_key]) + + for attempt in range(1, retries + 1): + try: + stdout_buffer = io.StringIO() + + with contextlib.redirect_stdout(stdout_buffer): + result = func(*args, **kwargs) + + printed_output = stdout_buffer.getvalue() + + if printed_output: + print(printed_output, end="") + + if "Status code returned a value of 404" in printed_output or "Not Found" in printed_output: + print(f"[INFO] Reactome returned 404 for {func_key}, args={args}, kwargs={kwargs}. Skipping.") + + ReactomeAnalysisFunctions.reactome_call_cache[cache_key] = copy.deepcopy(default) + return default + + if "Status code returned a value of 429" in printed_output or "Too Many Requests" in printed_output: + wait_time = max(60, sleep * attempt * 10) + + print( + f"[WARNING] Reactome rate limit reached for {func_key}, " + f"args={args}, kwargs={kwargs} ({attempt}/{retries}). " + f"Waiting {wait_time} seconds before retrying." + ) + + if attempt < retries: + time.sleep(wait_time) + continue + + print( + f"[WARNING] Reactome rate limit persisted after {retries} attempts. " + "Returning default." + ) + return default + + + if ( + "Status code returned a value of 500" in printed_output + or "Internal Server Error" in printed_output + or "Status code returned a value of 521" in printed_output + ): + print( + f"[WARNING] Reactome server error for {func_key}, args={args}, kwargs={kwargs} " + f"({attempt}/{retries}). Retrying." + ) + + if attempt < retries: + time.sleep(sleep * attempt) + continue + + print(f"[WARNING] Reactome server error persisted after {retries} attempts. Returning default.") + return default + + if result is None: + print( + f"[INFO] Empty response for {func_key}, args={args}, kwargs={kwargs}. " + f"No status code detected. Skipping." + ) + return default + + time.sleep(0.5) + ReactomeAnalysisFunctions.reactome_call_cache[cache_key] = copy.deepcopy(result) + return result + + except ( + requests.exceptions.ConnectionError, + requests.exceptions.Timeout, + requests.exceptions.ReadTimeout, + ) as e: + wait_time = max(30, sleep * attempt * 5) + print( + f"[WARNING] Connection problem for {func_key}, " + f"args={args}, kwargs={kwargs} ({attempt}/{retries}): {e}. " + f"Waiting {wait_time} seconds before retrying." + ) + + if attempt < retries: + time.sleep(wait_time) + continue + + return default + + except Exception as e: + error_msg = str(e) + + if "404" in error_msg or "Not Found" in error_msg: + print(f"[INFO] Reactome returned 404 for {func_key}, args={args}, kwargs={kwargs}. Skipping.") + + ReactomeAnalysisFunctions.reactome_call_cache[cache_key] = copy.deepcopy(default) + return default + + if "429" in error_msg or "Too Many Requests" in error_msg: + wait_time = max(60, sleep * attempt * 10) + + print( + f"[WARNING] Reactome rate limit reached for {func_key}, " + f"args={args}, kwargs={kwargs} ({attempt}/{retries}). " + f"Waiting {wait_time} seconds before retrying." + ) + + if attempt < retries: + time.sleep(wait_time) + continue + + print( + f"[WARNING] Reactome rate limit persisted after {retries} attempts. " + "Returning default." + ) + return default + + if "500" in error_msg or "Internal Server Error" in error_msg or "521" in error_msg: + print( + f"[WARNING] Reactome server error for {func_key}, args={args}, kwargs={kwargs} " + f"({attempt}/{retries}). Retrying." + ) + + if attempt < retries: + time.sleep(sleep * attempt) + continue + + print(f"[WARNING] Reactome server error persisted after {retries} attempts. Returning default.") + return default + + print(f"[WARNING] Reactome call failed for {func_key}, args={args}, kwargs={kwargs} ({attempt}/{retries}): {e}") + + if attempt < retries: + time.sleep(sleep * attempt) + continue + + return default + + return default + + @staticmethod + def check_attributes(pathway_step_obj: Any) -> None: + for attr in dir(pathway_step_obj): + if not attr.startswith("_"): + try: + getattr(pathway_step_obj, attr) + except Exception: + pass + @staticmethod + def top_level_pathways(pathways_list: List[Dict[str, Any]]) -> List[Dict[int, Dict[str, str]]]: + pathways_dic: List[Dict[int, Dict[str, str]]] = [] + for pathway in pathways_list: + pathway_code = pathway['stId'] + top_level = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.event_ancestors, + id=pathway_code, + default=[] + ) + dic: Dict[int, Dict[str, str]] = {} + for level in top_level: + for index, item in enumerate(level[::-1]): + dic[index] = {'name': item['displayName'], 'id': item['stId']} + pathways_dic.append(dic) + unique_dicts: List[Dict[int, Dict[str, str]]] = [] + seen: Set[frozenset] = set() + for item in pathways_dic: + frozen_item = frozenset((k, frozenset(v.items())) for k, v in item.items()) + if frozen_item not in seen: + seen.add(frozen_item) + unique_dicts.append(item) + return unique_dicts + + @staticmethod + + def collect_disease_link(stId: str) -> str: + if stId in ReactomeAnalysisFunctions.disease_link_cache: + return ReactomeAnalysisFunctions.disease_link_cache[stId] + + links: List[str] = [] + + disease = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.query_id, + id=stId, + enhanced=True, + default=None + ) + + if not disease or not disease.get("isInDisease"): + ReactomeAnalysisFunctions.disease_link_cache[stId] = "" + return "" + + for x in disease.get("crossReference", []): + if x.get("databaseName") == "Mondo" and x.get("url"): + links.append(x.get("url")) + + result = "_".join(links) + + ReactomeAnalysisFunctions.disease_link_cache[stId] = result + return result + + @staticmethod + def get_pathway_order(biopax_obj: Any) -> Dict[str, List[Dict[str, List[str]]]]: + step_processes: Dict[str, List[Dict[str, List[str]]]] = {} + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.BiochemicalReaction): + Step_Process_Of = pathway_step_obj.step_process_of + single_step_processes: List[Dict[str, List[str]]] = [] + if Step_Process_Of: + for process in list(Step_Process_Of): + single_step_process: Dict[str, List[str]] = {} + next_step = process.next_step + next_step_of = process.next_step_of + if next_step: + next_steps: List[str] = [] + for step in list(next_step): + next_steps.extend(i.name for i in step.step_process) + single_step_process['next_step'] = next_steps + if next_step_of: + next_steps_of: List[str] = [] + for step in list(next_step_of): + next_steps_of.extend(i.name for i in step.step_process) + single_step_process['previous'] = next_steps_of + single_step_processes.append(single_step_process) + reaction_name = ( + getattr(pathway_step_obj, "display_name", None) + or " ".join(getattr(pathway_step_obj, "name", []) or []) + ) + + if not reaction_name: + print("[WARNING] Reaction without display_name/name") + continue + + step_processes[reaction_name] = single_step_processes + return step_processes + + @staticmethod + def biopax_model_contains_uniprot(biopax_model: Any, uniprot_ac: str) -> bool: + for protein in biopax_model.get_objects_by_type(pybiopax.biopax.Protein): + + for xref in getattr(protein, "xref", []): + if str(getattr(xref, "id", "")) == uniprot_ac: + return True + + entity_reference = getattr(protein, "entity_reference", None) + if entity_reference is not None: + for xref in getattr(entity_reference, "xref", []): + if str(getattr(xref, "id", "")) == uniprot_ac: + return True + + return False + + @staticmethod + def build_reaction_information( + reaction: Dict[str, Any], + reaction_biopax_model: Any, + pathways_to_not_analyze: List[str] + ) -> List[Dict[str, Any]]: + reaction_stID = reaction["stId"] + reaction_name = reaction["displayName"] + + biochemical_reactions = list( + reaction_biopax_model.get_objects_by_type(pybiopax.biopax.BiochemicalReaction) + ) + + reaction_components: List[Any] = [ + biochemical_reaction.name + for biochemical_reaction in biochemical_reactions + ] + + should_analyze = True + + if len(reaction_components) > 1: + pathways_component: List[str] = [] + + for biochemical_reaction in biochemical_reactions: + for path_name in list(getattr(biochemical_reaction, "pathway_component_of", [])): + pathways_component.append(path_name.display_name) + + if any( + any(pathway_to_skip in pathway_component for pathway_to_skip in pathways_to_not_analyze) + for pathway_component in pathways_component + ): + should_analyze = False + + if not should_analyze: + print( + f"[INFO] Skipping {reaction_stID}: contains pathway components marked as not to analyze." + ) + return [] + + biochemical_information = ReactomeAnalysisFunctions.get_biochemical( + reaction_biopax_model + ) + complexes_information = ReactomeAnalysisFunctions.get_complex( + reaction_biopax_model + ) + control_information = ReactomeAnalysisFunctions.get_control_information( + reaction_biopax_model + ) + pathway_information = ReactomeAnalysisFunctions.get_pathway( + reaction_biopax_model + ) + proteins_information = ReactomeAnalysisFunctions.get_protein( + reaction_biopax_model + ) + disease = ReactomeAnalysisFunctions.collect_disease_link( + reaction_stID + ) + + if len(reaction_components) > 1: + display_names = reaction_components + else: + display_names = [reaction_name] + + reaction_information_list: List[Dict[str, Any]] = [] + + for reaction_name_component in display_names: + if isinstance(reaction_name_component, list): + display_name = " ".join(map(str, reaction_name_component)) + else: + display_name = str(reaction_name_component) + + reaction_information: Dict[str, Any] = { + "stID": reaction_stID, + "Display_Name": display_name, + "biochemical": biochemical_information, + "control_information": control_information, + #"catalytic": catalysis_information, + "complex": complexes_information, + "pathway": pathway_information, + "proteins": proteins_information, + "disease": disease + } + + reaction_information_list.append(reaction_information) + + return reaction_information_list + + @staticmethod + def get_complex(biopax_obj: Any) -> List[Dict[str, Any]]: + reaction_complexes: List[Dict[str, Any]] = [] + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.Complex): + complexes: Dict[str, Any] = {} + Name = pathway_step_obj.name + complex_name = "/".join(Name) + Component_stoichiometry = pathway_step_obj.component_stoichiometry + components: List[Dict[str, Any]] = [] + complex_ids: List[str] = [] + for reac_id in pathway_step_obj.xref: + if reac_id.id.startswith("R-HSA"): + complex_ids.append(reac_id.id) + complex_ids = list(set(complex_ids)) + for stoc in Component_stoichiometry: + st_id = None + stoc_coefficient = stoc.stoichiometric_coefficient + physical_entity = stoc.physical_entity + + for x in getattr(stoc.physical_entity, "xref", []): + xref_id = str(getattr(x, "id", "")) + + if ( + xref_id.startswith("R-HSA-") + or xref_id.startswith("R-ALL-") + or xref_id.startswith("R-COV-") + or xref_id.startswith("R-HIV-") + or xref_id.startswith("R-FLU-") + or xref_id.startswith("R-ECO-") + or xref_id.startswith("R-MTU-") + ): + st_id = xref_id + break + + components.append({ + "physical_entity": physical_entity, + "stoc_coefficient": stoc_coefficient, + "stId": st_id + }) + complexes["complex_name"] = complex_name + complexes["complexes_ids"] = complex_ids + complexes["complex_components"] = components + Participant = pathway_step_obj.participant_of + complexes['Participant'] = [i.display_name for i in Participant] + Evidence = pathway_step_obj.evidence + complexes['evidence'] = Evidence + reaction_complexes.append(complexes) + return reaction_complexes + + @staticmethod + def get_protein(biopax_obj: Any) -> List[Dict[str, Any]]: + proteins: List[Dict[str, Any]] = [] + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.Protein): + protein_features: Dict[str, Any] = {} + cellular_location = pathway_step_obj.cellular_location + entity_reference = pathway_step_obj.entity_reference + component_of_info = [] + + for comp in pathway_step_obj.component_of: + comp_id = None + + for ref in getattr(comp, "xref", []): + ref_id = getattr(ref, "id", None) + if ref_id and ref_id.startswith("R-HSA"): + comp_id = ref_id + break # Use the first Reactome identifier found for this component. + + component_of_info.append({ + "name": getattr(comp, "display_name", str(comp)), + "stId": comp_id + }) + controller_of = pathway_step_obj.controller_of + display_name = pathway_step_obj.display_name + evidence = pathway_step_obj.evidence + feature = pathway_step_obj.feature + get_plain_names = pathway_step_obj.get_plain_names + list_types = pathway_step_obj.list_types + member_physical_entity = pathway_step_obj.member_physical_entity + member_physical_entity_of = pathway_step_obj.member_physical_entity_of + names = pathway_step_obj.name + refs = pathway_step_obj.xref + stId = None + for ref in refs: + if ref.id.startswith("R-HSA"): + stId = ref.id + uniprot_acs: List[str] = [] + if '-' not in display_name: + if entity_reference: + uniprot_acs.append(entity_reference.name[0].split(" ")[0].split(":")[1]) + else: + splitted_names = display_name.split("-") + for splitted_name in splitted_names: + if splitted_name not in ["p", "S", "T", "Y"]: + uniprot_ac = UniprotFunctions.uniprot_gene_to_uniprot_ac(splitted_name) + if uniprot_ac and uniprot_ac not in uniprot_acs: + uniprot_acs.append(uniprot_ac) + if len(uniprot_acs) == 0: + uniprot_acs = [] + if entity_reference: + uniprot_acs.append(entity_reference.name[0].split(" ")[0].split(":")[1]) + not_feature = pathway_step_obj.not_feature + participant_of = pathway_step_obj.participant_of + standard_name = pathway_step_obj.standard_name + uid = pathway_step_obj.uid + uid_index = re.search(r'\d+$', uid).group() + if feature: + features: List[Dict[str, Any]] = [] + for i in feature: + features_dic: Dict[str, Any] = {} + features_dic["feature_location"] = i.feature_location + features_dic["feature_location_type"] = i.feature_location_type + features_dic["feature_of"] = i.feature_of + features_dic["uid"] = i.uid + features.append(features_dic) + if "modification_type" in dir(i): + features_dic["modification_type"] = i.modification_type + protein_features["feature"] = features + protein_features["cellular_location"] = cellular_location + protein_features["component_of"] = component_of_info + protein_features["controller_of"] = controller_of + protein_features["member_physical_entity"] = member_physical_entity + protein_features["member_physical_entity_of"] = member_physical_entity_of + protein_features["display_name"] = display_name + protein_features["available_names"] = "_".join(names) + protein_features["stId"] = stId + protein_features["uniprot_ac"] = "_".join(uniprot_acs) + proteins.append(protein_features) + return proteins + + @staticmethod + def get_control_information(biopax_obj: Any) -> List[Dict[str, Any]]: + reaction_controllers: List[Dict[str, Any]] = [] + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.Control): + reaction_controller: Dict[str, Any] = {} + control_type = pathway_step_obj.control_type + controlled = pathway_step_obj.controlled.display_name + controlled_of = pathway_step_obj.controlled_of + controller = pathway_step_obj.controller + controller_stid = [] + for controller_obj in pathway_step_obj.controller: + for x in controller_obj.xref: + if str(x.id).startswith("R-HSA-"): + controller_stid.append(x.id) + display_name = pathway_step_obj.display_name + get_plain_names = pathway_step_obj.get_plain_names() + interaction_type = pathway_step_obj.interaction_type + participant = pathway_step_obj.participant + participant_of = pathway_step_obj.participant_of + pathway_component_of = pathway_step_obj.pathway_component_of + standard_name = pathway_step_obj.standard_name + evidence = pathway_step_obj.evidence + step_process_of = pathway_step_obj.step_process_of + reaction_controller['control_type'] = control_type + reaction_controller['controlled'] = controlled + reaction_controller['controlled_of'] = controlled_of + reaction_controller['controller'] = controller + reaction_controller['controller_stid'] = controller_stid + reaction_controller['participant'] = participant + reaction_controller['pathway_component_of'] = pathway_component_of + reaction_controller['evidence'] = evidence + reaction_controller['control_type'] = control_type + reaction_controller['interaction_type'] = interaction_type + reaction_controllers.append(reaction_controller) + return reaction_controllers + + @staticmethod + def get_catalysis(biopax_obj: Any) -> Dict[str, List[Any]]: + all_catalysis: Dict[str, List[Any]] = defaultdict(list) + + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.Catalysis): + catalysis_direction = pathway_step_obj.catalysis_direction + cofactor = pathway_step_obj.cofactor + control_type = pathway_step_obj.control_type + controlled = pathway_step_obj.controlled.display_name + controlled_of = pathway_step_obj.controlled_of + controller = pathway_step_obj.controller + display_name = pathway_step_obj.display_name + interaction_type = pathway_step_obj.interaction_type + participant = pathway_step_obj.participant + pathway_component_of = pathway_step_obj.pathway_component_of + standard_name = pathway_step_obj.standard_name + evidence = pathway_step_obj.evidence + + all_catalysis["display_name"].append(display_name) + all_catalysis["standard_name"].append(standard_name) + all_catalysis["catalysis_direction"].append(catalysis_direction) + all_catalysis["cofactor"].append(cofactor) + all_catalysis["control_type"].append(control_type) + all_catalysis["controlled"].append(controlled) + all_catalysis["controlled_of"].append(controlled_of) + all_catalysis["controller"].append(controller) + all_catalysis["participant"].append(participant) + all_catalysis["pathway_component_of"].append(pathway_component_of) + all_catalysis["evidence"].append(evidence) + all_catalysis["interaction_type"].append(interaction_type) + + return dict(all_catalysis) + + @staticmethod + def get_pathway(biopax_obj: Any) -> None: + # This function in the original script prints debug information and + # returns nothing. It has been retained for completeness but does not + # contribute to the data collection. + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.Pathway): + dir(pathway_step_obj) + return None + + @staticmethod + def get_biochemical(biopax_obj: Any) -> List[Dict[str, Any]]: + all_reactions_parameters: List[Dict[str, Any]] = [] + for pathway_step_obj in biopax_obj.get_objects_by_type(pybiopax.biopax.BiochemicalReaction): + reaction_parameters: Dict[str, Any] = {} + Name = pathway_step_obj.name + reaction_name = " ".join(Name) + Display_Name = pathway_step_obj.display_name + Left = pathway_step_obj.left + Right = pathway_step_obj.right + Participant = pathway_step_obj.participant + Evidence = pathway_step_obj.evidence + Standard_Name = pathway_step_obj.standard_name + Conversion_Direction = pathway_step_obj.conversion_direction + Delta_G = pathway_step_obj.delta_g + Delta_H = pathway_step_obj.delta_h + Delta_S = pathway_step_obj.delta_s + K_eq = pathway_step_obj.k_e_q + Pathway_Component_Of = pathway_step_obj.pathway_component_of + pathways_component: List[str] = [] + if Pathway_Component_Of: + for path_name in list(Pathway_Component_Of): + pathways_component.append(path_name.display_name) + reaction_parameters['Left'] = Left + reaction_parameters['Right'] = Right + reaction_parameters['Participant'] = Participant + reaction_parameters['Pathway_Component_Of'] = pathways_component + reaction_parameters['Pathways_component'] = pathways_component + reaction_parameters['Evidence'] = Evidence + reaction_parameters['Standard_Name'] = Standard_Name + reaction_parameters['Conversion_Direction'] = Conversion_Direction + reaction_parameters['K_eq'] = K_eq + reaction_parameters['Delta_G'] = Delta_G + reaction_parameters['Delta_H'] = Delta_H + reaction_parameters['Delta_S'] = Delta_S + all_reactions_parameters.append(reaction_parameters) + return all_reactions_parameters + diff --git a/tools/reactome_to_mavisp/reactome_pipeline/uniprot_utils.py b/tools/reactome_to_mavisp/reactome_pipeline/uniprot_utils.py new file mode 100644 index 0000000..e31ee38 --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_pipeline/uniprot_utils.py @@ -0,0 +1,101 @@ +import re +from typing import List, Optional + +import pandas as pd +import requests + +class UniprotFunctions: + """Helper functions for UniProt lookups.""" + + @staticmethod + def uniprot_gene_to_uniprot_ac(gene_name: str) -> Optional[str]: + gene_name = re.sub(r"\(.*", "", str(gene_name)).strip() + + params = { + "query": f"gene:{gene_name}", + "fields": "accession,gene_names,protein_name,organism_name", + "format": "json" + } + + try: + response = requests.get( + "https://rest.uniprot.org/uniprotkb/search", + params=params, + timeout=20 + ) + + if response.status_code != 200: + return None + + data = response.json() + + for entry in data.get("results", []): + if ( + entry.get("organism", {}).get("scientificName") == "Homo sapiens" + and entry.get("entryType") == "UniProtKB reviewed (Swiss-Prot)" + ): + return entry.get("primaryAccession") + + return None + + except requests.exceptions.RequestException: + return None + + @staticmethod + def uniprot_id_to_uniprot_ac(uniprot_id: str) -> Optional[str]: + uniprot_url = f"https://rest.uniprot.org/uniprotkb/{uniprot_id}.json" + uniprot_response = requests.get(uniprot_url) + if uniprot_response.status_code == 200: + gene_information = uniprot_response.json() + if gene_information['organism']['scientificName'] == 'Homo sapiens' and gene_information['entryType'] == "UniProtKB reviewed (Swiss-Prot)": + return gene_information['primaryAccession'] + return None + else: + return None + @staticmethod + def uniprot_ac_to_protein_name(uniprot_ac: str) -> Optional[str]: + uniprot_url = f"https://rest.uniprot.org/uniprotkb/{uniprot_ac}.json" + + try: + response = requests.get(uniprot_url, timeout=20) + + if response.status_code != 200: + return None + + data = response.json() + + protein_description = data.get("proteinDescription", {}) + + recommended_name = protein_description.get("recommendedName", {}) + full_name = recommended_name.get("fullName", {}) + + protein_name = full_name.get("value") + + if protein_name: + return protein_name + + submission_names = protein_description.get("submissionNames", []) + if submission_names: + return submission_names[0].get("fullName", {}).get("value") + + return None + + except requests.exceptions.RequestException: + return None + + @staticmethod + def process_gene_names(gene_names: str) -> str: + if pd.isna(gene_names) or not isinstance(gene_names, str): + return "None_uniprot_ac" + genes = gene_names.split('_') + accessions: List[str] = [] + for gene in genes: + single_proteins = gene.split("-") + for protein in single_proteins: + try: + accession = UniprotFunctions.uniprot_gene_to_uniprot_ac(protein.strip().upper()) + if accession and accession not in accessions: + accessions.append(str(accession)) + except Exception: + pass + return "_".join(accessions) \ No newline at end of file diff --git a/tools/reactome_to_mavisp/reactome_pipeline/workflow.py b/tools/reactome_to_mavisp/reactome_pipeline/workflow.py new file mode 100644 index 0000000..275e6fb --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_pipeline/workflow.py @@ -0,0 +1,660 @@ +import os +from collections import defaultdict +from typing import Any, Dict, List, Set, Tuple + + +import pandas as pd +import networkx as nx + +try: + import reactome2py as rc + from reactome2py import analysis, content, utils +except ImportError: + rc = None # type: ignore + +try: + import pybiopax +except ImportError: + pybiopax = None # type: ignore + +from reactome_pipeline.graph_utils import PathwayGraphFunctions +from reactome_pipeline.data_processing import DataProcessingFunctions +from reactome_pipeline.reactome_analysis import ReactomeAnalysisFunctions +from reactome_pipeline.uniprot_utils import UniprotFunctions + +class ReactomeScript: + """Encapsulate the original script's workflow in an object oriented form.""" + + def __init__( + self, + uniprot_ac: str, + skip_pathway_order: bool = False, + output_dir: str = "." + ) -> None: + self.uniprot_ac = uniprot_ac + self.skip_pathway_order = skip_pathway_order + self.output_dir = output_dir + os.makedirs(self.output_dir, exist_ok=True) + + def output_path(self, *parts: str) -> str: + return os.path.join(self.output_dir, *parts) + + def retrieve_reactome_pathways(self) -> List[Dict[str, Any]]: + return ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.mapping, + id=self.uniprot_ac, + resource="UniProt", + species="9606", + by="pathways", + default=[] + ) + def resolve_target_information(self) -> Dict[str, Any]: + """ + Retrieve target-related information. + + search_fireworks is used only as an optional extra check. + If it fails, the workflow continues using: + 1. UniProt-derived protein name. + 2. Reactome complexes associated with the UniProt accession. + 3. Direct UniProt accession matching inside BioPAX reaction models. + """ + target_protein_stId: List[str] = [] + target_protein_stId_other: List[str] = [] + + # Primary target name source: UniProt. + uniprot_target_name = UniprotFunctions.uniprot_ac_to_protein_name( + self.uniprot_ac + ) + + # Optional target name source: Reactome search_fireworks. + reactome_target_name = None + + target_protein_info = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.search_fireworks, + query=self.uniprot_ac, + default={"entries": []} + ) + + entries = target_protein_info.get("entries", []) if target_protein_info else [] + + if not entries: + print( + f"[WARNING] search_fireworks failed or returned no entries for " + f"{self.uniprot_ac}. Continuing without Reactome protein-form IDs." + ) + + for entry in entries: + if entry.get("stId"): + target_protein_stId.append(entry["stId"]) + + if reactome_target_name is None: + reactome_target_name = ( + entry.get("displayName") + or entry.get("name") + ) + + # Final target name priority: + # 1. UniProt protein name + # 2. Reactome display/name + # 3. UniProt accession fallback + target_name = ( + uniprot_target_name + or reactome_target_name + or self.uniprot_ac + ) + + # Optional extra check: alternative Reactome forms of the target protein. + for prot_id in target_protein_stId: + other_forms = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.entity_other_form, + id=prot_id, + default=[] + ) + + if other_forms: + for form in other_forms: + if form.get("stId"): + target_protein_stId_other.append(form["stId"]) + + all_target_protein_stId = list( + set(target_protein_stId + target_protein_stId_other) + ) + + # Main useful information: complexes associated with the UniProt accession. + target_protein_complexes_names: List[str] = [] + + complexes_list = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.entities_complexes, + id=self.uniprot_ac, + resource="UniProt", + default=[] + ) + + for complex_entry in complexes_list: + if complex_entry.get("stId"): + target_protein_complexes_names.append(complex_entry["stId"]) + + return { + "target_name": target_name, + "target_protein_stId": target_protein_stId, + "target_protein_stId_other": target_protein_stId_other, + "all_target_protein_stId": all_target_protein_stId, + "target_protein_complexes_names": target_protein_complexes_names, + } + def write_pathway_ordering(self, reactome_pathways: List[Dict[str, Any]]) -> None: + """ + Build and write reaction ordering information for each Reactome pathway. + """ + reactions_order: Dict[str, Dict[str, List[Dict[str, List[str]]]]] = {} + + for i in reactome_pathways: + biopax_model = ReactomeAnalysisFunctions.safe_model_from_reactome_cached(i['stId']) + + if biopax_model is None: + print(f"[WARNING] Skipping pathway order for {i['stId']}: BioPAX unavailable.") + continue + + reaction_order = ReactomeAnalysisFunctions.get_pathway_order(biopax_model) + reactions_order[i['stId']] = reaction_order + + for hrsa, reac in reactions_order.items(): + reactions_lists: List[Tuple[str, str, str]] = [] + + for key, value_list in reac.items(): + for value in value_list: + next_steps = value.get('next_step', []) + previous_steps = value.get('previous', []) + + if next_steps and previous_steps: + for next_step in next_steps: + for next_item in next_step: + for previous_step in previous_steps: + for prev_item in previous_step: + reactions_lists.append((key, next_item, prev_item)) + + if not next_steps and previous_steps: + for previous_step in previous_steps: + for prev_item in previous_step: + reactions_lists.append((key, "", prev_item)) + + if next_steps and not previous_steps: + for next_step in next_steps: + for next_item in next_step: + reactions_lists.append((key, next_item, "")) + + G, starting_nodes, ending_nodes, edge_rows, node_rows = PathwayGraphFunctions.make_pathway_graph( + reactions_lists, + hrsa + ) + + pathway_output_dir = self.output_path("pathways_order", hrsa) + os.makedirs(pathway_output_dir, exist_ok=True) + + pd.DataFrame(edge_rows).to_csv( + os.path.join(pathway_output_dir, "graph_edges.csv"), + index=False + ) + + pd.DataFrame(node_rows).to_csv( + os.path.join(pathway_output_dir, "graph_nodes.csv"), + index=False + ) + + simple_paths: List[List[str]] = [] + + for start_node in starting_nodes: + for end_node in ending_nodes: + if start_node != end_node and start_node in G.nodes and end_node in G.nodes: + paths = list(nx.all_simple_paths(G, source=start_node, target=end_node)) + simple_paths.extend(paths) + + paths_to_process = PathwayGraphFunctions.remove_duplicates_order(simple_paths) + filtered_paths = PathwayGraphFunctions.find_non_subsequences(paths_to_process) + + if filtered_paths: + ordered_paths_rows = [] + + for path_id, path in enumerate(filtered_paths): + for step, reaction in enumerate(path): + ordered_paths_rows.append({ + "pathway_id": hrsa, + "path_id": path_id, + "step": step, + "reaction": reaction + }) + + pd.DataFrame(ordered_paths_rows).to_csv( + os.path.join(pathway_output_dir, "ordered_paths.csv"), + index=False + ) + + def collect_candidate_reactions( + self, + reactome_pathways: List[Dict[str, Any]], + skipped_reactions: List[Dict[str, str]] + ) -> Tuple[List[Dict[int, Any]], List[str]]: + """ + Retrieve candidate reactions from the lowest-level Reactome pathways. + + Returns: + pathways_with_reactions: + Pathway hierarchy dictionaries with an added 'reactions' field. + pathways_to_not_analyze: + Lowest-level pathway names used to avoid treating pathways as reactions. + """ + hierarchy_to_low_level = ReactomeAnalysisFunctions.top_level_pathways(reactome_pathways) + + # ------------------------------ Filter steps ----------------------------- # + + # Lowest-level pathway names are excluded later to avoid treating pathway + # containers as actual reactions. + + pathways_to_not_analyze: List[str] = [] + for dic in hierarchy_to_low_level: + last_pathway = list(dic.keys())[-1] + pathways_to_not_analyze.append(dic[last_pathway]['name']) + + # Retrieve and keep only true reaction events from each lowest-level pathway. + + pathways_with_reactions: List[Dict[int, Any]] = [] + for dic in hierarchy_to_low_level: + last_pathway = list(dic.keys())[-1] + std_id = dic[last_pathway]['id'] + reactions = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.pathway_contained_event, + id=std_id, + default=[] + ) + filtered_reactions: List[Dict[str, Any]] = [] + for reaction in reactions: + + if not isinstance(reaction, int) and reaction['displayName'] not in pathways_to_not_analyze: + filtered_reactions.append(reaction) + elif isinstance(reaction, int): + rhsa = "R-HSA-" + str(reaction) + reaction = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.query_id, + id=rhsa, + enhanced=False, + default=None + ) + if reaction is None: + print(f"[WARNING] Skipping reaction {rhsa}: Reactome returned no usable data.") + skipped_reactions.append({ + "uniprot_ac": self.uniprot_ac, + "reaction_id": rhsa, + "pathway_id": std_id, + "pathway_name": dic[last_pathway].get("name", ""), + "reason": "query_id_returned_none_after_retries" + }) + + continue + if reaction.get('displayName') not in pathways_to_not_analyze and reaction.get('className') == "Reaction": + filtered_reactions.append(reaction) + dic['reactions'] = filtered_reactions + pathways_with_reactions.append(dic) + return pathways_with_reactions, pathways_to_not_analyze + + def is_target_reaction( + self, + reaction: Dict[str, Any], + target_info: Dict[str, Any] + ) -> bool: + """ + Return True if the reaction contains the target protein, one of its + Reactome protein forms, or one of its associated complexes. + """ + rsa = reaction["stId"] + + reaction_cache_key = (self.uniprot_ac, rsa) + + if reaction_cache_key in ReactomeAnalysisFunctions.reaction_target_cache: + return ReactomeAnalysisFunctions.reaction_target_cache[ + reaction_cache_key + ] + + is_target_reaction = False + + # Retrieve physical entities participating in the reaction. + components = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.participants_physical_entities, + rsa, + default=[] + ) + + components_name: List[str] = [] + for component in components: + if component and component.get("stId"): + components_name.append(component["stId"]) + + # Detect reactions containing any complex associated with the target protein. + if any( + target in component_id + for component_id in components_name + for target in target_info["target_protein_complexes_names"] + ): + is_target_reaction = True + + # Detect reactions containing any Reactome form of the target protein. + reaction_proteins = ReactomeAnalysisFunctions.safe_model_from_reactome_cached(rsa) + + if reaction_proteins is not None: + protein_stid: List[str] = [] + + for protein in reaction_proteins.get_objects_by_type(pybiopax.biopax.Protein): + for protein_ref in getattr(protein, "xref", []): + protein_ref_id = getattr(protein_ref, "id", "") + + if protein_ref_id.startswith("R-HSA"): + protein_stid.append(protein_ref_id) + + if any( + target in protein_id + for protein_id in protein_stid + for target in target_info["all_target_protein_stId"] + ): + is_target_reaction = True + + if ReactomeAnalysisFunctions.biopax_model_contains_uniprot( + reaction_proteins, + self.uniprot_ac + ): + is_target_reaction = True + + ReactomeAnalysisFunctions.reaction_target_cache[ + reaction_cache_key + ] = is_target_reaction + + return is_target_reaction + + + + def extract_target_reaction_information( + self, + pathways_with_reactions: List[Dict[int, Any]], + pathways_to_not_analyze: List[str], + target_info: Dict[str, Any], + skipped_reactions: List[Dict[str, str]] + ) -> List[Dict[Any, Any]]: + """ + Filter candidate reactions for the target protein and extract BioPAX-level + reaction information. + """ + pathways_with_reactions_information: List[Dict[Any, Any]] = [] + for pathway in pathways_with_reactions: + pathways_with_reaction_information: Dict[Any, Any] = {} + pathways_keys = [key for key in pathway.keys() if key != 'reactions'] + reactions_to_analyze: List[Dict[str, Any]] = [] + unique_reactions = { + reaction["stId"]: reaction + for reaction in pathway["reactions"] + if reaction and reaction.get("stId") + } + + reactions = list(unique_reactions.values()) + + for reaction in reactions: + if self.is_target_reaction( + reaction, + target_info + ): + reactions_to_analyze.append(reaction) + + # Extract BioPAX-level annotations for target reactions. + + reactions_information: List[Dict[str, Any]] = [] + + for reaction in reactions_to_analyze: + reaction_stID = reaction["stId"] + + reaction_info_cache_key = ( + reaction_stID, + tuple(sorted(pathways_to_not_analyze)) + ) + + if reaction_info_cache_key in ReactomeAnalysisFunctions.reaction_information_cache: + reactions_information.extend( + ReactomeAnalysisFunctions.reaction_information_cache[reaction_info_cache_key] + ) + continue + + reaction_biopax_model = ReactomeAnalysisFunctions.safe_model_from_reactome_cached( + reaction_stID + ) + + if reaction_biopax_model is None: + print(f"[WARNING] Skipping reaction {reaction_stID}: BioPAX unavailable.") + skipped_reactions.append({ + "uniprot_ac": self.uniprot_ac, + "reaction_id": reaction_stID, + "pathway_id": pathway[pathways_keys[-1]].get("id", ""), + "pathway_name": pathway[pathways_keys[-1]].get("name", ""), + "reason": "biopax_unavailable_after_retries" + }) + ReactomeAnalysisFunctions.reaction_information_cache[reaction_info_cache_key] = [] + continue + + reaction_information_list = ReactomeAnalysisFunctions.build_reaction_information( + reaction, + reaction_biopax_model, + pathways_to_not_analyze + ) + + ReactomeAnalysisFunctions.reaction_information_cache[reaction_info_cache_key] = reaction_information_list + reactions_information.extend(reaction_information_list) + for key in pathways_keys: + pathways_with_reaction_information[key] = pathway[key] + pathways_with_reaction_information["reaction"] = reactions_information + pathways_with_reactions_information.append(pathways_with_reaction_information) + return pathways_with_reactions_information + + def remove_duplicate_pathway_reactions( + self, + result_df: pd.DataFrame + ) -> pd.DataFrame: + """ + Remove duplicated reaction rows caused by overlapping Reactome pathway + hierarchies. + + When the same reaction appears in multiple nested pathways, the row + associated with the shorter, less specific pathway path is removed. + """ + grouped = result_df.groupby("reaction_name") + paths_to_remove: List[List[str]] = [] + for reaction in grouped.groups.keys(): + reaction_df = grouped.get_group(reaction) + # list of R-HSA id of the lowest-level pathway + lowest_pathway_id = list(set(reaction_df["lowest_pathway_id"].to_list())) + # list of R-HSA id of reaction + reaction_id = list(set(reaction_df["reaction_id"].to_list())) + # find the duplicated reactions + reaction_dict: Dict[str, Set[str]] = defaultdict(set) + if len(lowest_pathway_id) != len(reaction_id): + # create a dictionary with the reaction name and the R-HSA id of + # the lowest pathway that contains the same reaction + for q in zip(reaction_df["lowest_pathway_id"].to_list(), reaction_df["reaction_id"].to_list()): + reaction_dict[q[1]].add(q[0]) + reaction_with_multiple_pathways = {k: list(v) for k, v in reaction_dict.items() if len(v) > 1} + + # for each double reaction and R-HSA id of the lowest pathway to compare: + # extract the df and build a string with all the R-HSA ids of the pathways from that reaction, + # put the two strings to be compared in a list + + for reaction_id_val, pathways in reaction_with_multiple_pathways.items(): + paths_to_check: List[str] = [] + for pathway in pathways: + df_for_comparison = result_df[(result_df["reaction_name"] == reaction) & (result_df["reaction_id"] == reaction_id_val) & (result_df["lowest_pathway_id"] == pathway)] + df_for_comparison = df_for_comparison.astype(str) + pathway_cols = [col for col in result_df.columns if "pathway" in col and "id" in col] + pathway_cols.append("reaction_id") + for _, row in df_for_comparison.iterrows(): + combined_pathway = "_".join(row[pathway_cols].dropna()) + paths_to_check.append(combined_pathway) + paths_to_remove.append(list(set(paths_to_check))) + + # Identify the shorter duplicated pathway path and mark its + # lowest_pathway_id/reaction_id pair for removal. + + lines_to_remove: List[Tuple[str, str]] = [] + for paths in paths_to_remove: + compare_pathways = [path.split("_") for path in paths] + processed_pathways: List[List[str]] = [] + for pathway in compare_pathways: + pathway_no_nan: List[str] = [] + for reaction_elem in pathway: + if reaction_elem != "nan": + pathway_no_nan.append(reaction_elem) + processed_pathways.append(pathway_no_nan) + to_remove = min(processed_pathways, key=len) + lines_to_remove.append((to_remove[-2], to_remove[-1])) + + result_df_filtered = result_df[~result_df.apply(lambda row: (row["lowest_pathway_id"], row["reaction_id"]) in set(lines_to_remove), axis=1)] + + return result_df_filtered + + def order_reactions_by_pathway_files( + self, + result_df_filtered: pd.DataFrame + ) -> pd.DataFrame: + """ + Reorder reactions using pathway ordering files when available. + + If no ordering file exists for a pathway, rows are kept in their original + order and marked as unordered. + """ + ordered_result_df_filtered = pd.DataFrame() + + for stid in set(result_df_filtered["lowest_pathway_id"].dropna().tolist()): + + ordered_paths_file = self.output_path( + "pathways_order", + stid, + "ordered_paths.csv" + ) + + if os.path.exists(ordered_paths_file): + ordered_paths_df = pd.read_csv(ordered_paths_file) + + df_sorted = DataProcessingFunctions.reorder_dataframe_from_ordered_paths( + result_df_filtered, + stid, + ordered_paths_df + ) + + else: + df_sorted = result_df_filtered[ + result_df_filtered["lowest_pathway_id"] == stid + ].copy() + df_sorted["ordered"] = False + + ordered_result_df_filtered = pd.concat( + [ordered_result_df_filtered, df_sorted], + ignore_index=True + ) + return ordered_result_df_filtered + + + def build_and_write_result( + self, + pathways_with_reactions_information: List[Dict[Any, Any]], + target_info: Dict[str, Any], + skipped_reactions: List[Dict[str, str]] + ) -> bool: + """ + Build the final result DataFrame, clean duplicated pathway rows, + order reactions when pathway ordering is available, and write outputs. + """ + if skipped_reactions: + pd.DataFrame(skipped_reactions).drop_duplicates().to_csv( + self.output_path("skipped_reactions.csv"), + index=False + ) + result_df = DataProcessingFunctions.process_data(pathways_with_reactions_information) + if result_df.empty: + print(f"[INFO] No Reactome reactions found for {self.uniprot_ac}.") + return False + + sorted_columns = sorted(result_df.columns, key=DataProcessingFunctions.column_sort_key) + result_df = result_df[sorted_columns] + + # Propagate complex and stoichiometry annotations within each reaction group. + if "complex_of" in result_df.columns: + result_df = result_df.groupby("reaction_name", group_keys=False).apply(DataProcessingFunctions.fill_complex_info) + + # Remove rows representing protein families rather than individual proteins. + result_df = result_df.loc[result_df['is_a_protein_family'] == False] + if result_df.empty: + print( + f"[INFO] No Reactome reactions left for {self.uniprot_ac} " + "after removing protein families." + ) + return False + + result_df_filtered = self.remove_duplicate_pathway_reactions(result_df) + + ordered_result_df_filtered = self.order_reactions_by_pathway_files( + result_df_filtered + ) + + if ordered_result_df_filtered.empty: + print(f"[INFO] Final Reactome output is empty for {self.uniprot_ac}.") + return False + + ordered_result_df_filtered.insert( + 0, + "target_uniprot_ac", + self.uniprot_ac + ) + + ordered_result_df_filtered.insert( + 1, + "target_name", + target_info["target_name"] + ) + + ordered_result_df_filtered.to_csv( + self.output_path("result.csv"), + sep=",", + index=False + ) + + return True + + def run(self) -> bool: + """Execute the full Reactome analysis workflow.""" + if rc is None or pybiopax is None: + raise ImportError( + "reactome2py and pybiopax must be installed to run this script" + ) + + skipped_reactions: List[Dict[str, str]] = [] + + reactome_pathways = self.retrieve_reactome_pathways() + + if not self.skip_pathway_order: + self.write_pathway_ordering(reactome_pathways) + + target_info = self.resolve_target_information() + + pathways_with_reactions, pathways_to_not_analyze = ( + self.collect_candidate_reactions( + reactome_pathways, + skipped_reactions + ) + ) + + pathways_with_reactions_information = ( + self.extract_target_reaction_information( + pathways_with_reactions, + pathways_to_not_analyze, + target_info, + skipped_reactions + ) + ) + + return self.build_and_write_result( + pathways_with_reactions_information, + target_info, + skipped_reactions + ) + \ No newline at end of file diff --git a/tools/reactome_to_mavisp/reactome_post_process.py b/tools/reactome_to_mavisp/reactome_post_process.py new file mode 100644 index 0000000..daaadc2 --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_post_process.py @@ -0,0 +1,356 @@ +#!/usr/bin/env python3 + +from __future__ import annotations + +import argparse +import re +from pathlib import Path +from typing import Dict, List, Optional + +import pandas as pd + + +MERGED_REACTION_BASE_COLUMNS = [ + "target_uniprot_ac", + "target_name", + "highest_pathway", + "highest_pathway_id", + "disease_name", +] + +MERGED_REACTION_TAIL_COLUMNS = [ + "lowest_pathway", + "lowest_pathway_id", + "reaction_name", + "reaction_id", + "reaction_Left", + "reaction_Right", + "reaction_Conversion_Direction", +] + +HIGHEST_PATHWAY_COLUMNS = [ + "highest_pathway", + "highest_pathway_id", + "target_uniprot_ac", + "target_name", +] + + +def parse_arguments() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Post-process Reactome result.csv files generated inside " + "UniProt-specific folders. The script creates: " + "merged_reaction.csv, merged_highest_pathways.csv, and " + "disease_single_sequence_site.csv." + ) + ) + + parser.add_argument( + "-i", + "--input_dir", + required=True, + type=str, + help="Main Reactome output directory containing UniProt-specific folders." + ) + + parser.add_argument( + "-o", + "--output_dir", + required=False, + default=None, + type=str, + help=( + "Directory where merged output files will be written. " + "Default: same as --input_dir." + ) + ) + + parser.add_argument( + "--result_filename", + required=False, + default="result.csv", + type=str, + help="Name of the result file inside each UniProt folder. Default: result.csv" + ) + + return parser.parse_args() + + +def find_result_files(input_dir: Path, result_filename: str) -> List[Path]: + """ + Find result files directly inside first-level UniProt subfolders. + + Example: + input_dir/P04637/result.csv + input_dir/Q8N726/result.csv + """ + return sorted(input_dir.glob(f"*/{result_filename}")) + + +def infer_target_uniprot_ac(result_file: Path) -> str: + """Infer the target UniProt accession from the parent folder name.""" + return result_file.parent.name + + +def harmonize_target_columns(df: pd.DataFrame, result_file: Path) -> pd.DataFrame: + """ + Ensure target_uniprot_ac and target_name are present. + + Priority for target_uniprot_ac: + 1. existing target_uniprot_ac column + 2. existing uniprot_ac column + 3. parent folder name + + Priority for target_name: + 1. existing target_name column + 2. existing protein column + 3. empty string + """ + df = df.copy() + + if "target_uniprot_ac" not in df.columns: + if "uniprot_ac" in df.columns: + df["target_uniprot_ac"] = df["uniprot_ac"] + else: + df["target_uniprot_ac"] = infer_target_uniprot_ac(result_file) + + if "target_name" not in df.columns: + if "protein" in df.columns: + df["target_name"] = df["protein"] + else: + df["target_name"] = "" + + return df + + +def load_result_files(result_files: List[Path]) -> pd.DataFrame: + dfs: List[pd.DataFrame] = [] + + for result_file in result_files: + try: + df = pd.read_csv(result_file) + + if df.empty: + print(f"[INFO] Skipping empty file: {result_file}") + continue + + df = harmonize_target_columns(df, result_file) + df["source_result_file"] = str(result_file) + df["source_folder"] = result_file.parent.name + dfs.append(df) + + except Exception as error: + print(f"[WARNING] Could not read {result_file}: {error}") + + if not dfs: + return pd.DataFrame() + + return pd.concat(dfs, ignore_index=True) + + +def get_pathway_columns(df: pd.DataFrame) -> List[str]: + """ + Return pathway_1/pathway_1_id ... pathway_n/pathway_n_id in numeric order. + """ + pathway_numbers = set() + + for column in df.columns: + match = re.fullmatch(r"pathway_(\d+)(_id)?", str(column)) + if match: + pathway_numbers.add(int(match.group(1))) + + pathway_columns: List[str] = [] + + for number in sorted(pathway_numbers): + pathway_col = f"pathway_{number}" + pathway_id_col = f"pathway_{number}_id" + + if pathway_col in df.columns: + pathway_columns.append(pathway_col) + + if pathway_id_col in df.columns: + pathway_columns.append(pathway_id_col) + + return pathway_columns + + +def select_existing_columns(df: pd.DataFrame, columns: List[str]) -> pd.DataFrame: + """ + Select requested columns, creating missing columns as empty strings. + + This makes the script robust if one result.csv lacks a specific column. + """ + df = df.copy() + + for column in columns: + if column not in df.columns: + df[column] = "" + + return df[columns] + + +def find_sequence_site_column(df: pd.DataFrame) -> Optional[str]: + """ + Find sequence_site column independently of case/style. + + It accepts: + sequence_site + SequenceSite + """ + candidates: Dict[str, str] = { + str(col).lower(): str(col) + for col in df.columns + } + + for key in ["sequence_site", "sequencesite"]: + if key in candidates: + return candidates[key] + + return None + + +def is_single_numeric_sequence_site(value: object) -> bool: + """ + True only for a single integer-like residue position. + + Accepted: + 330 + "330" + + Rejected: + "330_331" + "330-331" + "330;331" + empty / NaN + """ + if pd.isna(value): + return False + + value_str = str(value).strip() + + return bool(re.fullmatch(r"\d+", value_str)) + + +def write_outputs(concatenated_df: pd.DataFrame, output_dir: Path) -> None: + output_dir.mkdir(parents=True, exist_ok=True) + + pathway_columns = get_pathway_columns(concatenated_df) + + merged_reaction_columns = ( + MERGED_REACTION_BASE_COLUMNS + + pathway_columns + + MERGED_REACTION_TAIL_COLUMNS + ) + + merged_reaction_df = select_existing_columns( + concatenated_df, + merged_reaction_columns + ).drop_duplicates() + + merged_reaction_file = output_dir / "merged_reaction.csv" + + merged_reaction_df.to_csv( + merged_reaction_file, + index=False + ) + + merged_highest_pathways_df = select_existing_columns( + concatenated_df, + HIGHEST_PATHWAY_COLUMNS + ).drop_duplicates() + + merged_highest_pathways_file = output_dir / "merged_highest_pathways.csv" + + merged_highest_pathways_df.to_csv( + merged_highest_pathways_file, + index=False + ) + + sequence_site_column = find_sequence_site_column(concatenated_df) + + if sequence_site_column is None: + print( + "[WARNING] No sequence_site/SequenceSite column found. " + "Writing empty disease_single_sequence_site.csv." + ) + disease_single_sequence_site_df = concatenated_df.iloc[0:0].copy() + + else: + disease_single_sequence_site_df = concatenated_df[ + ( + concatenated_df["highest_pathway"] + .astype(str) + .str.strip() + .eq("Disease") + ) + & concatenated_df[sequence_site_column].apply( + is_single_numeric_sequence_site + ) + ].copy().drop_duplicates() + + disease_single_sequence_site_file = output_dir / "disease_single_sequence_site.csv" + + disease_single_sequence_site_df.to_csv( + disease_single_sequence_site_file, + index=False + ) + + print(f"[INFO] Output written: {merged_reaction_file}") + print(f"[INFO] Rows: {len(merged_reaction_df)}") + + print(f"[INFO] Output written: {merged_highest_pathways_file}") + print(f"[INFO] Rows: {len(merged_highest_pathways_df)}") + + print(f"[INFO] Output written: {disease_single_sequence_site_file}") + print(f"[INFO] Rows: {len(disease_single_sequence_site_df)}") + + +def main() -> None: + args = parse_arguments() + + input_dir = Path(args.input_dir).resolve() + + if args.output_dir: + output_dir = Path(args.output_dir).resolve() + else: + output_dir = input_dir + + if not input_dir.exists(): + raise FileNotFoundError(f"Input directory does not exist: {input_dir}") + + if not input_dir.is_dir(): + raise NotADirectoryError(f"Input path is not a directory: {input_dir}") + + result_files = find_result_files( + input_dir=input_dir, + result_filename=args.result_filename + ) + + if not result_files: + print(f"[WARNING] No {args.result_filename} files found inside {input_dir}") + return + + concatenated_df = load_result_files(result_files) + + if concatenated_df.empty: + print("[WARNING] No valid result files could be loaded.") + return + + if "highest_pathway" not in concatenated_df.columns: + raise ValueError("Missing required column: highest_pathway") + + write_outputs( + concatenated_df=concatenated_df, + output_dir=output_dir + ) + + print(f"[INFO] Found result files: {len(result_files)}") + print( + "[INFO] Concatenated rows before filtering/deduplication: " + f"{len(concatenated_df)}" + ) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/tools/reactome_to_mavisp/reactome_to_mavisp.py b/tools/reactome_to_mavisp/reactome_to_mavisp.py new file mode 100644 index 0000000..1037aa4 --- /dev/null +++ b/tools/reactome_to_mavisp/reactome_to_mavisp.py @@ -0,0 +1,212 @@ +""" +Reactome analysis workflow for UniProt accessions. + +Given one UniProt accession, or a file containing multiple UniProt accessions, +the script queries Reactome to retrieve associated pathways, identifies +candidate reactions involving the target protein, extracts BioPAX-level +protein/complex/reaction annotations, and writes cleaned CSV outputs. + +Main outputs +------------ +For each UniProt accession, the script creates a dedicated output folder +containing: + +- result.csv + Final cleaned table of Reactome reactions involving the target protein. + +- skipped_reactions.csv + Reactions that could not be processed because Reactome or BioPAX returned + no usable data after retries. + +- pathways_order/ + Optional pathway-ordering files generated when pathway ordering is not + skipped. + +At the global output level, the script can also write: + +- entries_not_in_reactome.csv + UniProt accessions that were not found in Reactome or did not produce a + valid final output. + +Usage +----- +Run the workflow for one UniProt accession: + + python reac_classified.py -u P04637 + +Run the workflow for multiple UniProt accessions: + + python reac_classified.py -uf uniprot_list.txt + +Skip pathway ordering: + + python reac_classified.py -u P04637 -s +""" +from __future__ import annotations + +import os +import argparse +import shutil +from typing import Dict, List + +import pandas as pd + +try: + import reactome2py as rc + from reactome2py import analysis, content, utils +except ImportError: + rc = None # type: ignore + +from reactome_pipeline.reactome_analysis import ReactomeAnalysisFunctions +from reactome_pipeline.workflow import ReactomeScript + +def parse_arguments() -> argparse.Namespace: + """Parse command line arguments.""" + parser = argparse.ArgumentParser( + description=( + 'This script retrieves pathways associated with a UniProt accession from Reactome, ' + 'orders the reactions, collects protein level data and writes several CSV outputs. ' + 'It is a class‑based refactoring of the original reac.py script.' + ) + ) + parser.add_argument( + '-uf', '--uniprot_file', + dest='uniprot_file', + required=False, + type=str, + help='Text file containing one UniProt accession per line.' + ) + + parser.add_argument( + '-o', '--output_dir', + dest='output_dir', + default='reactome_outputs', + type=str, + help='Main output directory.' + ) + parser.add_argument( + '-u', '--uniprot_ac', dest='uniprot_ac', default='Q8N726', type=str, + help='UniProt accession of the protein of interest (default: Q8N726).' + ) + parser.add_argument( + '-s', '--skip_pathway_order', dest='skip_pathway_order', action='store_true', + help='Skip the pathway ordering phase (analogous to supplying --pdb_file in the original script).' + ) + return parser.parse_args() + +def read_uniprot_accessions(uniprot_file: str) -> List[str]: + accessions: List[str] = [] + + with open(uniprot_file) as handle: + for line in handle: + accession = line.strip() + + if not accession: + continue + + if accession.startswith("#"): + continue + + accessions.append(accession) + + return accessions + +def uniprot_has_reactome_entry(uniprot_ac: str) -> bool: + """ + Return True if Reactome contains at least one mapped pathway for the UniProt accession. + """ + if rc is None: + raise ImportError("reactome2py is required to query Reactome.") + + reactome_pathways = ReactomeAnalysisFunctions.safe_reactome_call( + rc.content.mapping, + id=uniprot_ac, + resource="UniProt", + species="9606", + by="pathways", + default=[] + ) + + return bool(reactome_pathways) + + +def write_missing_reactome_entries( + output_dir: str, + missing_entries: List[Dict[str, str]] +) -> None: + """ + Write UniProt accessions not producing a valid Reactome output. + """ + if not missing_entries: + return + + os.makedirs(output_dir, exist_ok=True) + + output_file = os.path.join(output_dir, "entries_not_in_reactome.csv") + + new_df = pd.DataFrame(missing_entries) + + if os.path.exists(output_file): + old_df = pd.read_csv(output_file) + final_df = pd.concat([old_df, new_df], ignore_index=True) + final_df = final_df.drop_duplicates(subset=["uniprot_ac", "status"]) + else: + final_df = new_df.drop_duplicates(subset=["uniprot_ac", "status"]) + + final_df.to_csv(output_file, index=False) + +def main() -> None: + args = parse_arguments() + + if args.uniprot_file: + uniprot_accessions = read_uniprot_accessions(args.uniprot_file) + else: + uniprot_accessions = [args.uniprot_ac] + + missing_reactome_entries: List[Dict[str, str]] = [] + + for uniprot_ac in uniprot_accessions: + print(f"[INFO] Running Reactome analysis for {uniprot_ac}") + + protein_output_dir = os.path.join(args.output_dir, uniprot_ac) + + # First check: is this UniProt accession known by Reactome? + if not uniprot_has_reactome_entry(uniprot_ac): + print(f"[INFO] {uniprot_ac} is not present in Reactome.") + + missing_reactome_entries.append({ + "uniprot_ac": uniprot_ac, + "status": "not_found_in_reactome", + "reason": "Reactome search_fireworks returned no entries" + }) + + continue + + workflow = ReactomeScript( + uniprot_ac=uniprot_ac, + skip_pathway_order=args.skip_pathway_order, + output_dir=protein_output_dir + ) + + has_valid_output = workflow.run() + + if not has_valid_output: + print(f"[INFO] Removing empty output folder for {uniprot_ac}") + + if os.path.isdir(protein_output_dir): + shutil.rmtree(protein_output_dir) + + missing_reactome_entries.append({ + "uniprot_ac": uniprot_ac, + "status": "no_valid_reactome_output", + "reason": "Reactome entry found, but no valid reactions remained after filtering" + }) + + write_missing_reactome_entries( + output_dir=args.output_dir, + missing_entries=missing_reactome_entries + ) + + +if __name__ == '__main__': + main() \ No newline at end of file