Source code for camparitraj.ctprotein

"""
CTProtein is the main protein trajectory class in CAMPARITRAJ

"""

##
##                                       _ _              _ 
##   ___ __ _ _ __ ___  _ __   __ _ _ __(_) |_ _ __ __ _ (_)
##  / __/ _` | '_ ` _ \| '_ \ / _` | '__| | __| '__/ _` || |
## | (_| (_| | | | | | | |_) | (_| | |  | | |_| | | (_| || |
##  \___\__,_|_| |_| |_| .__/ \__,_|_|  |_|\__|_|  \__,_|/ |
##                     |_|                             |__/ 
##
##
## Alex Holehouse (Pappu Lab and Holehouse Lab)
## Simulation analysis package
## Copyright 2014 - 2021
##

import mdtraj as md
import numpy as np
from numpy import linalg as LA
from itertools import combinations
from scipy import stats
import scipy.optimize as SPO
from numpy.random import choice

from .configs import DEBUGGING
from .ctdata import THREE_TO_ONE, DEFAULT_SIDECHAIN_VECTOR_ATOMS, ALL_VALID_RESIDUE_NAMES
from .ctexceptions import CTException
from . import ctmutualinformation, ctio, cttools, ctpolymer, ctutils

from . _internal_data import BBSEG2

import scipy.cluster.hierarchy


## Order of standard args:
## 1. correctOffset
## 2. stride
## 3. weights
#  4. verbose
##
##


[docs]class CTProtein: """ 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 `CTProtein.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. """ ## >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ## ## A note for the code: ## The residue offset operation is performed by the __get_offset_residue fuction. ## For ANY function where a specific residue is supplied there must also be the ## option to perform this offset or not (i.e. each public facicing function ## should be able to stand alone and not rely on an offset being performed ## by another function) BUT the option to perform the offseting provided so it ## in functions that call other functions which can perform the offset it will ## only need to be performed once. ## ## >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> # ........................................................................ #
[docs] def __init__(self, traj, residue_offset): """ Initialize a CTProtein object instance using trajectory information, and information about offsets. Parameters ---------- traj: `cttrajectory.CTTrajectory` An instance of a system's trajectory populated via `cttrajectory.CTTrajectory`. residue_offset: int Similar to `atom_offset`, this is the residue index which will serve as the marker for residue index 0. The residue offset operation is performed by the `__get_offset_residue` fuction. For ANY function where a specific residue is supplied there must also be the option to perform this offset or not (i.e. each public facing function should be able to stand alone and not rely on an offset being performed by another function) BUT the option to perform the offseting provided so it in functions that call other functions which can perform the offset it will only need to be performed once. """ # set the trajectory object for easy access self.traj = traj self.topology = traj.topology # WARNING - at the moment it seems that while the trajectory # is a subset of full trajectory, the topology retains # its full content of RESIDUES (though not atoms) # this is somewhat annoying # same for residue o self.residue_offset = residue_offset if DEBUGGING: ctio.debug_message("Creating protein") ctio.debug_message("Residue offset : " + str(residue_offset)) r_string = '' for r in self.topology.chain(0).residues: r_string = r_string + (str) ctio.debug_message("Residue string from residues in self.topology.chain(0).residues: %s" %(r_string)) # delete the vaiable to avoid any possible introduction of this var into the namespace del r_string # initialze various protein-centric data self.__num_residues = sum( 1 for _ in self.topology.residues) # initialize some empty values that are populated on demand by functions that drive local # memoization. self.__amino_acids_3LTR = None self.__amino_acids_1LTR = None self.__residue_index_list = None self.__CA_residue_atom = {} self.__residue_atom_table = {} self.__residue_COM = {} (self.__resid_with_CA, self.__idx_with_CA) = self.__get_resid_with_CA() # define if caps are present or not - specifically, if the resid 0 is in the CA-containing if 0 in self.resid_with_CA: self.__ncap = False else: self.__ncap = True if (self.n_residues - 1) in self.idx_with_CA: self.__ccap = False else: self.__ccap = True
# <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> # # Properties @property def resid_with_CA(self): """ Return a list of resids that have CA atoms. Note these residues have already have their offset correction applied, and so are valid to be used for selection via the underlying mdtraj topology object (self.topology). The values associated with this property are built using the internal function `__get_resid_with_CA()` Returns -------- list of integers See Also ------------ idx_with_CA """ return self.__resid_with_CA @property def idx_with_CA(self): """ Return a list of zero-indexed IDs that have CA atoms. Note these indices start at zero regardless of what the internal .topology resid numbering is. Recall that for for mdtraj 1.9.5 or greater this should be the same as the values returned by `resid_with_CA`. The values associated with this property are built using the internal function `__get_resid_with_CA()` Returns -------- list of integers See Also ------------ resid_with_CA """ return self.__idx_with_CA @property def ncap(self): """ Flag that returns if an N-terminal capping residue is present (or not). Returns ---------- bool True if N-teminal cap is present, False if not """ return self.__ncap @property def ccap(self): """ Flag that returns if a C-terminal capping residue is present (or not). Returns ---------- bool True if C-teminal cap is present, False if not """ return self.__ccap @property def n_frames(self): """ Returns the number of frames in the trajectory Returns ---------- int Returns the number of frames in the simulation trajectory """ return self.traj.n_frames @property def n_residues(self): """ Returns the number of residues in the protein (including caps) Returns ---------- int Returns the number of frames in the protein """ return self.__num_residues @property def residue_index_list(self): """ 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. For backend implementation details, this returns the res.index values from all res in topology.residues (mdtraj.mdtrajectory) """ if self.__residue_index_list == None: reslist = [] for res in self.topology.residues: reslist.append(res.index) self.__residue_index_list = reslist return self.__residue_index_list def __repr__(self): return "CTProtein (%s): %i res and %i frames" % (hex(id(self)), self.n_residues, self.n_frames) def __len__(self): return (self.n_residues, self.n_frames) # ........................................................................ # def __check_weights(self, weights, stride=None): """ Function that checks a passed weights-array is usable and matches the number of frames (avoids a lot of heartache when something breaks deep inside the code). NOTE: This also typecasts weights to a `numpy.array` which allows them to be indexed directly using a list of values. Parameters ---------- weights : array_like An `numpy.array` object that corresponds to the number of frames within an input trajectory. stride : {None}, optional The stepsize used when iterating across the frames. A value of `None` will set the step size to 0. Returns ------- numpy.array An `np.array` object containing trajectory frames selected per `stride` number of frames. """ if weights is not False: if len(weights) != self.n_frames: raise CTException('Passed frame weights array is %i in length, while there are actually %i frames - these must match' % (len(weights), self.n_frames)) if stride > 1: ctio.warning_message("WARNING: Using stride with weights is ALMOST certainly not a good idea unless the weights are\ncalculated for every stride-th residue", with_frills=True) return np.array(weights[list(range(0,self.n_frames,stride))]) else: return np.array(weights) return False # ........................................................................ # def __get_first_and_last(self, R1, R2, withCA=False): """ Function which returns first and last residue for a range, correcting for the residue offset problem. The returned tup Parameters ---------- R1 : int or False/None {False} First residue in range - can be an integer (assumes first residue in chain indexed at 0). If `False` assume we start at 0. R2 : int or False/None {False} Last residue in range - can be an integer (assumes first residue in chain indexed at 0). If `False` assumes we're using the whole chain. withCA : bool {False} Flag which, if `True` and R1 or R2 are `False`, selects R1/R2 values that contain a CA, which basically means caps are dealt with here if present. Returns ------- tuple: Returns a tuple with three positions: - [0] = R1 with relevant offsets applid (`int`). - [1] = R2 with relevant offsets applid (`int`). - [2] = String that can be passed directly to topology select to extract the atoms associated with these positions. """ # first define as if we're starting from first and last residue with/without caps if R1 == False or R1 == None: if withCA: if self.ncap: R1 = 1 else: R1 = 0 else: R1 = 0 if R2 == False or R2 == None: if withCA: if self.ccap: R2 = self.n_residues - 2 else: R2 = self.n_residues - 1 else: R2 = self.n_residues - 1 # then apply the systematic offset R1 = R1 + self.residue_offset R2 = R2 + self.residue_offset # finally flip around if R1 is larger than R2 if R1 > R2: tmp = R2 R2 = R1 R1 = tmp return (R1, R2, "resid %i to %i" %(R1,R2)) # ........................................................................ #
[docs] def get_offset_residue(self, R1): """ Returns the true residue index (TRI) for this protein by taking into account the residue offset. Checks the residue first to ensure a valid residue has been passed. NOTE: As of MDTraj 1.9.5 this is ACTUALLY overkill as when new trajectories are created the resid numbering resets to starting at 0 (as opposed to retaining the old trajectory numbering). However, for 1.9.4 or lower this did not happen, such that usin the get_offset_residue() function provides long-lasting backwards compatibility. Parameters ----------- R1: int The zero-indexed input residue offset (`int`). This value is examined by `__check_single_residue()`. Returns -------- TRI: int The updated index which reflects the input offset to determine the **true** starting residue index. See Also -------- __check_single_residue """ self.__check_single_residue(R1) return R1 + self.residue_offset
# ........................................................................ # def __check_stride(self, stride): """ Checks that a passed stride value doesn't break everything. Returns `None` or raises a `CTException`. Parameters ---------- stride: int The non-zero number of steps to perform while iterating across a trajectory. Raises ------ CTException When the stride is larger than the number of available frames in the trajectory, or less than 1. """ if stride > self.n_frames: raise CTException('stride (%i) is larger than the number of frames (%i)' %(stride, self.n_frames)) if stride < 1: raise CTException('stride (%i) is less than 1' %(stride)) # ........................................................................ # def __check_single_residue(self, R1): """ Internal function that checks that a single residue provided makes sense in the context of this protein. NOTE this checks BEFORE an offset is applied (i.e. we assume once a residue offset has been made that the resid's validity has been established). Returns `None` or raises a `CTException`. Parameters ---------- R1: int The zero-indexed residue index (`int`) whose index is checked and validated. Raises ------ CTException When the residue ID is greater than the chain length, or when the distances explored are greater than the chain size. """ if R1 < 0: raise CTException("Trying to use a negative residue index [residue index = %i]"%R1) if R1 >= self.n_residues: raise CTException("Trying to use a residue ID greater than the chain length [residue index = %i, chain length = %i] " % (R1, self.n_residues)) if (R1 - self.residue_offset) >= self.n_residues: raise CTException("Trying to explore distances greater than chain size [residue index = %i, chain length = %i] " % (R1, self.n_residues)) # ........................................................................ # def __check_contains_CA(self, R1, R1_org=None, has_correctOffset=True): """ Function which checks if residue R1 (which is the residue AFTER an offset correction has been applied) contains a C-alpha atom. Returns `None` or raises a `CTException`. Parameters ---------- R1: int The zero-indexed residue index (`int`) whose index is checked and validated. R1_org: int or None {None} Raises ------ CTException When the `CTProtein` has an uncorrected offset and residue `R1` lacks a C-alpha atom. Or, when the offset residue `R1_org` which maps to `R1` lacks a C-alpha. And, when the offset converted residue lacks a C-alpha atom. """ exception_message = '' if has_correctOffset is False: if R1 not in self.idx_with_CA: raise CTException("Residue index %i (where first residue = 0) lacks a C-alpha carbon" % (R1)) return if R1 not in self.resid_with_CA: if R1_org: raise CTException("Residue %i (which maps to %i) lacks a C-alpha carbon [%s]" % (R1_org, R1, self.get_amino_acid_sequence()[R1])) else: raise CTException("Residue %i (offset converted residue ID) lacks a C-alpha carbon [%s]" %(R1, self.get_amino_acid_sequence()[R1])) # ........................................................................ # def __get_subtrajectory(self, traj, stride): """ Internal function which returns a subtrajectory. Expects `traj` to be an `mdtraj` trajectory object and `stride` to be an `int`. Parameters ---------- traj: mdtraj.Trajectory An instance of an `mdtraj.Trajectory` which is non-empty - i.e. contains at least 1 frame. stride: int The non-zero number of steps to perform while iterating across the input trajectory, `traj`. Returns ----------- mdtraj.Trajectory A sliced trajectory which contains the frames selected every `stride` step from the input trajectory, `traj`. """ stride = int(stride) self.__check_stride(stride) if stride == 1: return traj else: return traj.slice(list(range(0, self.n_frames, stride))) # ........................................................................ # def __get_resid_with_CA(self): """ Internal function which should only be needed during initialization. Defines the list of residues where CA atoms are present, and the list of zero-indexed residue indices which contain CA atoms. In the case that the first resid in the self.topology object is 0 then these two lists are the same (which for MDTraj 1.9.5 or higher should always be the same) but for systems with multiple protein chains in MDTraj 1.9.4 or lower the residuesWithCA and idxWithCA willn differe for the second chain and onwards. This list is then assigned to the property variable `self.resid_with_CA` and `self.idx_with_CA`. Note this is quite slow for large trajectories Returns ----------- tuple A 2-tuple that is comprised of lists: - [0] := The list of residue indices which contain C-alpha atoms selected from the topology. - [1] := The list of zero-indexed residue indices which contain C-alpha atoms. See also ------------ resid_with_CA idx_with_CA """ # initialize empty lists residuesWithCA = [] idxWithCA = [] idx = -1 # for each residue in the topology for res in self.topology.residues: idx = idx + 1 # try and get the CA atomic index try: self.get_CA_index(res.index, correctOffset=False) # if we could get a CA then append this residue index residuesWithCA.append(res.index) idxWithCA.append(idx) except CTException: continue return (residuesWithCA, idxWithCA) # ........................................................................ # def __residue_atom_lookup(self, resid, atomname=None): """ Memoisation function to lookup the atomic index of a specific residues atom. Originally I'd assumed the underlying MDTraj `topology.select()` operation was basically a lookup, BUT it turns out it's actually *really* expensive, so this method converts atom/residue lookup information into a dynamic O(1) operation, greatly improving the performance of a number of different methods in the processes. Parameters ---------- resid: int The residue index to lookup. This index must correspond to the range derived from the corrected offset. If the residue has not been cached, it will be added to the lookup table for later reuse. atomname: str or None {None} The name of the atom to lookup which will return the corresponding residue ID. Like the previous parameter, if that residue does not exist in the lookup table it will be added for later reuse. Returns ------- list A list containing all the atoms corresponding to a given residue id that match the input residue id (`resid`) or, the residue corresponding to the atom name (`atomname`). """ # if resid is not yet in table, create an empty dicitionary if resid not in self.__residue_atom_table: self.__residue_atom_table[resid] = {} # if we haven't specified an atom look up ALL the atoms associated with this residue if atomname is None: # if all_atoms not yet associated with this residue if 'all_atoms' not in self.__residue_atom_table[resid]: self.__residue_atom_table[resid]['all_atoms'] = self.topology.select('resid %i'%(resid)) # return set of all atoms return self.__residue_atom_table[resid]['all_atoms'] # if atom-name not yet associated with this resid lookup # the atomname from the underlying topology if atomname not in self.__residue_atom_table[resid]: self.__residue_atom_table[resid][atomname] = self.topology.select('resid %i and name %s'%(resid, atomname)) # at this point we know the resid-atomname pair is in the table # so goahead and look it up! return self.__residue_atom_table[resid][atomname] # ........................................................................ # def __get_selection_atoms(self, region=None, backbone=True, heavy=False, correctOffset=True): """ Function which returns a list of atoms associated with the residues defined by the region keyword. If no region is supplied this returns the entire region (NME/ACE caps included). Parameters ---------- region : `np.array`, `list`, or `tuple` {None} An array_like object of size 2 which defines the first and last residue (INCLUSIVE) for a region to be examined. backbone: bool {True} Boolean flag to determine if only the backbone atoms should be returned, or if all the full chain's atoms should be included (i.e. including sidechain). heavy: bool {False} Boolean flag to determine if we should only select heavy atoms or not (i.e. not H). 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 Returns ------- selectionatoms A `numpy.array` comprised of atom indices corresponding to the residues in a given region. Raises ------ CTException When the input region is larger than 2. """ # perform offset if necessary if correctOffset and not region == None: tmp = [] for i in region: tmp.append(self.get_offset_residue(i)) region = tmp if region is None: pass elif not len(region) == 2: CTException("Trying to select a subsection of atoms, but the provided 'region' tuple/list is not of exactly length two [region=%s].\nCould indicate a problem, so be safe raising an exception" % (str(region))) if not region == None and len(region) == 2: if backbone: if heavy: selectionatoms = self.topology.select('backbone and resid %i to %i and not type H)' % (region[0], region[1])) else: selectionatoms = self.topology.select('backbone and resid %i to %i' % (region[0], region[1])) else: if heavy: selectionatoms = self.topology.select('resid %i to %i and not type H' % (region[0], region[1])) else: selectionatoms = self.topology.select('resid %i to %i' % (region[0], region[1])) else: if backbone: if heavy: selectionatoms = self.topology.select('backbone and resid %i to %i and not type H' % ( self.residue_offset, self.residue_offset + self.n_residues)) else: selectionatoms = self.topology.select('backbone and resid %i to %i' % ( self.residue_offset, self.residue_offset + self.n_residues)) else: if heavy: selectionatoms = self.topology.select('resid %i to %i and not type H' % (self.residue_offset, self.residue_offset + self.n_residues)) else: selectionatoms = self.topology.select('resid %i to %i' % (self.residue_offset, self.residue_offset + self.n_residues)) return selectionatoms # ........................................................................ #
[docs] def print_residues(self, verbose=True): """ 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. Returns a list of lists, where each list element is itself a list of two elements, index position and the resname-resid from the PDB file. Parameters ---------- verbose : bool {True} If set to `True`, `print_residues()` will print out to screen and also return a list. If set to False, means nothing is printed to the screen. Returns ------- return_list List containing a mapping of the zero-indexed residues and their names. """ AA = self.get_amino_acid_sequence() return_list = [] for i in range(0, len(AA)): if verbose is True: print("%i --> %s" %(i, AA[i])) return_list.append([i,AA[i]]) return return_list
# ........................................................................ #
[docs] def get_amino_acid_sequence(self, oneletter=False, numbered=True): """ Returns the protein's amino acid sequence. Parameters ---------- 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. Returns ------- list A list comprised of the 1-letter or 3-letter names of the amino acid sequence. """ if oneletter: if self.__amino_acids_1LTR == None: res = [] for i in self.topology.residues: res.append(THREE_TO_ONE[str(i)[0:3].upper()]) self.__amino_acids_1LTR = "".join(res) else: if self.__amino_acids_3LTR == None: res = [] for i in self.topology.residues: res.append(str(i)[0:3]+"-"+str(i)[3:]) self.__amino_acids_3LTR = res # if numbered requsted if numbered: if oneletter: return self.__amino_acids_1LTR else: return self.__amino_acids_3LTR # else strip out the numbering with list comprehension else: if oneletter: return [x.split('-')[0] for x in self.__amino_acids_1LTR] else: return [x.split('-')[0] for x in self.__amino_acids_3LTR]
# ........................................................................ #
[docs] def get_CA_index(self, residueIndex, correctOffset=True): """ Get the CA atom index for the residue defined by residueIndex. Again does this via memoization - i.e. the first time a specific residue is requested the function looks up the information and then stores it locally in case its needed again. Defensivly checks for errors. Parameters ---------- 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. Returns ------- list A list of size 1 containing the CA atom index for the residue index, `residueIndex`. Raises ------ CTException When the number of CA atoms do not equal 1. """ # perform offset if necessary if correctOffset: residueIndex = self.get_offset_residue(residueIndex) if residueIndex not in self.__CA_residue_atom: return_val = self.__residue_atom_lookup(residueIndex, 'CA') if len(return_val) == 1: self.__CA_residue_atom[residueIndex] = return_val[0] else: raise CTException("get_CA_index - unable to find residue %i [corrected resid]" % residueIndex) return self.__CA_residue_atom[residueIndex]
# ........................................................................ #
[docs] def get_multiple_CA_index(self, resID_list=None, correctOffset=True): """ Returns the atom indices associated with the C-alpha (CA) atom for the residues defined in the resID_list OR for all residues, if no list is provided. Parameters ---------- 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 select the list of residues with C-alphas, whose indices have been corrected. 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. Returns ------- CAlist: list The list (`int`) of C-Alpha indices of the input list of residue IDs. """ # if we've just passed a single unlisted integer # then just return the single residue associated # with if type(resID_list) == int: return ([self.get_CA_index(resID_list, correctOffset)]) # if no value passed grab all the residues if resID_list == None: # this list uses corrected residue index positions. Note now the # correctOffset flag is irrelevant because we didn't pass any # positions to offset! resID_list = self.resid_with_CA else: if correctOffset: tmp = [] for resid in resID_list: tmp.append(self.get_offset_residue(resid)) resID_list = tmp CAlist = [] for res in resID_list: try: CAlist.append(self.get_CA_index(res, correctOffset=False)) except CTException as e: print(e) continue return CAlist
# ........................................................................ #
[docs] def calculate_all_CA_distances(self, residueIndex, mode='CA', onlyCterminalResidues=True, correctOffset=True, stride=1): """ 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. Parameters ---------- residueIndex: int Defines the residue index to select the CA from. mode: str {'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. 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. Returns ------- numpy.array Array containing the end-to-end distance measures based on the input mode. Raises ------ CTException If the input mode is nether 'CA' or 'COM'. """ # validate input ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') # first check the residue we passed is valid and offset if required if correctOffset: residueIndex = self.get_offset_residue(residueIndex) # determine atomic index of CA atom for the residue you passed in try: CA_base = self.get_CA_index(residueIndex, correctOffset=False) except CTException: # if we couldln't find a C-alpha for this residue then nothing # makes sense so return -1 return -1 # list of atomic indices for C-alpha atoms we care about CAlist = [] ## ## Block to use if using CA as base for computing distances ## if mode == 'CA': for residue in self.resid_with_CA: # only compute distances bigger than the residueIndex # such that by default full iteration calculates the # non-redundant map (i.e. the half diagonal) if residue <= residueIndex and onlyCterminalResidues: continue # note we have the correctOffset set the false because the residue # index here is the true residue index as was pre-calculated CAlist.append(self.get_CA_index(residue, correctOffset=False)) # now we construct a nested list of lists of pairs to compute distances # between pairs=[] for CA in CAlist: pairs.append([CA_base, CA]) if len(pairs) == 0: return np.array([]) local_traj = self.__get_subtrajectory(self.traj, stride) # 10* for angstroms return 10*md.compute_distances(local_traj, np.array(pairs)) if mode == 'COM': ## ## Block to use if using COM as base for computing distances ## return_distances = [] for residue in self.resid_with_CA: # only compute distances bigger than the residueIndex # such that by default full iteration calculates the # non-redundant map (i.e. the half diagonal) if residue <= residueIndex and onlyCterminalResidues: continue return_distances.append(self.get_inter_residue_COM_distance(residueIndex, residue)) # finally convert list to numpy array and flip so returns in same format as CA mode return np.array(return_distances).transpose()
# ........................................................................ #
[docs] def get_distance_map(self, mode='CA', RMS=False, stride=1, weights=False, verbose=True): """ 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. Parameters ---------- mode : str {'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. 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 Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns ------- tuple A 2-tuple containing: - [0] := The distance map derived from the measurements between CA atoms. - [1] := The standard deviation corresponding to the distance map. """ ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') weights = self.__check_weights(weights, stride) # get the list of residues which have CA (typically this means we exclude # ACE and NME if they're present, but they may not be residuesWithCA = self.resid_with_CA # initialize empty matrices that we're gonna fill up distanceMap = np.zeros([len(residuesWithCA),len(residuesWithCA),]) stdMap = np.zeros([len(residuesWithCA),len(residuesWithCA),]) SM_index=0 for resIndex in residuesWithCA[0:-1]: ctio.status_message("On protein residue %i (overall residue index = %i) of %i [distance calculations]"% (SM_index, resIndex, int(len(residuesWithCA))), verbose) # get all CA-CA distances between the residue of index resIndex and every other residue. Note this gives the non-redudant upper # triangle. No need to correct for offset because this was done when we retrived the set of residues with CA full_data = self.calculate_all_CA_distances(resIndex, mode=mode, stride=stride, correctOffset=False) # if we want root mean square then NOW square each distances if RMS: full_data = np.power(full_data,2) # calculate mean and standard deviation if weights is not False: mean_data = np.average(full_data,0,weights=weights) # if we want RMS then NOW take square root of <rij^2> if RMS: mean_data = np.sqrt(mean_data) std_data = None else: mean_data = np.mean(full_data,0) # if we want RMS then NOW take square root of <rij^2> if RMS: mean_data = np.sqrt(mean_data) std_data = np.std(full_data,0) # update the maps appropriately and increment the counter distanceMap[SM_index][1+SM_index:len(residuesWithCA)] = mean_data stdMap[SM_index][1+SM_index:len(residuesWithCA)] = std_data SM_index=SM_index+1 return (distanceMap, stdMap)
# ........................................................................ #
[docs] def get_polymer_scaled_distance_map(self, nu=None, A0=None, min_separation=10, mode='fractional-change', stride=1, weights=False, verbose=True): """ 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 a homopolymer model fit to the data. 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 function then assess how well this relationship holds for ALL possible 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 and 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). NOTE that the reduced chi-squared value will be -1 if nu and A0 are provided manually. 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() function, and then uses this model to determine pairwise deviations. Parameters ---------- 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 : str {'fractional-change'} Defines the mode in which deviation from a homopolymer model is calculated. Options are: `'fractional-change', 'signed-fractional-change', 'signed-absolute-change', 'scaled'`. *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. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns ------- tuple (len = 4) return[0] : 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. return[1] : float defining the scaling exponent nu return[2] : float defining the A0 prefactor return[3] : reduced chi-squared fitting to the polymer model (goodness of fit) Raises ------ CTException """ # First validate keyword ctutils.validate_keyword_option(mode, ['fractional-change','scaled', 'signed-fractional-change', 'signed-absolute-change'], 'mode') # next check that the minimum separation requested makes sense... (this is only a partial check) if min_separation < 1: raise CTException("Minimum separation to be used must be greater than 0") ## next see if nu or A0 have been provided... # If NEITHER provided then do fitting here and now! if nu == None and A0 == None: ctio.status_message("Fitting data to homopolymer mode...", verbose) # remind that this is the old get_scaling_exponent_v2() # ALSO note this uses COM which is also how the distance map is computed SE = self.get_scaling_exponent(verbose=False, weights = weights, stride = stride) nu = SE[0] A0 = SE[1] REDCHI = SE[7] elif nu is None: raise CTException("A0 parameter provided [%1.5f] but nu was not. Must provide BOTH or neither (in which case fitting is done)" %( A0)) elif A0 is None: raise CTException("nu parameter provided [%1.5f] but A0 was not. Must provide BOTH or neither (in which case fitting is done)" %( nu)) # else both were provided so we double check they're valid.. else: # sanity check nu if nu <= 0 or nu > 1: raise CTException("Nu parameter must be in interval 0 < nu <= 1 (and probably should be between 0.33 and 1.0...)") # sanity check A0 if A0 <= 0: raise CTException("A0 paameter must be greater than 0") # not computing a reduced chi REDCHI = -1 # We now define a function which will evaluate how we assess the deviation (or lack thereof) from the # traditional polymer scaling behaviour. This just saves us having if/then statements that get evaluated # on each loop of the analysis script below, and also makes it easy to add in additional 'modes' if mode == 'fractional-change': def d_funct(dMap_val, p_val): return abs(dMap_val - p_val)/p_val default_val = 0 elif mode == 'signed-fractional-change': def d_funct(dMap_val, p_val): return (dMap_val - p_val)/p_val default_val = 0 elif mode == 'signed-absolute-change': def d_funct(dMap_val, p_val): return (dMap_val - p_val) default_val = 0 elif mode == 'scaled': def d_funct(dMap_val, p_val): return dMap_val/p_val default_val = 1 ## Now we've set everything up we can actually compute some numbers # first compute and get the distance map (note [0] means we get first element # which is mean distance ([1] is STDEV) distance_map = self.get_distance_map(mode='COM', RMS=True, stride=stride, weights=weights, verbose=False)[0] # get distance map dimensions (will be a square so just take X-dim) dimensions = distance_map.shape[0] if dimensions <= min_separation: raise CTException('The minimum separation is shorter than the chain length') # compute expected distance given the standard polymer scaling model expected_distances = cttools.powermodel(list(range(0,dimensions)), nu, A0) # initialize the return matrix and then populate for distances that are # above the minimum threshold. We only populate upper right triangle return_matrix = np.zeros((dimensions, dimensions)) for i in range(0,dimensions): for j in range(0,dimensions): if j-i < min_separation: return_matrix[i,j] = default_val else: return_matrix[i,j] = d_funct(distance_map[i,j], expected_distances[(j-i)]) return (return_matrix, nu, A0, REDCHI)
# ........................................................................ #
[docs] def get_local_heterogeneity(self, fragment_size=10, bins=None, stride=20, verbose=True): """ 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. Parameters ----------- fragment_size : int {10} Size of local region that is considered to be a single unit over which structural heterogeneity is examined. Should be between 2 and the length of the sequence. bins : np.ndarray {np.arange(0,1,0.01)} Bins used to capture the heterogeneity at each position. If default a set of bins from 0 to 1 with an interval of 0.01 is used 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. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns ------- tuple (len = 4) return[0] : list of floats of len *n*, where each float reports on the mean RMSD-deviation at a specific position along the sequence as defined by the fragment_size return[1] : list of floats of len *n* , where each float reports on the standard deviation of the RMSD-deviation at a specific position along the sequence as defined by the fragment_size return[2] : List of np.ndarrays of len *n*, where each sub-array reports on the histogram values associated with the full RMSD distribution at a given position along the sequence return[3] : np.ndarray which corresponds to bin values for each of the histograms in return[2] Raises ------ CTException """ # validate bins if bins is None: bins = np.arange(0,10,0.01) else: try: if len(bins) < 2: raise CTException('Bins should be a numpy defined vector of values - e.g. np.arange(0,1,0.01)') except TypeError: raise CTException('Bins should be a list, vector, or numpy array of evenly spaced values') try: bins = np.array(bins, dtype=float) except ValueError: raise CTException('Passed bins could not be converted to a numpy array of floats') # check stride is ok self.__check_stride(stride) # get the residue IDXs were going to use res_idx_list = self.residue_index_list n_frames = self.n_frames # check the fragment_size is appropriate if fragment_size > len(res_idx_list): raise CTException('fragment_size is larger than the number of residues') if fragment_size < 2: raise CTException('fragment_size must be 2 or larger') meanData = [] stdData = [] histo = [] # cycle over each sub-region in the sequence for frag_idx in res_idx_list[0:-fragment_size]: tmp = [] ctio.status_message("On range %i" % frag_idx, verbose) # for each frame in ensemble, calculate RMSD for that sub-region compared to # all other sub-regions (i.e. we're doing a 1-vs-all RMSD calculation for EACH # frame (after adjusting for stride) for a subregion of the protein for j in range(0, n_frames, stride): tmp.extend(self.get_RMSD(j ,-1, region=[frag_idx, frag_idx+fragment_size])) # compute a histogram for this large dataset (b,c) = np.histogram(tmp,bins) histo.append(b) meanData.append(np.mean(tmp)) stdData.append(np.std(tmp)) return (meanData, stdData, histo, bins)
# ........................................................................ #
[docs] def get_D_vector(self, stride=20, verbose=True): """ Function to calculate the vector of D values used to calculate the Phi parameter from Lyle et al[1]. The stride parameter 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. Parameters ------------ 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 verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns --------- np.ndarray Returns a numpy array of D values (i.e. the D_vector) References -------------- [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 the list of residues which have CA (typically this means we exlcude # ACE and NME) residuesWithCA = self.resid_with_CA # compute number of frames exactly (this is empyrical but ensures we're always consistent with the projected # dimensions in the first for-loop) tmp = self.calculate_all_CA_distances(residuesWithCA[0], stride=stride, onlyCterminalResidues=True, correctOffset=False) n_frames = np.shape(tmp)[0] all_distances = np.zeros([len(residuesWithCA),len(residuesWithCA), n_frames]) # first compute upper triangle only (lower traingle is identical and doesn't change the answer # so we stick with the upper traingle only) SM_index=0 for resIndex in residuesWithCA[0:-1]: ctio.status_message("Calculating non redundant distance for res. %i " % resIndex, verbose) vals = self.calculate_all_CA_distances(resIndex, stride=stride, onlyCterminalResidues=True, correctOffset=False) # have to include a -1 here because we don't have a self:self distance all_distances[SM_index][0:(len(residuesWithCA)-1)-SM_index] = vals.transpose() SM_index=SM_index+1 # number of residues we're calculating distances between n_res = np.shape(all_distances)[0] # reshape to n_frame x n_rij 2D array all_distance_tmp = all_distances.transpose().reshape(n_frames,n_res*n_res) # find the idx of non-zero elements in the first frame (will be true for all frames - zeros # originate because we only computed the upper triangle, so 1/2 of elements in each row of # all_distance_tmp are zero non_zero_idx = np.nonzero(all_distance_tmp[0])[0] # finally extract out the positions that are nonzero, setting us up for the phi analysis non_zero_distance = all_distance_tmp[:,non_zero_idx] # calculate the D-vector of all frames D_vector = [] for A in range(0, n_frames): ctio.status_message("Running PHI calculation on frame %i of %i" % (A, n_frames), verbose) for B in range(A+1, n_frames): """ # This is the old less efficient implementation of the algorithm, kept in case, # for some reason, the new implementation has issues... # get the vector of distances for frame A and frame B - this extracts all the # inter-residue distances for frame A and frame B VA = all_distances.transpose()[A].flatten() VB = all_distances.transpose()[B].flatten() # remove zero entries (all real distances MUST be greater # than zero because there's a hardsphere potential stopping things actually # reaching zero distance (distance is center-of-mass calculated for Ca) # we need to calculate the allowed set NZ_A = np.nonzero(VA)[0] NZ_B = np.nonzero(VB)[0] # this finds the positions that in both vectors were non-zero all_zero = np.intersect1d(NZ_A, NZ_B) # now for those positions extract a new length-match vector for frame # A and B that contains only non-zero positions VA = VA[all_zero] VB = VB[all_zero] """ VA = non_zero_distance[A] VB = non_zero_distance[B] # and compute the D value for comparing these two frames D_vector.append(1 - np.dot(VA,VB)/(np.linalg.norm(VA)*np.linalg.norm(VB))) return np.array(D_vector)
# ........................................................................ #
[docs] def get_RMSD(self, frame1, frame2=-1, region=None, backbone=True, correctOffset=True, stride=1): """ 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. Parameter -------------- 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 region : list/tuple of length 2 {None} Defines the first and last residue (INCLUSIVE) for a region to be examined. By default is set to None which means the entire protein is used 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 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 Returns ---------- np.ndarray Returns a numpy array of either 1 (if two frames are passed) or nframes length which corresponds to the RMSD difference either between two frames or between one frame and ALL other frames """ # get the selection atoms (perform correction if required) selectionatoms = self.__get_selection_atoms(region=region, backbone=backbone, correctOffset=correctOffset) # set the reference trajectory we're working with ref = self.traj # if a second frame number was provided with which we're going to work with if frame2 > -1 and isinstance( frame2, int ): # our target is now a single (i.e. doing RMSD of two structures) target = self.traj.slice(frame2) else: # else we're going to carry out an RMSD comparison between *every* stride-th # frame and frame 1 target = self.__get_subtrajectory(self.traj, stride) # return the RMSD comparison in Angstroms return 10*md.rmsd(target, ref, frame1, atom_indices=selectionatoms)
# ........................................................................ #
[docs] def get_Q(self, protein_average = True, region = None, beta_const = 50.0, lambda_const = 1.8, native_contact_threshold = 4.5, correctOffset = True, stride = 1, 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 Parameters ----------- 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. By default is set to None which means the entire protein is used beta_const : float {50} Constant used for computing Q in reciprocal nanometers. Default is 50 and probably should not be changed without good reason. lambda_const : float {1.8} Constant value is 1.8 for all-atom simulations. Probably should not be changed without good reason native_contact_threshold : float {4.5} Threshold in Angstroms typically used for all-atom simulations and again probably should not be changed without good reason 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 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 {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. Returns ----------- If protein_average=True a single vector is returned with the overall protein average fraction of native contacts associated with each frame for each residue. 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 """ # SET native_state_frame = 0 n_res = self.n_residues # less stringent weights test cos trajectory is one frame too long because we probably loaded the PDB file # as a frame #weights = self.__check_weights(weights, stride) if weights is not False and stride != 1: raise CTException('For get_Q() weights must be set for EACH frame and stride=1') # if we're using a subregion # NOTE this is WAY more elegant than the previous way of doing this but there *used* to be problems with MDTraj doing # things like this... selectionatoms = self.__get_selection_atoms(region, backbone=False, heavy=True, correctOffset=correctOffset) # extract out the native state frame native = self.traj.slice(native_state_frame) # get the sub-trajectory to be used target = self.__get_subtrajectory(self.traj, stride) # now align the entire trajectory to the 'native' frame target.superpose(target, frame=native_state_frame, atom_indices=selectionatoms) try: BETA_CONST = float(beta_const) # in reciprocal nm (1/nm) LAMBDA_CONST = float(lambda_const) # For all-atom simulations # Native contact threshold distance in nm (not param is passed in Angstroms but the calculation # happens expecting nanometers so have to update (hence divide by 10) NATIVE_CUTOFF = float(native_contact_threshold)/10 except ValueError as e: raise CTException('Could not convert constant into float for setting constants in get_Q().\nSee below:\n\n%s' % (str(e))) # use all pairs of atoms that are over 3 away in sequence space heavy_pairs = np.array( [(i,j) for (i,j) in combinations(selectionatoms, 2) if abs(native.topology.atom(i).residue.index - \ native.topology.atom(j).residue.index) > 3]) # compute the distances between these pairs in the native state heavy_pairs_distances = md.compute_distances(native[0], heavy_pairs)[0] # and get the pairs s.t. the distance is less than NATIVE_CUTOFF. This returns the # set of interatomic residues that define the native contacts native_contacts = heavy_pairs[heavy_pairs_distances < NATIVE_CUTOFF] # now compute these distances for the whole trajectory r = md.compute_distances(target, native_contacts) # and recompute them for just the native state r0 = md.compute_distances(native[0], native_contacts) # If we're just computing the protein average then this returns the Q value for the whole protein on a per-frame basis if protein_average: if weights is not False: raise CTException('Reweighting for frame averaged should be done with trajectory weights OUTSIDE of CTraj') q = np.mean(1.0 / (1 + np.exp(BETA_CONST * (r - LAMBDA_CONST * r0))), axis=1) return q else: # if the analysis is to be re-weighted uses the weights here on a per-frame basis if weights is not False: q_full = (1.0 / (1 + np.exp(BETA_CONST * (r - LAMBDA_CONST * r0)))).transpose() q = [] for i in q_full: # note i[1:] means we ignore the first (native) frame q.append(np.average(i[1:], 0, weights)) q = np.array(q) else: # check this makes sense - aboe we do i[1:] should probably correct this to remove # the native state structure? Anyway, here we're averaging over every from for each # residue q = np.mean(1.0 / (1 + np.exp(BETA_CONST * (r - LAMBDA_CONST * r0))), axis=0) # get the set of unqiue atoms which are involved in native contacts unique_native_contact_atoms = np.unique(np.hstack((np.transpose(native_contacts)[0],np.transpose(native_contacts)[1]))) res2at = {} res2res = {} # for each unqiue atom for atom in unique_native_contact_atoms: # determine the name of the residue it's from and update the res2at dictionary # if needed local_res = str(self.topology.atom(atom).residue) if local_res not in res2at: res2at[local_res] = [] res2res[int(local_res[3:])] = local_res # now for every native contact pair that atom is involved in, # associate the fraction of the time it's native with the # residue in question. for pair_idx in range(0, len(native_contacts)): if atom in native_contacts[pair_idx]: res2at[local_res].append(q[pair_idx]) # we now have a dictionary where, for each residue, we have the # ensemble average fraction of the simulation each atom was making # native contacts. Note different residues have different numbers # of atoms (obviously...) SO each entry in res2at is going to be # a variable length # now construct an n-res by n-rex empty matrix res_res_matrix = np.zeros((self.n_residues, self.n_residues)) res_res_matrix_count = np.zeros((n_res, n_res)) # and for each pair of atoms as identified previously, look up # which residue they're from, determine the q score for that pairwise # interaction, and increment the associated positions on the nres by # nres matrix, keeping count of how many such increments we make in the # res_res_matrix_count matrix for pair_idx in range(0,len(native_contacts)): pair = native_contacts[pair_idx] R1 = self.topology.atom(pair[0]).residue.index R2 = self.topology.atom(pair[1]).residue.index res_res_matrix[R1,R2] = q[pair_idx] + res_res_matrix[R1,R2] res_res_matrix[R2,R1] = q[pair_idx] + res_res_matrix[R2,R1] res_res_matrix_count[R1,R2] = 1 + res_res_matrix_count[R1,R2] res_res_matrix_count[R2,R1] = 1 + res_res_matrix_count[R2,R1] # finaly, pairwise division accounts for the fact that some residues # have more atoms than others normalized_res_matrix = res_res_matrix / res_res_matrix_count # just as a convenience, build a sorted list of the residues which makes # the data a bit easier to play with going forward. res2res_keys = list(res2res.keys()) np.sort(res2res_keys) sorted_residues = [] for lk in res2res_keys: sorted_residues.append(res2res[lk]) # finally, return all the things mentioned in the function description (note we nan-to-num # to remove all the NaNs from the norlaized res matrix (generated by dividiving by zeor) return (q, native_contacts, res2at, sorted_residues, np.nan_to_num(normalized_res_matrix,0))
# ........................................................................ # #
[docs] def get_contact_map(self, distance_thresh=5.0, mode='closest-heavy', 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). Parameters ----------------- distance_thresh : float {5.0} Distance threshold used to define a 'contact' in Angstroms. Contacts are taken as frames in which the atoms defined by the scheme are within $distance_thresh angstroms of one another mode : string {'closest-heavy'} Mode allows the user to define differnet modes for computing contacts. The default value is 'closest-heavy'. 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. 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. 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. Returns: --------------- tuple of size 2 Returns a tuple where: 0 - contact map 1 - contact order """ ctutils.validate_keyword_option(mode, ['closest-heavy', 'ca', 'closest', 'sidechain', 'sidechain-heavy'] , 'mode') if weights is not False: if int(stride) != 1: raise CTException("Cannot accomodate weighst and non-one stride") # check weights are correct weights = self.__check_weights(weights, stride) # set the distance threshold to a value in nm (we use A by default) distance_thresh_in_nm = float(distance_thresh/10.0) # build a substractectory based on the stride argument subtraj = self.__get_subtrajectory(self.traj, stride) # ensure we only select main chain atoms (no termini) - NOTE, this is a REALLY useful design pattern - # should consider re-writing the code to use this... mainchain_atoms = self.topology.select('(not resname NME) and (not resname ACE)') # compute the contactmap and square-form it (map per frame) # CMAP is a [N_FRAMES x N_RES x N_RES] array CMAP_nonsquare = md.compute_contacts(subtraj.atom_slice(mainchain_atoms), scheme=mode) CMAP = md.geometry.squareform(CMAP_nonsquare[0], CMAP_nonsquare[1]) # extract the normalization factor used to compute fractional # contacts normalization_factor = np.shape(CMAP)[0] # build a MASK where distance is not zero (i.e. where distances were calculated MASK = (CMAP[0] != 0)*1 # if no weights... if weights is False: # for each frame set true/false if less than threshold, # then convert bools to ints and sum over all frames # and finally normalize by the normalization factor. This # gives us the _normalized_ contact map (i.e. each element is between 0 and 1) normalized_contact_map = (np.sum(1*(CMAP < distance_thresh_in_nm),0)*MASK) / float(normalization_factor) # else, if weights... else: # if we use weights then we multiply each frame's contact map by the weight and sum n_frames = CMAP.shape[0] normalized_contact_map = np.zeros((CMAP.shape[1],CMAP.shape[1])) for fid in range(0,n_frames): normalized_contact_map = normalized_contact_map + (np.ndarray.astype((CMAP[fid]<distance_thresh_in_nm),int)*MASK)*weights[fid] # we can further reduce the dimensionality to ask which residues are most involved in contacts with outher # residues in general (i.e. without caring about what those residues are). This gives us a normalized # contact order. n_res = int(np.shape(normalized_contact_map)[0]) # So this is a kind of funky line, but basically because we don't calculate over distances that are # that are # i to i # i to i-1 # i to i-2 # i to i+1 # i to i+2 # which means for MOST residues the max possible is n_res-5, but for those at the end and start # there is no i-1/i-2 for the 0th residues, so the line below builds a vector that for each residues # calculates the TRUE max fractional contacts contact_order_normalization_vector = n_res - np.hstack((np.hstack(([3,4],np.repeat(5,n_res-4))),[4,3])) normalized_contact_order = np.sum(normalized_contact_map,0)/contact_order_normalization_vector return (normalized_contact_map, normalized_contact_order)
# ........................................................................ # #
[docs] def get_clusters(self, region=None, n_clusters=10, backbone=True, correctOffset=True, stride=20): """ 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 ........................................ OPTIONS ........................................ region [list/tuple of length 2] {[]} Defines the first and last residue (INCLUSIVE) for a region to be examined 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. 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. """ # build an empty distance matrix if self.n_frames % stride == 0: distance_dims = int(self.n_frames/stride) else: distance_dims = int((self.n_frames/stride))+1 distances = np.zeros((distance_dims, distance_dims)) idx=0 # build an all vs. all RMSD matrix based on the parameters provided for every # stride-th frame for i in range(0,self.n_frames,stride): distances[idx] = self.get_RMSD(i, stride=stride, region=region, backbone=backbone, correctOffset=correctOffset) idx=idx+1 # CLUSTERING # having computed the RMSD distance matrix we do Ward based hierachical clustering # and then separate out into n_clusters # we feed ward a redundant distance matrix # See: http://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.ward.html#scipy.cluster.hierarchy.ward linkage = scipy.cluster.hierarchy.ward(distances) # linkage is the hierachical clustering encoded as a linkage matrix labels = scipy.cluster.hierarchy.fcluster(linkage, t=n_clusters, criterion='maxclust') # get a subtrajectory which corresponds to the trajectory examined # in the all vs. all comparison (i.e. a trajectory made of every stride-th # frame subtraj = self.__get_subtrajectory(self.traj, stride) # if we're looking at a region further extract out ONLY the atoms # associated with that subregion if region is not None: selectionatoms = self.__get_selection_atoms(region, backbone) subtraj = subtraj.atom_slice(selectionatoms) # we now build n_cluster separate trajectories contaning conformations from the clustering cluster_trajs = [] cluster_distance_matricies = [] cluster_centroids = [] cluster_members = [] cluster_frames = [] # regardless of how many clusters we *think* we should have, extract the number of labeles # we'll actually have... final_labels = list(set(labels)) # for each cluster for i in final_labels: # determine the indices associated with frames which are associated # with the i-th cluster IDXs=np.where(labels == i) IDXs=IDXs[0] cluster_frames.append(IDXs) # record how many frames are associated with the i-th cluster cluster_members.append(len(IDXs)) # create the trajectory and append to the cluster_trajectory list cluster_trajs.append(subtraj.slice(IDXs)) # create the appropriate submatrix cluster_distances=np.zeros((len(IDXs), len(IDXs))) # initial distances is a subset of the all vs all RMSD distance # matrix which now only includes rows associated with the frames # from the i-th cluster initial_dist=distances[IDXs] # for each frame associated with the i-th cluster for k in range(0,len(IDXs)): # add the full all vs. all set of distances between each of the frames from # the i-th cluster and all the other frames from the i-th clster cluster_distances[k] = initial_dist[k][IDXs] # at this point the cluster_distances matrix is an all vs. all RMSD distance # matrix for all the frames in i-th cluster - this effectivly gives you a way # to think about how well an RMSD cluster represents those structures cluster_distance_matricies.append(cluster_distances) # we determine the frame closest to the centroid of the cluster cluster_centroids.append(np.exp(-1*cluster_distances / cluster_distances.std()).sum(axis=1).argmax()) return (cluster_members, cluster_trajs, cluster_distance_matricies, cluster_centroids, cluster_frames)
# ........................................................................ # #
[docs] def get_inter_residue_COM_distance(self, R1, R2, correctOffset=True, stride=1): """ 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. ........................................ OPTIONS ........................................ R1 [int] Residue index of first residue R2 [int] Residue index of second residue 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. 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. """ # first check that the residues we're working with make sense if correctOffset: R1 = self.get_offset_residue(R1) R2 = self.get_offset_residue(R2) else: R1 = int(R1) R2 = int(R2) # get COM of the two residues for every stride-th frame (only split if R1 not in self.__residue_COM: atoms1 = self.__residue_atom_lookup(R1) TRJ_1 = self.traj.atom_slice(atoms1) TRJ_1 = self.__get_subtrajectory(TRJ_1, stride) self.__residue_COM[R1] = md.compute_center_of_mass(TRJ_1) if R2 not in self.__residue_COM: atoms2 = self.__residue_atom_lookup(R2) TRJ_2 = self.traj.atom_slice(atoms2) TRJ_2 = self.__get_subtrajectory(TRJ_2, stride) self.__residue_COM[R2] = md.compute_center_of_mass(TRJ_2) COM_1 = self.__residue_COM[R1] COM_2 = self.__residue_COM[R2] # calculate distance # note 10* to get angstroms d = 10*np.sqrt(np.square(np.transpose(COM_1)[0] - np.transpose(COM_2)[0]) + np.square(np.transpose(COM_1)[1] - np.transpose(COM_2)[1])+np.square(np.transpose(COM_1)[2] - np.transpose(COM_2)[2])) # finally fill in the table return d
# ........................................................................ # #
[docs] def get_inter_residue_COM_vector(self, R1, R2, correctOffset=True, stride=1): """ 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_inter_residue_COM_distance) ........................................ OPTIONS ........................................ R1 [int] Residue index of first residue R2 [int] Residue index of second residue 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. 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. """ if correctOffset: R1 = self.get_offset_residue(R1) R2 = self.get_offset_residue(R2) else: R1 = int(R1) R2 = int(R2) TRJ_1 = self.traj.atom_slice(self.topology.select('resid %i' % R1 )) TRJ_1 = TRJ_1.slice(list(range(0, self.n_frames, stride))) TRJ_2 = self.traj.atom_slice(self.topology.select('resid %i'% R2 )) TRJ_2 = TRJ_2.slice(list(range(0, self.n_frames, stride))) COM_1 = md.compute_center_of_mass(TRJ_1) COM_2 = md.compute_center_of_mass(TRJ_2) # note 10* to get Angstroms return (COM_1 - COM_2)
# ........................................................................ # #
[docs] def get_inter_residue_atomic_distance(self, R1, R2, A1='CA', A2='CA', mode='atom', correctOffset=True, stride=1): """ 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. ........................................ OPTIONS ........................................ 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 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. 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. """ # check mode keyword is valid ctutils.validate_keyword_option(mode, ['atom', 'ca', 'closest-heavy', 'closest', 'sidechain', 'sidechain-heavy'] , 'mode') # define R1 and R2 - if offset needed perform, else # define the first and last residue INCLUDING caps if correctOffset: R1 = self.get_offset_residue(R1) R2 = self.get_offset_residue(R2) else: # cast... R1 = int(R1) R2 = int(R2) try: # if atom mode was used if mode == 'atom': if A1 == 'CA' and A2 == 'CA': subtraj = self.__get_subtrajectory(self.traj, stride) distances = 10*md.compute_contacts(subtraj, [[R1,R2]],scheme='ca')[0].ravel() else: atom1 = self.__residue_atom_lookup(R1,A1) if len(atom1) == 0: raise CTException('Unable to find atom [%s] in residue R1 (%i)' % (A1, R1)) TRJ_1 = self.traj.atom_slice(atom1) TRJ_1 = self.__get_subtrajectory(TRJ_1, stride) atom2 = self.__residue_atom_lookup(R2,A2) if len(atom2) == 0: raise CTException('Unable to find atom [%s] in residue R1 (%i)' % (A2, R2)) TRJ_2 = self.traj.atom_slice(atom2) TRJ_2 = self.__get_subtrajectory(TRJ_2, stride) COM_1 = md.compute_center_of_mass(TRJ_1) COM_2 = md.compute_center_of_mass(TRJ_2) distances = 10*np.sqrt(np.square(np.transpose(COM_1)[0] - np.transpose(COM_2)[0]) + np.square(np.transpose(COM_1)[1] - np.transpose(COM_2)[1])+np.square(np.transpose(COM_1)[2] - np.transpose(COM_2)[2])) except IndexError as e: ctio.exception_message("This is likely because one of [%s] or [%s] is not a valid atom type for the residue in question. Full error printed below\n%s" %( A1,A2, str(e)), e, with_frills=True) # parse any of the allowed modes in compute_contacts (see http://mdtraj.org/1.8.0/api/generated/mdtraj.compute_contacts.html # for more details!) if mode == 'closest' or mode == 'ca' or mode == 'closest-heavy' or mode == 'sidechain' or mode == 'sidechain-heavy': subtraj = self.__get_subtrajectory(self.traj, stride) try: distances = 10*md.compute_contacts(subtraj, [[R1,R2]], scheme=mode)[0].ravel() except ValueError as e: raise CTException('Your current version of mdtraj does not support [%s] - please update mdtraj to 1.8.0 or later to facilitate support. Alternatively this may be because residue %i or %i is not parsed correctly by MDTraj. Original error:\n%s' % (mode(), R1, R2, str(e))) # note 10* to get Angstroms return distances
# ........................................................................ # #
[docs] def get_residue_mass(self, R1, correctOffset=True): """ Returns the mass associated with a specific residue. ........................................ OPTIONS ........................................ 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. """ if correctOffset: R1 = self.get_offset_residue(R1) else: R1 = int(R1) # get the atoms associated with the resite of interest res_atoms = self.topology.select('resid %i'%(R1)) totalMass=0 for atom in res_atoms: totalMass = totalMass + self.topology.atom(atom).element.mass return totalMass
# ........................................................................ # #
[docs] def get_asphericity(self, R1=None, R2=None, correctOffset=True, verbose=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). ........................................ OPTIONS ........................................ R1 [int] {None} Index value for first residue in the region of interest. If not provided (False) then first residue is used. R1 [int] {None} 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. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! """ # define R1 and R2 - if offset needed perform, else # define the first and last residue INCLUDING caps if R1 is None: R1 = 0 else: if correctOffset: R1 = self.get_offset_residue(R1) else: # cast... R1 = int(R1) if R2 is None: R2 = self.n_residues - 1 else: if correctOffset: R2 = self.get_offset_residue(R2) else: R2 = int(R2) # compute the gyration tensor gyration_tensor_vector = self.get_gyration_tensor(R1, R2, correctOffset=correctOffset, verbose=verbose) asph_vector = [] for gyr in gyration_tensor_vector: # calculate the eigenvalues of the gyration tensor! (EIG, norm) = LA.eig(gyr) # finally calculate the instantanous asphericity and append to the growing vector asph = 1 - 3*((EIG[0]*EIG[1] + EIG[1]*EIG[2] + EIG[2]*EIG[0])/np.power(EIG[0]+EIG[1]+EIG[2],2)) asph_vector.append(asph) return np.array(asph_vector)
# ........................................................................ # #
[docs] def get_gyration_tensor(self, R1=None, R2=None, correctOffset=True, verbose=True): """ Returns the instantaneous gyration tensor associated with each frame. Parameters --------------- R1 : int {None} Index value for first residue in the region of interest. If not provided (False) then first residue is used. R1 : int {None} 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 verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns ----------- np.ndarray Returns a numpy array where each position is the frame-specific gyration tensor value """ # define R1 and R2 - if offset needed perform, else # define the first and last residue INCLUDING caps if R1 is None: R1 = 0 else: if correctOffset: R1 = self.get_offset_residue(R1) else: # cast... R1 = int(R1) if R2 is None: R2 = self.n_residues - 1 else: if correctOffset: R2 = self.get_offset_residue(R2) else: R2 = int(R2) # switch em around so the A to B syntax makes sense if R1 > R2: tmp = R2 R2 = R1 R1 = tmp all_positions_all_frames = self.traj.atom_slice(self.topology.select('resid %i to %i'%(R1,R2))) gyration_tensor_vector = [] count = 1 for frame in all_positions_all_frames: # quick status update... if count % 500 == 0: ctio.status_message("On frame %i of %i [computing gyration tensor]" % (count, len(all_positions_all_frames)), verbose) count = count + 1 # compute the center of mass for the relevant atoms COM = md.compute_center_of_mass(frame) # calculate the COM to position difference matrix DIF = frame.xyz[0] - COM # old slow method """ T_PRE = 0.0 # for each atom in the frame calculate the outer product (dyadic product) of # the difference between the position at the overall center of mass. This creates a # 3x3 gyration tensor for pos in frame.xyz[0]: T_PRE = T_PRE + np.outer(pos - COM, pos - COM) T = T_PRE/len(frame.xyz[0]) """ # compute the gyration tensor - this syntax is WAY faster than the above T_new = np.sum(np.einsum('ij...,i...->ij...',DIF,DIF),axis=0)/len(frame.xyz[0]) gyration_tensor_vector.append(T_new) return np.array(gyration_tensor_vector)
# ........................................................................ # #
[docs] def get_end_to_end_distance(self, 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 Parameters --------------- mode [string] {'CA'} String, must be one of either 'CA' or 'COM'. 'CA' = alpha carbon 'COM' = center of mass (associated withe the residue) Returns --------- """ ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') # extract first and last residue AFTER offset correction (i.e. residue number from original PDB file) start = self.resid_with_CA[0] end = self.resid_with_CA[-1] if mode == 'CA': distance = self.get_inter_residue_atomic_distance(start, end, stride=1, correctOffset=False) elif mode == 'COM': distance = self.get_inter_residue_COM_distance(start, end, stride=1, correctOffset=False) return distance
# ........................................................................ # #
[docs] def get_radius_of_gyration(self, R1=None, R2=None, 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. Parameters --------------- R1 : int {None} Index value for first residue in the region of interest. If not provided (None) then first residue is used. R2 : int {None} 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. Returns ----------- np.ndarray Returns a numpy array with per-frame instantaneous radius of gyration """ # define R1 and R2 - if offset needed perform, else # define the first and last residue INCLUDING caps if R1 is None: R1 = 0 else: R1 = int(R1) if R2 is None: R2 = self.n_residues-1 else: R2 = int(R2) # note we correct offset regardless of if R1/R2 are passed or not if correctOffset: R1 = self.get_offset_residue(R1) R2 = self.get_offset_residue(R2) # switch em around so the A to B syntax makes sense if R1 > R2: tmp = R2 R2 = R1 R1 = tmp # in angstroms return 10*md.compute_rg(self.traj.atom_slice(self.topology.select('resid %i to %i'%(R1, R2))))
# ........................................................................ # #
[docs] def get_hydrodynamic_radius(self, R1=None, R2=None, alpha1=0.216, alpha2=4.06, alpha3=0.821, correctOffset=True): """ Returns the apparent hydrodynamic radius as calculated based on the approximation derived by Nygaard et al. [1]. Returns a hydrodynamic radius in Angstroms. Parameters (alpha1/2/3 should not be altered to recapitulate behaviour defined by Nygaard et al. References: ------------- [1] Nygaard M, Kragelund BB, Papaleo E, Lindorff-Larsen K. An Efficient Method for Estimating the Hydrodynamic Radius of Disordered Protein Conformations. Biophys J. 2017;113: 550–557. Radius of gyration is returned in Angstroms. Parameters --------------- R1 : int {False} Index value for first residue in the region of interest. If not provided (False) then first residue is used. R2 : int {False} Index value for last residue in the region of interest. If not provided (False) then last residue is used. alpha1 : float {0.216} First parameter in equation (7) from Nygaard et al. alpha2 : float {4.06} Second parameter in equation (7) from Nygaard et al. alpha3 : float {0.821} Third parameter in equation (7) from Nygaard et al. 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. Returns ----------- np.ndarray Returns a numpy array with per-frame instantaneous hydrodynamic radii """ # first compute the rg rg = self.get_radius_of_gyration(R1, R2, correctOffset) # precompute N_033 = np.power(self.n_residues, 0.33) N_060 = np.power(self.n_residues, 0.60) Rg_over_Rh = ((alpha1*(rg - alpha2*N_033)) / (N_060 - N_033)) + alpha3 return (1/Rg_over_Rh)*rg
# ........................................................................ # #
[docs] def get_t(self, R1=None, R2=None, correctOffset=True): """ Returns the <t>, a dimensionless parameter which describes the size of the ensemble. Parameters --------------- R1 : int {None} Index value for first residue in the region of interest. If not provided (False) then first residue is used. R2 : int {None} 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. Returns ----------- np.ndarray Returns a numpy array with per-frame instantaneous t-values """ # define R1 and R2 - if offset needed perform, else # define the first and last residue INCLUDING caps if R1 is None: R1 = 0 else: R1 = int(R1) if R2 is None: R2 = self.n_residues-1 else: R2 = int(R2) # first get the instantanoues RG rg = self.get_radius_of_gyration(R1, R2, correctOffset) n_res = self.n_residues c_length = n_res * 3.6 # next define the exponent exponent = 4.0/(np.power(n_res,0.3333)) # define a function which returns the instantaneous t value for # a given rg def inst_t(i): return 2.5*np.power((1.75*(i/(c_length))),exponent) # compile into a vectorized version K = np.vectorize(inst_t) # run over all rg return K(rg)
# ........................................................................ # #
[docs] def get_internal_scaling(self, R1=None, R2=None, mode='COM', mean_vals=False, correctOffset=True, stride=1, weights=False, verbose=True): """ 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 how this is calculated in CAMPARI, which includes the caps, if present. ........................................ OPTIONS ........................................ 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. 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... 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. 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. 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. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! """ if weights is not False: if int(stride) != 1: raise CTException("For get_scaling_exponent with weights stride MUST be set to 1. If this is a HUGE deal for you please contact alex and he'll try and update the code to accomodate this, but for now we suggest creating a sub-sampled trajectory and loading that") # check weights are correct weights = self.__check_weights(weights, stride) # check stride is ok self.__check_stride(stride) # check mode is OK ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') # process the R1/R2 to set the position after offset correction out = self.__get_first_and_last(R1, R2, withCA = True) R1 = out[0] R2 = out[1] # check the R1/R2 positions contina a CA (not R1 and R2 here are after # offset has been applied. This throws an exception if no CA self.__check_contains_CA(R1) self.__check_contains_CA(R2) max_seq_sep = (R2 - R1) + 1 # if chain is too short... if max_seq_sep < 1: return ([], []) seq_sep_distances = [] seq_sep_vals = [] for seq_sep in range(0, max_seq_sep): ctio.status_message("Internal Scaling - on sequence separation %i of %i" %(seq_sep, max_seq_sep-1), verbose) tmp = [] seq_sep_vals.append(seq_sep) for pos in range(0, max_seq_sep-seq_sep): # define the two positions A = R1 + pos B = R1 + pos+seq_sep # get the distance for every stride-th frame between those two positions using either the CA # mode or the COM mode if mode == 'CA': distance = self.get_inter_residue_atomic_distance(A, B, stride=stride, correctOffset=False) elif mode == 'COM': distance = self.get_inter_residue_COM_distance(A, B, stride=stride, correctOffset=False) # if weights were provided subsample from the set of distances using the weights vector if weights is not False: distance = choice(distance, len(distance), p=weights) tmp = np.concatenate((tmp,distance)) seq_sep_distances.append(tmp) if mean_vals: mean_is = [np.mean(i) for i in seq_sep_distances] return (seq_sep_vals, mean_is) else: return (seq_sep_vals, seq_sep_distances)
# ........................................................................ # #
[docs] def get_internal_scaling_RMS(self, R1=None, R2=None, mode='COM', stride=1, correctOffset=True, weights=False, verbose=True): """ If :math:`r_{i,j} = \langle \langle \sum \sigma_{1}` equals :math:`\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. ........................................ OPTIONS ........................................ R1 [int] {None} Index value for first residue in the region of interest. If not provided (False) then first residue is used. R1 [int] {None} 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. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! """ # compute the non RMS internal scaling behaviour (seq_sep_vals, seq_sep_distances) = self.get_internal_scaling(R1=R1, R2=R2, mode=mode, stride=stride, correctOffset=correctOffset, weights=weights, verbose=verbose) # calculate RMS for each distance mean_is = [np.sqrt(np.mean(i*i)) for i in seq_sep_distances] return (seq_sep_vals, mean_is)
# ........................................................................ #
[docs] def get_scaling_exponent(self, inter_residue_min=15, end_effect=5, correctOffset=True, subdivision_batch_size=20, mode='COM', num_fitting_points=40, fraction_of_points=0.5, fraction_override=False, stride=1, weights=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() ........................................ OPTIONS ........................................ 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. 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. 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. 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. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! """ # check weights are OK weights = self.__check_weights(weights, stride) # check mode is OK ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') # compute max |i-j| distance being used... first = self.resid_with_CA[0] last= self.resid_with_CA[-1] max_separation = (last-first)+1 # if we're not using fraction override check the number of points requested makes sense given sequence length if not fraction_override and (max_separation - (end_effect+inter_residue_min)) < num_fitting_points: fraction_of_points=1.0 ctio.warning_message("For scaling exponent calculation, sequence not long enough to use %i points (only %i valid positions once end effects and low |i-j| are accounted for), switching to using the fraction of points mode (will use %i points instead)" % (num_fitting_points, (max_separation - (end_effect+inter_residue_min)), int(fraction_of_points*(max_separation - (end_effect+inter_residue_min))))) fraction_override = True if fraction_override: if fraction_of_points > 1.0: raise CTException("Using fraction_overide to define the number of points to fit in the linear loglog analysis, but requested over 1.0 fraction (fraction_of_points must lie between >=0 and 1.0") # again note int to round down here num_fitting_points = int(fraction_of_points*(max_separation - (end_effect+inter_residue_min))) if num_fitting_points < 3: raise CTException("Less than three points - cannot fit a straight line") if num_fitting_points < 10: ctio.warning_message("Warning: Scaling fit has only %i points - likely finite size effects!" % (num_fitting_points)) # check weights make sense... weights = self.__check_weights(weights, stride) # check that the number of points being used makes sense.. # This section determines the number of subdivisions performed for error # bootstrapping. If we have fewer frames than we can divide the data into # then we just use each frame individually (although now error bootstrapping # is probably meaningless! # note integer math used here to round down - also set if int(self.n_frames/stride) < int(subdivision_batch_size): num_subdivisions_for_error = int(self.n_frames/stride) else: num_subdivisions_for_error = int(int(self.n_frames/stride) / subdivision_batch_size) seq_sep_vals = [] seq_sep_RMS_distance = [] seq_sep_RMS_var_distance = [] seq_sep_RSTDS_distance = [] seq_sep_subsampled_distances = [] # for each possible sequence separation (|i-j| value) for seq_sep in range(1, max_separation): ctio.status_message("Internal Scaling - on sequence separation %i of %i" %(seq_sep, max_separation), verbose) tmp = [] seq_sep_vals.append(seq_sep) # collect all possible average seq-sep values weighted by the weights for pos in range(0, max_separation-seq_sep): # define the two positions A = first + pos B = first + pos + seq_sep if mode == 'CA': distance = self.get_inter_residue_atomic_distance(A, B, stride=stride, correctOffset=False) elif mode == 'COM': distance = self.get_inter_residue_COM_distance(A, B, stride=stride, correctOffset=False) # compute the ensemble average of the distances if weights is not False: distance = choice(distance, len(distance), p=weights) tmp.extend(distance) tmp = np.array(tmp) # add mean and std vals for this sequence sep seq_sep_RMS_distance.append(np.sqrt(np.mean(tmp*tmp))) seq_sep_RSTDS_distance.append(np.sqrt(np.std(tmp*tmp))) seq_sep_RMS_var_distance.append(np.sqrt(np.power(np.var(tmp),2))) if num_subdivisions_for_error > 0: # note we cast this to an int to ensure subdivision_size is always # the value added to RMS_local_append is always the same, because # len(tmp) will vary with sequence separation. Basically this means # we take ALL the data and subidivided it into num_subdivisions_for_error # chunks and then use this for error calculations subdivision_size = int(len(tmp)/num_subdivisions_for_error) # get shuffled indices idx = np.random.permutation(list(range(0,len(tmp)))) # split shuffled indices into $num_subdivisions_for_error sized chunks subdivided_idx = cttools.chunks(idx, subdivision_size) # finally subselect each of the randomly selected indicies RMS_local = [] for idx_set in subdivided_idx: # subselect a random set of distances and compute RMS RMS_local.append(np.sqrt(np.mean(tmp[idx_set]*tmp[idx_set]))) # add distribution of values for this sequence sep seq_sep_subsampled_distances.append(RMS_local) # now sub-select the bit of the curve we actually want for the separation, distance, and distance variance data # note we are RE DEFINING these three variables here seq_sep_vals = seq_sep_vals[inter_residue_min:-end_effect] seq_sep_RMS_distance = seq_sep_RMS_distance[inter_residue_min:-end_effect] seq_sep_RMS_var_distance = seq_sep_RMS_var_distance[inter_residue_min:-end_effect] ## next find indices for evenly spaced points in logspace. This whole sectino # leads to the identification of the indices in logspaced_idx, which are the # list indices that will given evenly spaced points when plotted in log space y_data = np.log(seq_sep_vals) y_data_offset = y_data - y_data[0] interval = y_data_offset[-1]/num_fitting_points integer_vals = y_data_offset/interval logspaced_idx = [] for i in range(0,num_fitting_points): [local_ix,_] = cttools.find_nearest(integer_vals, i) if local_ix in logspaced_idx: continue else: logspaced_idx.append(local_ix) # finally using those evenly-spaced log indices we extract out new lists # that have values which will be evenly spaced in logspace. Cool. fitting_separation = [seq_sep_vals[i] for i in logspaced_idx] fitting_distances = [seq_sep_RMS_distance[i] for i in logspaced_idx] fitting_variance = [seq_sep_RMS_var_distance[i] for i in logspaced_idx] # fit to a log/log model and extract params out = np.polyfit(np.log(fitting_separation), np.log(fitting_distances), 1) nu_best = out[0] R0_best = np.exp(out[1]) ## >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ### next calculated reduced chi-squared n_points = len(fitting_distances) chi2=0 for i in range(0, n_points): chi2 = chi2 + (np.power(np.log(fitting_distances[i]) - nu_best*np.log(fitting_separation[i])+R0_best,2))/fitting_variance[i] # finally calculated reduced chi squared correcting for 2 model parameters reduced_chi_squared_fitting = chi2 / (n_points-2) full_n_points = len(seq_sep_vals) chi2=0 for i in range(0, full_n_points): chi2 = chi2 + (np.power(np.log(seq_sep_RMS_distance[i]) - nu_best*np.log(seq_sep_vals[i])+R0_best,2))/seq_sep_RMS_var_distance[i] # finally calculated reduced chi squared correcting for 2 model parameters reduced_chi_squared_all = chi2 / (full_n_points-2) ## >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ### Finally run the subselection protocol to subsampled subselected = np.array(seq_sep_subsampled_distances).transpose() nu_sub = [] R0_sub = [] for i in range(0, num_subdivisions_for_error): local_distances = subselected[i][inter_residue_min:-end_effect] OF = np.polyfit(np.log(fitting_separation), np.log([local_distances[i] for i in logspaced_idx]), 1) nu_sub.append(OF[0]) R0_sub.append(np.exp(OF[1])) if num_subdivisions_for_error < 1: nu_sub.append(np.nan) R0_sub.append(np.nan) return [nu_best, R0_best, min(nu_sub), max(nu_sub), min(R0_sub), max(R0_sub), reduced_chi_squared_fitting, reduced_chi_squared_all, np.vstack((np.array(fitting_separation),np.array(fitting_distances))), np.vstack((seq_sep_vals, seq_sep_RMS_distance, cttools.powermodel(seq_sep_vals, nu_best, R0_best)))]
# ........................................................................ # #
[docs] def get_residue_COM(self, R1, correctOffset=True): """ Returns the COM vector for the residue across the trajectory. ........................................ OPTIONS ........................................ 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. """ # correct the offset if necessary if correctOffset: R1 = self.get_offset_residue(R1) else: R1 = int(R1) # get the atoms associated with the resite of interest return md.compute_center_of_mass(self.traj.atom_slice(self.topology.select('resid %i'%R1)))
# ........................................................................ # #
[docs] def get_all_SASA(self, probe_radius=0.14, mode='residue', stride=20): """ 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! ........................................ OPTIONS ........................................ 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. 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 """ ## >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ## internal function def get_sasa_based_on_type(basis, passed_mode): # get all reas all_areas = basis[0]*100 all_atoms = list(basis[1]) # note these CA resids are after offset correction has been applied. NOTE we need to # check if the atom selection here works for multichain protein systems # TODO CA_res = self.resid_with_CA # for each residue with a CA atom... residue_SASA=[] for i in CA_res: # get the atomic indices if passed_mode == 'sidechain': # for some reason 'sidechain' selection includes the backbone hydrogen atoms??!?! relevant_atom_idx = self.topology.select('resid %i and %s and (not name H HA HA2 HA3)' % (i,passed_mode)) if passed_mode == 'backbone': # for some reason 'backbone' ignores the backbone hydrogen atoms ?!?!? relevant_atom_idx = self.topology.select('(resid %i and %s) or (resid %i and name H HA HA2 HA3)' % (i,passed_mode, i)) # no atoms so create an empty list if len(relevant_atom_idx) == 0: per_res = list(np.zeros(len(all_areas.transpose()[0]))) else: # get SASA index of first atom index... per_res = all_areas.transpose()[all_atoms.index(relevant_atom_idx[0])] for atom in relevant_atom_idx[1:]: per_res = per_res + all_areas.transpose()[all_atoms.index(atom)] residue_SASA.append(per_res) return np.array(residue_SASA).transpose() ## >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> # validate input mode ctutils.validate_keyword_option(mode, ['residue', 'atom','backbone','sidechain','all'], 'mode') # downsample based on the stride target = self.__get_subtrajectory(self.traj, stride) # 100* to convert from nm^2 to A^2 if mode == 'residue': return 100*md.shrake_rupley(target, mode='residue', probe_radius=probe_radius) if mode == 'atom': return 100*md.shrake_rupley(target, mode='atom', probe_radius=probe_radius) if mode == 'sidechain' or mode == 'backbone' or mode == 'all': print("WARNING: Not tested on multiprotein systems") # run calc basis = md.shrake_rupley(target, mode='atom', probe_radius=probe_radius, get_mapping=True) # extract sidechains if mode == 'sidechain' or mode == 'all': SC_SASA = get_sasa_based_on_type(basis, 'sidechain') # extract backbone if mode == 'backbone' or mode == 'all': BB_SASA = get_sasa_based_on_type(basis, 'backbone') if mode == 'all': ALL_SASA = 100*md.shrake_rupley(target, mode='residue', probe_radius=probe_radius) if mode == 'sidechain': return SC_SASA if mode == 'backbone': return SC_SASA if mode == 'all': return (ALL_SASA, SC_SASA, BB_SASA)
# ........................................................................ # #
[docs] def get_site_accessibility(self, input_list, probe_radius=0.14, mode='residue_type', stride=20): """ 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. ........................................ OPTIONS ........................................ 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. 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' 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 """ ## First check mode is valid and then sanity check input ctutils.validate_keyword_option(mode, ['residue_type', 'resid'], 'mode') # if we're in 'residue_type' mode if mode == 'residue_type': # set the list of valid residues to the 20 AAs + the caps for res in input_list: if res not in ALL_VALID_RESIDUE_NAMES: raise CTException("Error: Tried to use get_site_accessibility in 'residue_type' mode but residue %s not found in list of valid residues %s" % (res, str(ALL_VALID_RESIDUE_NAMES))) # if we're in resid mode elif mode.lower() == 'resid': offset_input_idx=[] for res in input_list: R1 = int(res) # check residue makes sense and then calculate the offset self.__check_single_residue(R1) offset_input_idx.append(self.get_offset_residue(R1)) ## if we get here all input has been validated if mode == 'residue_type': # get AA list and convert into an ordered list of each # AA residue type SEQS = self.get_amino_acid_sequence(numbered=False) # initialze the empty idx v offset_input_idx = [] idx=0 for res in SEQS: if res in input_list: offset_input_idx.append(idx) idx=idx+1 # next compute ALL SASA for all residues (need the full protein based SASA ALL_SASA = np.transpose(self.get_all_SASA(stride=stride, probe_radius=probe_radius)) lookup = self.get_amino_acid_sequence() return_data={} for i in offset_input_idx: return_data[lookup[i]] = [np.mean(ALL_SASA[i]), np.std(ALL_SASA[i])] return return_data
# ........................................................................ # #
[docs] def get_regional_SASA(self, R1, R2, probe_radius=0.14, correctOffset=True, stride=20): """ 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! ........................................ OPTIONS ........................................ R1 [int] Index value for the first residue in the region R2 [int] Index value for the last residue in the region 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. 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 """ if correctOffset: R1 = self.get_offset_residue(R1) R2 = self.get_offset_residue(R2) # NOTE - we HAVE to compute SASA over the full ensemble to take into acount # atoms OUTSIDE the region getting in the way of the regional SASA total = self.get_all_SASA(stride=stride, probe_radius=probe_radius) total = np.transpose(total) regional_SASA = 0 for i in range(R1, R2): regional_SASA = regional_SASA + np.mean(total[i]) # return the mean sum of SASA for all atoms return regional_SASA
# ........................................................................ # #
[docs] def get_sidechain_alignment_angle(self, 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. """ # The initial section is responsible for selecting/defining the sidechain atom and catching any # bad inputs. It's a litte drawn out but useful for code clarity reasons. # note need to do this with the un-corrected residue index values, hence why it comes first if sidechain_atom_1 == 'default': resname_1 = self.get_amino_acid_sequence(numbered=False)[R1] resname_1 = cttools.fix_histadine_name(resname_1) try: sidechain_atom_1 = DEFAULT_SIDECHAIN_VECTOR_ATOMS[resname_1] if sidechain_atom_1 == 'ERROR': raise CTException('Residue lacks a valid sidechain (%s)' % resname_1) except KeyError: raise CTException('Cannot parse residue at position %i (residue name = %s) ' % (R1, resname_1)) if sidechain_atom_2 == 'default': resname_2 = self.get_amino_acid_sequence(numbered=False)[R2] resname_2 = cttools.fix_histadine_name(resname_2) try: sidechain_atom_2 = DEFAULT_SIDECHAIN_VECTOR_ATOMS[resname_2] if sidechain_atom_2 == 'ERROR': raise CTException('Residue lacks a valid sidechain (%s)' % resname_2) except KeyError: raise CTException('Cannot parse residue at position %i (residue name = %s) ' % (R2, resname_2)) ### At this point we have reasonable atom names defined! # if correctOffset: R1 = self.get_offset_residue(R1) R2 = self.get_offset_residue(R2) TRJ_1_SC = self.traj.atom_slice(self.topology.select('resid %i and name %s' % (R1, sidechain_atom_1) )) TRJ_1_CA = self.traj.atom_slice(self.topology.select('resid %i and name CA' % (R1) )) TRJ_2_SC = self.traj.atom_slice(self.topology.select('resid %i and name %s' % (R2, sidechain_atom_2) )) TRJ_2_CA = self.traj.atom_slice(self.topology.select('resid %i and name CA' % (R2) )) # compute CA-SC vector R1_vector = TRJ_1_SC.xyz - TRJ_1_CA.xyz R2_vector = TRJ_2_SC.xyz - TRJ_2_CA.xyz # finally compute the alignment for each frame and return # a vector of alignments nframes = R1_vector.shape[0] alignment=[] for i in range(0, nframes): # convert to unit vector V1 = R1_vector[i][0] / np.linalg.norm(R1_vector[i][0]) V2 = R2_vector[i][0] / np.linalg.norm(R2_vector[i][0]) # compute dot product and take arccos and do rad->deg to get angle # between the unit vectors alignment.append(np.rad2deg(np.arccos(np.dot(V1,V2)))) return np.array(alignment)
# ........................................................................ # #
[docs] def get_dihedral_mutual_information(self, angle_name='psi', bwidth = np.pi/5.0, stride=1, weights=False): """ 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. ........................................ OPTIONS ........................................ angle_name String, must be one of the following options 'chi1','phi, 'psi', 'omega'. 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... 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. 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. """ ## .................................................. ## SAFETY FIRST! ## # verify binwidth input values if bwidth > 2*np.pi or not (bwidth > 0): raise CTException('The bwidth parameter must be between 2*pi and greater than 0') # if stride was passed make sure it's ok self.__check_stride(stride) # if weights were passed make sure they're LEGIT! weights = self.__check_weights(weights, stride) # check ctutils.validate_keyword_option(angle_name, ['chi1', 'phi', 'psi', 'omega'], 'angle_name') ## .................................................. # define histogram bins based on passed bin width bins=np.arange(-np.pi, np.pi+bwidth, bwidth) # construct the selector dictionary selector = {"phi":md.compute_phi, "omega":md.compute_omega, "psi":md.compute_psi, "chi1":md.compute_chi1} # check the angle_name is an allowed name if angle_name not in list(selector.keys()): raise CTException('The variable angle_name was set to %s, which is not one of phi, omega, psi, chi1' % angle_name) # select and compute the relevant angles of the subtrajectroy fx = selector[angle_name] angles = fx(self.traj[0::stride]) # construct empty matrices SIZE = len(angles[0]) MI_mat = np.zeros((SIZE,SIZE)) # populate matrices note we only compute the upper right triangle # but can populate the lower half because it's a symmetrical matrix for i in range(0,SIZE): for j in range(i,SIZE): X = np.transpose(angles[1])[j] Y = np.transpose(angles[1])[i] MI = ctmutualinformation.calc_MI(X,Y, bins, weights) MI_mat[i,j] = MI MI_mat[j,i] = MI return MI_mat
# ........................................................................ # #
[docs] def get_local_to_global_correlation(self, mode='COM', n_cycles=100, max_num_pairs=10, stride=20, weights=False, verbose=True): """ 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. Parameters ---------- mode : str {'COM'} Must be one of either 'CA' or 'COM'. - 'CA' = alpha carbon. - 'COM' = center of mass (associated withe the residue). 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. 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. weights : bool {False} Flag that indicates if frame weights should be used or not. verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns ------- tuple 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. Raises ------ CTException Raised when the mode is not 'CA' or 'COM'; or, when the lengths of the Rg-stride derived calculations did not match the lengths of the internal distances; or, when the installed numpy version doesn't support the `aweights` keyword for the `numpy.cov` function. """ weights = self.__check_weights(weights, stride) ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') # note start and end here are after offset correction start = self.resid_with_CA[0] end = self.resid_with_CA[-1] FIRST_CHECK = True # start with a sequence separation of 1 for seq_sep in range(1, self.n_residues): ctio.status_message("On sequence separation %i" % (seq_sep), verbose) for pos in range(start, end - seq_sep): # define the two positions A = pos B = pos + seq_sep # get the distance for every stride-th frame between those two positions using either the CA # mode or the COM mode if mode == 'CA': distance = self.get_inter_residue_atomic_distance(A, B, stride=stride, correctOffset=False) elif mode == 'COM': distance = self.get_inter_residue_COM_distance(A, B, stride=stride, correctOffset=False) if FIRST_CHECK: all_distances = distance FIRST_CHECK = False else: all_distances = np.vstack((all_distances, distance)) full_rg = self.get_radius_of_gyration() stride_rg = full_rg[0::stride] if len(stride_rg) != len(all_distances[0]): raise CTException('Something when wrong when comparing stride-derived Rg and internal distances, this is a bug in the code...)') # total number of distance pairs n_pairs = len(all_distances) pair_selection_vector = np.arange(1,max_num_pairs,1) return_data = np.zeros((len(pair_selection_vector)*n_cycles,2)) weights = False weights = np.repeat(1.0/len(stride_rg),len(stride_rg)) idx=0 for n_selected in pair_selection_vector: ctio.status_message("On %i pairs selected" % n_selected, verbose) for i in range(0,n_cycles): # select n_select different values between 0 and n_pairs idx_selection = np.random.randint(0, n_pairs, n_selected) # this next line means that for each idx_selection value (i.e. for each pair of distances) # we calculate the squared value of each distance, then for each conformation SUM all the pairs # and then for each of those n_frame summs divide by 2 * n_selected squared # this gives a series of values for EACH conformation (so local mean is a 1xn_configurations vector) local_mean_square = np.sum(np.power(all_distances[idx_selection],2),0)/(2*(n_selected*n_selected)) # correlation between local mean square and rg_square - returns a SINGLE value # THIS is the covariance matrix, but we can pass a weighting factor to this making it better suited for # weighted ensembles if weights is not False: try: cov_matrix = np.cov(np.vstack((local_mean_square, np.power(stride_rg,2))), ddof=0, aweights=weights) except TypeError: # this probaby happend because the version of numpy doesn't support aweights. If this is the case try again # without weights assuming weights are equal if len(set(weights)) == 1: cov_matrix = np.cov(np.vstack((local_mean_square, np.power(stride_rg,2)))) else: raise CTException('Weights being passed to get_global_from_local but the current version of numpy (%s) may not support weights via the "aweights" keyword...' % np.version.full_version) # this computes the upper right square from the correlation matrix from the covariance matrix using (Rij = (Cij/(sqrt(Cii - Cjj))) where i and j are # 0 and 1 respectively c = cov_matrix[0,1]/(np.sqrt(cov_matrix[0,0]*cov_matrix[1,1])) else: # this is the correlation matrix and was the old way of doing this c = np.corrcoef(local_mean_square, np.power(stride_rg,2))[0][1] return_data[idx] = [n_selected, c] idx=idx+1 # leaving this here incase we want to re-introduce the 2D histogram information in later versions... # np.histogram2d(np.transpose(return_data)[0],np.transpose(return_data)[1], bins=[np.arange(1,n_pairs), np.arange(0,1,0.01)])) return (return_data, pair_selection_vector, np.mean(np.reshape(np.transpose(return_data)[1], (len(pair_selection_vector),n_cycles)),1), np.std(np.reshape(np.transpose(return_data)[1], (len(pair_selection_vector),n_cycles)),1))
# ........................................................................ # #
[docs] def get_end_to_end_vs_rg_correlation(self, mode='COM'): """ Computes the correlation between Rg^2 and end-to-end^2. Parameters --------- mode: str {'CA'} String, must be one of either 'CA' or 'COM'. - 'CA' = alpha carbon. - 'COM' = center of mass (associated withe the residue). Returns ------- float A single float describing the correlation (as calculated by np.corrcoef). """ # validate the keyword ctutils.validate_keyword_option(mode, ['CA', 'COM'], 'mode') # get end-to-end distance distance = self.get_end_to_end_distance(mode) # get radius of gyration full_rg = self.get_radius_of_gyration() # NOTE this is the same as Re^2 = RG^2/ 6 (standard Debye result for a Gaussian chain) - basically asking # about fractality more than a specific chain model, because examining correlation the scalar factor doesn't # matter , ie really this is Rg^2 vs Re^2 - scalar is irrelevant. But, this approach is consistent with # the approach in the get_local_to_global_correlation() function local_mean_square = np.power(distance,2)/(2) c = np.corrcoef(local_mean_square, np.power(full_rg,2))[0][1] return c
# ........................................................................ # #
[docs] def get_secondary_structure_DSSP(self, R1=None, R2=None, 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. Parameters ---------- R1 : int {None} Default value is None. Defines the value for first residue in the region of interest. If not provided (False) then first residue is used. R2 : int {None} Default value is None. Defines the value for last residue in the region of interest. If not provided (False) then last residue is used. correctOffset : Bool 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. If you're calling this function you can probably ignore this variable. Returns ------- ddsp_vector : np.array A 4xn numpy array (where n is the number of residues) in which column 1 defines the residue index and column 2-4 defines the fractional occupancy of helical (H), extended (E) (beta-like) and coil (C) states. Note the three classifications will sum to 1 (within numerical precision). """ # build R1/R2 values out = self.__get_first_and_last(R1, R2, withCA=True) R1_real = out[0] R2_real = out[1] # select the relevant subtrajectory (out[2] is the 'resid %i to %i' where %i and %i are R1 and R2) ats = self.traj.atom_slice(self.topology.select('%s' % out[2])) # compute DSSP over the selected subtrajectory dssp_data = md.compute_dssp(ats) C_vector = [] E_vector = [] H_vector = [] # note the + 1 because the R1 and R2 positions are INCLUSIVE whereas reslist = list(range(R1_real, R2_real+1)) n_frames = self.n_frames for i in range(len(reslist)): C_vector.append(float(sum(dssp_data.transpose()[i] == 'C'))/n_frames) E_vector.append(float(sum(dssp_data.transpose()[i] == 'E'))/n_frames) H_vector.append(float(sum(dssp_data.transpose()[i] == 'H'))/n_frames) return np.array((reslist, H_vector, E_vector, C_vector))
# ........................................................................ # #
[docs] def get_secondary_structure_BBSEG(self, R1=None, R2=None, correctOffset=True): """ Returns a dictionary where eack key-value pair is keyed by a BBSEG classification type (0-9) and each value is a vector showing the fraction of time each residue is in that particular BBSEG type. BBSEG classification types are listed below 0 - unclassified 1 - beta (turn/sheet) 2 - PII (polyproline type II helix) 3 - Unusual region 4 - Right handed alpha helix 5 - Inverse C7 Equatorial (gamma-prime turn) 6 - Classic C7 Equatorial (gamma turn) 7 - Helix with 7 Residues per Turn 8 - Left handed alpha helix Parameters R1 and R2 are optional and allow a local sub-region to be defined. The return dictionary provides a classification vector for each of the 9 possible types of classification (note 0 = unclassified). Parameters ---------- R1 : int Default value is False. Defines the value for first residue in the region of interest. If not provided (False) then first residue is used. R2 : int Default value is False. Defines the value for last residue in the region of interest. If not provided (False) then last residue is used. correctOffset : Bool 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. If you're calling this function you can probably ignore this variable. Returns ------- return_bbseg : dict Dictionary of 9 key-value pairs where keys are integers 0-8 and values are numpy arrays showing the fractional occupancy of each of the distinct types of defined secondary structure. Note the three classifications will sum to 1 (within numerical precision). """ # build R1/R2 values out = self.__get_first_and_last(R1, R2, withCA=True) # extract the phi/psi angles in degrees phi_data = np.degrees(md.compute_phi(self.traj.atom_slice(self.topology.select('%s'%(out[2]))))[1]) psi_data = np.degrees(md.compute_psi(self.traj.atom_slice(self.topology.select('%s'%(out[2]))))[1]) # extract the relevant information (note shape of phi_data and psi_data will be identical) # shape info here is (number_of_frames, number_of_residues) sized shape_info = np.shape(phi_data) all_classes = [] # for each frame iterate through and classify each residue, building a shape_info # sized matrix where each elements reflects the BBSEG classification of that residue # in a given frame for f in range(0, shape_info[0]): # so each step through the loop we're passing two vectors, each of which # is nres residues long all_classes.append(self.__phi_psi_bbseg(phi_data[f], psi_data[f])) # convert to a numpy array all_classes = np.array(all_classes) # finally cycle through each BBSEG classification type and average # over each frame (shape_info[0] is number of frames) return_bbseg = {} for c in range(0,9): return_bbseg[c] = list(np.sum((all_classes == c)*1,0)/shape_info[0]) return return_bbseg
# ........................................................................ # def __phi_psi_bbseg(self, phi_vector, psi_vector): """ Internal function that takes two equally-matched phi and psi angle vectors and based on the pairwise combination classified each pair of elements using the BBSEG2 definition. Definition was generated from the BBSEG2 file distributed with CAMPARI, and is encoded and stored in the _internal_data module. NOTE that because this is an internal function we do not double check that the phi_vector and psi_vectors are of the same length, but this is critical, so if this function is being called make sure this is true! Parameters ---------- phi_vector : iterable (list or numpy vector) ordered list of phi angles for a specific residue psi_vector : iterable (list or numpy vector) ordered list of psu angles for a specific residue Returns ------- classes : list A list of length equal to phi_vector and psi_vector that classifies each pair of phi/psi angles using the BBSEG2 definition. """ classes = [] for i in range(len(phi_vector)): phi = phi_vector[i] psi = psi_vector[i] fixed_phi = phi-(phi%10) fixed_psi = psi-(psi%10) # following corrections for edge cases if we hit if fixed_phi == 180.0: fixed_phi = 170.0 if fixed_psi == 180.0: fixed_psi = 170.0 # classify the phi/psi values in terms of BBSEG values classes.append(BBSEG2[fixed_phi][fixed_psi]) return classes # ........................................................................ #
[docs] def get_overlap_concentration(self): """ Returns the overlap concentration for the chain. The overlap concentration reflects the concentration at which a flexible polymer begins to 'collide' in trans with other polymers - i.e. the concentration at which the chains begin to overlap. Returns ------- float Molar concentration for the overlap concentration. """ return ctpolymer.get_overlap_concentration(np.mean(self.get_radius_of_gyration()))
# ........................................................................ # #
[docs] def get_angle_decay(self, atom1='C', atom2='N', return_full_matrix=False): """ Returns the a 4 by n numpy array in which column 1 gives residue number, column 2 is local helicity, No checking of atom1 and atom2... Parameters ---------- atom1: str {C} The first atom to use when calculating the angle decay. atom2: str {N} The second atom to use when calculating the angle decay. return_full_matrix: bool {False} Whether or not to return the full matrix along with the angle decay calculation. Returns ------- array_like, or 2-tuple If `array_like`, the matrix returned is comprised of only the angle decay. If a 2-tuple, both the angle decay matrix (index 0) and the full matrix is returned (index 1). """ # first compute all the C-N vector for each residue CN_vectors = [] CN_lengths = [] for i in self.resid_with_CA: # this extracts the C->N vector for each frame for each residue value = np.squeeze(self.traj.atom_slice(self.__residue_atom_lookup(i, atom1)).xyz) - np.squeeze(self.traj.atom_slice(self.__residue_atom_lookup(i, atom2)).xyz) # CN_vectors becomes a list where each element is [3 x nframes] array where 3 is the x/y/z vector coordinates CN_vectors.append(value) # CN_lengths extracts the ||v|| length of each vector (should be basically the same) CN_lengths.append(np.linalg.norm(value,axis=1)) # calculate the number of residues for which we have C->N vectors npos = len(CN_vectors) # initialize an empty dictionary - the keys for this are i-j sequence separation, # and values are the cos(theta) angle between pairs of vectors all_vals={} for i in range(1, npos): all_vals[i] = [] # precompute length_multiplier = {} for i1 in range(0, npos-1): length_multiplier[i1] = {} for j1 in range(i1+1, npos): length_multiplier[i1][j1] = CN_lengths[i1]*CN_lengths[j1] # we then cycle over the non-redudant set of pairwise residues in the protein for i1 in range(0, npos-1): for j1 in range(i1+1, npos): # for each frame calculate (u . v) / (||u|| * ||v||) # where u and v are vectors and "." is the dot product between each pair. We're only calculating PAIR-WISE dot product # of each [x,y,z] with [x,y,z] vector, so doing np.sum(Matrix*matrix) is SO SO SO much faster than anything else. We # also take the average to avoid storing a ton of numbers and generating these giant vectors all_vals[j1-i1].append(np.mean(np.sum(CN_vectors[i1]*CN_vectors[j1],axis=1)/length_multiplier[i1][j1])) return_matrix = [] return_matrix.append([0,1.0,0.0]) for k in all_vals: return_matrix.append([k, np.mean(all_vals[k]), np.std(all_vals[k])]) # if we want the nres by nres matrix with specific decay <cos(omega)> for each specific pairwise # residue-residue set if return_full_matrix: full_matrix = np.zeros((len(return_matrix),len(return_matrix))) # for 0 to the number of residues (i.e. each row in the [nres x nres] matrix for i in range(0, len(return_matrix)): # set the column selector (c) to zero c = 0 # iterate through for j in range(0, len(all_vals[i])): full_matrix[i,c] = all_vals[i][j] for j in range(len(all_vals[i]), len(return_matrix)): full_matrix[i,c] = 0.0 return (return_matrix, full_matrix) else: return return_matrix
#oxoxoxoxoxooxoxoxoxoxoxoxoxoxoxoxooxoxoxoxoxoxoxoxoxoxoxooxoxoxoxoxoxoxoxoxoxoxooxoxo # #
[docs] def get_local_collapse(self, window_size=10, bins=None, verbose=True): """ local collapse calculates a vectorial representation of the radius of gyration along a polypeptide chain. This makes it very easy to determine where you see local collapse vs. local expansion. Parameters ---------- window_size : int, default=10 Size of the window over which conformations are examined. Default is 10. bins : np.arange or list A range of values (np.arange or list) spanning histogram bins. Default is np.arange(0, 10, 0.1). verbose : bool Flag that by default is True determines if the function prints status updates. This is relevant because this function can be computationally expensive, so having some report on status can be comforting! Returns ------- tuple (len = 4) return[0] : list of floats of len *n*, where each float reports on the mean local collapse a specific position along the sequence as defined by the fragment_size return[1] : list of floats of len *n* , where each float reports on the standard deviation of the local collapse at a specific position along the sequence as defined by the fragment_size return[2] : List of np.ndarrays of len *n*, where each sub-array reports on the histogram values associated with the local collapse at a given position along the sequence return[3] : np.ndarray which corresponds to bin values for each of the histograms in return[2] """ # validate bins if bins is None: bins = np.arange(0,10,0.01) else: try: if len(bins) < 2: raise CTException('Bins should be a numpy defined vector of values - e.g. np.arange(0,1,0.01)') except TypeError: raise CTException('Bins should be a list, vector, or numpy array of evenly spaced values') try: bins = np.array(bins, dtype=float) except ValueError: raise CTException('Passed bins could not be converted to a numpy array of floats') n_residues = self.n_residues n_frames = self.n_frames # check the window is an appropriate size if window_size > n_residues: raise CTException('window_size is larger than the number of residues') meanData = [] stdData = [] histo = [] for i in range(window_size - 1, n_residues): ctutils.status_message("On range %i" % i, verbose) # get radius of gyration (now by default is in Angstroms # - in previous versions we performed a conversion here) tmp = self.get_radius_of_gyration(i - (window_size-1), i) (b, c) = np.histogram(tmp, bins) histo.append(b) meanData.append(np.mean(tmp)) stdData.append(np.std(tmp)) return (meanData, stdData, histo, bins)