On these lines I have been playing with text mining and modeling tools in R, mainly via the excellent tm package. One of the techniques I have been playing around with is Latent Dirichlet Allocation. This is a generative modeling approach, allowing one to associate a document (composed of a set of words) with a “topic”. Here, a topic is a group of words that have a higher probability of being generated from that topic than another topic. The technique assumes that a document is comprised of a mixture of topics – as a result, one can assign a document to different topics with different probabilities. There have been a number of applications of LDA in bioinformatics with some applications focusing on topic models as way to cluster objects such as genes [1, 2], whereas others have used it in the more traditional document grouping context [3].

In text mining scenario, developing an LDA model for a set of documents is relatively straightforward (in R) – perform a series of pre-processing steps (mainly to standardize the text) such as converting everything to lower case, removing stopwords and so on. At the end of this one has a series of documents, each one being represented as a bag of words. The collection of words across all documents can be converted to a document-term matrix (documents in the rows, words in the columns) which is then used as input to the LDA routine.

Those familiar with building predictive models with keyed fingerprints will find this quite familiar – the individual bit positions represent structural fragments, thus are the chemical analogs of words. Based on this observation I wondered what I would get (and what it would mean) by applying a technique like LDA to a collection structures and their fragments.

My initial thought is that the use of LDA to determine a set of topics for a collection of chemical structures is essentially a clustering of the molecules, with the terms associated with the topics being representative substructures for that “cluster”. With these topics in hand, it wil be interesting to see what (or whether) properties (physical, chemical , biological) may be correlated with the clusters/topics identified. The rest of this post describes a quick first look at this, using ChEMBL as the source of structures and R for performing pre-processing and modeling.

Structures & fragments

We had previously fragmented ChEMBL (v8) in house, so obtaining the data was just a matter of running an SQL query to identify all fragments that occured in 50 or molecules and retrieving their structures and the molecules they were associated with. This gives us 190,252 molecules covered by 6,110 fragments. While a traditional text document-based modeling project would involved a series of pre-processing steps, the only one I need to perform in this scenario is the removal of small (and thus likely very common) fragments such as benzene – the cheminformatics equivalent of removing stopwords. (Ideally I would also remove fragments that already occur in other fragments – the cheminformatics equivalent of stemming)

The data file I have is of the form

where natom is the number of atoms in the fragment. The R code to generate (relatively) clean data, read to feed to the LDA function looks like:

In the code above, we rearrange the data to create “documents” – identified by a title (the molecule identifier) with the body of the document being the space concatenated SMILES for the fragments associated with that molecule. In other words, a molecule (document) is constructed from a set of fragments (words). With the data arranged in this form we can go ahead and reuse code from the tm and topicmodels packages.

Finally, we’re ready to develop some models, starting of with 6 topics.

So, what are the topics that have been identified? As I noted above, each topic is really a set of “words” that have a higher probability of being generated by that topic. In the case of this model we obtain the following top 4 fragments associated with each topic (most likely fragments are at the top of the table):

Visual inspection clearly suggests distinct differences in the topics – topic 1 appears to be characterized primarily by the lack of aromaticity, whereas topic 2 appears to be characterized by quinoline and indole type structures. This is just a rough inspection of the most likely “terms” for each topic. It’s also interesting to look at how the molecules (a.k.a., documents) are assigned to the topics. The barchart indicates the distribution of molecules amongst the 6 topics.

As with other unsupervised clustering methods, the choice of k (i.e., the number of topics) is tricky. A priori there is no reason to choose one over the other. Blei in his original paper used “perplexity” as a measure of the models generalizability (smaller values are better). In this case, we can vary k and evaluate the perplexity: with 6 topics the perplexity is 1122, with 12 topics it drops to 786 and with 100 topics it drops to 308 – you can see that it seems to continuously decrease with increase in number of topics (which has been observed elsewhere, though in my case, the hyperparameters are kept constant). Wallach et al have discussed various approaches to evaluating topic models.

Numerical evaluation of these models is useful, but we’re more interested in how these assignments correlate with chemical or biological features. First, one could look at the structural homogenity of the molecules assigned to topics. For k = 6, this is probably not useful, as the individual groups are very large. With k = 100 one obtains a much more sensible estimate of homogeneity (but this is to be expected). Another way to evaluate the topics from chemical point of view is to look at some property or activity. Given that ChEMBL provides assay and target information for the molecules, we have many ways to perform this evaluation. As a brief example, we can consider activity distrbutions derived from the molecules associated with each topic. Most ChEMBL molecules have multiple activities associated with them as many are tested in multiple assays. To allow comparison we converted activities in a given assay to Z-scores, allow comparison of activitives across assays. Then for each molecule, we identified the minimum activity, only considering those activities that were annotated as IC₅₀ and as exact (i.e., not < or >). After removal of a few extreme outliers we obtain:

Clearly, within each group, the Z-scores cluster tightly around 0. It appears that the groups differentiate from each other in terms of the extreme values. Indeed plotting summary statistics for each group confirms this – in fact the median Z-score has a range of 0.05 and the mean Z-score a range of 0.11 across the six groups. In other words, the bulk of the groups are quite similar.

Other possibilities

The example shown here is rather simplistic and is the equivalent of unsupervised clustering. One obvious next step is to search the parameter space of the LDA model, evaluate different approaches to estimating the posterior distribution (EM or Gibbs sampling) and so on. A number of extensions to the basic LDA technique have been proposed, one of them being a supervised form of LDA.

It’d also be useful to look at this method on a slightly smaller, labeled dataset – I’ve run some preliminary experiments on the Bursi AMES but those results need a little more work. More generally, smaller datasets can be problematic as the number of unique fragments can be low. In addition fewer observations means that the estimates of the posterior distribution becomes fuzzier. One way around this is to develop a model on something like the ChEMBL dataset I used here and then apply that to smaller datasets. Obviously, this goes towards ideas of applicability – but given the size of ChEMBL, it may indeed “cover” many smaller datasets.

Is this useful?

At first sight, it’s an interesting method that identifies groupings in an unsupervised manner. Of course, one could easily run k-means or any of the hierarchical clustering methods to achieve the same result. However, the generative aspect of LDA models is what is of interest to me, but also seems the part that is difficult to map to a chemical setting – unlike topics in a document, which one can (usually) understand based on the likely terms for that topic, it’s not clear what a topic is for a collection of molecules in an unsupervised setting. And then, how does one infer the meaning of a topic from fragments? While it’s certainly true that certain fragments are associated with specific properties/activities, this is certainly not a given (unlike words, where each one does have an individual meaning). Furthermore, in an unsupervised setting like the one I’ve described here, fishing for a correlation between (some set of) properties and groupings of molecules is probably not the way to go.

Google searches are screamingly fast, so fast that the type-ahead feature is doing the search as you key characters in.  Why are all chemical searches so sloooow? … Ideally, as you sketch your mol in, the searches should be happening at the same pace, like the typeahead feature.

Now, he doesn’t specifically mention what type of chemical search – it could be exact matches, similarity searches, substructure or pharmacophore searches. The first two can be done very quickly and lend themselves easily to type ahead type search interfaces. In light of the work my colleague has been doing, the substructure searches are now also amenable to a type ahead interface.

So I quickly put together a simple web page that lets you type in a SMILES (or SMARTS) and as you type it retrieves the results of a substructure search via the NCTT Search Server REST API. (In some cases the depiction is broken – that’s a bug on my side). Of course, typing in SMILES is not the most intuitive of interfaces. Since Trung employs the ChemDoodle sketcher, an ideal interface would respond to drawing events (say drawing a bond or adding atoms etc) and pull up matches on the fly. Another obvious extension is to rank (or filter) the results – all the while, maintaining the near real time speed of the application.

As I said before, seriously fast substructure searches. It also helps that I can build these examples via a public REST API. I’m sure there are reasons for SOAP, XML and so on. But it’s 2011. So lets help make extensions and mashups easier.

UPDATE: Yes, it’s easy to create patterns (especially with SMARTS) that DoS the server. We have some filters for excessively generic patterns; so some queries may not behave in the expected manner

It is blazingly fast.

I decided to run some benchmarks via the REST interface that he provided, using a set of 1000 SMILES derived from an in-house fragmentation of the MLSMR. The 1000 structure subset is available here. For each query structure I record the number of hits, time required for the query and the number of atoms in the query structure. The number of atoms in the query structures ranged from 8 to 132, with a median of 16 atoms.

The figure below shows the distribution of hits matching the query and the time required to perform the query (on the server) for the 1000 substructures. Clearly, the bulk of the queries take less than 1 sec, even though the result set can contain more than 10,000 hits.

The figures below provide another look. On the left, I plot the number of hits versus the size of the query. As expected, the number of matches drops of with the size of the query. We also observe the expected trend between query times and the size of the result sets. Interestingly, while not a fully linear relationship, the slope of the curve is quite low. Of course, these times do not include retrieval times (the structures themselves are stored in an Oracle database and must be retrieved from there) and network transfer times.

Finally, I was also interested in getting an idea of the number of hits returned for a given size of query structure. The figure below summarizes this data, highlighting the variation in result set size for a given number of query atoms. Some of these are not valid (e.g., query structures with 35, 36, … atoms) as there were just a single query structure with that number of atoms.

Overall, very impressive. And it’s something you can play with yourself.

Unfortunately, what would have been interesting but was not provided was a comparison of the performance at the Oracle query level with other products such as JChem Cartridge and OrChem. Furthermore, the test case is just under a million molecules from Golovin & Henrick – the entire dataset (not just the keys) could probably reside in-memory on todays servers. How does the system perform when say faced with PubChem (34 million molecules)? The paper mentions a command line implementation of their search procedure, but as far as I can tell, the Oracle cartridge is not available.

The ABCD system has many useful and interesting features. But as with the other publications on this system, this paper is one more in the line of “Papers About Systems You Can’t Use or Buy“. Unfortunate.

Being new to the area, it’d be great to meet up over a beer, with people in the surrounding areas (NY/CT/RI) doing cheminformatics, predictive modeling and other life science related topics (any R user groups in the area?). If anybody’s interested, drop me a line (comment, mail or @rguha).

The full announcement is below:

For the 2011 Fall meeting in Denver (Aug 28 – Sep 1), CINF will be running an experimental session of lightning talks – short, strictly timed talks. The session does not have a specific topic, however, all talks should be related to cheminformatics and chemical information. One of the key features of this session is that we will not be using the traditional ACS abstract submission system, since that system precludes the inclusion of recent work in the program.

So, since we will be accepting abstracts directly, the expectation is that they be about recent work and developments, rather than rehashes of year-old work. In addition, talks should not be verbal versions of posters submitted for this meeting. Given the short time limits we don’t expect great detail – but we are expecting compact and informative presentations.

That’s the challenge.

What

• Talks should be no longer than 8 minutes in length. At 8 minutes, you will be asked to stop.
• Use as many slides as you want, as long as you can finish in 8 minutes
• Talks should not be rehashes of poster presentations
• Talks will run back to back, and questions & discussion will be held of off until the end

If you haven’t participated in these types of talks before here are some suggestions:

• No more than three slides for a 5 minute talk (but if you can pull of 20 slides in 8 minutes, more power to you)
• Avoid slides with too much text (and don’t paste PDF’s of papers!)
• A single chart per slide and make sure labels are readable at a distance

When

1:30pm, Wednesday, August 31st, 2011

Submissions run from July 14 to Aug 14

Where

Room 112, Colorado Convention Center

How

• Send in an abstract of about 100 – 120 words to cinf.flash@gmail.com
• We will let you know if you will be speaking by Aug 21 and we will need slide decks by Aug 24
• You must be registered for the meeting
• Note that the usual publication/copyright rules apply
• We will encourage live blogging and tweets (if we have net access)

• Directly evaluate molecular volume (based on group contributions) using get.volume
• Generate fingerprints using the hybridization state
• get.total.charge and get.total.formal.charge work sensibly
• Added a function (copy.image.to.clipboard) that copies the 2D depiction of a molecule to the system clipboard in PNG format
• Now, OS X users can view and copy molecule depictions. This is slower compared to the same operation on Windows or Linux since it involves shell’ing out via system. But it is better than not being able to view anything.

The two main additions to this version are

• Read support for the new FPS fingerprint format described by Andrew Dalke at the chemfp project. Note, it currently discards some of header information
• The fingerprint class now has a field, misc, (a list) that allows one to read in extra, arbitrary data that might be provided along with a fingerprint. Exactly what gets stored in this field depends on the line function used to read in the fingerprint data. Currently only the FPS parser returns extra data (when available) in this field.

As always, you can get the package source directly from the Github repository.

One of our platforms is ImageXpress from Molecular Devices, which stores images in a file-based image store and meta data and numerical image features in an Oracle database. While they do provide an API to interact with the database it’s a Windows only DLL. But since much of modeling requires I access the data from R, I needed a more flexible solution.

So, I’ve put together an R package that allows one to access numeric image data (i.e., descriptors) and images themselves. It depends on the ROracle package (which in turns requires an Oracle client installation).

Currently the functionality is relatively limited, focusing on my common tasks. Thus for example, given assay plate barcodes, we can retrieve the assay ids that the plate is associated with and then for a given assay, obtain the cell-level image parameter data (or optionally, aggregate it to well-level data). This task is easily parallelizable – in fact when processing a high content RNAi screen, I make use of snow to speed up the data access and processing of 50 plates.

Alternatively, given a plate id (this is the internal MetaXpress plate id) and a well location, one can obtain the path to the relevant image(s). With the images in hand, you could use EBImage to perform image processing entirely in R.

Currently, you cannot get the internal plate id based on the user assigned plate name (which is usually different from the barcode). Also the documentation is non-existant, so you need to explore the package to learn the functions. If there’s interest I’ll put in Rd pages down the line. As a side note, we also have a Java interface to the MetaXpress database that is being used to drive a REST interface to make our imaging data accessible via the web.

Of course, this is all specific to the ImageXpress platform – we have others such as InCell and Acumen. To have a comprehensive solution for all our imaging, I’m looking at the OME infrastructure as a means of, at the very least, have a unified interface to the images and their meta data.