Merge pull request #54 from haddocking/interaction-partners

Interaction partners

Author: sverhoeven

README.md

@@ -10,7 +10,7 @@ Python package to detect proteins in EM density maps.
 
 It uses
 
-- [Uniprot Sparql endpoint](https://sparql.uniprot.org/) to search for proteins and their measured or predicted 3D structures.
+- [protein-quest](https://github.com/haddocking/protein-quest) to search, retrieve and filter protein structures from Uniprot, PDBe and AlphaFold DB.
 - [powerfit](https://pypi.org/project/powerfit-em/) to fit protein structure in a Electron Microscopy (EM) density map.
 
 An example workflow:
@@ -60,6 +60,26 @@ protein-detective search \
 
 In `./mysession` directory, you will find session.db file, which is a [DuckDB](https://duckdb.org/) database with search results.
 
+<details>
+<summary>You can also include interaction partners in the search</summary>
+
+```shell
+protein-detective --log-level INFO search \
+    --taxon-id 9606 \
+    --reviewed \
+    --subcellular-location-uniprot nucleus \
+    --subcellular-location-go GO:0005634 \
+    --molecular-function-go GO:0003677 \
+    --interaction-partner-seed A8MT69 \
+    --interaction-partner-exclude B1APH4 \
+    --limit 100 \
+    ./mysession2
+```
+
+Which will add `Q96H22` which is an interaction partner of `A8MT69` in a macromolecular complex.
+
+</details>
+
 ### To retrieve a bunch of structures
 
 ```shell

pyproject.toml

@@ -17,7 +17,7 @@ dependencies = [
     "molviewspec>=1.6.0",
     "pandas>=2.3.0",
     "powerfit-em[opencl]>=3.0.6",
-    "protein-quest>=0.3.3",
+    "protein-quest>=0.4.0",
     "psutil>=7.0.0",
     "rich>=14.0.0",
     "rich-argparse>=1.7.1",

src/protein_detective/__version__.py

@@ -1 +1 @@
-__version__ = "0.4.0"
+__version__ = "0.5.0"

src/protein_detective/cli.py

@@ -8,7 +8,6 @@
 from protein_quest.alphafold.fetch import downloadable_formats
 from protein_quest.converter import converter
 from protein_quest.ss import SecondaryStructureFilterQuery
-from protein_quest.uniprot import Query
 from rich import print as rprint
 from rich.logging import RichHandler
 from rich_argparse import RawDescriptionRichHelpFormatter, RichHelpFormatter
@@ -19,6 +18,7 @@
     add_powerfit_parser,
     handle_powerfit,
 )
+from protein_detective.search import UniprotQuery
 from protein_detective.workflow import (
     filter_structures,
     retrieve_structures,
@@ -50,6 +50,23 @@ def add_search_parser(subparsers):
         action="append",
         help="Molecular function (GO term, e.g. GO:0003677). Can be specified multiple times.",
     )
+    parser.add_argument(
+        "--interaction-partner-seed",
+        type=str,
+        action="append",
+        help=dedent("""\
+            UniProt ID to use as interaction partner seed.
+            The search will be expanded to include structures identifiers of the found interaction partners.
+            Can be specified multiple times.
+        """),
+    )
+    parser.add_argument(
+        "--interaction-partner-exclude",
+        type=str,
+        action="append",
+        help="UniProt ID to exclude as found interaction partners. Can be specified multiple times.",
+    )
+
     parser.add_argument("--limit", type=int, default=10_000, help="Limit number of results")
 
 
@@ -131,20 +148,27 @@ def add_filter_parser(subparsers: argparse._SubParsersAction):
 
 
 def handle_search(args):
-    query = Query(
-        taxon_id=args.taxon_id,
-        reviewed=args.reviewed,
-        subcellular_location_uniprot=args.subcellular_location_uniprot,
-        subcellular_location_go=args.subcellular_location_go,
-        molecular_function_go=args.molecular_function_go,
+    query = converter.structure(
+        {
+            "taxon_id": args.taxon_id,
+            "reviewed": args.reviewed,
+            "subcellular_location_uniprot": args.subcellular_location_uniprot,
+            "subcellular_location_go": args.subcellular_location_go,
+            "molecular_function_go": args.molecular_function_go,
+            "interaction_partner_seeds": args.interaction_partner_seed or [],
+            "interaction_partner_excludes": args.interaction_partner_exclude or [],
+        },
+        UniprotQuery,
     )
     session_dir = Path(args.session_dir)
-    nr_uniprot, nr_pdbes, nr_prot2pdbes, nr_afs = search_structures_in_uniprot(query, session_dir, limit=args.limit)
+    result = search_structures_in_uniprot(query, session_dir, limit=args.limit)
     rprint(
-        f"Search completed: {nr_uniprot} UniProt entries found, "
-        f"{nr_pdbes} PDBe structures, {nr_prot2pdbes} UniProt to PDB mappings, "
-        f"{nr_afs} AlphaFold structures."
+        f"Search completed: {result.nr_uniprot_accessions} UniProt entries found, "
+        f"{result.nr_pdbs} PDBe structures, {result.nr_prot2pdb} UniProt to PDB mappings, "
+        f"{result.nr_afs} AlphaFold structures."
     )
+    if query.interaction_partner_seeds:
+        rprint(f"Included {result.nr_interaction_partners} Uniprot entries found as interaction partners.")
 
 
 def handle_retrieve(args):

src/protein_detective/db.py

@@ -7,17 +7,17 @@
 from importlib.resources import read_text
 from pathlib import Path
 
-from cattrs import unstructure
 from duckdb import ConstraintException, DuckDBPyConnection, InvalidInputException
 from duckdb import connect as duckdb_connect
 from pandas import DataFrame
 from protein_quest.alphafold.entry_summary import EntrySummary
 from protein_quest.alphafold.fetch import AlphaFoldEntry
 from protein_quest.converter import converter
-from protein_quest.uniprot import PdbResult, Query
+from protein_quest.uniprot import PdbResult
 
 from protein_detective.filter import FilteredStructure, FilterOptions
 from protein_detective.powerfit.options import PowerfitOptions
+from protein_detective.search import UniprotQuery
 
 logger = logging.getLogger(__name__)
 
@@ -111,14 +111,15 @@ def connect(session_dir: Path, read_only: bool = False) -> Iterator[DuckDBPyConn
     con.close()
 
 
-def save_query(query: Query, con: DuckDBPyConnection):
+def save_query(query: UniprotQuery, con: DuckDBPyConnection):
     """Save a UniProt search query to the database.
 
     Args:
         query: The UniProt search query to save.
         con: The DuckDB connection to use for saving the data.
     """
-    con.execute("INSERT INTO uniprot_searches (query) VALUES (?)", (unstructure(query),))
+    value = converter.dumps(query).decode()
+    con.execute("INSERT INTO uniprot_searches (query) VALUES (?)", (value,))
 
 
 def save_uniprot_accessions(uniprot_accessions: Iterable[str], con: DuckDBPyConnection) -> int:

src/protein_detective/search.py

@@ -0,0 +1,54 @@
+"""Module with search logic."""
+
+import logging
+from dataclasses import dataclass, field
+
+from protein_quest.uniprot import Query, search4macromolecular_complexes
+
+logger = logging.getLogger(__name__)
+
+
+@dataclass
+class UniprotQuery(Query):
+    """A UniProt search query with interaction partner options.
+
+    Parameters:
+        interaction_partner_seeds: A set of UniProt accessions to search for interaction partners.
+        interaction_partners_excludes: A set of UniProt accessions to exclude from interaction partner results.
+    """
+
+    interaction_partner_seeds: set[str] = field(default_factory=set)
+    interaction_partner_excludes: set[str] = field(default_factory=set)
+
+
+def search_for_interaction_partners(query: UniprotQuery, limit: int) -> set[str]:
+    """Searches for interaction partners in UniProt database and ComplexPortal.
+
+    Args:
+        query: The search query containing seeds and excludes.
+        limit: The maximum number of results to return from the database query.
+
+    Returns:
+        A set of unique UniProt accessions of interaction partners found.
+    """
+    if not query.interaction_partner_seeds:
+        logger.info("No interaction partner seeds provided; skipping search for interaction partners.")
+        return set()
+    logger.info("Searching for interaction partners of seeds %s", query.interaction_partner_seeds)
+    uniprot_accessions_of_partners: set[str] = set()
+
+    complexes = search4macromolecular_complexes(query.interaction_partner_seeds, limit)
+    for complex_entry in complexes:
+        uniprot_accessions_of_partners.update(complex_entry.members)
+
+    # Exclude seeds and excludes from results
+    uniprot_accessions_of_partners.difference_update(query.interaction_partner_seeds)
+    uniprot_accessions_of_partners.difference_update(query.interaction_partner_excludes)
+
+    logger.info(
+        "Found %d unique interaction partners in %d macromolecular complexes after excluding %d accessions",
+        len(uniprot_accessions_of_partners),
+        len(complexes),
+        len(query.interaction_partner_excludes),
+    )
+    return uniprot_accessions_of_partners

src/protein_detective/workflow.py

@@ -2,6 +2,7 @@
 
 import asyncio
 import logging
+from dataclasses import dataclass
 from pathlib import Path
 from typing import Literal
 
@@ -10,7 +11,7 @@
 from protein_quest.alphafold.fetch import fetch_many_async as af_fetch
 from protein_quest.alphafold.fetch import relative_to as af_relative_to
 from protein_quest.pdbe.fetch import fetch as pdbe_fetch
-from protein_quest.uniprot import Query, search4af, search4pdb, search4uniprot
+from protein_quest.uniprot import search4af, search4pdb, search4uniprot
 
 from protein_detective.db import (
     connect,
@@ -33,11 +34,31 @@
     filter_alphafold_structures,
     filter_pdbe_structures,
 )
+from protein_detective.search import UniprotQuery, search_for_interaction_partners
 
 logger = logging.getLogger(__name__)
 
 
-def search_structures_in_uniprot(query: Query, session_dir: Path, limit: int = 10_000) -> tuple[int, int, int, int]:
+@dataclass
+class UniprotSearchResult:
+    """Result of a UniProt search.
+
+    Parameters:
+        nr_uniprot_accessions: Number of UniProt accessions found.
+        nr_pdbs: Number of PDB structures found.
+        nr_prot2pdb: Number of UniProt to PDB mappings found.
+        nr_afs: Number of AlphaFold structures found.
+        nr_interaction_partners: Number of interaction partners found.
+    """
+
+    nr_uniprot_accessions: int
+    nr_pdbs: int
+    nr_prot2pdb: int
+    nr_afs: int
+    nr_interaction_partners: int
+
+
+def search_structures_in_uniprot(query: UniprotQuery, session_dir: Path, limit: int = 10_000) -> UniprotSearchResult:
     """Searches for protein structures in UniProt database.
 
     Args:
@@ -51,8 +72,12 @@ def search_structures_in_uniprot(query: Query, session_dir: Path, limit: int = 1
         and the number of AlphaFold structures found.
     """
     session_dir.mkdir(parents=True, exist_ok=True)
-
     uniprot_accessions = search4uniprot(query, limit)
+    logger.debug(uniprot_accessions)
+    uniprot_accessions_of_partners = search_for_interaction_partners(query, limit)
+    logger.debug(uniprot_accessions_of_partners)
+    nr_interaction_partners = len(uniprot_accessions_of_partners)
+    uniprot_accessions.update(uniprot_accessions_of_partners)
     pdbs = search4pdb(uniprot_accessions, limit=limit)
     af_result = search4af(uniprot_accessions, limit=limit)
 
@@ -62,7 +87,13 @@ def search_structures_in_uniprot(query: Query, session_dir: Path, limit: int = 1
         nr_pdbs, nr_prot2pdb = save_pdbs(pdbs, con)
         nr_afs = save_alphafolds(af_result, con)
 
-    return len(uniprot_accessions), nr_pdbs, nr_prot2pdb, nr_afs
+    return UniprotSearchResult(
+        nr_uniprot_accessions=len(uniprot_accessions),
+        nr_pdbs=nr_pdbs,
+        nr_prot2pdb=nr_prot2pdb,
+        nr_afs=nr_afs,
+        nr_interaction_partners=nr_interaction_partners,
+    )
 
 
 WhatRetrieve = Literal["pdbe", "alphafold"]