Hi,
I have taken over a project involving a complex of a protein with a
peptide with sulfated tyrosine.
I have a few questions.
I noticed, that the masses for several of the atoms of the modified
tyrosine were missing from the topology file. I have added them now, but
what are the consequences of missing masses for the structure calculations?
With the refine.py script I am getting errors for duplicated entries like
these
%RTFRDR-ERR: duplicate (P-)RESIdue name CIPP
%PARRDR-ERROR: duplication of bond C C
It seems that the topology file is somehow read twice, but I can't figure
out, where that happens.
I don't get those errors with anneal.py (I copied anneal.py and refine.py
from the examples in gb1-rdc and commented out the nonrelevant parts, as we
dont have RDCs).
Any help would be much appreciated!
Thanks and a happy, healthy and successful New Year,
Sabine
Here is the script I am using
xplor.requireVersion("2.45")
#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2005/05/10
#
(opts,args) = xplor.parseArguments(["quick"]) # check for command-line typos
quick=False
for opt in opts:
if opt[0]=="quick": #specify -quick to just test that the script runs
quick=True
pass
pass
numberOfStructures=20
if quick:
numberOfStructures=3
pass
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed() #explicitly set random seed
#
# annealing settings
#
command = xplor.command
#protocol.initParams("protein")
import psfGen
psfGen.residueTypes['protein'].append("TYS")
protocol.initTopology(["data/topallhdgntccr5_M.pro"])#["protein","data/
tys.pro"])
protocol.initParams(["data/parallhdgntccr5_SA.pro"])#,"data/tys.par"])
command = xplor.command
#addDisulfideBond('resid 10 and name SG and segid A', 'resid 34 and name
SG and segid A')
#psfGen.addDisulfideBond('resid 11 and name SG and segid A', 'resid 50 and
name SG and segid A')
#assign ( resid 10 and name SG and segid A) ( resid 34 and name SG
and segid A) 2.05 0.05 0.05
#assign ( resid 11 and name SG and segid A) ( resid 50 and name SG
and segid A) 2.05 0.05 0.05
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#from psfGen import seqToPSF
#seqToPSF('psf/rantes.seq',segName='A',disulfide_bonds=[])
## generate a random extended structure with correct covalent geometry
## saves the generated structure in the indicated file for faster startup
## next time.
#psfGen.addDisulfideBond('resid 11 and name SG and segid A', 'resid 50 and
name SG and segid A')
#psfGen.addDisulfideBond('resid 10 and name SG and segid A', 'resid 34 and
name SG and segid A')
##protocol.initStruct("g_new.psf") # - or from file
## generate a random extended structure with correct covalent geometry
## saves the generated structure in the indicated file for faster startup
## next time.
##
#protocol.genExtendedStructure("re_rantes_%d.pdb" %
#protocol.initialRandomSeed())
# or read an existing model
#
protocol.loadPDB("complex.pdb",deleteUnknownAtoms=True)
protocol.fixupCovalentGeom(maxIters=100,useVDW=1)
#
# a PotList contains a list of potential terms. This is used to specify
which
# terms are active during refinement.
#
from potList import PotList
potList = PotList()
# parameters to ramp up during the simulated annealing protocol
#
from simulationTools import MultRamp, StaticRamp, InitialParams
rampedParams=[]
highTempParams=[]
# compare atomic Cartesian rmsd with a reference structure
# backbone and heavy atom RMSDs will be printed in the output
# structure files
#
#from posDiffPotTools import create_PosDiffPot
#refRMSD = create_PosDiffPot("refRMSD","name CA or name C or name N",
#pdbFile='sa_ra_22.pdb',
#cmpSel="not name H*")
# orientation Tensor - used with the dipolar coupling term
# one for each medium
# For each medium, specify a name, and initial values of Da, Rh.
#
#from varTensorTools import create_VarTensor
#media={}
## medium Da rhombicity
#for (medium,Da,Rh) in [ ('t', -6.5, 0.62),
#('b', -9.9, 0.23) ]:
#oTensor = create_VarTensor(medium)
#oTensor.setDa(Da)
#oTensor.setRh(Rh)
#media[medium] = oTensor
#pass
# dipolar coupling restraints for protein amide NH.
#
# collect all RDCs in the rdcs PotList
#
# RDC scaling. Three possible contributions.
# 1) gamma_A * gamma_B / r_AB^3 prefactor. So that the same Da can be used
# for different expts. in the same medium. Sometimes the data is
# prescaled so that this is not needed. scale_toNH() is used for this.
# Note that if the expt. data has been prescaled, the values for rdc
rmsd
# reported in the output will relative to the scaled values- not the
expt.
# values.
# 2) expt. error scaling. Used here. A scale factor equal to 1/err^2
# (relative to that for NH) is used.
# 3) sometimes the reciprocal of the Da^2 is used if there is a large
# spread in Da values. Not used here.
#
#from rdcPotTools import create_RDCPot, scale_toNH
#rdcs = PotList('rdc')
#for (medium,expt,file, scale) in \
#[('t','NH' ,'tmv107_nh.tbl' ,1),
#('t','NCO','tmv107_nc.tbl' ,.05),
#('t','HNC','tmv107_hnc.tbl' ,.108),
#('b','NH' ,'bicelles_new_nh.tbl' ,1),
#('b','NCO','bicelles_new_nc.tbl' ,.05),
#('b','HNC','bicelles_new_hnc.tbl',.108)
#]:
#rdc = create_RDCPot("%s_%s"%(medium,expt),file,media[medium])
##1) scale prefactor relative to NH
## see python/rdcPotTools.py for exact calculation
## scale_toNH(rdc) - not needed for these datasets -
## but non-NH reported rmsd values will be wrong.
##3) Da rescaling factor (separate multiplicative factor)
## scale *= ( 1. / rdc.oTensor.Da(0) )**2
#rdc.setScale(scale)
#rdc.setShowAllRestraints(1) #all restraints are printed during analysis
#rdc.setThreshold(1.5) # in Hz
#rdcs.append(rdc)
#pass
#potList.append(rdcs)
#rampedParams.append( MultRamp(0.05,5.0, "rdcs.setScale( VALUE )") )
## calc. initial tensor orientation
## and setup tensor calculation during simulated annealing
##
#from varTensorTools import calcTensorOrientation, calcTensor
#for medium in media.keys():
#calcTensorOrientation(media[medium])
#rampedParams.append( StaticRamp("calcTensor(media['%s'])" % medium) )
#pass
# set up NOE potential
noe=PotList('noe')
potList.append(noe)
from noePotTools import create_NOEPot
for (name,scale,file) in [('unambig',10,"data/noe_inter_unambig.tab"),
('ambig',10,"data/noe_inter_ambig.tab"),
('rantes',1,"data/noe_rantes_4.tbl"),
('ccr5',1,"data/noe_ccr5_pdb_3.tbl")]:
pot = create_NOEPot(name,file)
pot.setPotType("soft") # if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
rampedParams.append( MultRamp(2,30, "noe.setScale( VALUE )") )
## set up J coupling - with Karplus coefficients
#from jCoupPotTools import create_JCoupPot
#jCoup = create_JCoupPot("jcoup","jna_coup.tbl",
#A=6.98,B=-1.38,C=1.72,phase=-60.0)
#potList.append(jCoup)
# Set up dihedral angles
from xplorPot import XplorPot
protocol.initDihedrals("data/dihedral_rantes_pdb.tbl",
#useDefaults=False # by default, symmetric sidechain
# restraints are included
)
potList.append( XplorPot('CDIH') )
highTempParams.append( StaticRamp("potList['CDIH'].setScale(10)") )
rampedParams.append( StaticRamp("potList['CDIH'].setScale(200)") )
# set custom values of threshold values for violation calculation
#
potList['CDIH'].setThreshold( 5 ) #5 degrees is the default value, though
# gyration volume term
#
# gyration volume term
#
#from gyrPotTools import create_GyrPot
#gyr = create_GyrPot("Vgyr",
#"resid 1:56") # selection should exclude disordered
tails
#potList.append(gyr)
#rampedParams.append( MultRamp(.002,1,"gyr.setScale(VALUE)") )
# hbda - distance/angle bb hbond term
#
protocol.initHBDA('data/hb_rantes_ccr5.tbl')
potList.append( XplorPot('HBDA') )
# hbdb - knowledge-based backbone hydrogen bond term
#
protocol.initHBDB()
potList.append( XplorPot('HBDB') )
#New torsion angle database potential
#
from torsionDBPotTools import create_TorsionDBPot
torsionDB = create_TorsionDBPot('torsionDB')
potList.append( torsionDB )
rampedParams.append( MultRamp(.002,2,"torsionDB.setScale(VALUE)") )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
from repelPotTools import create_RepelPot,initRepel
repel = create_RepelPot('repel')
potList.append(repel)
rampedParams.append( StaticRamp("initRepel(repel,use14=False)") )
rampedParams.append( MultRamp(.004,4, "repel.setScale( VALUE)") )
# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""initRepel(repel,
use14=True,
scale=0.004,
repel=1.2,
moveTol=45,
interactingAtoms='name CA'
)""") )
# Selected 1-4 interactions.
import torsionDBPotTools
repel14 = torsionDBPotTools.create_Terminal14Pot('repel14')
potList.append(repel14)
highTempParams.append(StaticRamp("repel14.setScale(0)"))
rampedParams.append(MultRamp(0.004, 4, "repel14.setScale(VALUE)"))
potList.append( XplorPot("BOND") )
potList.append( XplorPot("ANGL") )
potList['ANGL'].setThreshold( 5 )
rampedParams.append( MultRamp(0.4,1,"potList['ANGL'].setScale(VALUE)") )
potList.append( XplorPot("IMPR") )
potList['IMPR'].setThreshold( 5 )
rampedParams.append( MultRamp(0.1,1,"potList['IMPR'].setScale(VALUE)") )
# Give atoms uniform weights, except for the anisotropy axis
#
protocol.massSetup()
# IVM setup
# the IVM is used for performing dynamics and minimization in
torsion-angle
# space, and in Cartesian space.
#
from ivm import IVM
dyn = IVM()
# initially minimize in Cartesian space with only the covalent constraints.
# Note that bonds, angles and many impropers can't change with the
# internal torsion-angle dynamics
# breaks bonds topologically - doesn't change force field
#
#dyn.potList().add( XplorPot("BOND") )
#dyn.potList().add( XplorPot("ANGL") )
#dyn.potList().add( XplorPot("IMPR") )
#
#dyn.breakAllBondsIn("not resname ANI")
#import varTensorTools
#for m in media.values():
# m.setFreedom("fix") #fix tensor parameters
# varTensorTools.topologySetup(dyn,m) #setup tensor topology
#
#protocol.initMinimize(dyn,numSteps=1000)
#dyn.run()
# reset ivm topology for torsion-angle dynamics
#
dyn.reset()
#for m in media.values():
## m.setFreedom("fixDa, fixRh") #fix tensor Rh, Da, vary
orientation
#m.setFreedom("varyDa, varyRh") #vary tensor Rh, Da, vary
orientation
protocol.torsionTopology(dyn)
# minc used for final cartesian minimization
#
minc = IVM()
protocol.initMinimize(minc)
#for m in media.values():
#m.setFreedom("varyDa, varyRh") #allow all tensor parameters float
here
#pass
protocol.cartesianTopology(minc)
# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t = 3000. # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =0.5,
ivm=dyn,
rampedParams = rampedParams)
def accept(potList):
"""
return True if current structure meets acceptance criteria
"""
if potList['noe'].violations()>0:
return False
#if potList['rdc'].rms()>1.2: #this might be tightened some
#return False
if potList['CDIH'].violations()>0:
return False
if potList['BOND'].violations()>0:
return False
if potList['ANGL'].violations()>0:
return False
if potList['IMPR'].violations()>1:
return False
return True
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on
the
structure, and then writes out a pdb file, with remarks.
"""
# initialize parameters for high temp dynamics.
InitialParams( rampedParams )
# high-temp dynamics setup - only need to specify parameters which
# differfrom initial values in rampedParams
InitialParams( highTempParams )
# high temp dynamics
#
protocol.initDynamics(dyn,
potList=potList, # potential terms to use
bathTemp=init_t,
initVelocities=1,
finalTime=10, # stops at 10ps or 5000 steps
numSteps=5000, # whichever comes first
printInterval=100)
dyn.setETolerance( init_t/100 ) #used to det. stepsize. default:
t/1000
dyn.run()
# initialize parameters for cooling loop
InitialParams( rampedParams )
# initialize integrator for simulated annealing
#
protocol.initDynamics(dyn,
potList=potList,
numSteps=100, #at each temp: 100 steps or
finalTime=.2 , # .2ps, whichever is less
printInterval=100)
# perform simulated annealing
#
cool.run()
# final torsion angle minimization
#
protocol.initMinimize(dyn,
printInterval=50)
dyn.run()
# final all- atom minimization
#
protocol.initMinimize(minc,
potList=potList,
dEPred=10)
minc.run()
#do analysis and write structure when this function returns
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
structLoopAction=calcOneStructure,
pdbTemplate="comp_ref_4_STRUCTURE.pdb",
calcMissingStructs=True, #calculate only missing structures
doWriteStructures=True, #analyze and write coords after calc
genViolationStats=True,
averagePotList=potList,
averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
noe,potList['CDIH']],
#averageCrossTerms=refRMSD,
averageTopFraction=0.5, #report only on best 50% of structs
averageAccept=accept, #only use structures which pass
accept()
averageContext=FinalParams(rampedParams),
averageFilename="SCRIPT_complex_ave.pdb", #generate
regularized ave structure
averageFitSel="name CA",
averageCompSel="not resname ANI and not name H*" ).run()
_______________________________________________
Xplor-nih mailing list
[email protected]
https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
