"""Core module for extracting Protein-nucleic acid interfaces."""
from Bio.PDB.Polypeptide import PPBuilder, Polypeptide
from Bio.PDB.NeighborSearch import NeighborSearch
from Bio.PDB.PDBExceptions import PDBConstructionException
from Bio.PDB.Structure import Structure
from Bio.PDB.Model import Model
from Bio.PDB.Chain import Chain
from Bio.PDB.Residue import Residue
from Bio.PDB.Atom import Atom
import pandas as pd
from PDBNucleicAcids.NucleicAcid import NABuilder, NucleicAcid
from PDBNucleicAcids.NucleicAcid import DSNABuilder, DoubleStrandNucleicAcid
from PDBNucleicAcids.BasePairRules import BasePairRules
from PDBNucleicAcids.BasePairRules import WatsonCrickBasePairRules
import copy
import warnings
[docs]
class Interface:
"""
Extract Protein-nucleic acid interface.
Parameters
----------
contacts : list[tuple[Atom, Atom]]
List of pairs of atoms, first one is from the nucleic acids,
second one is from the protein.
search_radius : float | int, optional
Search radius, measured in Angstrom, within which Protein-nucleic acid
interactions are found. Default is 4.0
"""
def __init__(
self,
contacts: list[tuple[Atom, Atom]],
search_radius: int | float = 4.0,
) -> None:
# save input
self._contacts = contacts
# save parameters
self.search_radius = search_radius
# initialize other important arguments
# they will be filled with correct data during building
# by the class InterfaceBuilder
self._binding_domain: Polypeptide | None = None
self._bound_na_list: list[NucleicAcid] | None = None
self._bound_dsna_list: list[DoubleStrandNucleicAcid] | None = None
def __repr__(self) -> str:
"""Return string representation of the nucleic acid."""
pp_chain_id = self.get_protein_atoms()[0].parent.parent.id # type: ignore
na_chain_ids = [
na.get_chain_id()
for na in self.get_bound_nucleic_acids() # type: ignore
]
return f"<Interface chains={pp_chain_id}:\
{''.join(na_chain_ids)} contacts={len(self._contacts)} search_radius=\
{self.search_radius}>"
[docs]
def get_protein_atoms(self) -> list[Atom]:
"""
Get only protein atoms in the protein-nucleic acid interface.
Returns
-------
list[Atom]
List of protein atoms in the interface.
"""
return list({atom_pair[1] for atom_pair in self._contacts})
[docs]
def get_nucleic_acid_atoms(self) -> list[Atom]:
"""
Get only nucleic acid atoms in the protein-nucleic acid interface.
Returns
-------
list[Atom]
List of nucleic acid atoms in the interface.
"""
return list({atom_pair[0] for atom_pair in self._contacts})
[docs]
def get_aminoacids(self) -> list[Residue]:
"""
Get only protein residues in the protein-nucleic acid interface.
Returns
-------
list[Residue]
List of protein residues in the interface.
"""
return list(
{
atom_pair[1].parent # type: ignore
for atom_pair in self._contacts
}
)
[docs]
def get_nucleotides(self) -> list[Residue]:
"""
Get only nucleic acid residues in the protein-nucleic acid interface.
Returns
-------
list[Residue]
List of nucleic acid residues in the interface.
"""
return list(
{
atom_pair[0].parent # type: ignore
for atom_pair in self._contacts
}
)
[docs]
def as_dataframe(self) -> pd.DataFrame:
"""
Get all data from the interface, as a dataframe.
Contains the following data fields:
Residue hetero field
Residue number
Residue insertion code
Residue name
Atom name
Atom alternate location
Atom element
Atomic coordinates (x, y, z)
From both protein and nucleic acid atoms
Euclidean distance between atom pair in contact
Returns
-------
df : pd.DataFrame
All data from the interface.
"""
protein_chain_id: str = self.get_aminoacids()[0].id
data = []
for na_atom, prot_atom in self._contacts:
prot_res_hetfield = prot_atom.parent.id[0] # type: ignore
prot_res_number = prot_atom.parent.id[1] # type: ignore
prot_res_icode = prot_atom.parent.id[2] # type: ignore
prot_res_name = prot_atom.parent.resname # type: ignore
prot_atom_name = prot_atom.name
prot_atom_altloc = prot_atom.altloc
prot_atom_element = prot_atom.element
prot_atom_coord_x = prot_atom.coord[0]
prot_atom_coord_y = prot_atom.coord[1]
prot_atom_coord_z = prot_atom.coord[2]
dna_chain_id = na_atom.parent.parent.id # type: ignore
dna_res_hetfield = na_atom.parent.id[0] # type: ignore
dna_res_number = na_atom.parent.id[1] # type: ignore
dna_res_icode = na_atom.parent.id[2] # type: ignore
dna_res_name = na_atom.parent.resname # type: ignore
dna_atom_name = na_atom.name
dna_atom_altloc = na_atom.altloc
dna_atom_element = na_atom.element
dna_atom_coord_x = na_atom.coord[0]
dna_atom_coord_y = na_atom.coord[1]
dna_atom_coord_z = na_atom.coord[2]
euclidean_distance = na_atom - prot_atom
row = (
protein_chain_id,
prot_res_hetfield,
prot_res_number,
prot_res_icode,
prot_res_name,
prot_atom_name,
prot_atom_altloc,
prot_atom_element,
prot_atom_coord_x,
prot_atom_coord_y,
prot_atom_coord_z,
dna_chain_id,
dna_res_hetfield,
dna_res_number,
dna_res_icode,
dna_res_name,
dna_atom_name,
dna_atom_altloc,
dna_atom_element,
dna_atom_coord_x,
dna_atom_coord_y,
dna_atom_coord_z,
euclidean_distance,
)
data.append(row)
df = pd.DataFrame(
data,
columns=[
"protein_chain_id",
"prot_res_hetfield",
"prot_res_number",
"prot_res_icode",
"prot_res_name",
"prot_atom_name",
"prot_atom_altloc",
"prot_atom_element",
"prot_atom_coord_x",
"prot_atom_coord_y",
"prot_atom_coord_z",
"dna_chain_id",
"dna_res_hetfield",
"dna_res_number",
"dna_res_icode",
"dna_res_name",
"dna_atom_name",
"dna_atom_altloc",
"dna_atom_element",
"dna_atom_coord_x",
"dna_atom_coord_y",
"dna_atom_coord_z",
"euclidean_distance",
],
)
return df
[docs]
def get_binding_domain(self) -> Polypeptide | None:
"""
Get nucleic acid binding domain from the protein.
The output is the binding "gapped" subsequence of the full protein
found in the structure.
This method allows for "gaps" of unbound aminoacids inside the
binding domain, only the aminoacids at the ends are trimmed according
to being bound to nucleic acids (NAs) or not.
A visual example of "gaps":
```
Input full protein: MQMLLNHKPTKFNGAIDERFHWKVIQRISGSEG
NA-bound: **** **
Output binding domain: FNGAIDER
```
This method is only an inference of the NA-binding domain:
while the output will likely align with the annotated true domain,
it'll likely not infer the whole domain. This is because a domain is
defined by folding properties, while this method is much more naive.
This is why I implemented some "padding" on both ends of the binding
domain, it allows to be more lenient of the extent of the binding
domain.
Returns
-------
binding_domain : Polypeptide | None
Nucleic acid binding domain.
"""
return self._binding_domain
[docs]
def get_bound_nucleic_acids(self) -> list[NucleicAcid] | None:
"""
Get all nucleic acids bound by the protein.
The output nucleic acids (NAs) are subsequences
of the full NAs found in the structure,
since proteins might not bind the whole NA.
This method allows for "gaps" of unbound nucleotides inside the
NA, only the nucleotides at the ends are trimmed according
to being protein-bound or not.
A visual example of "gaps":
```
Input full NA: GATATACAAGCCA
Protein-bound: **** **
Output protein-bound NA: TATACAAG
```
Returns
-------
list[NucleicAcid] | None
List of nucleic acids bound by the protein.
"""
return self._bound_na_list
[docs]
def get_bound_double_strands(self) -> list[DoubleStrandNucleicAcid] | None:
"""
Get all double strand nucleic acids bound by the protein.
The output double strand nucleic acids (NAs) are subsequences
of the full DSNAs found in the structure,
since proteins might not bind the whole DSNA.
This method allows for "gaps" of unbound base pairs inside the
DSNA, only the base pairs at the ends are trimmed according
to being protein-bound or not.
A visual example of "gaps":
```
Input full DSNA: GATATACAAGCCA
|||||||||||||
TGGCTTGTATATC
Protein-bound: **** **
Output protein-bound DSNA: TATACAAG
||||||||
CTTGTATA
```
Returns
-------
list[NucleicAcid] | None
List of double strand nucleic acids bound by the protein.
"""
return self._bound_dsna_list
# def fixed_protein_atoms_number(self, num_atoms) -> None:
# """Filter contacts by a fixed number of protein atoms."""
# # cast list into dataframe, ready to be sorted
# df = pd.DataFrame(self.contacts, columns=["na_atom", "protein_atom"])
# df["euclidean_distance"] = df.apply(
# lambda row: row["na_atom"] - row["protein_atom"], axis=1
# )
# # aggregate: for each atom, its minimum distance from DSNA
# agg = df.groupby(["protein_atom"]).min()
# agg = agg.reset_index()
# agg = agg.sort_values(by="euclidean_distance", ascending=True)
# # get closest n atoms to DSNA
# top_protein_atoms = agg.head(num_atoms)["protein_atom"].tolist()
# if len(top_protein_atoms) <= num_atoms:
# raise Exception("Not enough atoms.")
# # select contacts by top n atoms
# selected_contacts = [
# (na_atom, protein_atom)
# for na_atom, protein_atom in self.contacts
# if protein_atom in top_protein_atoms
# ]
# self.contacts = selected_contacts
# def fixed_na_atoms_number(self, num_atoms) -> None:
# """Filter contacts by a fixed number of nucleic acid atoms."""
# # cast list into dataframe, ready to be sorted
# df = pd.DataFrame(self.contacts, columns=["na_atom", "protein_atom"])
# df["euclidean_distance"] = df.apply(
# lambda row: row["na_atom"] - row["protein_atom"], axis=1
# )
# # aggregate: for each atom, its minimum distance from DSNA
# agg = df.groupby(["na_atom"]).min()
# agg = agg.reset_index()
# agg = agg.sort_values(by="euclidean_distance", ascending=True)
# # get closest n atoms to DSNA
# top_na_atoms = agg.head(num_atoms)["na_atom"].tolist()
# if len(top_na_atoms) <= num_atoms:
# raise Exception("Not enough atoms.")
# # select contacts by top n atoms
# selected_contacts = [
# (na_atom, protein_atom)
# for na_atom, protein_atom in self.contacts
# if na_atom in top_na_atoms
# ]
# self.contacts = selected_contacts
[docs]
class InterfaceBuilder:
"""
Use atomic distance to find Protein-Nucleic acid interfaces.
Assuming you *only* want standard nucleotides and amino acids.
Parameters
----------
search_radius : float | int, optional
Search radius, measured in Angstrom, within which
Protein-Nucleic acid interactions are found. Default is 4.0
"""
def __init__(self, search_radius: float | int = 4.0) -> None:
self.search_radius = search_radius
def _extract_contacts(
self, pp: Polypeptide, na_list: list[NucleicAcid]
) -> list[tuple[Atom, Atom]]:
"""
Extract interface contacts (PRIVATE).
Parameters
----------
Raises
------
PDBConstructionException
In case of no protein-nucleic acids interactions.
Returns
-------
list[tuple[Atom, Atom]]
List of pairs of atoms, first one is from the nucleic acids,
second one is from the protein.
"""
# get all the atoms from the proteins
pp_atoms = []
for res in pp:
pp_atoms.extend(res.get_atoms())
pp_atoms = list(set(pp_atoms))
# get chain id for protein
pp_chain_id = pp[0].parent.id
# get all the atoms from the nucleic acids
na_atoms = []
for na in na_list:
for res in na:
na_atoms.extend(res.get_atoms())
na_atoms = list(set(na_atoms))
# get chain ids for nucleic acids
na_chain_ids = [na.get_chain_id() for na in na_list]
# build list of all atoms, both nucleic acid and protein
all_atoms = na_atoms + pp_atoms
# filter out hydrogens
all_atoms = [atom for atom in all_atoms if atom.element != "H"]
# look for contacts between nucleic acid and protein
# within a given radius
ns = NeighborSearch(all_atoms)
all_contacts = ns.search_all(self.search_radius)
if not all_contacts:
return []
# TODO using chain ids ignores the case where a chain has both protein
# and nucleic
# examples: 7B5Y
# filter possible contacts, meaning the atom couples
# with at least one nucleic acid atom
temp = [
(atom1, atom2)
for atom1, atom2 in all_contacts
if (
atom1.parent.parent.id in na_chain_ids
or atom2.parent.parent.id in na_chain_ids
)
]
# filter possible contacts, meaning the atom couples
# with one nucleic acid atom and one protein atom
contacts: list[tuple[Atom, Atom]] = [
(atom1, atom2)
for atom1, atom2 in temp
if (
atom1.parent.parent.id in na_chain_ids
and atom2.parent.parent.id == pp_chain_id
)
] + [
(atom2, atom1)
for atom1, atom2 in temp
if (
atom1.parent.parent.id == pp_chain_id
and atom2.parent.parent.id in na_chain_ids
)
]
# temp รจ utile per greedyness, prima prendi il DNA, che ha meno
# atomi, escludendo quelli intra-proteina
# poi escludi anche quelli intra-DNA
# inoltre ci assicuriamo che la col 0 contenga gli atomi di DNA e
# che col 1 contenga gli atomi di proteine
# Ce ne assicuriamo invertendo atom1 e atom2 nella seconda lista
return contacts
def _extract_binding_domain(
self,
face: Interface,
pp: Polypeptide,
upstream_pad: int = 0,
downstream_pad: int = 0,
) -> Polypeptide:
"""
Get nucleic acid binding domain from the protein.
The output is the binding "gapped" subsequence of the full protein
found in the structure.
This method allows for "gaps" of unbound aminoacids inside the
binding domain, only the aminoacids at the ends are trimmed according
to being bound to nucleic acids (NAs) or not.
A visual example of "gaps":
```
Input full protein: MQMLLNHKPTKFNGAIDERFHWKVIQRISGSEG
NA-bound: **** **
Output binding domain: FNGAIDER
```
This method is only an inference of the NA-binding domain:
while the output will likely align with the annotated true domain,
it'll likely not infer the whole domain. This is because a domain is
defined by folding properties, while this method is much more naive.
This is why I implemented some "padding" on both ends of the binding
domain, it allows to be more lenient of the extent of the binding
domain.
:param face: Interface already initialized with contacts.
:type face: Interface
:param pp: Polypeptide that binds nucleic acids.
:type pp: Polypeptide
:param upstream_pad: Number of non-binding residues,
upstream of the first binding residue, to take inside the binding
domain. Allows some leniency on what is considered a binding
domain.
:type upstream_pad: int
:param downstream_pad: Number of non-binding residues,
downstream of the last binding residue, to take inside the binding
domain. Allows some leniency on what is considered a binding
domain.
:type downstream_pad: int
Returns
-------
Polypeptide
Nucleic acid binding domain.
"""
bound_aminoacids = face.get_aminoacids()
# check if given chain id is actually a protein
# if len(pp_list) > 1:
# warnings.warn(
# f"Warning: {len(pp_list)} peptides are found. \
# Meaning the protein is uncontiguous.\
# Uncontiguous polypeptides will be joined \
# and treated as a whole protein."
# )
# protein = Polypeptide()
# for pp in pp_list:
# protein.extend(pp)
# # check if the protein has sorted residues
# protein_nums: list[int] = [res.id[1] for res in pp]
# if protein_nums != sorted(protein_nums):
# raise PDBConstructionException(
# f"Residues are not in correct order in proteic \
# chain id: {self.protein_chain_id}"
# )
bound_protein_nums: list[int] = [res.id[1] for res in bound_aminoacids]
# find the start and the end of the domain
for res in pp:
if res.id[1] == min(bound_protein_nums):
start: int = pp.index(res)
elif res.id[1] == max(bound_protein_nums):
end: int = pp.index(res)
# add padding
# TODO add padding as init parameter, also for nucleic acids
start = start - upstream_pad # type: ignore
end = end + downstream_pad # type: ignore
if start < 0:
start = 0
if end > len(pp) - 1:
end = len(pp) - 1
# finally get binding domain
binding_domain = pp[start : end + 1] # type: ignore
binding_domain = Polypeptide(binding_domain)
return binding_domain
def _extract_bound_nucleic_acids(
self, face: Interface, na_list: list[NucleicAcid]
) -> list[NucleicAcid]:
"""
Get all nucleic acids bound by the protein.
The output nucleic acids (NAs) are subsequences
of the full NAs found in the structure,
since proteins might not bind the whole NA.
This method allows for "gaps" of unbound nucleotides inside the
NA, only the nucleotides at the ends are trimmed according
to being protein-bound or not.
A visual example of "gaps":
```
Input full NA: GATATACAAGCCA
Protein-bound: **** **
Output protein-bound NA: TATACAAG
```
:param face: Interface already initialized with contacts.
:type face: Interface
:param na_list: List of nucleic acids found in the entity given
to the builder.
:type na_list: list[NucleicAcid]
Returns
-------
list[NucleicAcid] | None
List of nucleic acids bound by the protein.
"""
bound_nucleotides = face.get_nucleotides()
bound_na_list = []
for na in na_list:
# find the start and the end of the bound NA
start, end = None, None
for res in na:
if res in bound_nucleotides and start is None:
start: int | None = na.index(res)
if res in bound_nucleotides:
end: int | None = na.index(res)
# skip non bound NAs
if not start or not end:
continue
# address padding
# start = start - upstream_pad # type: ignore
# end = end + downstream_pad # type: ignore
# if start < 0: # type: ignore
# start = 0
# if end > len(protein) - 1: # type: ignore
# end = len(protein) - 1
# finally get a substring of the nucleic acid
sub_na = na[start : end + 1]
sub_na = NucleicAcid(sub_na)
bound_na_list.append(sub_na)
return bound_na_list
def _extract_bound_double_strands(
self, face: Interface, dsna_list: list[DoubleStrandNucleicAcid]
) -> list[DoubleStrandNucleicAcid]:
"""
Get all double-strand nucleic acids bound by the protein.
The output double stranded nucleic acids (DSNAs) are subsequences
of the full DSNAs found in the structure,
since proteins usually do not bind the whole DSNA.
This method allows for "gaps" of unbound base-pairs inside the
DSNA, only the base pairs at the ends are trimmed according
to being protein-bound or not.
A visual example of "gaps":
```
Input full DSNA: GATATACAAGCCA
|||||||||||||
TGGCTTGTATATC
Protein-bound: **** **
Output protein-bound DSNA: TATACAAG
||||||||
CTTGTATA
```
:param face: Interface already initialized with contacts.
:type face: Interface
:param dsna_list: List of double strand nucleic acids found in the
entity given to the builder.
:type dsna_list: list[DoubleStrandNucleicAcid]
Returns
-------
list[DoubleStrandNucleicAcid]
List of double-strand nucleic acids bound by the protein.
"""
bound_nucleotides = face.get_nucleotides()
bound_dsna_list = []
for dsna in dsna_list:
bound_dsna = copy.copy(dsna)
while (
len(bound_dsna) > 0
and bound_dsna[0].i_res not in bound_nucleotides
and bound_dsna[0].j_res not in bound_nucleotides
):
# if the FIRST base pair isn't bound by protein
# then discard it and check the next FIRST base pair
bound_dsna.pop(0)
while (
len(bound_dsna) > 0
and bound_dsna[-1].i_res not in bound_nucleotides
and bound_dsna[-1].j_res not in bound_nucleotides
):
# if the LAST base pair isn't bound by protein
# then discard it and check the next LAST base pair
bound_dsna.pop(-1)
if len(bound_dsna) > 0:
# in this case, there is an actual bound DSNA
bound_dsna_list.append(bound_dsna)
unbound_bps = []
for bp in bound_dsna:
if (
bp.i_res not in bound_nucleotides
and bp.j_res not in bound_nucleotides
):
unbound_bps.append(bp)
if unbound_bps:
warnings.warn(
f"Warning: There are {len(unbound_bps)} unbound \
base-pairs inside {bound_dsna} - {unbound_bps}"
)
return bound_dsna_list
[docs]
def build_interfaces(
self,
entity: Structure | Model | Chain,
pp_builder: PPBuilder = PPBuilder(),
standard_aminoacids: bool = True,
na_builder: NABuilder = NABuilder(),
dsna_builder: DSNABuilder = DSNABuilder(),
standard_nucleotides: bool = True,
pairing_rules: BasePairRules = WatsonCrickBasePairRules(),
) -> list[Interface]:
"""
Extract all Protein-Nucleic acid interfaces found in a PDB entity.
Parameters
----------
entity : L{Structure}, L{Model} or L{Chain}
Protein-nucleic acid interfaces are searched for in this object.
L{Structure} is the suggested input.
pp_builder : PPBuilder, optional
Polypeptide builder class from Biopython. Default is ``PPBuilder``
with default parameters.
standard_aminoacids: bool, optional
Use only standard aminoacids. This is the `aa_only` parameter
in the ``PPBuilder.build_peptides()`` method. Default is True.
na_builder : NABuilder, optional
Polypeptide builder class from PDBNucleicAcids.
Default is ``NABuilder`` with default parameters.
dsna_builder : DSNABuilder, optional
Polypeptide builder class from PDBNucleicAcids.
Default is ``DSNABuilder`` with default parameters.
standard_nucleotides: bool, optional
Use only standard nucleotides. This parameter is used
in the ``NABuilder.build_nucleic_acids()`` method and
in the ``DSNABuilder.build_double_strands()`` method.
Default is True.
pairing_rules : optional
Rules for proper base pairing class instance from PDBNucleicAcids.
This parameter is used in the
``DSNABuilder.build_double_strands()`` method.
Default is ``WatsonCrickBasePairRules()`` with default parameters.
Raises
------
PDBConstructionException: In case there is no protein
in the input entity.
PDBConstructionException: In case there is no nucleic acid
in the input entity.
Returns
-------
list[Interface]
List of all Protein-Nucleic acid interfaces found in a PDB entity.
"""
# build nucleic acids
na_list = na_builder.build_nucleic_acids(
entity=entity, standard_nucleotides=standard_nucleotides
)
# build double stranded nucleic acids
dsna_list = dsna_builder.build_double_strands(
entity=entity,
standard_nucleotides=standard_nucleotides,
pairing_rules=pairing_rules,
)
# check if there are nucleic acids
if not na_list:
raise PDBConstructionException(
f"No nucleic acids found in the input entity {entity}"
)
# get all the atoms from the nucleic acids
na_atoms = []
for na in na_list:
na_atoms.extend(na.get_atoms())
na_atoms = list(set(na_atoms))
# build the proteins
pp_list = pp_builder.build_peptides(
entity=entity, aa_only=standard_aminoacids
)
# check if there are proteins
if not pp_list:
raise PDBConstructionException(
f"No polypeptides found in the input entity {entity}"
)
face_list = []
for pp in pp_list:
# extract contacts
contacts = self._extract_contacts(pp=pp, na_list=na_list)
# skip empty interfaces
if not contacts:
continue
# build interface
face = Interface(
contacts=contacts,
search_radius=self.search_radius,
)
# fill other properties
face._binding_domain = self._extract_binding_domain(
face=face, pp=pp
)
face._bound_na_list = self._extract_bound_nucleic_acids(
face=face, na_list=na_list
)
face._bound_dsna_list = self._extract_bound_double_strands(
face=face, dsna_list=dsna_list
)
# save interface for output
face_list.append(face)
if not face_list:
warnings.warn(
f"No Protein-Nucleic acids interfaces found \
from polypeptides {pp_list} and nucleic acids {na_list}, \
with search radius {self.search_radius}, \
in entity {entity}"
)
return face_list