Model.write_data() — write derivative model data

write_data(output, file=None, surftyp=1, neighbor_cutoff=6.0, accessibility_type=8, probe_radius=1.4, psa_integration_step=0.1, dnr_accpt_lib='$LIB/donor_acceptor.lib', edat=None)

Requirements:

topology file

This command calculates and writes out the selected type(s) of data (from the list below) about the model. If file is specified and is non-empty, each property is written out to a file using file as the root of the file name. The last such property is also assigned to the B_iso field of each atom of the model (in the case of residue properties, each atom in the residue gets the value of the property). This can be accessed as Atom.biso, or written out with Model.write() to a PDB file, where it appears as the temperature factor.

The data to be calculated are specified by concatenating the corresponding keywords in the output variable:

• 'ALL': All types of data are written.

• 'PSA': The atomic and residue solvent accessibilities are written to .sol and .psa files, respectively. The algorithm for the solvent contact areas is described in [Richmond & Richards, 1978], and can be tuned by changing the probe_radius and psa_integration_step parameters. The normalization for the fractional areas is carried out as described in [Hubbard & Blundell, 1987], with the normalization factors courtesy of Simon Hubbard (personal communication)^6.3. The single reference is [Šali & Overington, 1994]. Accessibilities are calculated with scaled radii from the $MODELS_LIB library, as specified by Topology.submodel. The radii are scaled by EnergyData.radii_factor, which should usually be set to 1. If output also contains ATOMIC_SOL, atomic accessibilities in Ų are assigned to B_iso, otherwise residue accessibility of type accessibility_type (from 1 to 10, for the columns in the .psa file) is assigned. (In either case, atomic accessibilities in Atom.accessibility are updated.) If surftyp is 1, contact accessibility is calculated; if 2, surface accessibility is returned. Atoms with undefined coordinates are assigned zero accessibility.

• 'NGH': Residue neighbors of each residue are listed to a .ngh file. The MODELLER definition of a residue-residue contact used in restraints derivation is applied [Šali & Blundell, 1993]: Any pair of residues that has any pair of atoms within neighbor_cutoff Å of each other are in contact. The number of neighbors for each residue is assigned to B_iso.

• 'DIH': All the dihedral angle types defined in the $RESDIH_LIB library (mainchain, sidechain, and the virtual dihedral between four successive C_α atoms, starting with the previous residue) are written to a .dih file. One column from this file, as selected by accessibility_type, is also assigned to B_iso.

• 'SSM': Secondary structure assignments are written to a .ssm file, and also to B_iso (0 for unknown, 1 for strand, 2 for helix, and -2 for kink). The algorithm for secondary structure assignment depends on the C_α positions only and is based on the distance matrix idea described in [Richards & Kundrot, 1988]. For each secondary structure type, a ‘library’ C_α distance matrix was calculated by averaging distance matrices for several secondary structure segments from a few high resolution protein structures. Program DSSP was used to assign these secondary structure segments [Kabsch & Sander, 1983]. Outlier distances were omitted from the averaging. Currently, there are only two matrices: one for the α-helix and one for the β-strand. The algorithm for secondary structure assignment is as follows:
  1. For each secondary structure type (begin with a helix, which can thus overwrite parts of strand if they overlap):
     • Define the degree of the current secondary structure fit for each C_α atom by DRMS deviation ($P_1$) and maximal distance difference ($P_2$) obtained by comparing the library distance matrix with the distance matrix for a segment starting at the given C_α position;

     • Assign the current secondary structure type to all C_α's in all segments whose DRMS deviation and maximal distance difference are less than some cutoffs ($P_1<cut_1$, $P_1<cut_2$) and are not already assigned to ‘earlier’ secondary structure types;

  2. Split kinked contiguous segments of the same type into separate segments:

     Kinking residues have both DRMS and maximal distance difference beyond their respective cutoffs ($P_1>cut_3$, $P_2>cut_4$). The actual single kink residue separating the two new segments of the same type is the central kinking residue. Note: we are assuming that there are no multiple kinks within one contiguous segment of residues of the same secondary structure type.

  3. If the current secondary structure type is β-strand: Eliminate those runs of strand residues that are not close enough to other strand residues separated by at least two other residues: $P_3$ is minimal distance to a non-neighboring residue of the strand type ( $P_3 < cut_3$). Currently, only one pass of this elimination is done, but could be repeated until self-consistency.

  4. Eliminate those segments that are shorter than the cutoff ($cut_6$) length (e.g., 5 or 6).

  5. Remove the isolated kinking residues (those that occur on their own or begin or end a segment).

• 'CRV': Local mainchain curvatures, multiplied by 0.1, are assigned to the B_iso field. Local mainchain curvature at residue $i$ is defined as the angle (expressed in degrees, between 0 and 180) between the least-squares lines through C_α atoms $i-3$ to $i$ and $i$ to $i+3$.

• 'HBONDS': Hydrogen bonds between amino acid residues are written to a .hbnds file. A list of donors and acceptors in the 20 naturally occurring amino acids, specified by dnr_accpt_lib, is utilized in the computation of H-bonds. Hydrogen bonds are reported if the donor - acceptor distance is between 2.5 and 3.5Å and the donor-acceptor-acceptor antecedent angle is larger than 120 degrees. The model's B_iso field is not changed by this property.

Example: examples/commands/write_data.py

# Example for: Model.write_data()

# This will calculate solvent accessibility, dihedral angles,
# residue-residue neighbors, secondary structure, and local mainchain
# curvature for a structure in the PDB file.

from modeller import *

log.verbose()

# Get topology library for radii and the model without waters and HETATMs:
env = Environ()
env.io.atom_files_directory = ['../atom_files']
env.io.hetatm = False
env.io.water = False

env.libs.topology.read(file='$(LIB)/top_heav.lib')
mdl = Model(env, file='1fas')

# Calculate residue solvent accessibilities, dihedral angles,
# residue neighbors, and secondary structure, and write to files:
myedat = EnergyData()
myedat.radii_factor = 1.0 # The default is 0.82 (for soft-sphere restraints)
mdl.write_data(file='1fas', edat=myedat, output='PSA DIH NGH SSM')

# Calculate local mainchain curvature
mdl.write_data(output='CRV')

# Use the calculated curvature data (in Biso)
print("The following residues have local mainchain curvature")
print("greater than 90 degrees:")
print([r for r in mdl.residues if r.atoms[0].biso * 10 > 90.0])