[ Superimpose faq | Align 3D faq | Ramachandran faq | Hbonds faq | Interface residues | Res selections | Interface torsions | Packing density | Principal component analysis | Dihedral angle calculation | All torsions table | Hydrophobicity profile | Cavity analysis ]
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:
|
|
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:
read pdb "3znf"
display a_1.1//n,ca,c magenta
make sequence a_1.1
read pdb "1ard"
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
See also: transform, transformation vector, Rmsd( ) and Srmsd( ).
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:
read pdb "4fxc"
read pdb "1ubq"
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
plotRamaEps rs_selection l_addLabel l_addBoundaries
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:
read pdb "1crn"
convert a_1. # Note, one more object appeared in addition
# to the original (PDB) object 1crn
l_addLabel = yes # add residue labels to the plot
l_shadedBoundaries = yes # add allowed regions to the plot
plotRamaEps a_2. l_addLabel l_shadedBoundaries
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 to show the deviation angle of
the hydroben bond from linearity (see the GRAPHICS.hbondStyle preference ).
Examples:
read object "crn" # already converted
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
See also:
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:
read object "complex"
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:
read pdb "1crn"
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.
read pdb "1qau"
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:
read pdb "1qau"
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.
read object "complex"
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.
read object s_icmhome+"crn"
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.
read sequences s_icmhome+"zincFing" # read sequences from file,
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:
read pdb "1crn" # read it in
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:
read object "crn"
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
See also:
Acos( ),
Length( ),
Sum( ),
Vector( ).
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:
read pdb "1crn"
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
i_windowSize = Ask("Enter window size",windowSize)
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[2] = Real(Length(seq_))
r_tic = 1./Sqrt(Real(i_windowSize))
r_tic = Integer(r_tic*100.0)/100.
R_ruler[5] = -7.*r_tic
R_ruler[6] = 7.*r_tic
R_ruler[7] = r_tic
R_ruler[8] = 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
read object "crn"
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