Hi all, I used an adapted version of saxs_EI/dock.py to setup an run with RDC and SAXS data. I am struggling to get the desired results. I have only recently started using xplor-nih and I am not sure whether my script is sound.
Basically, what I would like to achieve is to calculate structures of a two domain protein, whereby the domain-domain orientations are defined by RDCs and the overall shape information derived from SAXS. I should mention that the inter-domain linker comprises 180 residues and there are barely any experimental constraints for the linker except for a few RDCs (and SAXS). Crystal structures of the two domains have been solved previously - I have made use of them and modelled the inter-domain linker into one pdb to create a tandem domain construct. Briefly, the N-, linker and C-terminal domains comprise of residues 40:168, 169:348 and 349:460 respectively. I have stumbled across saxs_EI/dock.py and believed it to be a good starting point for my calculations. The flow I am trying to follow is described below and the script is attached at the end of this message: - Randomize N-terminal domain position (keep C-terminal domain fixed). - Calculate C-terminal alignment tensor and keep it fixed. - Rotate N-terminal domain until the RDCs are satisfied. - Keep domains fixed; only allow translation of the domains. - High-temperature dynamics on the linker (treat domains as rigid bodies) - Refine structure using SAXS. Is there a better way to do this? Or a better template than saxs_EI/dock.py? When I run dock.py, I have to set cutoffRfactor at a rather 'high' value (>32) for the simulation to pass the RDC step but then I end up with a structure for the linker that is not meaningful, i.e., numerous 'unusual bonds', steric clashes, etc. The domains haven't changed during the simulation, as anticipated, but they're not aligned according to the RDCs as I would have hoped. It also appears as though the structure is not refined against the SAXS data too well. In the case where cutoffRfactor is set at a more reasonable value of 21, the simulation is stuck in an endless cycle trying to find the 'right' tensor. I'd appreciate any feedback. Many thanks. Romel Postdoctoral Fellow, Protein Structure and Biophysics ___________________________________________ AstraZeneca R&D | Discovery Sciences 50F23 Mereside, Alderley Park, Macclesfield, SK10 4TG, UK Tel: +44 (0)1625 518527 Fax: +44 (0) 1625 232693 [email protected]<mailto:[email protected]> import os print 'process id:', os.getpid() numberOfStructures=1 import protocol protocol.initRandomSeed(42) ensembleSize=1 cutoffRfactor=21.0 # need an R-factor better than this to accept a docked struct numDock=1 pdbTemplate="SCRIPT_STRUCTURE_MEMBER.fit" protocol.loadPDB("tandem_em.pdb",deleteUnknownAtoms=True) linker2Sel = "resid 169:348" #179 residue linker between the domains linker2Break = "(resid 348 and name c) or (resid 349 and name n)" nTerminus = "resid 40:168" # N-terminal domain cTerminus = "resid 349:460" # C-terminal domain 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") rdcs=PotList("rdc") from rdcPotTools import create_RDCPot rdc1 = create_RDCPot("rdc1","tandem_rdc_C.tbl",ten1) calcTensor(ten1) #this to set Da, Rh rdc1 = create_RDCPot("rdc1","tandem_rdc_N.tbl",ten1) rdc1.addRestraints(open("tandem_rdc_N.tbl").read()) rdc1.setShowAllRestraints(1) rdcs.append(rdc1) potList.append(rdcs) from avePot import AvePot from xplorPot import XplorPot from xplorPot import XplorPot protocol.initDihedrals("tandem_dihed_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)") ) nbond=AvePot(XplorPot,"VDW") potList.append(nbond) rampedParams.append( StaticRamp("""protocol.initNBond()""") ) 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)""") ) from solnXRayPotTools import create_solnXRayPot import solnXRayPotTools xray=create_solnXRayPot('xray',experiment='0691-2-1-zc.dat', numPoints=26, normalizeIndex=-3,preweighted=False) xrayCorrect=create_solnXRayPot('xray-c',experiment='0691-2-1-zc.dat', numPoints=26, normalizeIndex=-3,preweighted=False) solnXRayPotTools.useGlobs(xray) xray.addRigidRegion(cTerminus) xray.addRigidRegion(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 #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=70.3) # approximate value from Guinier analysis potList.append(Rgyr) from residueAffPotTools import create_ResidueAffPot contact = create_ResidueAffPot("contact") crossTerms.append(contact) # Statistical torsion angle potential # from torsionDBPotTools import create_TorsionDBPot torsionDBPot = create_TorsionDBPot('tDB') potList.append( torsionDBPot ) rampedParams.append( MultRamp(.002,2,"torsionDBPot.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 ivm import IVM dynRigid = IVM() dynRigid.fix(cTerminus) dynRigid.group(nTerminus) dynRigid.breakBond(linker2Break) protocol.torsionTopology(dynRigid) ten1.setFreedom("fix") dynTranslate = IVM() # this one fixes the orientation of the a and ab domains dynTranslate.fix(cTerminus) dynTranslate.group(nTerminus) dynTranslate.hinge("translate",nTerminus) dynTranslate.breakBond(linker2Break) protocol.cartesianTopology(dynTranslate) dynLinkerOnly = IVM() dynLinkerOnly.fix(cTerminus) dynLinkerOnly.fix(nTerminus) dynLinkerOnly.breakBond(linker2Break) protocol.cartesianTopology(dynLinkerOnly) 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 randomizeDomainPos( 'not (%s)' % cTerminus, deltaPos=45 ) from monteCarlo import randomizeTorsions randomizeTorsions(dynRigid) #fit RDCs protocol.initMinimize(dynRigid, potList=[rdcs], numSteps=1000) dynRigid.run() ; calcTensor(ten1); dynRigid.run(); calcTensor(ten1) #fix broken bond # if bond is excluded, linker2 can extend a long distance # to make an unphysical connection protocol.initMinimize(dynRigid, potList=[bond,rdcs], numSteps=1000) bond.setScale(1e-5) ; dynRigid.run() ; calcTensor(ten1) bond.setScale(1e-3) ; dynRigid.run() ; calcTensor(ten1) bond.setScale(1) ; dynRigid.run() ; 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,rdcs,Rgyr], numSteps=1000) dynTranslate.run() ; dynTranslate.setVerbose( 0 ) InitialParams(rampedParams) #done for xray-correct fitting protocol.initMinimize(dynTranslate, potList=[bond,rdcs,Rgyr,xray], numSteps=1000) dynTranslate.run() InitialParams(rampedParams) #done for xray-correct fitting #attempt to fix linker geom try: protocol.fixupCovalentGeom(sel=linker2Sel) except protocol.CovalentViolation: pass print 'number of bond violations:', bond.violations() if bond.violations()>refBondViolations: print 'too many bond violations. Discarding this configuration' continue # try fixing everything together protocol.initMinimize(dynTranslate, potList=[bond,angl,impr, rdcs,xray], numSteps=1000) dynTranslate.run() ; dynTranslate.run() InitialParams(rampedParams) # refit xray-correct before scoring scorePot = PotList() for p in (xray,rdcs,bond): 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=80, # stops at 800ps or 8000 steps numSteps=800, # 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() #do analysis and write structure loopInfo.writeStructure(potList) pass from simulationTools import StructureLoop StructureLoop(numStructures=numberOfStructures, pdbTemplate=pdbTemplate, structLoopAction=calcOneStructure, genViolationStats=True, averageRegularize=False, averagePotList=potList, averageCrossTerms=crossTerms, averageFitSel="name CA and (%s)" % cTerminus, averageSortPots=[xray,XplorPot('VDW')], averageTopFraction=0.5, #report on all structs averageContext=FinalParams(rampedParams)).run() ________________________________ AstraZeneca UK Limited is a company incorporated in England and Wales with registered number: 03674842 and a registered office at 2 Kingdom Street, London, W2 6BD. Confidentiality Notice: This message is private and may contain confidential, proprietary and legally privileged information. If you have received this message in error, please notify us and remove it from your system and note that you must not copy, distribute or take any action in reliance on it. Any unauthorised use or disclosure of the contents of this message is not permitted and may be unlawful. Disclaimer: Email messages may be subject to delays, interception, non-delivery and unauthorised alterations. Therefore, information expressed in this message is not given or endorsed by AstraZeneca UK Limited unless otherwise notified by an authorised representative independent of this message. No contractual relationship is created by this message by any person unless specifically indicated by agreement in writing other than email. Monitoring: AstraZeneca UK Limited may monitor email traffic data and content for the purposes of the prevention and detection of crime, ensuring the security of our computer systems and checking compliance with our Code of Conduct and policies.
_______________________________________________ Xplor-nih mailing list [email protected] https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
