Source code for schrodinger.ui.sequencealignment.structure_utils

from schrodinger import structure
from schrodinger.infra import mm
from schrodinger.structutils.analyze import calculate_sasa
from schrodinger.structutils.analyze import evaluate_asl

sasa_by_res_name = None


def calculate_sasa_dict(ignore_backbone=False, include_calpha=False):
    """
    This code comes from reactive_protein_residues.py script.

    Calculates the surface-accessable surface area (SASA) for each residue's
    side-chain, and stores them in a dictionary, where the res name is a key.

    :param ignore_backbone: boolean
    :param ignore_backbone: When True, backbone N and C atoms will be deleted
        from a fragment before surface area calculation.

    :param include_calpha: boolean
    :param include_calpha: When True C-alpha will be included in surface
        area calculation.
    """

    # Calculate the SASAs only the first time this function is called:
    global sasa_by_res_name
    if sasa_by_res_name is not None:
        return sasa_by_res_name
    mm.mmfrag_initialize(mm.error_handler)
    frag_handle = mm.mmfrag_new("peptide")
    res_name_st_list = []
    first = True
    while True:
        if first:
            which = mm.MM_FIRST
            first = False
        else:
            which = mm.MM_NEXT
        try:
            name, ct_handle = mm.mmfrag_get_fragment(frag_handle, which)
        except mm.MmException as e:
            if e.rc == mm.MMFRAG_DONE:
                break
        copy_ct = mm.mmct_ct_duplicate(ct_handle)
        res_st = structure.Structure(copy_ct, True)
        if ignore_backbone:
            backbone_atoms = evaluate_asl(
                res_st, "atom.ptype \" N\" or atom.ptype \" C\"")
            res_st.deleteAtoms(backbone_atoms)
        res_name_st_list.append((name, res_st))
    mm.mmfrag_delete(frag_handle)
    mm.mmfrag_terminate()
    sasas_by_res_name = {}
    asl = "sidechain"
    if include_calpha:
        asl += " OR atom.ptype \"CA\""
    for name, res_st in res_name_st_list:
        side_chain_atoms = evaluate_asl(res_st, asl)
        res_sasa = calculate_sasa(res_st, side_chain_atoms)
        # BIOLUM-1715 Prevent possible normalization problems
        if res_sasa < 1e-3:
            res_sasa = 1.0
        sasas_by_res_name[name] = res_sasa
        # Ev:126808 Add an entry for HID (same SASA as HIS):
        if name == "HIS":
            sasas_by_res_name["HID"] = res_sasa
    return sasas_by_res_name