| hi Guillermo, yes, i have received your last mails. And the attachment is my scripts.
|
#dock EIN to EIC
# repeat numDock times:
# randomize domain positions
# enforce symmetry
# satisfy rdcs
# dock using xray term (possibly including Rgyr to collapse) -
# allowing only translation
#
# keep coordinates from best run
# 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 sys
import os
os.system(sys.argv[1])
print 'process id:', os.getpid()
numberOfStructures=100
import protocol
protocol.initRandomSeed(42)
ensembleSize=1
cutoffRfactor=18.22 # need an R-factor better than this to accept a docked struct
numDock=2
pdbTemplate="SCRIPT_STRUCTURE_MEMBER.pdb"
protocol.loadPDB(sys.argv[1],deleteUnknownAtoms=True)
groupedResids = [(1,71), ( 85,152), ( 167,306), ( 316,412) ]
breakAtoms = ["C5'", "O5'"]
linkerResids = sum( [[( groupedResids[i][1],groupedResids[i+1][0])] for
i in range(len(groupedResids)-1)] , [] )
breakResids = [int(0.5*(pair[0]+pair[1])) for pair in linkerResids]
linkerSel = " or ".join(["resid %d:%d" % pair for pair in linkerResids])
import regularize
regularize.fixupCovalentGeomIVM("not resid %d:%d" % groupedResids[2],
rigidRegions=["resid %d:%d" % pair for
pair in groupedResids[:2]+
groupedResids[3:]],
)
from ensembleSimulation import EnsembleSimulation
esim=EnsembleSimulation('ens',ensembleSize)
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("medium1",axis=axisAtoms)
#media.append(ten1)
#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)
#rampedParams.append( StaticRamp('for t in media: calcTensor(t)') )
from avePot import AvePot
from xplorPot import XplorPot
#
#from xplorPot import XplorPot
#protocol.initDihedrals("torsion_linker.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)") )
#
selPairs=[(linkerSel,"not PSEUdo")]
for i,pair1 in enumerate(groupedResids):
for pair2 in groupedResids[i+1:]:
selPairs.append( ("resid %d:%d"%pair1,"resid %d:%d"%pair2) )
pass
pass
#selPairs=[(linkerSel,"not PSEUdo")]
for i,pair3 in enumerate(groupedResids):
for j,pair4 in enumerate(linkerResids):
selPairs.append( ("resid %d:%d"%pair3,"resid %d:%d"%pair4) )
pass
pass
from repelPotTools import create_RepelPot,initRepel
repel = create_RepelPot('repel',
selPairs=selPairs,
selection=AtomSel('not PSEUDO'))
potList.append(repel)
#from repelPotTools import initRepel
#highTempParams.append( StaticRamp("""protocol.initRepel(repel,
# use14=True,
# scale=0.004,
# repel=1.2,
# moveTol=45,
# onlyCA=0)""") )
#
#from repelPotTools import create_RepelPot,initRepel
#rampedParams.append( StaticRamp("initRepel(repel,use14=False)") )
#rampedParams.append( MultRamp(.004,4, "repel.setScale(VALUE)") )
#
#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.addRigidRegion(cTerminus)
#xray.addRigidRegion("segid A and (%s)" % nTerminus)
#xray.addRigidRegion("segid B and (%s)" % nTerminus)
#
#
#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())") )
#
#from residueAffPotTools import create_ResidueAffPot
#contact = create_ResidueAffPot("contact")
#crossTerms.append(contact)
#
#protocol.initRamaDatabase()
#rama= AvePot(XplorPot,'RAMA')
#potList.append(rama)
#rampedParams.append( MultRamp(.002,1,"potList['RAMA'].setScale(VALUE)") )
import protocol
from xplorPot import XplorPot
protocol.initCollapse(sel='all',Rtarget=75, scale=0.005)
rGyr = XplorPot('COLL')
potList.append(rGyr)
# Statistical torsion angle potential
#
from torsionDBPotTools import create_TorsionDBPot
torsionDBPot = create_TorsionDBPot('tDB',system="rna")
potList.append( torsionDBPot )
rampedParams.append( MultRamp(.002,2,"torsionDBPot.setScale(VALUE)") )
#1-4repulsions to prevent eclipsed conformations
#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)"))
bond= AvePot(XplorPot,"BOND")
refBondViolations=bond.violations() # don't wan't more violations than this
print 'refBondViolations:', refBondViolations
print 'initial bond energy:', bond.calcEnergy()
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 ivm import IVM
dynRigid = IVM()
#dynRigid.fix("resid %d:%d" % groupedResids[2])
for pair in groupedResids:
dynRigid.fix("resid %d:%d" % pair)
pass
for resid in breakResids:
dynRigid.breakBond("resid %d and (name %s or name %s)" %
(resid,breakAtoms[0],breakAtoms[1]))
pass
protocol.cartesianTopology(dynRigid)
#ten1.setFreedom("fix")
dynTranslate = IVM() # this one fixes the orientation of the a and ab domains
#dynTranslate.fix("resid %d:%d" % groupedResids[2])
for pair in groupedResids:
dynTranslate.hinge("full","resid %d:%d" % pair)
dynTranslate.group("resid %d:%d" % pair)
pass
for resid in breakResids:
dynTranslate.breakBond("resid %d and (name %s or name %s)" %
(resid,breakAtoms[0],breakAtoms[1]))
pass
protocol.cartesianTopology(dynTranslate)
from simulationTools import AnnealIVM
init_t = 3000
#cool = AnnealIVM(initTemp =init_t,
# finalTemp=25,
# tempStep =25,
# ivm=dynTranslate,
# 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()
InitialParams(rampedParams)
k=0
while k < numDock:
xplor.simulation.setAtomPosArr(startPos)
from atomAction import randomizeDomainPos
for pair in groupedResids:
# for pair in linkerResids:
randomizeDomainPos( "resid %d:%d" % pair,
deltaPos=150 )
pass
# generate a new structure with randomized torsion angles
# from monteCarlo import randomizeTorsions
# randomizeTorsions(dynRigid)
# protocol.writePDB('init_%d.pdb'% loopInfo.structNum)
# #fit RDCs
# protocol.initMinimize(dynRigid,
# potList=[dSymm,ncs,rdcs],
# numSteps=1000)
#
# dynRigid.run() ; calcTensor(ten1); dynRigid.run(); calcTensor(ten1)
# for pair in linkerResids:
# randomizeDomainPos( "resid %d:%d" % pair,
# deltaPos=40 )
# pass
from monteCarlo import randomizeTorsions
randomizeTorsions(dynRigid)
# protocol.writePDB('init1_%d.pdb'% loopInfo.structNum)
#fix broken bond
# if bond is excluded, linker2 can extend a long distance
# to make an unphysical connection
protocol.initMinimize(dynRigid,
potList=[bond,],
numSteps=1000)
# protocol.writePDB('init2_%d.pdb'% loopInfo.structNum)
# bond.setScale(1e-5) ; dynRigid.run() ;
# calcTensor(ten1)
# bond.setScale(1e-3) ; dynRigid.run() ;
# calcTensor(ten1)
bond.setScale(1) ; dynRigid.run() ;
# protocol.writePDB('init2_%d.pdb'% loopInfo.structNum)
# from simulationTools import analyze
# print analyze(rdcs)
# from rdcPotTools import Rfactor_infinite
# rFac = Rfactor_infinite(rdcs[0])
# print 'rFactor: %.2f' % rFac
# if rFac>cutoffRfactor:
# print 'bad R factor. Discarding this configuration'
# continue
# #dock using Rgyr, saxs, vdw, contact?
# protocol.initMinimize(dynTranslate,
# potList=[bond,dSymm,ncs,rdcs,Rgyr],
# numSteps=1000)
#
## xray.setVerbose(2)
## dynTranslate.setVerbose( dynTranslate.printNodeDef )
# dynTranslate.run() ;
# dynTranslate.setVerbose( 0 )
#
# InitialParams(rampedParams) #done for xray-correct fitting
# protocol.initMinimize(dynTranslate,
# potList=[bond,],
# numSteps=1000)
# bond.setScale(1) ; dynTranslate.run()
#
# protocol.writePDB('init3_%d.pdb'% loopInfo.structNum)
InitialParams(rampedParams) #done for xray-correct fitting
#attempt to fix linker geom
# try:
# protocol.fixupCovalentGeom(sel=linkerSel)
# except protocol.CovalentViolation:
# pass
# print 'number of bond violations:', bond.violations()
# if bond.violations()>refBondViolations:
# print 'too many bond violations. Discarding this configuration'
# continue
# protocol.writePDB('init4_%d.pdb'% loopInfo.structNum)
# try fixing everything together
protocol.initMinimize(dynTranslate,
potList=[bond,angl,impr,repel,rGyr],
numSteps=1000)
repel.setScale(100); dynTranslate.run() ; dynTranslate.run()
# protocol.writePDB('thir_%d.pdb'% loopInfo.structNum)
# InitialParams(rampedParams) # refit xray-correct before scoring
scorePot = PotList()
for p in (bond,repel):
scorePot.append(p)
pass
scoreEnergy = scorePot.calcEnergy()
print 'dock iteration: %d score energy: %.2f' %(k,scoreEnergy)
if scoreEnergy < minEnergy :
minPos = xplor.simulation.atomPosArr()
minEnergy = scoreEnergy
pass
k += 1
pass
xplor.simulation.setAtomPosArr(minPos)
InitialParams( rampedParams )
InitialParams( highTempParams )
#high-temp dynamics
protocol.initDynamics(dynTranslate,
potList=potList,
bathTemp=init_t,
initVelocities=1,
finalTime=1600, # stops at 1600ps or 16000 steps
numSteps=16000, # whichever comes first
printInterval=100)
dynTranslate.setETolerance( init_t/100 ) #used to det. stepsize. default: t/1000
dynTranslate.run()
# # initialize parameters for cooling loop
# InitialParams( rampedParams )
#
#
# # initialize integrator for simulated annealing
# #
# protocol.initDynamics(dynTranslate,
# potList=potList,
# numSteps=100, #at each temp: 100 steps or
# finalTime=.2 , # .2ps, whichever is less
# printInterval=100)
#
# # perform simulated annealing
# #
# cool.run()
#
# # refit SAXS for final minimization
# FinalParams( rampedParams )
# final torsion angle minimization
#
protocol.initMinimize(dynTranslate)
dynTranslate.run()
pass
from simulationTools import StructureLoop
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=pdbTemplate,
doWriteStructures=True,
structLoopAction=calcOneStructure,
genViolationStats=True,
averageRegularize=False,
averagePotList=potList,
averageCrossTerms=crossTerms,
averageFitSel="resid %d:%d" % groupedResids[2],
averageTopFraction=0.5, #report on all structs
averageContext=FinalParams(rampedParams)).run()
Yikan Ph. D candidate TsingHua University School of life Sciences Beijing, China
Hi Yikan,
I realized that my original response was only sent to you and not included in the mailing list. Here it is.
One possibility is that there's something wrong with your script(s).
I recommend basing your calculations on the scripts found in eginput/rna directory (within your Xplor-NIH directory). There you'll find:
* fold.py, for initial calculations starting from an extended conformation. * refine.py, for subsequent refinement.
In addition, you'll find all the restraints for an example system, in case you want to play with the scripts before venturing into your data. The scripts are highly commented, which should help in customizing them for your needs. There's also a paper describing them [Bermejo et al., 2016, Structure 24, 806-815].
Best,
Guillermo
hello all,
Recently, I am doing conformational sampling with Xplor-NIH. And I got many pdbs, but when I check them I found that some pdbs have wrong topology structures like the attachments(in the pink cycle ): The RNA linker did not run with the right way and It went through the helix. I have tried to pick the structures by the “total energy” or “repel” terms, but it make no sense. Are there someone know how to kick out these wrong structures or to void these errors in the sampling process? Thank you! Yikan Ph. D candidate TsingHua University School of life Sciences Beijing, China
|
_______________________________________________
Xplor-nih mailing list
[email protected]
https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
