 ICM Manual v.3.8
by Ruben Abagyan,Eugene Raush and Max Totrov
Jun 9 2022

 Prev ICM Language Reference Structure analysis Next

# How to optimally superimpose two 3D structures

 Optimal superposition implies optimization of the Ca-RMSD upon rigid body superposition of the equivalent residues/atoms. This set of equivalent positions can be predefined, or determined by sequence alignment, or automatically derived from structure. In the latter case the optimized value is the RMSD of a trial alignment is corrected by the alignment length to reward longer alignment with slightly worse RMSD. The resulting transformation vector is returned in R_out (see the transform command) . There are several different algorithms which can be applied: • superimpose by known residue or atom equivalences (see the superimpose command)
• superimpose by sequence alignment which is calculated on the fly:
```
superimpose a_1.1 a_2.1 align  # the sequences are generated on the fly
superimpose a_1.1 a_2.1 align minimize  # even better method
```
This procedure fails if two structures have no significant sequence similarity
• align .. command to find global structural alignment (returned by the ali_out variable) and superimpose accordingly. Here we do not rely on sequence alignement (although it can be added to the optimized scoring function with a certain weight).
align a_1.1 a_2.1
• Align(.. distance) for local structural alignment

Use superimpose command. It performs an optimal rotation and translation of one structure onto the other. If necessary, a sequence alignment may be done prior to superposition by specifying align option in the command line.
Example:
```
display a_1.1//n,ca,c magenta
make sequence a_1.1
display a_2.1//n,ca,c blue
make sequence a_2.1
show sequence
# somewhat different sequences of two Zn-fingers
# sequence alignment is required
superimpose a_1. a_2. align
```
Note, in this particular example, the whole structural similarity is not so high. However, better fit may be obtained if only portions of the structures are superimposed, for example:
```
superimpose a_1.1/16:27/n,ca,c a_2.1/116:127/n,ca,c align
transform g_2_skinGrob R_out  # rearrange a related grob accordingly
```

# How to optimally superimpose without the residue alignment

The core of this procedure is the
align ms_molecule1 ms_molecule2
command. There are two variants: a fast superposition using dynamic programming algorithm align [distance] ms_1 ms_2 or a more rigorous, but somewhat less stable and slow align heavy ms_1 ms_2 ... command. This first command is well described above and identifies only the best superposition. The initial superposition is then refined similarly to the find alignment command.
The second algorithm (option heavy ) identifies a number of possible superpositions (solutions) based on the Ca atom coordinates only. The first solution is the best hit. See also load solution command.
Examples:
```
display a_*.//ca,c,n
color molecule a_*.
align a_1.1 a_2.1
center
color red  as_out
color blue as2_out
show ali_out
```

# How to make a Ramachandran plot

Use macro plotRamaEps . The macro is invoked by
Important: if a PDB structure is analyzed, convert it first to get a proper ICM-object (true ICM-molecular object does not require prior preparation for building Ramachandran plot).
Example:
```
convert a_1.               # Note, one more object appeared in addition
# to the original (PDB) object 1crn

quit
```

# How to display hydrogen bonds

A list of hydrogen bonds may be calculated and displayed for an ICM-objects or non-ICM object with hydrogens or made into a separate object of atom pair distances with labels. If you are dealing with a PDB structure without hydrogens, convert it first. The command show hbond prints list of hydrogen bonds in the text window. After that they can be displayed. (Hydrogen bonds can also be calculated by minimize and show energy commands provided that the hydrogen bond energy term is switched on.)
The display hbond command allows one to show the deviation angle of the hydroben bond from linearity (see the GRAPHICS.hbondStyle preference ).

Examples:
```
show energy
display
show hbond 2.5 a_/1:15  # list of H-bonds with H-X distance < 2.5 A
# appears in the text window
display hbond 1.9       # H-bonds shorter than 1.9 A are shown
GRAPHICS.hbondStyle = 3
display hbond
```

# How to identify atoms or residues at the molecular interface

When two or more parts of the polypeptide chain(s) are near each other in space, they are referred to as a molecular interface. What "interface" is can be defined more specifically in the context of a particular study, so here only an example is given to illustrate how interface may be identified and displayed. Two ICM functions are to be used for that, viz. Sphere and Acc . Suppose, you analyze 3D structure of a complex of two molecules, and would like to see what residues are at the interface. It can be done by the following:
```
display a_1,2//!h*         # display both molecules
# without hydrogens

color a_1 yellow
color a_2 green

show area surface a_1//!h* a_1//!h*      # calculate surface of
# the 1st molecule only

color red Sphere(a_2//* a_1//* 4.) & Acc(a_1/*)# interface residue atoms
# of the 1st molecule
# in 4 A radius vicinity
# of the 2nd molecule

show area surface a_2//!h* a_2//!h*      # calculate surface of
# the 2nd molecule only

color blue Sphere(a_1//* a_2/* 4.)&Acc(a_2//*) # interface residue atoms
# of the 2nd molecule
# in 4 A radius vicinity
# of the 1st molecule
```
If, instead, you need to mark residues, convert the selection of the interface atoms to residues with the Res () command:
```
....                                         # same as above.
....
resAtInterf =  Res( Sphere(a_2//* a_1//* 4.)) & Acc(a_1/*)
display residue label resAtInterf
```
Note that in the Sphere command it does not matter if you specify the atom selection or a residue selection as an argument, since the function operates at the atom level anyway. The difference in the specification of the ICM selection in these two examples (usage of two slashes for atom selection, and one slash for residue selection):
Sphere(a_1.1/* 4.) versus Sphere(a_1.2//* 4.) and also Acc(a_1.1/*) versus Acc(a_1.2//*) for specifying residues and atoms, respectively.
Important: when calculating surface, be sure that you properly specify the selection arguments in the show area surface command.

# How to select accessible, buried, hydrophobic, residues.

Selecting exposed residues. Here is one way to compile a list of exposed residues. Use the Acc ( rs_ [ r_threshold ] ) function. It will return all residues for which relative residue solvent accessibility is larger than certain limit ( by default it is 25% of its fully accessible surface ). To use the function you need to get rid of water molecules and use the show surface are command. Follow this example:
```
delete a_W                # delete water molecules
show surface area         # compute exposed surface areas for each residue
show Acc( a_/* , 0.25 )   # show all residues exposed more than 25%)
show !Acc(a_/* , 0.25 )   # show buried residues
```

Converting the selection into other formats
You can also show the selection in a different format, e.g.
```
String( Acc( a_/* ))  # or
a_1crn.m/1,5:8,11:12,14:15,17:20,22:25,28:29,31,33:34,36:46
Label(  Acc(a_/*))
#>S string_array
T1
P5
S6
I7
V8
S11
...
```

Identifying buried polar residues
About 50% of all residues have relative accessibility less than 25%. Polar residues typically do not like to be buried. The charged residues like it even less, and if they are buried they usually form a salt bridge. Example:
```
read pdb "1qau" # neuronal nitric oxide synthase
display
display xstick a_/62,121 # buried asp and arg form a bridge
```

To identify buried charged residues, use combine the previous selection of the buried residues with specific residue type selection, e.g.
```
delete a_W
show surface area
show a_/asp,glu,arg,lys & ! Acc( a_/* 0.15 )
```
Here we used a more conservative threshold of 0.15 as the burial threshold. Feel free to modify the selection of residue types above to find other buried residues.
Identifying exposed hydrophobic patches
A similar technique can be used to identify hydrophobic patches:
```
delete a_W
show surface area
show a_/leu,ile,val,met,trp & ! Acc( a_/* )  # buried ones
show Acc( a_/leu,ile,val,met,trp )           # exposed hydrophobs
```

To find clusters of exposed hydrophobic residues, use the Sphere function. The Sphere function returns atoms, and they need to be converted to residues with the Res function.
```
exposed_hres =  Acc( a_/leu,ile,val,met,trp )
for i=1,Nof(exposed_hres )
nbrs = Res(Sphere( exposed_hres[i] exposed_hres , 5.0 ))
if Nof( nbrs ) >= 2 show nbrs   # show pairs of exposed hydrophobs
endfor
```

# How to identify torsions at the molecular interface

Identification of the torsions belonging to residues at the molecular interface is a necessary and non-trivial step in many tasks of the molecular modeling. An example below shows how this identification may be done in ICM. The same as above ICM-object "complex" is considered.
```
display a_1.1,2//!h*                  # display both molecules
# of the complex w/o hydrogens
color a_1.1 yellow
color a_1.2 green

show area surface a_1.1//!h* a_1.1//!h*   # calculate surface of
# the 1st molecule only
a1=Sphere(a_1.2//* 7.) & Acc(a_1.1//*) # define the 1st molecule atoms
# belonging to the interface
# of the 1st molecule
# in 7 A radius vicinity

v1=V_1.1//S & a1                      # identify standard geometry
# torsions of the 1st molecule
# belonging to the interface

color red Atom(v1)                    # color atoms which torsions
# belong to

# similar for the 2nd molecule
show area surface a_1.2//!h* a_1.2//!h*

a2=Sphere(a_1.1//* 7.) & Acc(a_1.2//*)
v2=V_1.2//S & a2
color blue Atom(v2)
```

# How to calculate packing density

The packing density analysis requires understanding of two types of surfaces: the skin (molecular surface) and solvent-accessible surface of water probe centers (which is one water radius away from the skin ). The following is an example of how it may be done for a fragment of a protein.
```
asel = a_/5:15
show volume skin asel
rskin = r_out
vwExpand = 0.
show volume surface asel
rsurf = r_out
print "skin volume = ", rskin, "; vw volume = ", rsurf
print "packing density = ", rskin/rsurf
```

# How to perform a principal component analysis

For a set of objects with given measure of similarity between each two of them, one can easily perform the principal component analysis or to solve a distance geometry problem by using Disgeo function. The following example shows how to get a two-dimensional distribution of amino acid sequences of a series of Zn-fingers given the distance between sequences is defined by Distance( sequence1 sequence2 ) function. This distance is essentially a measure of sequence similarity: the distance is 0. for two identical sequences, it is proportional to percent identity divided by 100. for very similar sequences and goes above one at about 30% sequence identity, tending to infinity for very small seq. identity numbers.
```
list sequences              # see them, then ...
group sequence alZnFing     # group them, then ...
align alZnFing              # align them, then ...
a = Distance(alZnFing)      # calculate a matrix of pairwise
# distances among them
n=Nof(a)                    # number of points
b=Disgeo(a)                 # calculate principal components
corMat = b[1:n,1:n-1]       # coordinate matrix [n,n-1] of n points
eigenV = b[1:n,n]           # vector with n sorted eigenvalues
xplot  = corMat[1:n,1]
yplot  = corMat[1:n,2]
# plot 2D distribution
plot xplot yplot CIRCLE display
```

# How to calculate a dihedral angle

Normal dihedral angles like torsion angles describing conformation are directly available through the Value( ) function. You need to convert an object into the ICM type if necessary.
Example:
```
convert                # quickly convert into ICM-object
show v_//phi,psi       # just show all phi,psi's
show v_/3:17/xi*       # show chi angles of residues from 3 to 17
a = Value(v_//phi,psi) # create a real array of values of spec. torsions
```
You can calculate a dihedral angle between any two planes defined by three atoms (for example two Phe rings) with the calcDihedralAngle macro. If this macro is loaded, you can do the following:
```
as_1 = a_/1/n,ca,c
as_2 = a_/3/n,ca,c
display a_//n,ca,c blue
color as_1 magenta
color as_2 green
calcDihedralAngle as_1 as_2
print "dihedral angle 1(n,ca,c) and 3(n,ca,c) (deg) = ", r_out
```
Note that the order of atoms in the selections as_1 and as_2 is determined only by the ICM-molecular tree (will be the same for a_//n,c or a_//c,n). Thus, any changes in the selections as_1 and as_2 not changing their content has no effect on the resulting dihedral angle:
```
as_1 = a_/1/ca,c,n
as_2 = a_/3/ca,c,n
calcDihedralAngle as_1 as_2
print "dihedral angle 1(n,ca,c) and 3(n,ca,c) (deg) = ", r_out
```

# How to print a table of the torsion angles

The simplest list can be generated by:
```
show V_//*               # or
show V_//phi*,psi*,omg*  # or
show V_//xi*             # side chain torsions
```
A nicer formatted output (one line per residue) may be generated with macro printTorsions, for example:
```
convert
printTorsions a_/2:15
```
Note, that you do not need to convert your molecular object if it is an ICM-object.

# How to build a hydrophobicity profile

First, define a hydrophobicity scale, for example that from Kyte and Doolittle, 1982 or use your favorite one. (Note, there should be 26 entries in the hydrophobicity parameters list hPhobInd corresponding to the 26 letters of the alphabet. Non-participating letters B,J,O,U,X,Z are marked by zero values.)
```
# define a hydrophobicity scale
hPhobInd = { 1.8,  0.0,  2.5, -3.5, -3.5,  2.8, -0.4, \
-3.2,  4.5,  0.0, -3.9,  3.8,  1.9, -3.5, \
0.0, -1.6, -3.5, -4.5, -0.8, -0.7,  0.0, \
4.2, -0.9,  0.0, -1.3,  0.0}

# make a macro
macro hPhobProfile seq_ i_windowSize
if (Type(i_windowSize)=="unknown") then
endif
R_window = Rarray(i_windowSize,1./i_windowSize)
R_hphob  = Smooth (Rarray(seq_,hPhobInd),  R_window)
R_ruler  = {0.,0.,10.,10.,0.,0.,0.,0.}
R_ruler = Real(Length(seq_))
r_tic = 1./Sqrt(Real(i_windowSize))
r_tic = Integer(r_tic*100.0)/100.
R_ruler = -7.*r_tic
R_ruler =  7.*r_tic
R_ruler = r_tic
R_ruler = r_tic
print R_ruler
s_legend = {"Hydrophobicity plot","Sequence","Hydrophobicity"}
xplot = Count(1 Length(seq_))
yplot = R_hphob
psfilename = "hphob"
plot xplot yplot R_ruler s_legend grid="xy" display psfilename
delete R_window R_ruler r_tic
endmacro
# now, an example
s = Sequence(a_)
hPhobProfile s 7
```

# How to display and characterize protein cavities

ICM offers fast, elegant and mathematically accurate way to identify, display, and measure protein cavities. An example session which displays all cavities with their surroundings, and calculates their volumes and surface areas:
Examples:
```
read object "crn"  # or whatever
make grob skin "g_skin"
split g_skin
nShells = i_out
display wire residue labels
for i=1,nShells
v = Volume(g_skin\$i)   # actually its surface is returned in r_out
s = r_out              # there is no need to use the explicit Area(g_skin\$i)
if(v > 0.) then        # note that cavities have negative volume!
display transparent smooth g_skin\$i
printf "Shell  %d: V=%f A=%f\n", i, v, s
else
display reverse smooth yellow g_skin\$i
center g_skin\$i
printf "CAVITY %d: V=%f A=%f R~%f\n", i, -v, s, -3.*v/s
endif
# add pause here for an interactive session
endfor
```

 PrevColor by charge HomeUp NextSequence