-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
Hello Paul--
> I am a newcomer to Xplor-NIH (but not to Xplor) so apologies is this has
> come up before. We are gathering solution scattering and NMR data
> together on a multi-domain protein and would like to model its (dynamic)
> structure on the basis similar to that reported for the 128 kD EIN dimer
> and it 146 kD complex with HPr (JACS 2010 132: 13026-45), i.e. including
> the option of ensemble averaging. This report would appear to add
> additional features to the example script provided with the Xplor-NIH
> distribution that address the Dickerson dodecamer DNA-duplex
> (xplor-nih-2.26/eginput/dna-refi/ensemble.py). Is there available any
> more up-to-date or illustrative example script than this one (have I
> missed something?), or is this a good enough model to base our effort on?
>
I'm attaching a script which we used for this purpose, although it
turned out that an ensemble size of one was the largest we could
justify. Specifically, this is the script for the EI-free
calculations. I intend to work this into an eginput example to go with
the distribution at some point.
The attached script does only the refinement stage- getting good
structures from random coords will take multiple runs.
Please let me know if you have questions.
best regards--
Charles
# ensemble free EI refinement
# the RDC alignment tensor symmetric
# The free-EI half of these structures will be fixed in space, but
# contribute to the SAXS and RDC calcs.
# then: attempt to make broken bond by moving only the linker region
# if we can't, print fail and go to next structure
# if we can, do refinement
import os
print 'process id:', os.getpid()
numberOfStructures=120
import protocol
protocol.initRandomSeed(42)
ensembleSize=1
# file index description
startStructure="genXrayDimer6_newc.pdb"
pdbTemplate="SCRIPT_STRUCTURE_MEMBER.fit"
bboneAtoms = "name C or name CA or name N or name O or name HN"
linker1Sel = "resid 22:24 or resid 144:146" #between alpha and alphaBeta
linker1Break1 = "(resid 24 and name c) or (resid 25 and name n)"
linker1Break2 = "(resid 146 and name c) or (resid 147 and name n)"
linker2Sel = "resid 255:261" #linker between N and C termini
linker2Break = "(resid 261 and name c) or (resid 262 and name n)"
nTerminus = "resid 1:254 or resid 601:685"
cTerminus = "resid 262:573"
alphaBetaSel = "resid 1:21 or resid 147:254"
alphaSel = "resid 25:143"
protocol.loadPDB(startStructure)
contactNterm=AtomSel(
"tag and segid A and (%s) and ( ((%s) or (%s)) around 14.5)" %
(nTerminus,cTerminus,linker2Sel))
contactCterm=AtomSel(
"tag and segid A and (%s) and ( ((%s) or (%s)) around 14.5)" %
(cTerminus,nTerminus,linker2Sel))
rstring=''
for r in contactNterm:
rstring += 'resid %d or ' % r.residueNum()
pass
rstring += "not all"
print map(lambda a:a.residueNum(),contactNterm)
print map(lambda a:a.residueNum(),contactCterm)
nonContactNtermA=AtomSel("segid A and (%s) and not (%s)" % (nTerminus,
rstring))
nonContactNtermB=AtomSel("segid B and (%s) and not (%s)" % (nTerminus,
rstring))
print len(nonContactNtermA)
print len(nonContactNtermB)
rstring=''
for r in contactCterm:
rstring += 'resid %d or ' % r.residueNum()
pass
rstring += "not all"
nonContactCterm=AtomSel("(%s) and not (%s)" % (cTerminus,rstring))
print len(nonContactCterm)
#protocol.addUnknownAtoms() - FIX: this doesn't work yet.
xplor.simulation.deleteAtoms("not known")
from ensembleSimulation import EnsembleSimulation
esim=EnsembleSimulation('ens',ensembleSize)
esim.setAveType('sum')
# note equal weights for saxs members- these are changed for rdc and saxs
# terms
ensIndex=esim.member().memberIndex()
import trace
#trace.addModuleName("varTensorTools")
#trace.addModuleName("solnXRayPotTools")
#trace.addModuleName("psfGen")
#from atomSelAction import SetProperty
#from vec3 import Vec3
#AtomSel("segid B").apply(SetProperty("pos",Vec3(0,0,0)))
# fit the c-termini of these structures
from potList import PotList
potList = PotList()
crossTerms=PotList('cross terms')
crossTerms=PotList('cross terms')
from simulationTools import StaticRamp, MultRamp
from simulationTools import InitialParams, FinalParams
rampedParams=[]
highTempParams=[]
media=[]
from varTensorTools import create_VarTensor, addAxisAtoms
from varTensorTools import calcTensor
axisAtoms=addAxisAtoms(resid=999)
ten1=create_VarTensor("ten-double",axis=axisAtoms)
media.append(ten1)
#FIX: should there be two alignment tensors (one each for singly/doubly bound)
ten1.setFreedom("varyDa, varyRh")
#ten1.setFreedom("fixDa, fixRh")
rdcs=PotList("rdc")
from rdcPotTools import create_RDCPot
rdc1 = create_RDCPot("rdc1","free_Mg_A.tbl",ten1)
calcTensor(ten1) #this to set Da, Rh
rdc1.addRestraints(open("free_Mg_B.tbl").read())
rdc1.setShowAllRestraints(1)
rdcs.append(rdc1)
potList.append(rdcs)
from simulationTools import analyze
print analyze(rdcs)
#import sys; sys.exit(0)
#rampedParams.append( StaticRamp('for t in media: calcTensor(t)') )
from avePot import AvePot
from xplorPot import XplorPot
# ``NCS'' term - keep dimer subunits identical
from posDiffPotTools import create_PosDiffPot
ncs=PotList("ncs")
ncs1 = create_PosDiffPot("ncs1","segid A","segid B")
ncs.append( ncs1 )
ncs.setScale(1000)
#from simulationTools import analyze
#print analyze(ncs)
#import sys; sys.exit()
xplor.command(r"""
ncs restraints
initialize
group
!equi ((segid A and resid 253:261) )
!equi ((segid B and resid 253:261) )
equi ((segid A) )
equi ((segid B) )
weight 1.
end
?
end
""")
ncs = AvePot(XplorPot,"NCS")
ncs.setScale(10)
potList.append( ncs )
#rampedParams.append( StaticRamp('ncs.setScale(10)') )
#highTempParams.append( StaticRamp('ncs.setScale(10)') )
#this applies to the doubly bound form.
from distSymmTools import create_DistSymmPot, genDimerRestraints
dSymm=PotList("dsymm")
dSymmP = create_DistSymmPot('dSymmP')
for r in genDimerRestraints(segids=['A','B'],
resids=range(10,380,10)):
# resids=range(340,380,10)):
dSymmP.addRestraint(r)
pass
dSymmP.setShowAllRestraints(True)
dSymmP.setEnsembleAveType("average")
dSymm.append(dSymmP)
dSymm.setScale(100)
potList.append(dSymm)
from xplorPot import XplorPot
protocol.initDihedrals("torsion_linker-20100524.tbl",
useDefaults=False)
dihed=AvePot(XplorPot,"CDIH")
potList.append(dihed)
highTempParams.append( StaticRamp("potList['CDIH'].setScale(0.1)") )
rampedParams.append( MultRamp(0.1,200,"potList['CDIH'].setScale(VALUE)") )
#
#enabledInteractions=[
# ("not (%s)" % cTerminus ,"not (%s)" % cTerminus ) ,
# ("segid A and not (%s)"%cTerminus ,"segid A and bondedto (%s)" % cTerminus ),
# ("segid B and not (%s)"%cTerminus ,"segid B and bondedto (%s)" % cTerminus ),
# ]
#
nbond=AvePot(XplorPot,"VDW")
potList.append(nbond)
#
#cmd ="constraints\n"
#for pair in enabledInteractions:
# cmd+="""inte (name CA and (%s))
# (name CA and (%s)) weights * 0 vdw 1 end\n""" % pair
# cmd+="""inte ( %s) ( %s) weights * 1 vdw 0 end\n""" % pair
# pass
#
#
#cmd += "end\n"
#
#xplor.command(cmd)
#print cmd
#htConsInteCMD=cmd
#
#cmd = "constraints\n"
#for pair in enabledInteractions:
# cmd += " inte (%s) (%s)\n" % pair
#cmd += "end\n"
#print cmd
#consInteCMD=cmd
#
rampedParams.append( StaticRamp("""protocol.initNBond()""") )
#rampedParams.append( StaticRamp("xplor.command(consInteCMD)") )
rampedParams.append( MultRamp(0.9,0.8,
"xplor.command('param nbonds repel VALUE end end')") )
rampedParams.append( MultRamp(.004,4,
"xplor.command('param nbonds rcon VALUE end end')") )
highTempParams.append( StaticRamp("""protocol.initNBond(rcon=0.1,
repel=1.2,
onlyCA=1)""") )
#highTempParams.append( StaticRamp("xplor.command(htConsInteCMD)") )
from solnXRayPotTools import create_solnXRayPot
import solnXRayPotTools
xray=create_solnXRayPot('xray',experiment='EI-Mg-saxswaxs-trunc.dat',
numPoints=26,
normalizeIndex=-3,preweighted=False)
xrayCorrect=create_solnXRayPot('xray-c',experiment='EI-Mg-saxswaxs-trunc.dat',
numPoints=26,
normalizeIndex=-3,preweighted=False)
solnXRayPotTools.useGlobs(xray)
xray.setNumAngles(50)
xrayCorrect.setNumAngles(500)
xray.setScale(400)
xray.setCmpType("plain")
potList.append(xray)
crossTerms.append(xrayCorrect)
print xray.calcEnergy()
from solnScatPotTools import fitParams
#fitParams(xrayCorrect)
#xray.calcGlobCorrect(xrayCorrect.calcd())
#corrects I(q) to higher accuracy, and include solvent contribution corrections
rampedParams.append( StaticRamp("fitParams(xrayCorrect)") )
rampedParams.append( StaticRamp("xray.calcGlobCorrect(xrayCorrect.calcd())") )
Rgyr=AvePot(XplorPot,'COLL')
protocol.initCollapse("not resname ANI",
Rtarget=46.4) # approximate value from Guinier analysis
crossTerms.append(Rgyr)
#rap term to keep members of ensemble in similar orientations
#from posRMSDPotTools import RAPPot
#rap = RAPPot("rap","((name CA and resid 2:250 and segid A) or (name ca and
resid 2:250 and segid B))")
#rap.setScale( 100.0 )
#rap.setPotType( "square" )
#rap.setTol( 25.0 )
#potList.append(rap)
#crossTerms.append(rap)
#measure of difference between ensemble members
#posRMSD = RAPPot("posRMSD","""((name ca and resid 2:250 and segid A) or
# (name ca and resid 2:250 and segid B) )""")
#crossTerms.append(posRMSD)
from orderPotTools import create_OrderPot
orderPot = create_OrderPot("s2","order.tbl")
orderPot.setScale(10000)
orderPot.setPotType('square')
orderPot.setShowAllRestraints(1)
#potList.add( orderPot )
crossTerms.append(orderPot)
from residueAffPotTools import create_ResidueAffPot
contact = create_ResidueAffPot("contact")
potList.append(contact)
protocol.initRamaDatabase()
rama= AvePot(XplorPot,'RAMA')
potList.append(rama)
rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
bond= AvePot(XplorPot,"BOND")
refBondViolations=bond.violations() # don't wan't more violations than this
print 'refBondViolations:', refBondViolations
potList.append( bond)
angl= AvePot(XplorPot,"ANGL")
potList.append( angl)
rampedParams.append( MultRamp(0.4,1,"potList['ANGL'].setScale(VALUE)") )
impr= AvePot(XplorPot,"IMPR")
potList.append( impr )
rampedParams.append( MultRamp(0.1,1,"potList['IMPR'].setScale(VALUE)") )
protocol.massSetup()
from atomSelAction import SetProperty
nonContactNtermA.apply( SetProperty("mass",
100./len(nonContactNtermA)) )
nonContactNtermB.apply( SetProperty("mass",
100./len(nonContactNtermB)) )
from ivm import IVM
dynRigid = IVM()
from atomSel import union
dynRigid.fix(union("(%s) and (%s)" % (cTerminus,bboneAtoms),
nonContactCterm))
dynRigid.group(union("segid A and (%s) and (%s)"% (nTerminus,bboneAtoms),
nonContactNtermA))
dynRigid.group(union("segid B and (%s) and (%s)"% (nTerminus,bboneAtoms),
nonContactNtermB))
dynRigid.breakBond("segid A and (%s)" % linker2Break)
dynRigid.breakBond("segid B and (%s)" % linker2Break)
protocol.cartesianTopology(dynRigid,sel=linker2Sel)
protocol.torsionTopology(dynRigid)
from simulationTools import AnnealIVM
init_t = 3000
cool = AnnealIVM(initTemp =init_t,
finalTemp=25,
tempStep =25,
ivm=dynRigid,
rampedParams = rampedParams)
def calcOneStructure(loopInfo):
""" this function calculates a single structure, performs analysis on the
structure, and then writes out a pdb file, with remarks.
"""
minEnergy = 1e30 # big number
startPos = xplor.simulation.atomPosArr()
dockFilename='dock_%d_0.fit' % loopInfo.structNum
protocol.initCoords(dockFilename)
InitialParams( rampedParams )
InitialParams( highTempParams )
#initial minimization
protocol.initMinimize(dynRigid,
potList=potList,
numSteps=1000, # whichever comes first
printInterval=100)
dynRigid.run()
#high-temp dynamics
protocol.initDynamics(dynRigid,
potList=potList,
bathTemp=init_t,
initVelocities=1,
finalTime=800, # stops at 800ps or 8000 steps
numSteps=8000, # whichever comes first
printInterval=100)
from cdsVector import sum
for i in range(esim.size()):
print "member velocity sum:",sum(esim.members(i).atomVelArr())
pass
dynRigid.setETolerance( init_t/100 ) #used to det. stepsize. default:
t/1000
dynRigid.run()
# initialize parameters for cooling loop
InitialParams( rampedParams )
# initialize integrator for simulated annealing
#
protocol.initDynamics(dynRigid,
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(dynRigid)
dynRigid.run()
#do analysis and write structure
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=pdbTemplate,
structLoopAction=calcOneStructure,
genViolationStats=True,
averageRegularize=False,
averagePotList=potList,
averageCrossTerms=crossTerms,
averageSortPots=[rdcs,xray],
averageTopNum=20, #report on the top 20 structs
averageFitSel="name CA and (%s)" % cTerminus,
averageContext=FinalParams(rampedParams)).run()
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Processed by Mailcrypt 3.5.9 <http://mailcrypt.sourceforge.net/>
iEYEARECAAYFAkzr6J0ACgkQPK2zrJwS/lZ0wgCdEWFBj01H4Hj9UfNWgEv8A215
NrMAn3kvrO+WK4AxBcQ1BkkaF4AptWrE
=A8WH
-----END PGP SIGNATURE-----
_______________________________________________
Xplor-nih mailing list
[email protected]
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
