ctprotein module

This is where the ctprotein shinanigans occurs. Maybe we reference [BiDB79] yeah

This is the docstring at the top of ctprotein

class camparitraj.ctprotein.CTProtein(traj, atom_offset, residue_offset)

CTProtein objects are initialized with a trajectory subset that contains only the atoms a specific, single protein. This means that a CTProtein object allows operations to performed on a single protein. A single trajectory may have multiple proteins in it. indexing with a protein object assumes that the protein is indexed from 0 to n, n is the number of residues. This is an important idea to emphasize - it means that you (the user) will need to determine the correct residue index for a region of interest being examined. The region will NOT (necessarily) correspond to the residue index in the PDB file used.

To make this easier the function .print_residues() will print the mapping of residue index to residue name and residue number.

To re-iterate:

Residue number is the number of the residue in the PDB file.

Residue index is the index value associated with a residue in a specific protein and will always begin from 0 - note this will include the peptide caps (ACE/NME) if present.

__init__(traj, atom_offset, residue_offset)

Initialize self. See help(type(self)) for accurate signature.

calculateAllCAdistances(residueIndex, stride=1, mode='CA', onlyCterminalResidues=True, correctOffset=True)

Calculate the full set of distances between C-alpha atoms. Note that by default this explicitly works in a way to avoid computing redundancy where we ONLY compute distances between residue i and residues greater than i up to the final residue. This behaviour is defined by the ‘onlyCterminalResidues’ flag.

Distance is returned in Angstroms.

Can be fed a mode ‘COM’ keyword to calcule center of mass distances instead of CA distances.

residueIndex [int] Defines the residue index to select the CA from

stride [int] {1} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

mode [string] {‘CA’} String, must be one of either ‘CA’ or ‘COM’. ‘CA’ = alpha carbon ‘COM’ = center of mass (associated withe the residue)

onlyCterminalResidues [Bool] (True) This variable means that only residues C-terminal of the residueIndex value will be considered. This is useful when performing an ALL vs. ALL matrix as it ensures that only the upper triangle is calculated if we iterate over all residues, but may not be deseriable in other contexts.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again

get_CAindex(residueIndex, correctOffset=True)

Get the CA atom index for the residue defined by residueIndex

Defensivly checks for errors.

residueIndex [int] Defines the residue index to select the CA from

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again

get_DVector(stride=20)

Function to calculate the vector of D values used to calculate the Phi parameter from Lyle et al[1]. The stride defines the spacing between frames which are analyzed. This is just for practical purposes. The Phi calulation computes a D value for each frame vs. frame comparison - for a 2000 frame simulation this would be 4 Million D values if every value was calculated which is a bit much, so the stride let’s you define how many frames you should skip.

Importantly, the DVector calculated here measures dij (see part III A of the paper) as the CA-CA distance and NOT the average inter-atomic distance. This has two effects:

  1. Heterogeneity is now, formally, a measure over backbone heterogeneity and not full protein heterogeneity - this may be desirable (arguably it’s a more interpratable measure of conformational change) but if the interatomic version is required this could be implemented.

  2. It is much more efficient than the original version

For a 2000 frame trajectory of a 80 residue protein with a stride=20 allows the calculation to take about 5 seconds. However, as protein size increases the computational cost of this process grows rapidly.

stride [int] {20} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

[1] Lyle, N., Das, R. K., & Pappu, R. V. (2013). A quantitative measure for protein conformational heterogeneity. The Journal of Chemical Physics, 139(12), 121907.

get_Q(stride=1, protein_average=True, region=None, correctOffset=True, weights=False)

Function which will calculate the fraction of native contacts in each frame of the trajectory, where the ‘native’ state is defined as a specific frame (1st frame by default - note this means the native state frame = 0 as we index from 0!). In earlier versions the ‘native state frame’ was a variable, but this ends up being extremely messy when weights are considered, so assume the native state frame is always frame 0.

Native contacts are defined using the definition from Best, Hummer, and Eaton, “Native contacts determine protein folding mechanisms in atomistic simulations” PNAS (2013) 10.1073/pnas.1311599110. The implementation is according to the code helpfully provided at http://mdtraj.org/latest/examples/native-contact.html

If protein_average=True a single vector is returned with the overall protein average fraction of native contacts associated with each frame. If protein_average is set to false a 4-position tuple is returned, where each of the four positions has the following identity

idx | 0 | The fraction of the time a contact is native for all each native contact (vector of length $N_NATIVE CONTACTS)

1 | The native contact definition (same length as 1) where each element is a pair of atoms which are considered

native

2 | The residue-by-residue dictionary-native contacts dictionary. Keys are residue name-number and each key-associated

value is the fractional native contacts for atoms associated with that residue. To get the residue-specific fraction of native contacts take the mean of the element values

3 | The ordered list of keys from 2 for easy plotting in a residue-residue manner

4 | A nres x nres array showing a 2D contact map defining inter-residue specific Q values

stride [int] {1} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

protein_average [bool] {True} This is the default mode, and means that the return vector is the AVERAGE fraction of native contacts over the protein surface for each frame (e.g. each value refers to a single frame). If set to false the simulation-average value at native-contact resolution is returned, instead returning $NATIVE_CONTACT number of values and an additional list of native contact pairs.

region [list/tuple of length 2] {None} Defines the first and last residue (INCLUSIVE) for a region to be examined

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

get_RMSD(frame1, frame2=-1, stride=1, region=None, backbone=True, correctOffset=True)

Function which will calculate the aligned RMSD between two frames, or between one frame and all frames in the trajectory. This can be done over the entire protein but we can also specificy a local region to perform this analysis over.

Units are Angstroms.

frame1 [int] Defines the frame to be used as a reference

frame2 [int] {-1} Defines the frame to be used as a comparison, OR if left blank or set to -1 means the entire trajectory

stride [int] {1} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

region [list/tuple of length 2] {[]} Defines the first and last residue (INCLUSIVE) for a region to be examined

backbone [bool] {True} Boolean flag for using either the full chain or just backbone. Generally backbone alone is fine so default to true.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again

get_all_SASA(stride=20, probe_radius=0.14, mode='residue')

Returns the Solvent Accessible Surface Area (SASA) for each residue from every stride-th frame. SASA is determined using shrake_rupley algorithm.

SASA is returned in Angstroms squared, BUT PROBE RADIUS is in nanometers!

stride [int] {20} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

probe_radius [float] {0.14} Radius of the solvent probe used in nm. Uses the Golden-Spiral algorithm. 0.14 nm is pretty standard. NOTE - the probe radius must be in nanometers

mode [string] {‘residue’,’atom’,’sidechain’,’backbone’, ‘all’} Defines the mode used to compute the SASA. Must be one of ‘residue’ or ‘atom’. For atom mode, extracted areas are resolved per-atom. For ‘residue’, this is computed instead on the per residue basis.

get_aminoAcidSequence(oneletter=False, numbered=True)

Returns the protein’s amino acid sequence

keyword [type] {default} Description ………………………………….

oneletter [bool] {False} If true returns a single sequence of one letter amino acid codes. If false get a list of 3 letter codes with residue number separated by a ‘-‘ character

numbered [bool] {True} If true the return value is a list of RESNAME-RESID strings, if false return value is a list of RESNAME in the correct order

get_asphericity(R1=False, R2=False, statusMessage=True, correctOffset=True)

Returns the asphericity associated with the region defined by the intervening stretch of residues between R1 and R2. This can be a somewhat slow operation, so a status message is printed for the impatient biophysicisit.

Asphericity is defined in many places - for my personal favourite explanation and definition see Page 65 of Andreas Vitalis’ thesis (Probing the Early Stages of Polyglutamine Aggregation with Computational Methods, 2009, Washington University in St. Louis).

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_clusters(region=None, stride=20, n_clusters=10, backbone=True, correctOffset=True)

Function to determine the structural clusters associated with a trajectory. This can be useful for identifying the most populated clusters. This approach uses Ward’s hiearchical clustering, which means we must define the number of clusters we want a-priori. Clustering is done using RMSD - BUT the approach taken here would be easy to re-implement in another function where you ‘simiarity’ metric was something else.

Returns a 4-place tuple with the following sub-elements:

[0] - cluster_members: A list of length n_clusters where each element corresponds to the number of frames in each of the 1-n_clusters cluster. i.e. if I had defined n_clusters=3 this would be a list of length 3

[1] - cluster_trajectories: A list of n_cluster mdtraj trajectory objects of the conformations in the cluster. This is particularly useful because it means you can perform any arbitrary analaysis on the cluster members

[2] - cluster distance_matrices: A list of n_clusters where each member is a square matrix that defines the structural distance between each of the members of each cluster. in other words, this quantifies how similar (in terms of RMSD, in units Angstroms).

[3] - cluster_centroids A list of n_clusters where each element is the index associated with each cluster trajectory that defines the cluster centroid (i.e. a ‘representative’ conformation). As an example - if we had 3 clusters this might look like [4,1,6], which means the 4th, 1st, and 6th frame from each of the respective mdtraj trajectories in the cluster_trajectories list would correspond to the centroid.

[4] - cluster frames: List of lists, where each sublist contains the frame indices associated with that cluster. Useful if clustering on a single chain and want to use that information over an entire trajectory

region [list/tuple of length 2] {[]} Defines the first and last residue (INCLUSIVE) for a region to be examined

stride [int] {20} Defines the spacing betwen frames to compare with - i.e. take every $stride-th frame. Setting stride=1 would mean every frame is used, which would mean you’re doing an all vs. all comparions, which would be ideal BUT may be slow.

n_clusters [int] {10} Number of clusters to be returned through Ward’s clustering algorithm.

backbone [bool] {True} Flag to determine if backbone atoms or full chain should be used

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_contact_map(distance_thresh=5.0, stride=1, weights=False)

get_contact_map() returns 2-position tuple with the contact map (N x N matrix) and a contact order vector (N x 1) that describe the contacts (heavy atom - heavy atom interactions) made by each of the residues over the simulation.

Each element is normalized such that it is between 0 and 1. i+1 and i+2 elements are excluded, and the minimum distances is the distance between the closest heavy atoms on each residue (backbone or sidechain).

stride [int] {1} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame. Note this operation may scale poorly as protein length increases at which point increasing the stride may become necessary.

distance_thresh [float] {5.0} Distance threshold used to define a ‘contact’ in Angstroms. Contacts are

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

get_dihedral_mutual_information(angle_name='psi', weights=False, stride=1, bwidth=0.6283185307179586)

Generate the full mutual information matrix for a specific diehdral type. The resulting matrix describes the mutual information between each dihedral angle as defined by the variable angle_name. A weights parameter can be passed if frames are to be re-weighted, but this requires that the (number of frames) / stride = the number of weights.

The mutual information for a pair of angles is determined by generating a histogram of each dihedral induvidually (p(phi1), p(phi2)) and the joint probability histogram (p(phi1,phi2)), and then computing the Shannon entropy associated with the single and joing probability histograms (H_phi1, H_phi2, H_phi1_phi2). The mutual information is then returned as

H_phi1 + H_phi2 - (H_phi1 * H_phi2 )

The easiest way to interpret these results is to normalize the inferred matrix using an equivalent matrix generated using a limiting polymer model (e.g. an EV or FRC simulation).

Return: Mutual information matrix ( n x n) where n is the number of that type of bonds in the protein.

angle_name String, must be one of the following options ‘chi1’,’phi, ‘psi’, ‘omega’.

stride [int] {20} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame.

bwidth [np.array} {np.pi/5.0) The width of the bins that will stretch from -pi to pi. np.pi/5.0 is probablty the smallest binsize you want - even np.pi/2 should work well. You may want to experiment with this parameter…

get_distanceMap(weights=False, verbose=True, mode='CA', RMS=False, stride=1)

Function to calculate the CA defined distance map for a protein of interest. Note this function doesn’t take any arguments and instead will just calculate the complete distance map.

NOTE that distance maps are calculated between all CA-CA distances and NOT center of mass positions. This also means ACE/NME caps are EXCLUDED from this anlysis.

Distance is described in Angstroms.

stride [int] {1} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

weights [list or array of floats] Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

verbose [bool] {True} Print messages (or not)

mode [string] {‘CA’} String, must be one of either ‘CA’ or ‘COM’. ‘CA’ = alpha carbon ‘COM’ = center of mass (associated withe the residue)

RMS [Bool] {False} If set to False, scaling map reports ensemble average distances (this is the standard and default behaviour). If True, then the distance reported is the root mean squared (RMS) = i.e. SQRT(<r_ij^2>), which is the formal order parameter that should be used for polymeric distance properties.

get_end_to_end_distance(mode='COM')

Returns the vector assocaiated with the end-to-end distance for the simulation.

The vector of End-to-end distances is returned in angstroms

mode [string] {‘CA’} String, must be one of either ‘CA’ or ‘COM’. ‘CA’ = alpha carbon ‘COM’ = center of mass (associated withe the residue)

get_end_to_end_vs_rg_correlation(mode='COM')

Computes the correlation between Rg^2 and end-to-end^2.

Return: A single float describing the correlation (as calculated by np.corrcoef)

mode [string] {‘CA’} String, must be one of either ‘CA’ or ‘COM’. ‘CA’ = alpha carbon ‘COM’ = center of mass (associated withe the residue)

get_gyration_tensor(R1=False, R2=False, statusMessage=True, correctOffset=True)

Returns the instantaneous gyration tensor associated with each frame.

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_interResidueCOMDistance(R1, R2, stride=1, correctOffset=True)

Function which calculates the complete set of distances between two residues’ centers of mass (COM) and returns a vector of the distances.

Distance is returned in Angstroms.

R1 [int] Residue index of first residue

R2 [int] Residue index of second residue

stride [int] {20} Defines the spacing betwen frames to compare with - i.e. take every $stride-th frame. Setting stride=1 would mean every frame is used, which would mean you’re doing an all vs. all comparions, which would be ideal BUT may be slow.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_interResidueCOMVector(R1, R2, stride=1, correctOffset=True)

Function which calculates the complete set of distances between two residues’ centers of mass (COM) and returns the inter-residue distance vector.

NOTE: This gives a VECTOR and not the distance between the two centers of mass (which is calculated by get_interResidueCOMDistance)

R1 [int] Residue index of first residue

R2 [int] Residue index of second residue

stride [int] {20} Defines the spacing betwen frames to compare with - i.e. take every $stride-th frame. Setting stride=1 would mean every frame is used, which would mean you’re doing an all vs. all comparions, which would be ideal BUT may be slow.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_interResidue_atomic_distance(R1, R2, A1='CA', A2='CA', stride=1, mode='atom', correctOffset=True)

Function which returns the distance between two specific atoms on two residues. The atoms selected are based on the ‘name’ field from the topology selection language. This defines a specific atom as defined by the PDB file. By default A1 and A2 are CA (C-alpha) but one can define any residue of interest.

We do not perform any sanity checking on the atom name - this gets really hard - so have an explicit try/except block which will warn you that you’ve probably selected an illegal atom name from the residues.

Distance is returned in Angstroms.

R1 [int] Residue index of first residue

R2 [int] Residue index of second residue

A1 [string] {CA} Atom name of the atom in R1 we’re looking at

A2 [string {CA} Atom name of the atom in R2 we’re looking at

stride [int] {1} Defines the spacing betwen frames to compare with - i.e. take every $stride-th frame. Setting stride=1 would mean every frame is used, which would mean you’re doing an all vs. all comparions, which would be ideal BUT may be slow.

mode [string] {‘atom’} Mode allows the user to define differnet modes for computing atomic distance. The default is ‘atom’ whereby a pair of atoms (A1 and A2) are provided. Other options are detailed below and are identical to those offered by mdtraj in compute_contacts

‘ca’ - same as setting ‘atom’ and A1=’CA’ and A2=’CA’, this uses the C-alpha atoms

‘closest’ - closest atom associated with each of the residues, i.e. the is the point

of closest approach between the two residues

‘closest-heavy’ - same as closest, except only non-hydrogen atoms are considered

‘sidechain’ - closest atom where that atom is in the sidechain. Note this requires

mdtraj version 1.8.0 or higher.

‘sidechain-heavy’ - closest atom where that atom is in the sidechain and is heavy.

Note this requires mdtraj version 1.8.0 or higher.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_internal_scaling(R1=False, R2=False, mode='COM', stride=1, correctOffset=True, mean_vals=False, weights=False)

Calculates the raw internal scaling info for the protein in the simulation. R1 and R2 define a sub-region to operate over if sub-regional analysis is required. When residues are not provided the full protein’s internal scaling (EXCLUDING* the caps, if present) is calculated.

Distance is measured in Angstroms.

Returns two lists of the same length

  1. List of arrays, where each array is the simulation average set of inter-residue distances for the primary sequence separation defined by the equivalent position in the second array. Each array in this list will be a different length as there are many more i to i+1 pairs of residues than i to i+10 (e.g. in a 15 residue sequence).

  2. The sequence separation associated with each set of measurements (i.e. a single list which will normally be 0,1,2,3,…n where n = number of residues in sequence.

The internal scaling profile is a plot of sequence separation vs. mean through-space distance for all pairs of residues at a given sequence separation. What this means is that if we had a 6 residue peptide the internal scaling profile would be calculated as follows.

sequence separation = 0 average distance(average distance of 1-to-1, 2-to-2, 3-to-3, etc.)

sequence separation = 1 average distance(average distance of 1-to-2, 2-to-3, 3-to-4, etc.)

sequence separation = 2 average distance(average distance of 1-to-3, 2-to-4, 3-to-5, etc.)

sequence separation = 3 average distance(average distance of 1-to-4, 2-to-5, 3-to-6.)

sequence separation = 4 average distance(average distance of 1-to-5, 2-to-6)

sequence separation = 5 average distance(average distance of 1-to-6)

The residues considered for internal scaling analysis DO NOT include the ACE/NME peptide caps if present. This differs from CAMPARI, which DOES include the peptide caps.

For more information on this and other ideas for how polymer-physics can be a useful way of thinking about proteins, take a look at

Mao, A.H., Lyle, N., and Pappu, R.V. (2013). Describing sequence-ensemble relationships for intrinsically disordered proteins. Biochem. J 449, 307-318.

and

Pappu, R.V., Wang, X., Vitalis, A., and Crick, S.L. (2008). A polymer physics perspective on driving forces and mechanisms for protein aggregation - Highlight Issue: Protein Folding. Arch. Biochem. Biophys. 469, 132-141.

  • The exclusion of the caps is different to CAMPARI, which includes the caps, if present.

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

mode [‘CA’ or ‘COM’] {‘CA’} Defines the mode used to define the residue position, either the residue center or mass or the residue CA atom. The provided mode must be one of these two options.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

mean_vals [Book] {False} This is False by default, but if true the mean IS is returned instead of the explicit values. In reality the non-default behaviour is probably preferable…

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

get_internal_scaling_RMS(R1=False, R2=False, mode='COM', stride=1, correctOffset=True, weights=False)

If \(r_{i,j} = \langle \langle \sum \sigma_{1}\) equals \(\sigma_{2}\) then etc, etc.

Calculates the averaged internal scaling info for the protein in the simulation in terms of root mean square (i.e. sqrt(<Rij^2>) vs | i - j |.

R1 and R2 define a sub-region to operate over if sub-regional analysis is required. When residues are not provided the full protein’s internal scaling (EXCLUDING* the caps, if present) is calculated.

Distance is measured in Angstroms.

Returns two lists of the same length:

  1. sequence separation (|i - j|)

  2. mean sqrt(<Rij^2>)

The internal scaling profile is a plot of sequence separation vs. mean through-space distance for all pairs of residues at a given sequence separation. What this means is that if we had a 6 residue peptide the internal scaling profile would be calculated as follows.

sequence separation = 0 average distance(average distance of 1-to-1, 2-to-2, 3-to-3, etc.)

sequence separation = 1 average distance(average distance of 1-to-2, 2-to-3, 3-to-4, etc.)

sequence separation = 2 average distance(average distance of 1-to-3, 2-to-4, 3-to-5, etc.)

sequence separation = 3 average distance(average distance of 1-to-4, 2-to-5, 3-to-6.)

sequence separation = 4 average distance(average distance of 1-to-5, 2-to-6)

sequence separation = 5 average distance(average distance of 1-to-6)

The residues considered for internal scaling analysis DO NOT include the ACE/NME peptide caps if present. This differs from CAMPARI, which DOES include the peptide caps.

For more information on this and other ideas for how polymer-physics can be a useful way of thinking about proteins, take a look at

Mao, A.H., Lyle, N., and Pappu, R.V. (2013). Describing sequence-ensemble relationships for intrinsically disordered proteins. Biochem. J 449, 307-318.

and

Pappu, R.V., Wang, X., Vitalis, A., and Crick, S.L. (2008). A polymer physics perspective on driving forces and mechanisms for protein aggregation - Highlight Issue: Protein Folding. Arch. Biochem. Biophys. 469, 132-141.

  • The exclusion of the caps is different to CAMPARI, which includes the caps, if present.

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

mode [‘CA’ or ‘COM’] {‘CA’} Defines the mode used to define the residue position, either the residue center or mass or the residue CA atom. The provided mode must be one of these two options.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

get_local_heterogeneity(fragmentSize=10, stride=20, bins=array([0., 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1., 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4, 1.41, 1.42, 1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5, 1.51, 1.52, 1.53, 1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6, 1.61, 1.62, 1.63, 1.64, 1.65, 1.66, 1.67, 1.68, 1.69, 1.7, 1.71, 1.72, 1.73, 1.74, 1.75, 1.76, 1.77, 1.78, 1.79, 1.8, 1.81, 1.82, 1.83, 1.84, 1.85, 1.86, 1.87, 1.88, 1.89, 1.9, 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97, 1.98, 1.99, 2., 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 2.09, 2.1, 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19, 2.2, 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3, 2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4, 2.41, 2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5, 2.51, 2.52, 2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6, 2.61, 2.62, 2.63, 2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7, 2.71, 2.72, 2.73, 2.74, 2.75, 2.76, 2.77, 2.78, 2.79, 2.8, 2.81, 2.82, 2.83, 2.84, 2.85, 2.86, 2.87, 2.88, 2.89, 2.9, 2.91, 2.92, 2.93, 2.94, 2.95, 2.96, 2.97, 2.98, 2.99, 3., 3.01, 3.02, 3.03, 3.04, 3.05, 3.06, 3.07, 3.08, 3.09, 3.1, 3.11, 3.12, 3.13, 3.14, 3.15, 3.16, 3.17, 3.18, 3.19, 3.2, 3.21, 3.22, 3.23, 3.24, 3.25, 3.26, 3.27, 3.28, 3.29, 3.3, 3.31, 3.32, 3.33, 3.34, 3.35, 3.36, 3.37, 3.38, 3.39, 3.4, 3.41, 3.42, 3.43, 3.44, 3.45, 3.46, 3.47, 3.48, 3.49, 3.5, 3.51, 3.52, 3.53, 3.54, 3.55, 3.56, 3.57, 3.58, 3.59, 3.6, 3.61, 3.62, 3.63, 3.64, 3.65, 3.66, 3.67, 3.68, 3.69, 3.7, 3.71, 3.72, 3.73, 3.74, 3.75, 3.76, 3.77, 3.78, 3.79, 3.8, 3.81, 3.82, 3.83, 3.84, 3.85, 3.86, 3.87, 3.88, 3.89, 3.9, 3.91, 3.92, 3.93, 3.94, 3.95, 3.96, 3.97, 3.98, 3.99, 4., 4.01, 4.02, 4.03, 4.04, 4.05, 4.06, 4.07, 4.08, 4.09, 4.1, 4.11, 4.12, 4.13, 4.14, 4.15, 4.16, 4.17, 4.18, 4.19, 4.2, 4.21, 4.22, 4.23, 4.24, 4.25, 4.26, 4.27, 4.28, 4.29, 4.3, 4.31, 4.32, 4.33, 4.34, 4.35, 4.36, 4.37, 4.38, 4.39, 4.4, 4.41, 4.42, 4.43, 4.44, 4.45, 4.46, 4.47, 4.48, 4.49, 4.5, 4.51, 4.52, 4.53, 4.54, 4.55, 4.56, 4.57, 4.58, 4.59, 4.6, 4.61, 4.62, 4.63, 4.64, 4.65, 4.66, 4.67, 4.68, 4.69, 4.7, 4.71, 4.72, 4.73, 4.74, 4.75, 4.76, 4.77, 4.78, 4.79, 4.8, 4.81, 4.82, 4.83, 4.84, 4.85, 4.86, 4.87, 4.88, 4.89, 4.9, 4.91, 4.92, 4.93, 4.94, 4.95, 4.96, 4.97, 4.98, 4.99, 5., 5.01, 5.02, 5.03, 5.04, 5.05, 5.06, 5.07, 5.08, 5.09, 5.1, 5.11, 5.12, 5.13, 5.14, 5.15, 5.16, 5.17, 5.18, 5.19, 5.2, 5.21, 5.22, 5.23, 5.24, 5.25, 5.26, 5.27, 5.28, 5.29, 5.3, 5.31, 5.32, 5.33, 5.34, 5.35, 5.36, 5.37, 5.38, 5.39, 5.4, 5.41, 5.42, 5.43, 5.44, 5.45, 5.46, 5.47, 5.48, 5.49, 5.5, 5.51, 5.52, 5.53, 5.54, 5.55, 5.56, 5.57, 5.58, 5.59, 5.6, 5.61, 5.62, 5.63, 5.64, 5.65, 5.66, 5.67, 5.68, 5.69, 5.7, 5.71, 5.72, 5.73, 5.74, 5.75, 5.76, 5.77, 5.78, 5.79, 5.8, 5.81, 5.82, 5.83, 5.84, 5.85, 5.86, 5.87, 5.88, 5.89, 5.9, 5.91, 5.92, 5.93, 5.94, 5.95, 5.96, 5.97, 5.98, 5.99, 6., 6.01, 6.02, 6.03, 6.04, 6.05, 6.06, 6.07, 6.08, 6.09, 6.1, 6.11, 6.12, 6.13, 6.14, 6.15, 6.16, 6.17, 6.18, 6.19, 6.2, 6.21, 6.22, 6.23, 6.24, 6.25, 6.26, 6.27, 6.28, 6.29, 6.3, 6.31, 6.32, 6.33, 6.34, 6.35, 6.36, 6.37, 6.38, 6.39, 6.4, 6.41, 6.42, 6.43, 6.44, 6.45, 6.46, 6.47, 6.48, 6.49, 6.5, 6.51, 6.52, 6.53, 6.54, 6.55, 6.56, 6.57, 6.58, 6.59, 6.6, 6.61, 6.62, 6.63, 6.64, 6.65, 6.66, 6.67, 6.68, 6.69, 6.7, 6.71, 6.72, 6.73, 6.74, 6.75, 6.76, 6.77, 6.78, 6.79, 6.8, 6.81, 6.82, 6.83, 6.84, 6.85, 6.86, 6.87, 6.88, 6.89, 6.9, 6.91, 6.92, 6.93, 6.94, 6.95, 6.96, 6.97, 6.98, 6.99, 7., 7.01, 7.02, 7.03, 7.04, 7.05, 7.06, 7.07, 7.08, 7.09, 7.1, 7.11, 7.12, 7.13, 7.14, 7.15, 7.16, 7.17, 7.18, 7.19, 7.2, 7.21, 7.22, 7.23, 7.24, 7.25, 7.26, 7.27, 7.28, 7.29, 7.3, 7.31, 7.32, 7.33, 7.34, 7.35, 7.36, 7.37, 7.38, 7.39, 7.4, 7.41, 7.42, 7.43, 7.44, 7.45, 7.46, 7.47, 7.48, 7.49, 7.5, 7.51, 7.52, 7.53, 7.54, 7.55, 7.56, 7.57, 7.58, 7.59, 7.6, 7.61, 7.62, 7.63, 7.64, 7.65, 7.66, 7.67, 7.68, 7.69, 7.7, 7.71, 7.72, 7.73, 7.74, 7.75, 7.76, 7.77, 7.78, 7.79, 7.8, 7.81, 7.82, 7.83, 7.84, 7.85, 7.86, 7.87, 7.88, 7.89, 7.9, 7.91, 7.92, 7.93, 7.94, 7.95, 7.96, 7.97, 7.98, 7.99, 8., 8.01, 8.02, 8.03, 8.04, 8.05, 8.06, 8.07, 8.08, 8.09, 8.1, 8.11, 8.12, 8.13, 8.14, 8.15, 8.16, 8.17, 8.18, 8.19, 8.2, 8.21, 8.22, 8.23, 8.24, 8.25, 8.26, 8.27, 8.28, 8.29, 8.3, 8.31, 8.32, 8.33, 8.34, 8.35, 8.36, 8.37, 8.38, 8.39, 8.4, 8.41, 8.42, 8.43, 8.44, 8.45, 8.46, 8.47, 8.48, 8.49, 8.5, 8.51, 8.52, 8.53, 8.54, 8.55, 8.56, 8.57, 8.58, 8.59, 8.6, 8.61, 8.62, 8.63, 8.64, 8.65, 8.66, 8.67, 8.68, 8.69, 8.7, 8.71, 8.72, 8.73, 8.74, 8.75, 8.76, 8.77, 8.78, 8.79, 8.8, 8.81, 8.82, 8.83, 8.84, 8.85, 8.86, 8.87, 8.88, 8.89, 8.9, 8.91, 8.92, 8.93, 8.94, 8.95, 8.96, 8.97, 8.98, 8.99, 9., 9.01, 9.02, 9.03, 9.04, 9.05, 9.06, 9.07, 9.08, 9.09, 9.1, 9.11, 9.12, 9.13, 9.14, 9.15, 9.16, 9.17, 9.18, 9.19, 9.2, 9.21, 9.22, 9.23, 9.24, 9.25, 9.26, 9.27, 9.28, 9.29, 9.3, 9.31, 9.32, 9.33, 9.34, 9.35, 9.36, 9.37, 9.38, 9.39, 9.4, 9.41, 9.42, 9.43, 9.44, 9.45, 9.46, 9.47, 9.48, 9.49, 9.5, 9.51, 9.52, 9.53, 9.54, 9.55, 9.56, 9.57, 9.58, 9.59, 9.6, 9.61, 9.62, 9.63, 9.64, 9.65, 9.66, 9.67, 9.68, 9.69, 9.7, 9.71, 9.72, 9.73, 9.74, 9.75, 9.76, 9.77, 9.78, 9.79, 9.8, 9.81, 9.82, 9.83, 9.84, 9.85, 9.86, 9.87, 9.88, 9.89, 9.9, 9.91, 9.92, 9.93, 9.94, 9.95, 9.96, 9.97, 9.98, 9.99]))

Function to calculate the vector of D values used to calculate the Phi parameter from Lyle et al[1]. The stride defines the spacing between frames which are analyzed. This is just for practical purposes. The Phi calulation computes a D value for each frame vs. frame comparison - for a 2000 frame simulation this would be 4 Million D values if every value was calculated which is a bit much, so the stride let’s you define how many frames you should skip.

For a 2000 frame trajectory of a 80 residue protein with a stride=20 allows the calculation to take about 5 seconds. However, as protein size increases the computational cost of this process grows rapidly.

fragmentSize [int] {10} Size of local region that is considered to be a single unit over which structural heterogeneity is examined.

stride [int] {20} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame.

bins [numpy array] {np.arange(0,1,0.01)} Bins used to capture the heterogeneity at each position

get_local_to_global_correlation(mode='COM', stride=20, n_cycles=100, max_num_pairs=10, weights=False)

Method to analyze how well ensemble average distances taken from a finite number of inter-residue distances correlate with global dimensions as measured by the radius of gyration. This is a new analysis, that is best explained through a formal write up.

Return: Returns a four-place tuple. [0] - This is an 2 by (n_cycles*max_num_pairs) array, where the first column is the number of pairs

and the second column is the Rg-lambda correlation for a specific set of pairs (the pairs in question are not included). This can be thought of as the ‘raw’ data, and may only be useful if distributions of the correlation are of interest (e.g. for generating 2D histograms).

[1] - Array with the number of pairs used (e.g. if max_num_pairs = 10 then this would be

[1,2,3,4,5,6,7,8,9]

[2] - Array with the mean correlation associated with the numbers of pairs in position 2

and the radius of gyration

[3] - Array with the standard deviation of the correlation associated with the number of

pairs in position 2 and the radius of gyration. Note the inclusion of the standard deviation makes the assumption that the distribution is Gaussian which may or may not be true.

mode [string] {‘CA’} String, must be one of either ‘CA’ or ‘COM’. ‘CA’ = alpha carbon ‘COM’ = center of mass (associated withe the residue)

stride [int] {20} Defines the spacing between frames for calculating the ensemble average. As stride gets larger this analysis gets slower. It is worth experimenting with to see how the results change as a function of stride, but in theory the accuracy should remain fixed but precision improved as stride is reduced.

n_cycles [int] {100} Number of times, for each number of pairs, we re-select a different set of paris to use. This depends (weakly) on the number of residues, but we do not recommend a value < 50. For larger proteins this number should be increased. Again, it’s worth examining how your results chain as a function of n_cycles to determine the optimal tradeoff between speed and precision.

max_num_pairs [int} {10} The maximum number of pairs to be consider for correlation analysis. In general we’ve found above 10-20 the average correlation tends to plateau towards 1 (note it will NEVER be 1) so 1-20 is a reasonable set to use.

get_multiple_CAindex(resID_list=None, correctOffset=True)

Returns the atom indices associated with the C-alpha atom for the residues defined in the resID_list OR for all residues if no list is provided.

resID_list [list of int] {None} Defines a list of residues for which the C-alpha atom index will be retrieved. If no list is provided we simply

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again

get_numberOfFrames()

Returns the number of frames associated with a trajectory

get_numberOfResidues()

Returns the number of residues in this protein object

get_persistence_length()

Returns the a 4 by n numpy array inwhich column 1 gives residue number, column 2 is local helicity,

get_polymer_scaled_distance_map(nu=None, A0=None, min_separation=10, mode='fractional-change', weights=False, stride=1)

Function that allows for a global assesment of how well all i-j distances conform to standard polymer scaling behaviour (i.e. r_ij = A0*|i-j|^{nu}).

Essentially, this generates a distance map (2D matrix of i vs. j distances) where that distance is either normalized by the expected distance for a provided homopolymer model, or quantifies the fractional deviation from that homopolymer mode. These two modes are explained in more detail below.

In this standard scaling relationship: r_ij : Average inter-residue distance of residue i and j A0 : Scaling prefactor. Note this is NOT the same numerical value as the R0 prefactor that

defines the relationship Rg = R0*N^{nu}.

|i-j| : Sequence separation between residues i and j nu : The intrinsic polymer scaling exponent

This is the scaling behaviour expected for a standard homopolymer. This fucnction then assess how well this relationship holds for ALL inter-residue distances.

This function returns a four position tuple. Position 1 is an n x n numpy matrix (where n = sequence length), where the element is either the default value OR quantifes the deviation from a polymer model in one of two ways. Positions two, three are the nu and A0 values used, respectively. Finally, position 4 will be the reduced chi-square fitting to the polymer model for the internal scaling profile (i.e. how A0 and nu are originally calculated).

If no options are provided, the function calculates the best fit to a homopolymer mode using the default parameters associated with the get_scaling_exponent_v2() function, and then uses this model to determine pairwise deviations.

nu [float] {None} Scaling exponent used (if provided). Note for a provided nu to be used, both nu and A0 must be provided.

A0 [float] {None} Scaling prefactor used (if provided). Note for a provided A0 to be used, both A0 and nu must be provided.

min_separation [int] {10} Minimum distance for which deviations are calculated. At close distances, we expect local steric effects to cause deviations from a polymer model, so this value defines the threshold minimum distance to be used.

mode [‘fractional-change’, ‘signed-fractional-change’, ‘signed-absolute-change’, or ‘scaled’] {‘fractional-change’} Defines the mode in which deviation from a homopolymer model is calculated.

fractional-change: each inter-residue deviation is calculated as:

d_ij = abs(r_ij - polymer_ij)/polymer_ij

Where r_ij is the mean distance from the simulation for residues i and j, and polymer_ij is the expected distance for any pair of residues that are separated by |i-j| distance in the polymer model

signed-fractional-change: each inter-residue deviation is calculated as:

d_ij = (r_ij - polymer_ij)/polymer_ij

i.e. the same as the fractional-change, except a sign is now also included. Positive values mean there is expansion with respect to the homopolymer behaviour, while negative values mean there is contraction with respect to the homopolymer model.

signed-absolute-change: each inter-residue deviation is calculated as:

d_ij = (r_ij - polymer_ij)

i.e. the same as the signed-fractional-change, except now it is no longer fraction but in absolute distance units. This can be useful for getting a sense of by how-much the real behaviour deviates from the model in terms of Angstroms

scaled: each inter-residue deviation is calculated as:

d_ij = r_ij/polymer_ij

Where r_ij is and polymer_ij are defined as above.

weights [list or array of floats] {None} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

get_radius_of_gyration(R1=False, R2=False, correctOffset=True)

Returns the radius of gyration associated with the region defined by the intervening stretch of residues between R1 and R2. When residues are not provided the full protein’s radius of gyration (INCLUDING the caps, if present) is calculated.

Radius of gyration is returned in Angstroms.

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_regional_SASA(R1, R2, stride=20, probe_radius=0.14, correctOffset=True)

Returns the Solvent Accessible Surface Area (SASA) for a local region in every stride-th frame. SASA is determined using shrake_rupley algorithm.

SASA is returned in Angstroms squared, BUT PROBE RADIUS is in nanometers!

R1 [int] Index value for the first residue in the region

R2 [int] Index value for the last residue in the region

stride [int] {20} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

probe_radius [float] {0.14} Radius of the solvent probe used in nm. Uses the Golden-Spiral algorithm. 0.14 nm is pretty standard. NOTE - the probe radius must be in nanometers

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_residue_COM(R1, correctOffset=True)

Returns the COM vector for the residue across the trajectory.

R1 [int] Index value for the residue of interest

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_residue_index_list()

Function that returns the list of residue index values for this protein. i.e. this is the correctly offset residue list and will include NME/ACE caps if present.

No arguments taken.

get_residue_mass(R1, correctOffset=True)

Returns the mass associated with a specific residue.

R1 [int] Residue index of residue to examine

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_scaling_exponent(inter_residue_min=15, end_effect=5, stride=1, correctOffset=True, weights=False, subdivision_batch_size=20, mode='COM')

Estimation for the A0 and nu exponents for the standard polymer relationship

sqrt(<Rij^2>) = A0|i-j|^(nu)

Nu reports on the solvent quality, while the prefactor (A0) reports on the average chain persistence length. For polymers that are above a 0.5 scaling exponent this works, but below this they deviate from fractal behaviour, so formally this relationship stops working. In practice, the best possible fit line does still track with relative compactness.

Returns a 9 position tuple with the following associated values: 0 - best nu 1 - best A0 2 - minimum nu identified in bootstrap fitting 3 - maximum nu identified in bootstrap fitting 4 - minimum A0 identified in bootstrap fitting 5 - maximum A0 identified in bootstrap fitting 6 - Root mean square error (fit vs. real) 7 - chi2 (fit vs. real) 8 - 3-column array, where col 1 is the sequence separation, col 2 is

the real spatial separation observed and col 3 is the best fit curve.

NOTE: Despite their precision nu and A0 should be treated as qualitative metrics, and are subject to finite

chain effects. The idea of a polymer scaling behaviour is only necessarily useful in the case of a homopolymer, whereas heterpolymers engender massive averaging that can mask underlying conformational complexity. We _strongly_ caution against over interpretation of the scaling exponent.

inter_residue_min [int] {15} Minimum distances used when selecting pairs of residues. This 25 threshold was determined previously, and essentially avoids scenarios where the two residues (i and j) are close to each other. The goal of this limit is to avoid finite chain size strongly influencing the scaling exponent limit.

end_effect [int] {5} Avoid pairs where one of the residues is $end_effect residues from the end. Helps mitigate end-effects. 5 chosen as it’s around above the blob-length in a polypeptide. Note that for homopolymers this is much less of an issue.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

get_scaling_exponent_v2(inter_residue_min=15, end_effect=5, stride=1, correctOffset=True, weights=False, subdivision_batch_size=20, mode='COM', num_fitting_points=40, fraction_of_points=0.5, fraction_override=False, verbose=True)

Estimation for the A0 and nu exponents for the standard polymer relationship

sqrt(<Rij^2>) = A0|i-j|^(nu)

Nu reports on the solvent quality, while the prefactor (A0) reports on the average chain persistence length. For polymers that are above a 0.5 scaling exponent this works, but below this they deviate from fractal behaviour, so formally this relationship stops working. In practice, the best possible fit line does still track with relative compactness.

Returns a 9 position tuple with the following associated values: 0 - best nu

1 - best A0

2 - minimum nu identified in bootstrap fitting

3 - maximum nu identified in bootstrap fitting

4 - minimum A0 identified in bootstrap fitting

5 - maximum A0 identified in bootstrap fitting

6 - reduced chi^2 for the fit region

7 - reduced chi^2 for ALL points

8 - 2-column array, where col 1 is the sequence separation and col 2

is the real spatila separation for the ACTUAL data used to fit to the polymer model (i.e. these points are uniformly spaced from one another on a log-log plot). Reduced chi^2 for the fit region is calculated using this dataset.

9 - 3-column array, where col 1 is the sequence separation, col 2 is

the real spatial separation observed and col 3 is the best fit curve, for ALL i-j distances. Reduced chi^2 for all points is calculated using this dataset.

NOTE: Despite their precision nu and A0 should be treated as qualitative metrics, and are subject to finite

chain effects. The idea of a polymer scaling behaviour is only necessarily useful in the case of a homopolymer, whereas heterpolymers engender massive averaging that can mask underlying conformational complexity. We _strongly_ caution against over interpretation of the scaling exponent. For a better assement of how your chain actually deviates from homopolymer behaviour, see the function get_polymer_scaled_distance_map()

inter_residue_min [int] {15} Minimum distances used when selecting pairs of residues. This 25 threshold was determined previously, and essentially avoids scenarios where the two residues (i and j) are close to each other. The goal of this limit is to avoid finite chain size strongly influencing the scaling exponent limit.

end_effect [int] {5} Avoid pairs where one of the residues is $end_effect residues from the end. Helps mitigate end-effects. 5 chosen as it’s around above the blob-length in a polypeptide. Note that for homopolymers this is much less of an issue.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

weights [list or array of floats] {False} Defines the frame-specific weights if re-weighted analysis is required. This can be useful if an ensemble has been re-weighted to better match experimental data, or in the case of analysing replica exchange data that is re-combined using T-WHAM.

mode [string, either ‘COM’ or ‘CA’] {‘COM’} Defines the mode in which the internal scaling profile is calculated, can use either COM (center of mass) of each residue or the CA carbon of each residue. COM is more appropriate as CA will inherently give a larger profile.

num_fitting_points [int] {40} Number of evenly spaced points to used to fit the scaling exponent in loglog space. 40 seems to be a decent number that scales well

fraction_of_points [float between 0 and 1] {0.5} This is only used if fraction_override is set to True OR the sequence has less than the num_of_fitting_points residues. Means that instead of using a an absolute number of points (e.g. 40) to fit the loglog data, we use this fraction of residues. i.e. if the protein had 20 residues and fraction_of_points = 0.5 we’d use 10 points

fraction_override [bool] {False} If set to False then fraction_of_points ONLY used if the length of the sequence is less than the num_fitting points. If true then we explicitly use fraction_of_points and ignore num_fitting_points.

get_secondary_structure_DSSP(R1=False, R2=False, correctOffset=True)

Returns the a 4 by n numpy array inwhich column 1 gives residue number, column 2 is local helicity, column 3 is local ‘extended’ (beta strand/sheet) and column 4 is local coil on a per-residue basis.

Parameter R1 and R2 define a local region if only a sub-region is required.

Return vector provides normalized secondary structure (between 0 and 1) which reflects the fraction of the simulation each residue is in that particular secondary structure type.

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

get_sidechain_alignment_angle(R1, R2, sidechain_atom_1='default', sidechain_atom_2='default', correctOffset=True)

Function that computes the angle alignment between two residue sidechains. Sidechain vectors are defined as the unit vector between the CA of the residue and a designated ‘sidechain’ atom on the sidechain. The default sidechain atoms are listed below, but custom atom names can also be provided using the sidechain_atom_1/2 variables.

get_site_accessibility(input_list, probe_radius=0.14, stride=20, mode='residue_type')

Function to compute site/residue type accessibility. This can be done using one of two modes. Under ‘residue_type’ mode, the input_list should be a list of canonical 3-letter amino acid names. Under ‘resid’ mode, the input list should be a list of residue id positions (recall that resid ALWAYS starts from 0 and will include the ACE cap if present).

Returns a dictionary with the key = residue “type-number” string and the value equal the average and standard devaition associated with the SASA. Useful for looking at comparative accessibility of the same residue accross the sequence.

input_list [list] List of either residue names (e.g. [‘TRP’,’TYR’,’GLN’] or resid values ([1,2,3,4]) which will be taken and the SASA calculated for.

stride [int] {20} Defines the spacing between frames to compare - i.e. if comparing frame1 to a trajectory we’d compare frame 1 and every stride-th frame

probe_radius [float] {0.14} Radius of the solvent probe used in nm. Uses the Golden-Spiral algorithm. 0.14 nm is pretty standard. NOTE - the probe radius must be in nanometers

mode [string] {‘residue_type’} Mode used to examine sites. MUST be one of ‘residue_type’ or ‘resid’

get_t(R1=False, R2=False, correctOffset=True)

Returns the <t>, a dimensionless parameter which describes the size of the ensemble.

R1 [int] {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used.

R1 [int] {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used.

correctOffset [Bool] {True} Defines if we perform local protein offset correction or not. By default we do, but some internal functions may have already performed the correction and so don’t need to perform it again.

print_residues()

Function to help determine the mapping of residue ID to PDB residue value. Prints the mapping between resid and PDB residue, and returns this information in a list

References

BiDB79

Biskup, J.; Dayal, U.; Bernstein, P.A..: Synthesizing independent database schemas. In: ACM SIGMOD 1979 Int. Conf. On Management of Data Proceedings, S. 143-151.