Documentation

Realign — superpose molecules

realign computes the rotation and translation that puts one molecule onto another, then applies that same transformation to any number of other files. The usual reason to want this: GRIM has matched two complexes and told you how their interaction patterns correspond, and you now want the actual protein and ligand moved into that frame so you can look at them together.

Run it

The reference (GRIM_ints.mol2) and the mobile file (4amj_INTS_M.mol2) are interaction pseudoatom files. The transformation that superposes them is then applied to the 4AMJ ligand and protein.

cd "$ICHEM_DIR/dataset/test/REALIGN"
IChem --wMob realign GRIM_ints.mol2 4amj_INTS_M.mol2 4amj_lig.mol2 4amj_prot.mol2

What you get

realign is quiet on success. Each input gets a rot_ copy in the reference frame:

ls rot_*
rot_4amj_INTS_M.mol2  rot_4amj_lig.mol2  rot_4amj_prot.mol2

These three files are byte-identical to the stored reference outputs in REALIGN/output/, which is a useful way to confirm your installation is sound.

Options worth knowing

IFP — interaction fingerprints

IFP turns a protein–ligand complex into a bit string laid out residue by residue: eleven bits per binding-site residue, one per interaction type. Because the layout follows the binding site, two fingerprints are only comparable when they come from the same site — which is exactly what you want when scoring docking poses.

Fingerprint one complex

cd "$ICHEM_DIR/dataset/test/IFP"
IChem --all IFP site.mol2 ligand.mol2

First the interaction table, one row per contact detected:

HBond_PROT        	|ND2  |421  |ASN 312-A	|O17  |1    |CAU   0-XX	|0	|2.76961	|159.236
HBond_LIG         	|OD1  |78   |ASP 113-A	|O17  |1    |CAU   0-XX	|1	|2.60795	|134.667
Hydrophobic       	|CZ   |349  |PHE 289-A	|C16  |2    |CAU   0-XX	|3	|3.74597	|/
Ionic_LIG         	|OD1  |78   |ASP 113-A	|N19  |4    |CAU   0-XX	|4	|3.60381	|/
Ionic_LIG         	|OD2  |79   |ASP 113-A	|N19  |4    |CAU   0-XX	|5	|2.94042	|/
Hydrophobic       	|CZ3  |46   |TRP 109-A	|C21  |6    |CAU   0-XX	|6	|4.43435	|/

In the first section, all intermolecular interactions between active site and ligand are tabulated.

Then the fingerprint itself:

|A      M82|A      V86|A     W109|A     T110|A     D113| ... |A     N312|A     Y316
000000000000000000000010000000000100000000000000101000010 ... 00010000000000

In the second section, the interaction fingerprint is displayed with 11 bits/residue in normal (basic) mode. The first line describes the chain identifier and residue number. The second line is the bit string (1: interaction, 0: no interaction). The order is the following:

Score docking poses against a reference

Give IFP a third file and it fingerprints the poses, then reports the Tanimoto similarity of each to the reference ligand's fingerprint.

IChem --all IFP site.mol2 docked.mol2 ligand.mol2 | tail -1
DOCKED	REF	0.7

This docked pose reproduces 70 % of the crystallographic interaction pattern — the kind of number you rank a docking run on.

Choosing a bit format

Options are all adjustable and listed in the IChem full menu page.

INTS — interactions as pseudoatoms

ints detects the same contacts as IFP but stores each one as a pseudoatom in a MOL2 file rather than as a bit. Because a pseudoatom has coordinates instead of a residue slot, the result no longer depends on the binding site — you can compare interaction patterns across unrelated proteins. This is the representation GRIM consumes.

Generate the pseudoatoms

cd "$ICHEM_DIR/dataset/test/TIFP"
IChem -logf 2rh1_ints.txt -type CENT ints site.mol2 ligand.mol2
CENT

Two files come out:

For this complex the 28 pseudoatoms break down as 22 hydrophobic contacts, 3 hydrogen bonds, 2 ionic bonds and 1 edge-to-face aromatic interaction. Each carries a residue name encoding its type — SEC, ALC, ASC, GLC, PHC:

      1 N         -33.9740   10.7140    9.4395 N.am      1 SEC1        0.0000
      2 O         -33.9650   11.1100    6.9705 O.2       2 ALC2        0.0000
      4 OD1       -34.0775    9.7565    6.6855 O.co2     3 ASC4        0.0000
      6 CA        -31.2970    9.5295    9.4820 C.3       4 GLC5        0.0000

Atomic description of the interaction pseudoatoms (MOL2 file):

Where the pseudoatom sits

-type decides the geometry, and the suffix on the output filename follows it:

MERG is what GRIM wants, since matching on protein-side, ligand-side and central points together is far more discriminating.

Triplet fingerprints (TIFPs)

Adding -fgps encodes every triplet of pseudoatoms as a fixed-length integer vector — alignment-free, and comparable across any two complexes.

IChem -fgps STD ints site.mol2 ligand.mol2 2rh1_full.fgp
# SAVING FINGERPRINT

IChem --small -fgps STD ints site.mol2 ligand.mol2 2rh1_small.fgp
# SAVING FINGERPRINT

Each line is the source name, a tab, then the vector: 12 509 integers at full length, or 210 with --small. Formats are STD (plain), SVM (sparse index:value) and CMP (compressed) — for long sparse vectors SVM is much the smaller. Point repeated runs at the same output file and the fingerprints accumulate, one row per complex, ready for sims.

GRIM — graph matching of interaction patterns

GRIM is a tool to match protein‐ligand complexes using a graph matching algorithm focusing on protein‐ligand interaction pseudoatoms. It thereby enables an interaction‐based alignment of different protein‐ligand complexes that can be quantified by an empirical scoring function (GrScore) and used to post‐process docking poses by similarity to known protein‐ligand interaction patterns. In a recent international contest (D3R Docking Challenge 2015), GRIM was ranked 2nd out of 44 scoring functions to predict the binding mode of 36 inhibitors prior to the release of protein‐bound X‐ray coordinates.

Compare two complexes

2RH1 (β2-adrenergic receptor) against 4AMJ (β1-adrenergic receptor). GRIM computes the interactions itself when given four MOL2 files.

cd "$ICHEM_DIR/dataset/test/GRIM"
IChem -sim 1 -rn 2rh1 -cn 4amj -score STD grim \
      2rh1_prot.mol2 2rh1_lig.mol2 4amj_prot.mol2 4amj_lig.mol2

7 files are outputted:

cat Grifp_res.csv
NCli	Ref	Comp	Simil	LIG	CENTER	PROT	SumCl	RMSD	RLig	RCent	RProt	CLig	CCent	CProt	GrSc	RMSDAl	NPol
1	2rh1	4amj	1	16	19	11	0.0587	0.1653	19	28	23	24	33	29	0.9740	0.3365	11

Alongside the scores, GRIM writes GRIM_prot.mol2, GRIM_lig.mol2 and GRIM_ints.mol2 — the comparison complex moved into the reference frame — plus GRIM.log, which lists every matched pair of interaction points. Those aligned files are the input to realign.

Reuse interactions you already computed

If you already ran ints -type MERG, hand GRIM the pseudoatom files directly and skip detection. --values reports scores without writing the aligned structures.

IChem -sim 1 -rn 2rh1 -cn 4amj --values grim 2rh1_INTS_M.mol2 4amj_INTS_M.mol2
NCli	Ref	Comp	Simil	LIG	CENTER	PROT	SumCl	RMSD	RLig	RCent	RProt	CLig	CCent	CProt	GrSc	RMSDAl	NPol
1	2rh1	4amj	1	16	19	11	0.0587	0.1653	19	28	23	24	33	29	0.9740	0.3365	11

Screen docking poses

Given a site (site.mol2), a reference ligand (ligand.mol2) and docking poses (docked.mol2), GRIM scores every pose against the reference binding mode.

The file docked.txt summarized the docking results: pose number, docking score (pKd), RMSD to the true X‐ray pose.

All docking poses (docked.mol2) are matched to the X‐ray pose (ligand.mol2) for rescoring based on the interaction pattern graph similarity score (GrScore).

cd screen
IChem -rn 2rh1 --values grim site.mol2 ligand.mol2 docked.mol2

cat GrScreen.csv
REF	1.05313	48	14	14	20	0.202949	0.0721811	0.00215763	11

GrScreen.csv keeps the best possible alignment in the following table format:

Grimscreen_res.tsv keeps them all. Ranking on the Grimscore column (GrSc) prioritises poses that best reproduce the reference interaction pattern, which is a considerably better filter than geometry alone.

The format of the Grimscreen_res.tsv file is explained below:

The pose with the highest Grimscore is not the top‐ranked one (ligand_000) but pose 002. Interestingly, it is the pose with the smallest rmsd to the true X‐ray pose (see docking.txt)

Options are all adjustable and listed in the IChem full menu page.

VolSite — cavity detection and druggability

volsite lays a grid over the protein, keeps the points that are buried enough to form a pocket, labels each with the pharmacophoric property a ligand atom would need there, and predicts how druggable the resulting cavity is. volsite can be run in two modes depending on whether coordinates of a bound ligand are given (ligand‐restricted mode) or not (unrestricted mode).

A known binding site

With a ligand supplied, detection is restricted to its site — fast, and the answer is about the pocket you care about.

cd "$ICHEM_DIR/dataset/test/VOLSITE/ligand"
IChem volsite protein.mol2 ligand.mol2
Ligand file : ligand.mol2
Protein name : Unknown
Druggability : 1.13
Prediction : Is Druggable

cat VolSite_Stat.csv
Name	NCav	Size	NPts	Volume	CA	CZ	O	OG	OD1	N	NZ	DU	Recovery	Drugg
CAU	1	4	115	388.125	31	37	3	12	0	22	6	4	53.913		/
CAU	1	6	160	540	44	44	6	20	0	30	10	6	38.75		/
CAU	1	8	205	691.875	58	56	6	22	1	37	15	10	30.2439		/
CAU	1	12	303	1022.62	83	74	9	26	8	56	30	17	20.462		/
CAU	1	ALL	321	1083.38	87	76	11	27	8	61	34	17	19.3146		1.12574

The VolSite_stat.csv summary recapitulates the properties of all detected cavities as follows:

Duggability estimation is performed on the full cavity (no truncation).

Each cut-off is also written as a MOL2 file (CAVITY_N1_4.mol2 … CAVITY_N1_ALL.mol2) that you can load into any viewer — the pocket appears as a cloud of property-coloured pseudoatoms.

The whole protein surface

Drop the ligand and VolSite scans everywhere, ranking what it finds. --desc adds a descriptor file for downstream modelling.

cd ../unrestricted
IChem --desc volsite protein.mol2
Run VolSite on : protein.mol2
Build the file descriptor.txt
6 Cavities found
End

cat VolSite_Stat.csv
Name	NCav	Size	Buriedness	NPts	Volume	CA	CZ	O	OG	OD1	N	NZ	DU	Drugg
/	1	ALL	84.9875	321	1083.38	87	76	11	27	8	61	34	17	1.12574
/	2	ALL	78.1878	149	502.875	63	15	4	15	8	12	13	19	0.71907
/	3	ALL	70.0841	133	448.875	30	9	6	5	0	32	28	23	-0.0102194
/	4	ALL	71.5261	80	270	36	11	3	5	2	7	0	16	0.543419
/	5	ALL	71.0452	78	263.25	53	9	0	0	0	2	0	14	1.05675
/	6	ALL	71.3473	75	253.125	18	5	5	0	18	1	18	10	-1.98361
END

In this mode, a new column (Buriedness) is added to VolSite_Stat.csv, describing the average buriedness of each cavity found.

The IChem ‐‐desc command outputs cavity properties as a vector of 89 reals that may be used for machine learning applications

The projection value is the number of regularly‐spaced 8 Å‐long vectors emitted from each cavity points intercepting the protein surface (maximum of 120).

Correspondence between projection value and buriedness

Cavity 1 is the orthosteric site found in the previous run — same 321 points, same 1.13 score — recovered here with no prior knowledge of where the ligand binds. Cavity 5 is small but scores 1.06, the sort of secondary pocket worth a look. Cavities 3 and 6 have negative scores and are not worth pursuing.

Cavity-based pharmacophores

--pharm converts the cavity into a pharmacophore you can screen with directly — no ligand required.

cd ../pharm
IChem --pharm volsite protein.mol2 ligand.mol2

ls Pharmacophore.*
Pharmacophore.chm  Pharmacophore.mol2  Pharmacophore.pml

30 features in total — 19 directional (H-bond donors and acceptors, written as vectors) and 11 positional (hydrophobic and aromatic points). The three files are the same model in three file formats: .chm for Discovery Studio and Pipeline Pilot, .pml for LigandScout, .mol2 for Shaper2. Add --outExclu for exclusion spheres.

Tuning detection

Options are all adjustable and listed in the IChem full menu page.

Utils — buried surface, fragments, waters

Three self-contained utilities that answer questions which come up constantly but do not warrant a module of their own.

bsa — buried surface area

How much of the ligand the protein actually wraps around. A good first-pass check on whether a pose is plausibly buried.

cd "$ICHEM_DIR/dataset/test/BSA"
IChem utils bsa protein.mol2 ligand.mol2
2rh1_protein	70.1934	347.125

70.19 % of carazolol's surface is buried, amounting to 347.1 Ų. For a well-enclosed orthosteric site, 70 % is about what you would expect; a solvent-exposed surface binder would come in far lower.

frag — RECAP fragmentation in context

Cuts the bound ligand at RECAP bonds and reports each fragment together with the interactions it makes — so you can see which piece is earning its keep.

cd ../Frag
IChem utils frag protein.mol2 ligand.mol2 | tail -2
2rh1_protein|FRAG	ID	NHeavy	NLink	NRing	NArRing	TOT	CA	CZ	O	OD1	N	NZ	ZN
2rh1_protein|FRAG	1	18	2	3	2	23	17	1	2	2	1	0	0

18 heavy atoms, 2 cut points, 3 rings of which 2 aromatic, making 23 interactions — 17 of them hydrophobic. Each fragment is written out with its interaction pseudoatoms (*_FRAG_1_MOLE.mol2 and *_FRAG_1_INTS.mol2), which is the raw material for bioisosteric replacement.