Copyright © 2020, Molsoft LLC
Nov 24 2024
| Prev | ICM Language Reference How to plot | Next |
[ 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 yIf 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) yUse on-line help:
help plotor to see the concise list of arguments
help command plot
Another way to make a file with a plot is to use the make plot command, e.g.
add column t Random(100,0.,1.) Random(100,0.,1.) Random(100,0.,1.)
make plot t "x=A;y=B;color=C;size=6;style=dots;;" output="aaa.png" size={500,500} # or pdf,eps,jpg
The type of file will be determined by the extension.
How to plot a histogram |
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).
Flexible peptide docking |
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 peptideNow 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 six grid potentials with the make map potential command:
make map potential "gc,gh,ge,gs,gb" Box() 0.5 # additional 'gl' grid is made for 'gc'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 fileThe 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 |
... __ __ 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 structuresThe 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 confsYou 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.
| Prev Density correlation | Home Up | Next Docking/VLS |