rdcPotTools
index
/home/schwitrs/xplor/python/rdcPotTools.py


 
Tools to aid in setup/analysis of dipolar coupling potential term.
 
This module provides functions to simplify the creation, manipulation and
analysis of rdcPot.RDCPot potential terms.
 
Note on dipolar coupling normalization and sign convention:
 
  Protein residual dipolars are usually normalized relative to
  measurements of N-HN - this is achieved by calling the function
  scale_toNH.  This normalization involves the gamma_A * gamma_B /
  rAB^3 prefactor (see the docs for rdcPot).  15N has a negative
  gyromagnetic ratio, but it was deemed simpler to treat as positive,
  so that if there is an experiment which does not include 15N, the
  dipolar coupling sign must be flipped.  If you would rather not use
  this convention, call the function correctGyromagneticSigns() at the
  beginning of your script, and then the correct signs will be used
  for the Da prefactors.

 
Functions
       
Da(potTerm)
DaMax(type, dist=0)
Return scale factor to scale other experiments to the NNH norm
convention.
DaScale(term, dist=0)
return scale factor to scale other experiments to the current
convention, set by calling either scale_toNH or scale_toCH.
RDC_analyze(potList)
Perform analysis of RDCPot terms and return nicely formatted
    summary.
 
    For each RDC term, the following is printed in the PDB header:
 
    
RMS
  root mean square difference between calculated and observed.
R-fac   
  Eq. 5 from Clore+Garrett, JACS 121, 9008 (1999).  
  Da is scaled using the current scaling convention.       
Devia
  Deviation of calculated values in an EnsembleSimulation calculation.
Da
  Alignment tensor Da.
Rh
  rhombicity of the alignment tensor
R-inf
  R-factor (in percent) for an infinite number of
  randomly distributed vectors.
  Eq. 3 from Clore+Garrett, JACS 121, 9008 (1999)
Viols
  number of RDCs with :math:`|calcd - obs|` values larger than threshold().
Chi^2
  A normalized :math:`\chi^2 =  1/N \sum_i (calcd_i-obs_i)^2 / err_i^2`. It
  is not meaningful if errors are not realistic.
corr
  Pearson's correlation coefficient between calculated and observed RDCs.
Rfactor(rdc, selection='all')
R-factor (in percent).
 
Eq. 5 from Clore+Garrett, JACS 121, 9008 (1999).
Da is scaled using the current scaling convention.
 
The selection argument can be used to choose a subset of
restraints whose atoms lie in selection.
Rfactor_infinite(rdc, selection='all')
R-factor (in percent) for an infinite number of
randomly distributed vectors.
Eq. 3 from Clore+Garrett, JACS 121, 9008 (1999).
 
The selection argument can be used to choose a subset of
restraints whose atoms lie in selection
Rh(potTerm)
chi2(pot, selection='all')
Compute the Chi^2 for the restraints associated with the
specified rdcPot.RDCPot term.
 
The selection argument can be used to choose a subset of
restraints whose atoms lie in selection
composite_Rfactor_infinite(rdcs, selection='all', normalize=True)
R-factor (in percent) for an infinite number of
randomly distributed vectors.
Eq. 3 from Clore+Garrett, JACS 121, 9008 (1999).
 
This version takes a list of rdcPot.RDCPot terms and produces a
composite R-factor by simply performing sums over all (mixed) RDCs. For
this to work the RDCs should be unscaled.
 
The selection argument can be used to choose a subset of
restraints whose atoms lie in selection.
 
The normalize argument specifies whether or not to rescale RDCs by
dividing by the Da prefactor to that they are on the same scale
prior to calculating the R-factor.
composite_chi2(rdcs, selection='all')
Compute the Chi^2 for the restraints associated with the
specified squence of rdcPot.RDCPot terms.
 
The selection argument can be used to choose a subset of
restraints whose atoms lie in selection
correctGyromagneticSigns()
Use the correct (negative) sign for 15N's gyromagnetic ratio.
 
Call this function before any other in this module.
correlation(rdc)
Return correlation of calculated to observed RDC values.
create_DaRatioPot(name, rdc1, rdc2, ratio)
NOTE: needs to be updated and transfered to varTensorTools.
 
create a potential term which restrains the ratio of
  rdc1.Da/rdc2.Da to be ratio
create_RDCPot(name, file=0, oTensor=0, esim=0, restraints='', subSel='')
Create an rdcPot.RDCPot1 with given name, given orientational
varTensor.VarTensor, the filename of an rdc assignment table and/or
a string of assignments, and an ensemble simulation.
 
The file argument can optionally be a sequence of filenames.
 
The subSel argument is passed to as the selectionFilter argument to
rdcPot.RDCPot.addRestraints so that only the subset of restraints
which match the selection will be loaded.
create_twoDaRatiosPot(name, rdc1, rdc2, rdc3, rdc4)
NOTE: needs to be updated and transferred to varTensorTools.
 
Create a potential term which restrains the ratio of
  rdc1.Da/rdc2.Da to be that of rdc3.Da/rdc4.Da
For EnsembleSimulations, three terms are returned in a PotList. These
are named 'inter', 'ratio1', and 'ratio2'. The first restrains the
ensemble average of the two ratios to be the same. The 'ratio' terms
restrain the ratio of each ensemble member to be the same. The terms
can be weighted separately using the scale member.
deviation_percent(rdc)
Deviation as a percent of total range of dipolar coupling values.
Da is scaled using the current scaling convention.
getRDCType(rdc)
Return the atoms involved in the RDC experiment from the atom selections.
makeTable(rdc)
Return the assignment table (a string) corresponding to the 
restraints associated with the specified rdcPot.RDCPot.
orderParameter(aSel, bSel, simulation=0)
Calculate the order parameter S^2 given two atom selections defining
a bond vector. The simulation parameter specifies an EnsembleSimulation
and defaults to the current simulation. It is assumed that the atomSels
correspond to single atoms: averaging is not performed
readNEF(potList, block, oTensor=None, subSel='all', minErr=0.1)
Given a RDCPot object, and a saveframe from a NEF record, read in
restraints. Each saveframe corresponds to a single varTensor.VarTensor 
object which is returned. The passed potList argument should be a
potList.PotList object. Restraints with with the same value of
scale, weight and distance_dependent are grouped into rdcPot.RDCPot
objects which are appended to the potList object.
minErr gives the minimum error value allowed.
 
oTensor is an optional varTensor.VarTensor object to use as an alignment
tensor. If not specified one is created.
 
The subSel argument is passed as the selectionFilter argument to
the RDCPot's addRestraints method.
recalculatePrefactors()
Recalculate gamma_A * gamma_B / rAB^3 prefactors.
 
Call this function after changing gyromagnetic ratios or bond distance
values.
scaleByErrs(rdc)
For the specified RDCPot term, set the per-restraint energy
scale factor to be 1/plusErr^2.
scaleDa(rdc, type, ref, custom)
Set the Da associated with input RDCs based on that for a reference pair.   
 
Given an instance of rdcPot.RDCPot1 (rdc argument) containing RDC
restraints exclusively involving nuclei of type i and j, and the
specification of this nucleus pair type (type argument, a string; see
below), this function sets the effective Da, Da_eff, for each restraint to:
 
Da_eff = {[gamma(i) gamma(j)/r_ij^3]/[gamma(m) gamma(n)/r_mn^3]} Da  [Eq. 1]
 
where gamma(x) is the gyromagnetic ratio of nucleus x, r_xy is the x-y
internuclear distance, and Da is the Da parameter of the corresponding
varTensor.VarTensor instance.  m and n indicate the reference nucleus
pair specified by the ref argument (a string; see below).
 
Valid values for the argument type/ref are:
 
Values        Pair      Distance (Angstroms)
'NHN','NH'    N-HN      1.041      [Ottiger & Bax 1998]
'NCO'         C-N       1.329      [Engh & Huber 1991]
'HNCO'        C-HN      2.105      [Cornilescu & Bax 2000] 
'CACO'        CA-C      1.525      [Engh & Huber 1991]
'CACB'        CA-CB     1.525      (as 'CACO')
'CAHA','CH'   CA-HA     1.104      [Cornilescu & Bax 2000]
'CAHN'        CA-HN     2.560 
'methyl'      MeC-MeH   1.517      (C-CH3 distance; see below)
'HH'          1H-1H     1          (variable-distance RDC)
'HP'          1H-31P    1          (variable-distance RDC)
 
Bibliography: Cornilescu & Bax (2000) JACS 122:10143
              Engh & Huber (1991) Acta Crystallogr A47:392
              Ottiger & Bax (1998) JACS 120:12334
              Ottiger & Bax (1999) JACS 121:4690
 
where the atom pair and its (RDC-effective) interatomic distance is
specified.  Named atoms in the pairs (e.g., CA), are assumed to belong to
protein.  All atoms in the pairs are assumed to belong to the same residue,
except for options 'NCO', 'HNCO' and 'CAHN' (where they are sequential), and
'HH' and 'HP' (which correspond to variable-distance RDCs; see below).  The
'methyl' option specifies the one-bond methyl C-H RDC, the Da_eff of which
is calculated as [Ottiger & Bax 1999]:
 
Da_eff = -3.17 {[gamma(c)^2/r_cc^3]/[gamma(m) gamma(n)/r_mn^3]} Da   [Eq. 2]
 
where gamma(c) is the 13C gyromagnetic ratio and r_cc is the C-CH3 bond
distance in proteins.  Options 'HH' and 'HP' are special cases that
correspond to variable-distance RDCs, for which Eq. 1 is used with unit-
valued r_ij, as the distance dependence is assumed explicitly accounted for
by appropriate configuration of rdc via rdc.setUseDistance(True) outside
this function (see rdcPot.RDCPot1 documentation).
 
The specific values for the [gamma(i) gamma(j)/r_ij^3] factors used for the
different type/ref options can be retrieved from the Da_prefactor dictionary
in this module.  For example,
 
print(Da_prefactor['CAHA'])
 
prints [gamma(CA) gamma(HA)/r_CAHA^3].
 
Alternatively, a custom prefactor, K, can be specified via the argument
custom (a number) that completely overrides Eqs. 1 and 2 above with:
 
Da_eff = K Da                                                        [Eq. 3]
 
Specifying any value for type/ref along with custom has no effect.
scaleRDC(rdc, type='', ref='NHN', custom=None)
Set the Da associated with input RDCs based on that for a reference pair.   
 
This function is useful for mixing RDCs from different nucleus pair types
without pre-normalization to a reference pair.
 
rdc is an instance of rdcPot.RDCPot1, corresponding to RDCs of the
nucleus pair type optionally specified via the argument type (if type is not
specified, it is guessed from the atom selections of the first restraint of
rdc with function getRDCType).  The type of the reference nucleus pair is
indicated with argument ref.
 
The function that does the actual work is scaleDa, which shares the same
arguments.  See scaleDa documentation for a more detailed description of all
the arguments and functionality.
scale_toCH(rdc, type='')
Set the Da associated with input RDCs based on that for the C-H pair.   
 
This function is useful for mixing RDCs from different nucleus pair types
without pre-normalization to the directly bonded 13C-1H pair.
 
rdc is an instance of rdcPot.RDCPot1, corresponding to RDCs of the
nucleus pair type optionally specified via argument type (if type is not
specified, it is guessed from the atom selections of the first restraint of
rdc with function getRDCType).  If type is set to 'HH' or 'HP', the distance
dependence of rdc is enabled (via rdc.setUseDistance(True)) as a side
effect.
 
The function that does most of the actual work is scaleDa (which shares the
rdc and type arguments), with its argument ref set to 'CH'.  See scaleDa
documentation for a more detailed description of all the arguments and
functionality.
 
This function is a less general version of scaleRDC.
scale_toNH(rdc, type='')
Set the Da associated with input RDCs based on that for the N-H pair.   
 
This function is useful for mixing RDCs from different nucleus pair types
without pre-normalization to the directly bonded 15N-1H pair.
 
rdc is an instance of rdcPot.RDCPot1, corresponding to RDCs of the
nucleus pair type optionally specified via argument type (if type is not
specified, it is guessed from the atom selections of the first restraint of
rdc with function getRDCType).  If type is set to 'HH' or 'HP', the distance
dependence of rdc is enabled (via rdc.setUseDistance(True)) as a side
effect.
 
The function that does most of the actual work is scaleDa (which shares the
rdc and type arguments), with its argument ref set to 'NH'.  See scaleDa
documentation for a more detailed description of all the arguments and
functionality.
 
This function is a less general version of scaleRDC.
spaceSeparatedToRestraint(inString, atom1Name='N', atom2Name='HN', residCol=1, rdcCol=2, errCol=None, defaultErr=0.1, segid=None)
Convert string restraint table (inString) consisting of columns
for resid and RDC values to an Xplor-NIH readable restraint
table. Column numbers start with 1. This function assumes that the RDCs are
intraresidue with atom names spacified by atom1Name and atom2Name.
 
For RDCs associated with the amide bond, the resid column in inString 
refers to the preceeding (N-terminal) residue.
writeNEF(rdc, name)
Return a formatted NEF record for the restraints in the given RDCPot
object as a string with the specified saveframe name.
writeNEF_VarTensor(vten, name)
Return a formatted NEF record for the restraints associated with the
varTensor.VarTensor object vten as a string with the specified saveframe
name.

 
Data
        Da_prefactor = {'CACB': 12764422.29508197, 'CACO': 12764422.29508197, 'CAHA': 133768561.96321149, 'CAHN': 10728558.721184729, 'CH': 133768561.96321149, 'HH': 715668433.9204, 'HNCO': 19297674.55632729, 'HP': 290018215.18, 'NCO': 7772415.869992006, 'NH': 64302723.53092173, ...}
gamma_C = 6728.3
gamma_H = 26751.98
gamma_N = 2711.6
gamma_P = 10841.0
headerHelpString = '\nRMS\n root mean square difference between calcu...oefficient between calculated and observed RDCs.\n'
r_cacb = 1.525
r_caco = 1.525
r_caha = 1.104
r_cahn = 2.56
r_ccme = 1.517
r_cohn = 2.105
r_nco = 1.329
r_nhn = 1.041