Hi Yikan,
I notice that your script seems to be more complicated than it has to be. I'm
not sure this will solve your problem, but I'm attaching a simpler script based
on eginput/rna/fold.py. It's basically the script used to fold an RNA
molecule from an extended conformation (described in the paper I mentioned
before). I modified it so that you load your input pdb structure and define
the different rigid bodies (look for the lines with an inline comment "#
COMPLETE !", where you must enter the right information). Then, it does
torsion angle dynamics using the remaining torsional degrees of freedom (i.e.,
at the linker(s)).
I commented all the experimental restraints. If you want to have the radius of
gyration term (which appears in your script), you can introduce it yourself.
Also, you might want to use it as a starting point for more complicated stuff,
like "breaking" the linkers and randomizing the position of the domains.
A suggestion to try to avoid (and detect) knotted structures is to increase the
atomic radii (see the 'repel' argument in initRepel).
Best,
Guillermo
________________________________
From: ZhangYikan <[email protected]>
Sent: Monday, December 11, 2017 10:43 AM
To: Bermejo, Guillermo (NIH/CIT) [E]; [email protected]
Subject: Re: [Xplor-nih] Wrong topology structures in conformational sampling
hello Guillermo,
I want these domains move (as rigid body)randomly and freely to sample all
possible conformations of the RNA.
Yikan
Ph. D candidate
TsingHua University
School of life Sciences
Beijing, China
tel:010-62773783
在 2017年12月11日，下午11:27，Bermejo, Guillermo (NIH/CIT) [E]
<[email protected]<mailto:[email protected]>> 写道：
Hi Yikan,
Can you tell me what you are trying to do?
Best,
Guillermo
________________________________
From: ZhangYikan <[email protected]<mailto:[email protected]>>
Sent: Friday, December 8, 2017 6:56 PM
To: Bermejo, Guillermo (NIH/CIT) [E];
[email protected]<mailto:[email protected]>
Subject: Re: [Xplor-nih] Wrong topology structures in conformational sampling
hi Guillermo,
yes, i have received your last mails. And the attachment is my scripts.
Yikan
Ph. D candidate
TsingHua University
School of life Sciences
Beijing, China
tel:010-62773783
在 2017年12月9日，上午12:43，Bermejo, Guillermo (NIH/CIT) [E]
<[email protected]<mailto:[email protected]>> 写道：
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
________________________________
From:
[email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
on behalf of ZhangYikan <[email protected]<mailto:[email protected]>>
Sent: Friday, December 8, 2017 6:53 AM
To: [email protected]<mailto:[email protected]>
Subject: [Xplor-nih] Wrong topology structures in conformational sampling
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
tel:010-62773783
#
# Based on eginput/rna/fold.py. 11 Dec. 2017.
#
#
# Total number of structures.
#
nstructures = 200
#
# Base name for output PDB files.
# ---
# This string must contain the "STRUCTURE" literal to be replaced by the
# structure number in the PDB filename. The (optional) "SCRIPT" literal is
# replaced by the name of this file (or stdin if redirected using <).
#
outfilename = 'SCRIPT_STRUCTURE.sa'
#
# Set random seed.
#
import protocol
protocol.initRandomSeed(3421) # by specific seed number
# Load coordinates.
protocol.loadPDB('initial.pdb' , deleteUnknownAtoms=True) # COMPLETE!!!
#
# Create a PotList() to contain the energy terms that will be active
# during structure calculations.
#
from potList import PotList
potList = PotList()
#
# Lists highTempParams and rampedParams will hold simulationTools.StaticRamp and
# simulationTools.MultRamp objects to handle parameter changes between the high
# termperature and annealing stage, and within the annealing stage (e.g., ramped
# force constants).
#
from simulationTools import StaticRamp, MultRamp, InitialParams, AnnealIVM
highTempParams = []
rampedParams = []
# Below, the entire setup of each energy term (including their treatment during
# the high temperature and annealing stages) is performed in self-contained
# sections, so that removal of a term or addition of a new one can be done
# simply by commenting out or adding the corresponding section, respcetively.
###
### Set up distance restraint potential (e.g., from NOEs).
###
##import noePotTools
##noe = PotList('noe')
##for (name, scale, table) in [('all', 1, 'noe_all.tbl'),
## # add entries for more restraint tables here
## ]:
## pot = noePotTools.create_NOEPot(name, table)
## # pot.setPotType("soft") # if you suspect there are bad NOEs
## pot.setScale(scale)
## noe.append(pot)
##potList.append(noe)
##rampedParams.append(MultRamp(2, 30, "noe.setScale( VALUE )"))
##
##
###
### Set up torsion angle restraint potential (e.g., from J-couplings).
###
##import xplorPot
##dihedralTables = ['dihed_new.tbl', 'gamma.tbl']
##protocol.initDihedrals(dihedralTables)
##potList.append(xplorPot.XplorPot('CDIH'))
##highTempParams.append(StaticRamp("potList['CDIH'].setScale(10)"))
##rampedParams.append(StaticRamp("potList['CDIH'].setScale(200)"))
##
##
###
### Set up potential for base-pair planarity restraints.
###
##import xplorPot
##protocol.initPlanarity('plane.tbl')
##potList.append(xplorPot.XplorPot('PLAN'))
### (The setup of this term remains unchanged throughout; no need to involve
### highTempParams and/or rampedParams.)
#
# Set up statistical torsion angle potential (torsionDB).
#
import torsionDBPotTools
torsiondb = torsionDBPotTools.create_TorsionDBPot(name='torsiondb',system='rna')
potList.append(torsiondb)
rampedParams.append(MultRamp(0.5, 4, "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 C1' atoms
highTempParams.append( StaticRamp("""initRepel(repel,
use14=True,
scale=0.004,
repel=1.2,
moveTol=45,
interactingAtoms="name C1'"
)""") )
# Selected 1-4 interactions.
# ---
# 1-4 interactions are ommited from the above term (repel) because they are
# mostly affected by torsionDB (torsiondb). However, torsionDB doesn't affect
# torsions of terminal, protonated groups (e.g., methyl rotation). Thus, 1-4
# interactions involved in such torsions are set up here. (Note: a more elegant
# solution might be implemented in the near future.)
#
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)"))
#
# Set up bond length potential.
# (Needed even if no Cartesian minimization is used, for "broken" rings.)
#
potList.append(xplorPot.XplorPot('BOND'))
# (The setup of this term remains unchanged throughout; no need to involve
# highTempParams and/or rampedParams.)
#
# Set up bond angle potential.
# (Needed even if no Cartesian minimization is used, for "broken" rings.)
#
potList.append(xplorPot.XplorPot('ANGL'))
rampedParams.append(MultRamp(0.4, 1.0, "potList['ANGL'].setScale(VALUE)"))
#
# Set up improper dihedral angle potential.
# (Needed even if no Cartesian minimization is used, for "broken" rings.)
#
potList.append(xplorPot.XplorPot('IMPR'))
rampedParams.append(MultRamp(0.1, 1.0, "potList['IMPR'].setScale(VALUE)"))
#
# Done with energy terms.
#
# Give atoms uniform weights, except for anisotropy axes (if any).
#
protocol.massSetup()
#
# Set up IVM object(s).
#
# IVM object for torsion-angle dynamics/minimization.
import ivm
dyn = ivm.IVM()
# Define rigid domains.
dyn.group('resid 1:15') # COMPLETE !!! domain i
dyn.group('resid 20:55') # COMPLETE !!! domain j
# Argument flexRiboseRing below is a string that selects residues whose ribose
# rings will have all endocyclic angles flexible. In general, all ribose rings
# should be selected. In this example, the non-RNA ligand (residue 34) has to
# be excluded.
protocol.torsionTopology(dyn, flexRiboseRing='all')
### Optional IVM object for final Cartesian minimization.
##minc = ivm.IVM()
##protocol.cartesianTopology(minc)
#
# Temperature set up.
#
temp_ini = 3500.0 # initial temperature
temp_fin = 25.0 # final temperature
def calcOneStructure(loopInfo):
"""Calculate a structure.
"""
# Generate initial structure by randomizing torsion angles.
import monteCarlo
monteCarlo.randomizeTorsions(dyn)
# Fix up covalent geometry.
# (The torsion restraints may include ring torsions and distort geometry.)
while True:
try:
protocol.fixupCovalentGeom(maxIters=100, useVDW=1)
break
except protocol.CovalentViolation:
pass
#
# High Temperature Dynamics Stage.
#
# Initialize parameters for high temperature dynamics.
InitialParams(rampedParams)
InitialParams(highTempParams) # purposedly overides some
# setups in rampedParams
# Set up IVM object and run.
protocol.initDynamics(dyn,
potList=potList,
bathTemp=temp_ini,
initVelocities=1,
finalTime=15, # run for finalTime or
numSteps=15001, # numSteps * 0.001, whichever is less
printInterval=100)
dyn.setETolerance(temp_ini/100) # used to find stepsize (default: temp/1000)
dyn.run()
#
# Simulated Annealing Stage.
#
# Initialize parameters for annealing.
InitialParams(rampedParams)
# Set up IVM object for annealing.
protocol.initDynamics(dyn,
potList=potList,
finalTime=0.2, # run for finalTime or
numSteps=201, # numSteps * 0.001, whichever is less
printInterval=100)
# Set up cooling loop and run.
AnnealIVM(initTemp=temp_ini,
finalTemp=temp_fin,
tempStep=12.5,
ivm=dyn,
rampedParams=rampedParams).run()
#
# Torsion angle minimization.
#
protocol.initMinimize(dyn,
potList=potList,
printInterval=50)
dyn.run()
## #
## # Cartesian minimization (optional).
## #
## protocol.initMinimize(minc,
## potList=potList,
## dEPred=10)
## minc.run()
from simulationTools import StructureLoop
StructureLoop(numStructures=nstructures,
pdbTemplate=outfilename,
doWriteStructures=True,
structLoopAction=calcOneStructure,
# Arguments for generating structure statistics:
genViolationStats=True,
averageSortPots=potList, # terms for structure sorting
averageTopFraction=0.1, # top fraction of structs. to report on
averagePotList=potList, # terms analyzed
averageFitSel='not (name H* or resname ANI)', # selection to fit
).run() # to average struct.
# and report precision
_______________________________________________
Xplor-nih mailing list
[email protected]
https://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
