Hello, I have some questions about refinement a protein-DNA complex structure combine NMR restraints and SAXS data. We first used the script sry_final.inp (attachment) which is from the eginput/deprecated/sry example in the xplor-NIH distribution to calculated complex structure. In this step, the DNA bases are grouped as rigid-bodies. Then we want to add SAXS data into this script to refine the complex structure. So I found the script refine.py in the eginput/saxsref example. But these two scripts were written in different languages, the SAXS data module could not directly used in the rigid-body refinement script sry_final.inp.
So I want to how to add the SAXS data in the script sry_final.inp to refine protein-DNA complex, because in this script we can fix the DNA bases as rigid-bodies. If it counld not work, how to fix the DNA-bases in the SAXS data script refine.py. Or if there are some better solutions by which we can both fix DNA bases and introduce SAXS data in the refinement. Thanks! Best regards. Qianwen Li ######################################################################## To unsubscribe from the XPLOR-NIH list, click the following link: http://list.nih.gov/cgi-bin/wa.exe?SUBED1=XPLOR-NIH&A=1
sry_final.inp
Description: Binary data
xplor.requireVersion("2.25")
#
# 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
outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=100
if quick:
numberOfStructures=3
pass
# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(3421) #explicitly set random seed
#
# annealing settings
#
command = xplor.command
protocol.initParams("protein")
# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF('protG.seq')
#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("gb1_extended_%d.pdb" %
# protocol.initialRandomSeed())
# or read an existing model
#
protocol.loadPDB("1E8L.pdb",
model=49,swapProtons=True,fixMethylImpropers=True)
protocol.addUnknownAtoms()
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()
crossTerms = 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='1E8L.pdb',
cmpSel="not name H*")
crossTerms.append(refRMSD)
# 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 [ ('m1', 16.11, 0.31),
('m2', 13.32, 0.16) ]:
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 \
[('m1','NH' ,'rdc1.tbl' ,1),
('m2','NH' ,'rdc2.tbl' ,1),
]:
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.01,0.3, "rdcs.setScale( VALUE )") )
# calc. initial tensor orientation
# and setup tensor calculation during simulated annealing
#
from varTensorTools import calcTensorOrientation, calcTensor
for medium in media.keys():
calcTensor(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 [('dist',1,"noe-noviol.tbl"),
('hb' ,1,"hbond.tbl"),
#add entries for additional tables
]:
pot = create_NOEPot(name,file)
#pot.setPotType("soft") # - if you think there may be bad NOEs
pot.setScale(scale)
noe.append(pot)
pass
#noe['dist'].setPotType("soft") #some of the distance bounds are too small
rampedParams.append( MultRamp(2,100, "noe.setScale( VALUE )") )
# Set up dihedral angles
from xplorPot import XplorPot
protocol.initDihedrals("dihedral.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)") )
# hbdb - bb hbonds from database
#
protocol.initHBDB()
potList.append( XplorPot('HBDB') )
# Statistical torsion angle potential
#
from torsionDBPotTools import create_TorsionDBPot
torsionDBPot = create_TorsionDBPot('tDB')
potList.append( torsionDBPot )
rampedParams.append( MultRamp(.002,2,"torsionDBPot.setScale(VALUE)") )
#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList.append( XplorPot('VDW') )
rampedParams.append( StaticRamp("protocol.initNBond()") )
rampedParams.append( MultRamp(0.9,0.8,
"command('param nbonds repel VALUE end end')") )
rampedParams.append( MultRamp(.004,4,
"command('param nbonds rcon VALUE end end')") )
# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""protocol.initNBond(cutnb=100,
rcon=0.004,
tolerance=45,
repel=1.2,
onlyCA=1)""") )
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)") )
#SAXS term
from solnXRayPotTools import create_solnXRayPot
import solnXRayPotTools
xray=create_solnXRayPot('xray',experiment='lsz_5.0mgmL_saxswaxs.dat',
numPoints=50,
useInternalSpline=True,
normalizeIndex=-3,preweighted=False)
xrayCorrect=create_solnXRayPot('xray-c',experiment='lsz_5.0mgmL_saxswaxs.dat',
numPoints=50,
useInternalSpline=True,
normalizeIndex=-3,preweighted=False)
solnXRayPotTools.useGlobs(xray)
xray.setNumAngles(50)
xrayCorrect.setNumAngles(500)
xray.setScale(40)
xray.setCmpType("plain")
potList.append(xray)
crossTerms.append(xrayCorrect)
print xray.calcEnergy()
from solnScatPotTools import fitParams
#correct I(q) to higher accuracy, and include solvent contribution corrections
#stride=10 specifies that this fit is performed every 100th temperature during
#simulated annealing. During refinement, infrequent solvent correction updates
#are ok because the structure doesn't change too much.
#
rampedParams.append(
StaticRamp(
"fitParams(xrayCorrect);xray.calcGlobCorrect(xrayCorrect.calcd())",
stride=100))
# 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 =12.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=150, # stops at 150ps or 15000 steps
numSteps=15000, # 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=400, #at each temp: 400 steps or
finalTime=1 , # 1 ps, 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
loopInfo.writeStructure(potList)
pass
from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
pdbTemplate=outFilename,
structLoopAction=calcOneStructure,
genViolationStats=1,
averagePotList=potList,
# averageSortPots=[potList['BOND'],potList['ANGL'],potList['IMPR'],
# noe,rdcs,potList['CDIH']],
averageCrossTerms=crossTerms,
averageTopFraction=0.2, #report only on best 50% of structs
averageContext=FinalParams(rampedParams),
averageFilename="SCRIPT_ave.pdb", #generate regularized ave structure
averageFitSel="name CA",
averageCompSel="not resname ANI and not name H*" ).run()
