model.write_data() -- write derivative MODEL data

edat = <energy_data>		objective function parameters
file = <str:1>	'default'	root of output filename(s)
accessibility_type = <int:1>	8	type of solvent accessibility: 1-10
surftyp = <int:1>	1	Surface Type for accessibility calculations 1= contact; 2=surface
psa_integration_step = <float:1>	0.1	integration step for WRITE_DATA
probe_radius = <float:1>	1.4	probe_radius for WRITE_DATA
neighbor_cutoff = <float:1>	6.0	for defining atom-atom contacts in WRITE_DATA
output = <str:1>	'LONG'	what to calculate and write out: 'ALL' | 'PSA' | 'ATOMIC_SOL' | 'NGH' | 'DIH' | 'SSM' | 'CRV' | 'CAV' | 'CROSS-SECTIONS'

Requirements:

topology file

Description:

This command writes the selected types of data about the MODEL to a corresponding file and to the `fourth' column of the model. The root of the output filenames is given by the file variable. In addition to the output files, the $B_{\mbox{iso}}$ field of the model (`fourth column' in the PDB file) will be assigned the last selected property from the following list: atomic or residue accessibility, dihedral type accessibility_type (from 1 to 9 for $\alpha$, $\Phi$, $\Psi$, $\omega$, $\chi_1$, $\chi_2$, $\chi_3$, $\chi_4$, and $\chi_5$; where $\alpha$ is the virtual dihedral angle between four successive ${C}_\alpha$ atoms), number of residue neighbors, the secondary structure type, and the local mainchain curvature. For accessibility, when output contains ATOMIC_SOL, atomic accessibilities in $\mbox{\AA}^2$ are assigned to $B_{\mbox{iso}}$, otherwise residue accessibility of type accessibility_type (from 1 to 10, for the columns in the .psa file) is assigned. If surftyp is 1, contact accessibility is calculated; if 2, surface accessibility is returned.

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

• 'ALL': All types of data are written to the corresponding files.

• 'PSA': The atomic and residue solvent accessibilities are written to the .sol and .psa files, respectively. The algorithm for the solvent contact areas is described in [Richmond & Richards, 1978]. 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). 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 energy_data.radii_factor, which should usually be set to 1.

• 'CAV': The protein and internal cavity volumes are written out. The calculation on a grid is used. The grid unit is specified by grid_unit in angstroms (say 1.4). The radii are scaled by energy_data.radii_factor, which should usually be set to 1. The cross-sections are written to file file.cav when output contains CROSS-SECTIONS. The number_of_steps is the number of small shifts along x, y, and z that are used in the averaging of the protein and cavity volumes with respect to small changes in the relative position of the protein and the grid; the total number of calculations is therefore equal to the third power of number_of_steps. If orient is True, the structure is oriented before the volume calculation such that the moment of inertia are parallel to the x, y, and z coordinate axes (this orientation minimizes the size of the grid). However, the coordinates of the MODEL are not changed upon exit from this routine (you need to use model.orient() to change the orientation of the MODEL).

• '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 6Å of each other are in contact.

• 'DIH': All the dihedral angle types defined in the $RESDIH_LIB library (virtual ${C}_\alpha$, mainchain, and sidechain dihedral angles) are written to a .dih file.

• 'SSM': Secondary structure assignments are written to a .ssm file. The algorithm for secondary structure assignment depends on the ${C}_\alpha$ positions only and is based on the distance matrix idea described in [Richards & Kundrot, 1988]. For each secondary structure type, a `library' ${C}_\alpha$ 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 $\alpha$-helix (secondary structure type 2) and one for the $\beta$-strand (type 1). 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}_\alpha$ 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}_\alpha$ position;

     • Assign the current secondary structure type to all ${C}_\alpha$'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. The kink residue type is $-2$.

  3. If the current secondary structure type is $\beta$-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 are written to a .crv file. Local mainchain curvature at residue $i$ is defined as the angle between the least-squares lines through ${C}_\alpha$ atoms $i-3$ to $i$ and $i$ to $i+3$.

Example: examples/commands/write_data.py

# Example for: model.write_data()

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

log.verbose()

# Get topology library for radii and the model without waters and HETATMs:
env = environ()
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, and 
# residue-residue neighbors:
myedat = energy_data()
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 CRV')