biopython/Bio/PDB/DSSP.py

# Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
#
# This file is part of the Biopython distribution and governed by your
# choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
# Please see the LICENSE file that should have been included as part of this
# package.

r"""Use the DSSP program to calculate secondary structure and accessibility.

You need to have a working version of DSSP (and a license, free for academic
use) in order to use this. For DSSP, see https://swift.cmbi.umcn.nl/gv/dssp/.

The following Accessible surface area (ASA) values can be used, defaulting
to the Sander and Rost values:

    Ahmad
        Ahmad et al. 2003 https://doi.org/10.1002/prot.10328
    Miller
        Miller et al. 1987 https://doi.org/10.1016/0022-2836(87)90038-6
    Sander
        Sander and Rost 1994 https://doi.org/10.1002/prot.340200303
    Wilke
        Tien et al. 2013 https://doi.org/10.1371/journal.pone.0080635

The DSSP codes for secondary structure used here are:

    =====     ====
    Code      Structure
    =====     ====
     H        Alpha helix (4-12)
     B        Isolated beta-bridge residue
     E        Strand
     G        3-10 helix
     I        Pi helix
     T        Turn
     S        Bend
     \-       None
    =====     ====

Usage
-----
The DSSP class can be used to run DSSP on a PDB or mmCIF file, and provides a
handle to the DSSP secondary structure and accessibility.

**Note** that DSSP can only handle one model, and will only run
calculations on the first model in the provided PDB file.

Examples
--------
Typical use::

    from Bio.PDB import PDBParser
    from Bio.PDB.DSSP import DSSP
    p = PDBParser()
    structure = p.get_structure("1MOT", "/local-pdb/1mot.pdb")
    model = structure[0]
    dssp = DSSP(model, "/local-pdb/1mot.pdb")

Note that the recent DSSP executable from the DSSP-2 package was
renamed from ``dssp`` to ``mkdssp``. If using a recent DSSP release,
you may need to provide the name of your DSSP executable::

    dssp = DSSP(model, '/local-pdb/1mot.pdb', dssp='mkdssp')

DSSP data is accessed by a tuple - (chain id, residue id)::

    a_key = list(dssp.keys())[2]
    dssp[a_key]

The dssp data returned for a single residue is a tuple in the form:

    ============ ===
    Tuple Index  Value
    ============ ===
    0            DSSP index
    1            Amino acid
    2            Secondary structure
    3            Relative ASA
    4            Phi
    5            Psi
    6            NH-->O_1_relidx
    7            NH-->O_1_energy
    8            O-->NH_1_relidx
    9            O-->NH_1_energy
    10           NH-->O_2_relidx
    11           NH-->O_2_energy
    12           O-->NH_2_relidx
    13           O-->NH_2_energy
    ============ ===

"""

import os
import re
import subprocess
import warnings
from io import StringIO

from Bio.Data.PDBData import protein_letters_3to1
from Bio.Data.PDBData import residue_sasa_scales
from Bio.PDB.AbstractPropertyMap import AbstractResiduePropertyMap
from Bio.PDB.MMCIF2Dict import MMCIF2Dict
from Bio.PDB.PDBExceptions import PDBException
from Bio.PDB.PDBParser import PDBParser

# Match C in DSSP
_dssp_cys = re.compile("[a-z]")

# Maximal ASA of amino acids
# Used for relative accessibility
residue_max_acc = residue_sasa_scales


def version(version_string):
    """Parse semantic version scheme for easy comparison."""
    return tuple(map(int, (version_string.split("."))))


def ss_to_index(ss):
    """Secondary structure symbol to index.

    H=0
    E=1
    C=2
    """
    if ss == "H":
        return 0
    if ss == "E":
        return 1
    if ss == "C":
        return 2
    assert 0


def dssp_dict_from_pdb_file(in_file, DSSP="dssp", dssp_version="3.9.9"):
    """Create a DSSP dictionary from a PDB file.

    Parameters
    ----------
    in_file : string
        pdb file

    DSSP : string
        DSSP executable (argument to subprocess)

    dssp_version : string
        Version of DSSP executable

    Returns
    -------
    (out_dict, keys) : tuple
        a dictionary that maps (chainid, resid) to
        amino acid type, secondary structure code and
        accessibility.

    Examples
    --------
    How dssp_dict_from_pdb_file could be used::

        from Bio.PDB.DSSP import dssp_dict_from_pdb_file
        dssp_tuple = dssp_dict_from_pdb_file("/local-pdb/1fat.pdb")
        dssp_dict = dssp_tuple[0]
        print(dssp_dict['A',(' ', 1, ' ')])

    """
    # Using universal newlines is important on Python 3, this
    # gives text handles rather than bytes handles.
    # Newer version of DSSP executable is named 'mkdssp',
    # and calling 'dssp' will hence not work in some operating systems
    # (Debian distribution of DSSP includes a symlink for 'dssp' argument)
    try:
        if version(dssp_version) < version("4.0.0"):
            DSSP_cmd = [DSSP, in_file]
        else:
            DSSP_cmd = [DSSP, "--output-format=dssp", in_file]
        p = subprocess.Popen(
            DSSP_cmd,
            universal_newlines=True,
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE,
        )
    except FileNotFoundError:
        if DSSP == "mkdssp":
            raise
        if version(dssp_version) < version("4.0.0"):
            DSSP_cmd = ["mkdssp", in_file]
        else:
            DSSP_cmd = ["mkdssp", "--output-format=dssp", in_file]
        p = subprocess.Popen(
            DSSP_cmd,
            universal_newlines=True,
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE,
        )
    out, err = p.communicate()

    # Alert user for errors
    if err.strip():
        warnings.warn(err)
        if not out.strip():
            raise Exception("DSSP failed to produce an output")

    out_dict, keys = _make_dssp_dict(StringIO(out))
    return out_dict, keys


def make_dssp_dict(filename):
    """DSSP dictionary mapping identifiers to properties.

    Return a DSSP dictionary that maps (chainid, resid) to
    aa, ss and accessibility, from a DSSP file.

    Parameters
    ----------
    filename : string
        the DSSP output file

    """
    with open(filename) as handle:
        return _make_dssp_dict(handle)


def _make_dssp_dict(handle):
    """Return a DSSP dictionary, used by mask_dssp_dict (PRIVATE).

    DSSP dictionary maps (chainid, resid) to an amino acid,
    secondary structure symbol, solvent accessibility value, and hydrogen bond
    information (relative dssp indices and hydrogen bond energies) from an open
    DSSP file object.

    Parameters
    ----------
    handle : file
        the open DSSP output file handle

    """
    dssp = {}
    start = 0
    keys = []
    for line in handle:
        sl = line.split()
        if len(sl) < 2:
            continue
        if sl[1] == "RESIDUE":
            # Start parsing from here
            start = 1
            continue
        if not start:
            continue
        if line[9] == " ":
            # Skip -- missing residue
            continue

        dssp_index = int(line[:5])
        resseq = int(line[5:10])
        icode = line[10]
        chainid = line[11]
        aa = line[13]
        ss = line[16]
        if ss == " ":
            ss = "-"
        try:
            NH_O_1_relidx = int(line[38:45])
            NH_O_1_energy = float(line[46:50])
            O_NH_1_relidx = int(line[50:56])
            O_NH_1_energy = float(line[57:61])
            NH_O_2_relidx = int(line[61:67])
            NH_O_2_energy = float(line[68:72])
            O_NH_2_relidx = int(line[72:78])
            O_NH_2_energy = float(line[79:83])

            acc = int(line[34:38])
            phi = float(line[103:109])
            psi = float(line[109:115])
        except ValueError as exc:
            # DSSP output breaks its own format when there are >9999
            # residues, since only 4 digits are allocated to the seq num
            # field.  See 3kic chain T res 321, 1vsy chain T res 6077.
            # Here, look for whitespace to figure out the number of extra
            # digits, and shift parsing the rest of the line by that amount.
            if line[34] != " ":
                shift = line[34:].find(" ")

                NH_O_1_relidx = int(line[38 + shift : 45 + shift])
                NH_O_1_energy = float(line[46 + shift : 50 + shift])
                O_NH_1_relidx = int(line[50 + shift : 56 + shift])
                O_NH_1_energy = float(line[57 + shift : 61 + shift])
                NH_O_2_relidx = int(line[61 + shift : 67 + shift])
                NH_O_2_energy = float(line[68 + shift : 72 + shift])
                O_NH_2_relidx = int(line[72 + shift : 78 + shift])
                O_NH_2_energy = float(line[79 + shift : 83 + shift])

                acc = int(line[34 + shift : 38 + shift])
                phi = float(line[103 + shift : 109 + shift])
                psi = float(line[109 + shift : 115 + shift])
            else:
                raise ValueError(exc) from None
        res_id = (" ", resseq, icode)
        dssp[(chainid, res_id)] = (
            aa,
            ss,
            acc,
            phi,
            psi,
            dssp_index,
            NH_O_1_relidx,
            NH_O_1_energy,
            O_NH_1_relidx,
            O_NH_1_energy,
            NH_O_2_relidx,
            NH_O_2_energy,
            O_NH_2_relidx,
            O_NH_2_energy,
        )
        keys.append((chainid, res_id))
    return dssp, keys


class DSSP(AbstractResiduePropertyMap):
    """Run DSSP and parse secondary structure and accessibility.

    Run DSSP on a PDB/mmCIF file, and provide a handle to the
    DSSP secondary structure and accessibility.

    **Note** that DSSP can only handle one model.

    Examples
    --------
    How DSSP could be used::

        from Bio.PDB import PDBParser
        from Bio.PDB.DSSP import DSSP
        p = PDBParser()
        structure = p.get_structure("1MOT", "/local-pdb/1mot.pdb")
        model = structure[0]
        dssp = DSSP(model, "/local-pdb/1mot.pdb")
        # DSSP data is accessed by a tuple (chain_id, res_id)
        a_key = list(dssp.keys())[2]
        # (dssp index, amino acid, secondary structure, relative ASA, phi, psi,
        # NH_O_1_relidx, NH_O_1_energy, O_NH_1_relidx, O_NH_1_energy,
        # NH_O_2_relidx, NH_O_2_energy, O_NH_2_relidx, O_NH_2_energy)
        dssp[a_key]

    """

    def __init__(self, model, in_file, dssp="dssp", acc_array="Sander", file_type=""):
        """Create a DSSP object.

        Parameters
        ----------
        model : Model
            The first model of the structure
        in_file : string
            Either a PDB file or a DSSP file.
        dssp : string
            The dssp executable (ie. the argument to subprocess)
        acc_array : string
            Accessible surface area (ASA) from either Miller et al. (1987),
            Sander & Rost (1994), Wilke: Tien et al. 2013, or Ahmad et al.
            (2003) as string Sander/Wilke/Miller/Ahmad. Defaults to Sander.
        file_type: string
            File type switch: either PDB, MMCIF or DSSP. Inferred from the
            file extension by default.

        """
        self.residue_max_acc = residue_max_acc[acc_array]

        # create DSSP dictionary
        if file_type == "":
            file_type = os.path.splitext(in_file)[1][1:]
        file_type = file_type.upper()
        if file_type == "CIF":
            file_type = "MMCIF"
        assert file_type in [
            "PDB",
            "MMCIF",
            "DSSP",
        ], "File type must be PDB, mmCIF or DSSP"
        # If the input file is a PDB or mmCIF file run DSSP and parse output:
        if file_type == "PDB" or file_type == "MMCIF":
            # Newer versions of DSSP program call the binary 'mkdssp', so
            # calling 'dssp' will not work in some operating systems
            # (Debian distribution of DSSP includes a symlink for 'dssp' argument)
            try:
                version_string = subprocess.check_output([dssp, "--version"], text=True)
                dssp_version = re.search(r"\s*([\d.]+)", version_string).group(1)
                dssp_dict, dssp_keys = dssp_dict_from_pdb_file(
                    in_file, dssp, dssp_version
                )
            except FileNotFoundError:
                if dssp == "dssp":
                    dssp = "mkdssp"
                elif dssp == "mkdssp":
                    dssp = "dssp"
                else:
                    raise
                version_string = subprocess.check_output([dssp, "--version"], text=True)
                dssp_version = re.search(r"\s*([\d.]+)", version_string).group(1)
                dssp_dict, dssp_keys = dssp_dict_from_pdb_file(
                    in_file, dssp, dssp_version
                )
        # If the input file is a DSSP file just parse it directly:
        elif file_type == "DSSP":
            dssp_dict, dssp_keys = make_dssp_dict(in_file)

        dssp_map = {}
        dssp_list = []

        def resid2code(res_id):
            """Serialize a residue's resseq and icode for easy comparison."""
            return f"{res_id[1]}{res_id[2]}"

        # DSSP outputs label_asym_id from the mmCIF file as the chain ID
        # But MMCIFParser reads in the auth_asym_id
        # Here we create a dictionary to map label_asym_id to auth_asym_id
        # using the mmCIF file
        if file_type == "MMCIF" and version(dssp_version) < version("4.0.0"):
            mmcif_dict = MMCIF2Dict(in_file)
            mmcif_chain_dict = {}
            for i, c in enumerate(mmcif_dict["_atom_site.label_asym_id"]):
                if c not in mmcif_chain_dict:
                    mmcif_chain_dict[c] = mmcif_dict["_atom_site.auth_asym_id"][i]
            dssp_mapped_keys = []

        # Now create a dictionary that maps Residue objects to
        # secondary structure and accessibility, and a list of
        # (residue, (secondary structure, accessibility)) tuples
        for key in dssp_keys:
            chain_id, res_id = key
            if file_type == "MMCIF" and version(dssp_version) < version("4.0.0"):
                chain_id = mmcif_chain_dict[chain_id]
                dssp_mapped_keys.append((chain_id, res_id))
            chain = model[chain_id]
            try:
                res = chain[res_id]
            except KeyError:
                # In DSSP, HET field is not considered in residue identifier.
                # Thus HETATM records may cause unnecessary exceptions.
                # (See 3jui chain A res 593.)
                # Try the lookup again with all HETATM other than water
                res_seq_icode = resid2code(res_id)
                for r in chain:
                    if r.id[0] not in (" ", "W"):
                        # Compare resseq + icode
                        if resid2code(r.id) == res_seq_icode:
                            # Found a matching residue
                            res = r
                            break
                else:
                    raise KeyError(res_id) from None

            # For disordered residues of point mutations, Biopython uses the
            # last one as default, But DSSP takes the first one (alternative
            # location is blank, A or 1). See 1h9h chain E resi 22.
            # Here we select the res in which all atoms have altloc blank, A or
            # 1. If no such residues are found, simply use the first one appears
            # (as DSSP does).
            if res.is_disordered() == 2:
                for rk in res.disordered_get_id_list():
                    # All atoms in the disordered residue should have the same
                    # altloc, so it suffices to check the altloc of the first
                    # atom.
                    altloc = res.child_dict[rk].get_list()[0].get_altloc()
                    if altloc in tuple("A1 "):
                        res.disordered_select(rk)
                        break
                else:
                    # Simply select the first one
                    res.disordered_select(res.disordered_get_id_list()[0])

            # Sometimes point mutations are put into HETATM and ATOM with altloc
            # 'A' and 'B'.
            # See 3piu chain A residue 273:
            #   <Residue LLP het=H_LLP resseq=273 icode= >
            #   <Residue LYS het=  resseq=273 icode= >
            # DSSP uses the HETATM LLP as it has altloc 'A'
            # We check the altloc code here.
            elif res.is_disordered() == 1:
                # Check altloc of all atoms in the DisorderedResidue. If it
                # contains blank, A or 1, then use it.  Otherwise, look for HET
                # residues of the same seq+icode.  If not such HET residues are
                # found, just accept the current one.
                altlocs = {a.get_altloc() for a in res.get_unpacked_list()}
                if altlocs.isdisjoint("A1 "):
                    # Try again with all HETATM other than water
                    res_seq_icode = resid2code(res_id)
                    for r in chain:
                        if r.id[0] not in (" ", "W"):
                            if resid2code(r.id) == res_seq_icode and r.get_list()[
                                0
                            ].get_altloc() in tuple("A1 "):
                                res = r
                                break

            (
                aa,
                ss,
                acc,
                phi,
                psi,
                dssp_index,
                NH_O_1_relidx,
                NH_O_1_energy,
                O_NH_1_relidx,
                O_NH_1_energy,
                NH_O_2_relidx,
                NH_O_2_energy,
                O_NH_2_relidx,
                O_NH_2_energy,
            ) = dssp_dict[key]

            res.xtra["SS_DSSP"] = ss
            res.xtra["EXP_DSSP_ASA"] = acc
            res.xtra["PHI_DSSP"] = phi
            res.xtra["PSI_DSSP"] = psi
            res.xtra["DSSP_INDEX"] = dssp_index
            res.xtra["NH_O_1_RELIDX_DSSP"] = NH_O_1_relidx
            res.xtra["NH_O_1_ENERGY_DSSP"] = NH_O_1_energy
            res.xtra["O_NH_1_RELIDX_DSSP"] = O_NH_1_relidx
            res.xtra["O_NH_1_ENERGY_DSSP"] = O_NH_1_energy
            res.xtra["NH_O_2_RELIDX_DSSP"] = NH_O_2_relidx
            res.xtra["NH_O_2_ENERGY_DSSP"] = NH_O_2_energy
            res.xtra["O_NH_2_RELIDX_DSSP"] = O_NH_2_relidx
            res.xtra["O_NH_2_ENERGY_DSSP"] = O_NH_2_energy

            # Relative accessibility
            resname = res.get_resname()
            try:
                rel_acc = acc / self.residue_max_acc[resname]
            except KeyError:
                # Invalid value for resname
                rel_acc = "NA"
            else:
                if rel_acc > 1.0:
                    rel_acc = 1.0
            res.xtra["EXP_DSSP_RASA"] = rel_acc
            # Verify if AA in DSSP == AA in Structure
            # Something went wrong if this is not true!
            # NB: DSSP uses X often
            resname = protein_letters_3to1.get(resname, "X")
            if resname == "C":
                # DSSP renames C in C-bridges to a,b,c,d,...
                # - we rename it back to 'C'
                if _dssp_cys.match(aa):
                    aa = "C"
            # Take care of HETATM again
            if (resname != aa) and (res.id[0] == " " or aa != "X"):
                raise PDBException(f"Structure/DSSP mismatch at {res}")

            dssp_vals = (
                dssp_index,
                aa,
                ss,
                rel_acc,
                phi,
                psi,
                NH_O_1_relidx,
                NH_O_1_energy,
                O_NH_1_relidx,
                O_NH_1_energy,
                NH_O_2_relidx,
                NH_O_2_energy,
                O_NH_2_relidx,
                O_NH_2_energy,
            )

            dssp_map[(chain_id, res_id)] = dssp_vals
            dssp_list.append(dssp_vals)

        if file_type == "MMCIF" and version(dssp_version) < version("4.0.0"):
            dssp_keys = dssp_mapped_keys
        AbstractResiduePropertyMap.__init__(self, dssp_map, dssp_keys, dssp_list)