[ Plot simple | Plot histogram | Plot 3D 2Dfunction | Plot 3D shape | Peptide docking ]
How to make a simple plot y=F( x) |
If your function is determined on an uneven set of abscissa values,
stored in array x, then use:
plot display x y
If your function is determined on a regular 1D-grid of values varying
from x0 to x1, and is stored in the array y, use:
plot display Rarray(Nof(y), x0, x1) y
Use on-line help:
help plot
or to see the concise list of arguments
help command plot
Function Histogram returns the 2*N matrix which may be used on the
fly with a command like
plot display BAR red Histogram(a 50)
More examples:
a=Random(0. 100. 10000)
b=.04*(Count(1 50)*Count(1 50))
annot = {"Histogram with quadratic ruler" "Random value" "N"}
plot display BAR red Histogram(a b) annot
b=Sqrt(100.*Count(1 100))
annot = {"Histogram with square root ruler" "Random value" "N"}
plot display green BAR Histogram(a b) annot
How to make a 3D-surface plot of a 2D-function |
In the most common case you may have a matrix N*M my_matrix
representing a two-dimensional function of x and y
real arrays.
First create a
graphics object,
then display it and save as a picture if you wish:
read matrix "ram"
# let's try a "wire" representation first
x=Rarray(Nof(ram),-180.,180.) #this scale may be nonlinear
y=x
make grob "matwire" ram x y # scales for x,y. Def. z-scale=1.
display matwire
# now let's have a look at a solid surface
make grob solid "matsolid" ram x y 0.9 # scales for x,y and z
display matsolid solid
# adjust the view by rotating/translating
# connect matsolid # you may move them with respect to each other
# and save image
write image rgb "img_matsolid"
How to create a new graphics object of a specific shape |
First, use one of the provided
graphics objects,
stored in
*.gro
files to try these kind of ICM non-molecular objects. You can scale
graphics objects by multiplying them by a real number and shifting (e.g. g1=10.*g + {1.,0.,0.}) and/or
adjust their positions either by
translate or
rotate
commands or by
connecting
to the selected
grob and dragging/rotating them with the mouse in the
graphics window, for example:
read grob "icos"
my_icos = g_icos*0.2
display my_icos g_icos
translate add my_icos {5. 15. 0.}
You may also create a new or edit an existing graphics objects (see
provided
*.gro
files to learn about format).
Let us consider a problem of flexible peptide docking to a receptor.
The peptide may be completely free or constrained (e.g. have a helical
fragment in the middle), and the receptor will be assumed to have
a fixed conformation. We will be going through a series of interactive
or non-interactive steps. Feel free to put the commands below into
several scripts and run them sequentially.
We will follow a particular example through all the steps. In this example
we will be docking a constrained 11-residue alpha-helical peptide to a
Bcl-XL molecule. It is a good idea to create a separate directory for the
project, say, projectX , to keep all the files in one place.
Preparing the receptor binding surface for peptide docking.
First, we need to get the receptor model as an ICM object.
It may come from several sources.
The most common one is a pdb-structure converted with the convertObject macro,
or the build model command which builds a model by homology, followed
by regularization and refinement with the refineModel macro.
We will create two objects:
a_1.|the original pdb object with receptor with some other peptide
|
a_2.|converted receptor (without the peptide)
|
Example:
% cd projectX
% icmgl -G
nice "1g5j" # an open form of the receptor
delete a_w* # delete waters but keep bound peptide
copy a_ "rec"
delete a_rec.2 # keep only the main chain
undisplay window
convertObject a_rec. yes yes yes no # also optimizes hydrogens
#
display box Box( a_1.2 4. ) # display a box around the peptide
Now you may rotate the scene and adjust the box size by pulling a box corner
with the leftMouseClick.
The box needs to cover the receptor surface you need to include in the simulation,
but needs to be minimized to reduce the computational burden and memory use.
Sometimes the orientation of the box with respect to the receptor surface is not
favorable. In that case, use connect a_*. and rotate the objects with
the respect to the origin ( use axisLength = 15. and display origin to
to see the axis in addition to seeing the box ). If you change the absolute position
of the receptor model, you may want to save the new orientation of the pdb object
(e.g. write pdb a_1. "template" )
The result
Now we have two objects: (i) a pdb object in which we deleted the unnecessary molecules,
but kept the peptide for comparison and convenience in the box definition;
and (ii) an ICM-object with a receptor for which we will be calculating the grid maps.
Create grid potential maps
Now calculate five grid potentials with the make map potential command:
make map potential "gc,gh,ge,gs,gb" Box() 0.5
This command uses all molecules in the current object as a source
for the map calculations.
The Box() function will return the box you see on your screen.
Alternatively you can just provide the 6 numbers defining that box.
0.5 is the grid size.
Even though the ICM grid potential is extremely fast compared to other programs,
for large boxes you may need to increase the grid cell size. We do not advise to
go beyond 0.7A, however.
The result of this step is five grid maps around the docking surface of interest.
Building the peptide
The peptide of interest can be build with the following simple command:
build string "ECLKRIGDELD" # peptide sequence from Bax
rename a_ "pep"
set a_/* "HHHHHHHHHHH"
superimpose a_1.2/310:320/ca a_/1:11/ca align
Saving the results
Let us now save both the objects and the maps for future use.
write m_gc m_gh m_ge m_gs m_gb # gc.map .. files are created
write object a_*. "all" # write 3 objects in one multi-object file
Running a simple simulation from one starting conformation
The actual simulation does not need to be run interactively.
An ICM script file _dockpep is described below. Run the script by
% _dockpep >! f1.out # the output file
The file contains the following commands:
#!/usr/local/bin/icmng
read libraries
read object "all" # contains a_
set object a_pep. # the peptide is the current object
fix v_/3:10/phi,psi,omg # fix some backbone angles
nvar = Nof( v_//phi,psi,H,P ) # number of essential variables
# let us choose iteration limits depending on nvar
mncallsMC = 10000 + Integer(0.008*nvar*nvar*nvar*nvar*nvar)
mncalls = 170+nvar*3
print mncallsMC mncalls
# make sure that these two numbers are not insane
# or set smaller mncallsMC and mncalls (e.g. 10000 and 100) as a test
temperature = 600 # optimal temperature for the simulation
tolGrad = 0.001 # exit minimization when gradient is < 0.001
mnhighEnergy = 15 + nvar/2 #
set vrestraints a_/* # load preferred backbone and side-chain angle zones
#
vicinity = 2.5
compare a_pep. static # use sRmsd of the peptide
set terms "vw,14,hb,el,to,en,gc,gh,ge,gb,gs"
minimizeMethod = "conjugate"
#
# maps
read map "gc,gh,ge,gb,gs"
GRID.margin = Distance( a_pep./1/ca a_pep./11/ca )
#
# run montecarlo from a single start
# the reverse option will make more reasonables moves
#
# impose tethers ?
montecarlo reverse trajectory
# run montecarlo from multiple starts and merge the stacks
for i=1,Nof(<starting conditions>)
montecarlo reverse trajectory append
endfor
This simulation created several files:
f1.cnf | a stack file with mnconf (or less) best non-redundant conformations
|
f1.trj | a trajectory file with all accepted conformations (if you used the trajectory option)
|
f1.ou|the file with the text output of the script
|
The output file f1.ou contains information about every accepted conformation of the simulation, e.g.
...
__ __ 300 5 asp BPMC db da da -32.42 -23.66 100 1.29 97810
DY Visi 300 1 phe BPMC fb fb fb -32.42 -32.42 101 0.57 97911
__ __ 300 4 ser BPMC sa sg sd -32.42 -24.45 88 0.06 97999
DY Visi 300 5 asp BPMC dmn dmn dmn -32.42 -32.43 49 0.08 98048
...
Analyzing the results of a single simulation
Energy profile: Is the simulation long enough
The energies can be plotted to analyze the progression of the global energy optimization.
plotEnergy "f1.ou" 50.
plotBestEnergy "f1.ou" 50. "display append"
The first macro plots energies of each conformation generated and accepted in
the course of the ICM stochastic global energy optimization. Due to the nature
of the procedure, the energy may go up and down. The second macro plots
(and appends the plot to the previously generated def.eps file) only the
lowest energy achieved by a given iteration. This plot shows if and when a
better conformation is found. Naturally, the energy can not be improved after
the global minimum is found.
However, having a long flat plateau does not guarantee that the global minimum is found
since a new drop of energy may come after a very long simulation.
The time required for finding the global minimum in general depends on
the number of variables, number of atoms (the cost of a single energy calculation)
and complexity of the potential energy surface.
The best empirical way to make sure that the minimum is found is to compare
the lowest energies found from completely different starting conformations
(read about multiple start simulations below).
We use a window of 50 kcals/mole from the lowest energy. It helps to avoid
the dependence of the scale on the initial energy which may be quite high.
To analyze the improvements only use plotBestEnergy .
Graphical analysis of the results
You can generate all stack conformations as independent ICM objects with the
mkStackConf macro,e.g.
read object "all" # three objects, the peptide is the current object
mkStackConf 1 10 # make ICM objects from 10 best energy structures
The objects will be called a_pep1. , a_pep2. , etc.
Now these objects can be displayed simultaneously or separately and
further analyzed, e.g.
display a_pep*.
color a_pep?*. Energy( stack ) # in case you extracted all confs
You can loop through the conformations interactively either with
the original stack, e.g.
# create some view you like
for i=1,10 # or Nof(conf)
load conf i
pause # rotate and zoom, hit Return
endfor
, or use the objects you created with the mkStackConf macro.
Running a multiple start simulation
You can figure it out yourself :-) .
Convergence analysis
The key question we would like to ask is if at least two independent
simulations from random starting conformations identified a nearly identical
conformation with a similar energy as the lowest energy conformation.
A low-tech example with just two simulations f1 and f2 :
read object "all"
load conf 0 "f1" # 0 is the lowest energy conformaiton.
load conf 0 "f2"
You can also plot the best energy (see above) and compare the plots, as well
as perform graphical analysis to see if the best conformations are similar.
Analyzing the results of a multiple start simulation, merging stacks
Use the read stack append command to
merge all stack conformations together. Now you can
redefine the compare command and the vicinity parameter
depending on how you want to further compare and filter out
the accumulated conformations.
The compress command performs the compression.