Circular Fingerprints with the CDK and Clojure

One of the things I’m liking about Clojure is that it can be used as a quick prototyping language, a lot like Python. This is quite handy when playing with the CDK, as I can avoid a lot of the verbosity of Java code. As an example, I put together a very simplistic circular fingerprint (exemplified by molprint2d and ECFP’s). In this case, the fingerprint is simply composed of atom type names and is quite similar to an implementation by Noel.

First, some imports and simple helper functions that allow me to work with native Clojure sequences rather than having to go via Java methods

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
(import '(org.openscience.cdk.smiles SmilesParser))
(import '(org.openscience.cdk DefaultChemObjectBuilder))

(use 'clojure.contrib.duck-streams)
(use 'clojure.contrib.str-utils)

;; surely there must be a core function for this!
(defn has
  "true if val exists in coll, nil otherwise. Uses ="
  [coll val]
  (some #(= % val) coll))

;; http://markmail.org/message/56r3eflx4a6tasoe
(defn flatten
  "Flatten a list which may contain lists"
  [x]
  (let [s? #(instance? clojure.lang.Sequential %)]
    (filter
     (complement s?)
     (tree-seq s? seq x))))

;; Maybe these could be converted to macro's?
(defn get-atom-list
  "Get the atoms of a molecule as a Clojure sequence"
  [mol]
  (map #(identity %) (. mol atoms)))

(defn get-connected-atoms-list
  "Get the atoms connected to a specified atom of a molecule
   as a Clojure sequence"

  [mol atom]
  (map #(identity %) (. mol (getConnectedAtomsList atom))))

(defn getmol
  "Parse a SMILES string and return the molecule object. This is
   quite inefficient since each time it is called it creates a new
   parser object. Handy for quick testing"

  [smiles]
  (. (new SmilesParser (DefaultChemObjectBuilder/getInstance))  
     (parseSmiles smiles)))

These go into my user.clj, so are available when I start Clojure. With these in hand, we need a function which will get us the atoms in successive shells around a specified atom. In other words, atoms in the first shell are those connected directly to the specified atom. Atoms in the second shell lie two bonds away and so on. This is achieved by the following function

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
(defn get-shells
  "Get the atoms in successive shells surrounding a specified atom.
   The first shell equals the immediate neighbors, the second shell
   the neighbors at bond length 2 and so on. The return value is a list
   in which each element is a list of atoms. The first element of the
   top level list is the n'th shell. Note that the return value does not
   contain the starting atom (i.e., zeroth shell)"

  [mol start depth]
  (loop [shells (list (list start))
        visited (list)]
    (if (= (count shells) (+ 1 depth))
      (drop-last shells)
      (recur (conj shells
           (filter #(not (has visited %))
               (flatten (map #(get-connected-atoms-list mol %) (first shells)))))
         (set (flatten (list visited shells)))))))

Here, depth indicates how many shells we want to go out to. The return value is a list of lists, each sublist representing the atoms in a shell. The first element corresponds to the n’th shell. Note that the zeroth shell (i.e., the start atom) is not contained in the return value.

Before applying this function to all the atoms, we need to consider how we’re going to represent each shell. For example ECFP’s consider various atom properties (degree, H-bonding ability etc). But a very simple procedure is to simply replace each atom in a shell with its atom type name, sort the names and join them into a single string. The end result would be, for a single atom, a list of strings corresponding to the shells around that atom.

First, we have a simple function to convert a list of atoms into a list of CDK atom type names

1
2
3
(defn atype [atoms]
  (map #(. % getAtomTypeName)
       atoms))

Using this, we can then use it to convert the lists of “atom shells” around a given atom to a list of atom type strings

1
2
3
4
5
(defn make-atom-fp
  [mol atom depth]
  (map #(str-join "-" (sort %))
       (map atype
           (get-shells mol atom depth))))

The return value is still a list of strings representing each shell around the specified atom. In practice this list of strings would likely be reduced to a single string or number. The last thing to do is to apply the function to all the atoms in a molecule (using a depth of 2 for example)

1
(map #(make-atom-fp mol % 2) (get-atom-list mol))

An example of using these functions is:

1
2
3
4
user=> (def mol (getmol "C(O)(F)CN"))
#'user/mol
user=> (map #(make-atom-fp mol % 2) (get-atom-list mol))
(("N.sp3" "C.sp3-F-O.sp3") ("C.sp3-F" "C.sp3") ("C.sp3-O.sp3" "C.sp3") ("F-O.sp3" "C.sp3-N.sp3") ("C.sp3" "C.sp3"))

Alternatively, we could reduce this to a list of strings, one for each atom, by simply joining the strings for each shell of a given atom

1
2
user=> (map #(str-join "-" %) (map #(make-atom-fp mol % 2) (get-atom-list mol)))
("N.sp3-C.sp3-F-O.sp3" "C.sp3-F-C.sp3" "C.sp3-O.sp3-C.sp3" "F-O.sp3-C.sp3-N.sp3" "C.sp3-C.sp3")

This is a very simple example, and as I’m still wrapping my head around Lisp, it took some time to write. But once you get the hang of loop and recur, things start flowing!

2 thoughts on “Circular Fingerprints with the CDK and Clojure

  1. […] while back I had started playing with Clojure. It’s always been a spare-time hobby and not having had much spare time […]

  2. […] that are represented as variable length vectors of numbers or strings. An example would be circular fingerprints. Now, when reading fingerprints you have to indicate whether you’re loading binary […]

Leave a Reply

Your email address will not be published. Required fields are marked *