Assessing open source IUPAC Name generators

(Updated on 26/10/2026 to ensure that the Orthonym experiment uses similar settings as for NISPO and OpenClatura (rather then just PIN). The coverage results are now comparable.)

Open source IUPAC Name generators have something in common with buses; you wait forever for one, and then three come at once.

I am talking about NISPO, OpenClatura and Orthonym, which were released in quick succession over the summer/autumn, together with preprints that are likely destined for a journal near you. These were developed by Nicolas Runcie/Charlotte Deane (Oxford), Rostislav Fedorov/Adrian Mirza/Kevin Maik Jablonka (Helmholtz Jena/Berlin), and Kohulan Rajan/Achim Zielesny/Christoph Steinbeck (Jena), and at their essence they take a SMILES string and produce an IUPAC Name.

Many readers will be familiar with OPSIN, maintained by Daniel Lowe (and hosted online by my team), which has been around for a while and does the opposite transformation from IUPAC Name to structure; in fact, "NISPO" is "OPSIN" backwards :-). There are a small number of commercial IUPAC Name generators (one less now, with Revvity's purchase of ACD), but until this year nothing usable in the open space; previous open approaches have used transformer models to directly translate between SMILES and IUPAC Name, but these require a large training corpus and no sufficiently large open corpus exists (though see this set we extracted from patents).

If it's not obvious, the reason that three independent programs have arisen this year to generate IUPAC names is because frontier AI models are only now able to tackle problems of this complexity. NISPO, OpenClatura and Orthonym were developed independently (as far as I know), and so present an interesting glimpse into how the same problem lends itself to different solutions.

I'm not going to go into the details of each approach, nor am I going to compare features and so forth. Instead, I will do my own biased comparision of the coverage of these tools, and make some judgements about the quality of the names. So recall vs precision in a sense, or quantity vs quality. Before we begin, I note that these tools are under active development and at least one has made several releases to address bugs/features since I did this work. With this in mind, let me state that I am testing NISPO 0.1.7, OpenClatura 0.3.1 and Orthonym 1.0.9.

I am interested in the IUPAC names that biologists might want to see. To this end, my starting point is a set of IUPAC names extracted from PubMed Abstracts by LeadMine. I retain only those that OPSIN can parse and that are not indicated as ambiguous by OPSIN. This latter is important because we are going to be comparing to these extracted names later, and we want to make sure that the names and the SMILES strings have a 1:1 correspondance. All names resulting in the same SMILES were collated yielding 341,943 SMILES strings (with associated names).

Coverage

Let's start with how many names were generated.

	NISPO	OpenClatura	Orthonym
For how many was an IUPAC name generated?	336167	313588	327227
...that was read by OPSIN?	333468	312068	325908
...considered unambiguous by OPSIN?	327458	311878	325883
...where InChIs agreed?	319327	292223	325869

So, comparing the last line, a slight difference between them in the order: Orthonym ≥ NISPO > OpenClatura.

Quality

For any non-trival molecule, many different IUPAC names can be written, some of which are more painful than others. For example, ethane could be written as 1-methylmethane [1]. This is definitely not a Preferred IUPAC Name (PIN), but it's an IUPAC name just the same. How can we compare the generated names along this axis? First of all, we'll assume that shorter is better. Secondly, we'll assume that fewer locants is better. Finally, we can check whether the generated names match any found in a PubMed Abstract for that SMILES string. Let's start with this final test, and extend the table above:

Does the name match one of the original names?

8974

28439

33567

Orthonym is generating names that best recapitulate the original name. OpenClatura is a bit behind, while NISPO names are least likely to match the original.

For the name length and locant count comparison, we will restrict ourselves to the subset of 187489 entries where all three generate a name. If we consider the net pairwise comparison (i.e. X > Y minus Y < X, ignoring X = Y), then the original name is shortest and has the fewest locants. In terms of length, NISPO > OpenClatura > Orthonym (shortest). For locants, the original name has the fewest; for the others it's NISPO = OpenClatura > Orthonym (fewest). [2]

Slow cases

I experienced a small number of very slow cases for each software (which I caught with timeouts of 10s): 8 for NISPO, 14 for OpenClatura and 3 for Orthonym. The NISPO cases were around highly-branched symmetric molecules, while Orthonyn and OpenClatura's cases were highly-conjugated systems. I have reported all of these to the developers.

Concluding remarks

It's great to see three different writers developed in different ways and having different characteristics. Each will be able to use the others to identify strategies to improve quality or coverage as the case may be. It's very much early days and these tools are still in active development with versions released since I started looking into the comparison in early Sept.

The development of these three tools has been enabled by OPSIN; at the same time they are limited (at least in terms of round trip evaluation) to whatever OPSIN supports which, while extensive, is not exhaustive. OPSIN itself could perhaps benefit from a 'strict mode' on reading so that writers are forced to constrain themselves to strict IUPAC names. I think we are about to see a period of improvements to OPSIN in tandem with the development of these tools as each challenges the other. This seems to be already happening with the contribution of support for phane nomenclature by the author of...wait a sec...yet another IUPAC name generator?!! [3]

Footnotes

[1] I'm proud of the "1-". The "methylmethane" is an okay example, but the "1-" is the cherry on top that really grinds the gears of a chemist.

[2] I have thought about providing the original matrix of values, but you'll just find it confusing - I do myself each time I look at it.

[3] I thought about including this one also in the comparison but there's no associated preprint, and so it might not be ready. Also, at this rate, I'll never get this blog post out if I keep stopping to include new IUPAC name generators!

ANNalog at ACS Fall 2026 Chicago

Earlier this week, I presented on Wei Dai's work on ANNalog at ACS Fall 2026 in Chicago. ANNalog is a generative model that takes an input SMILES string and generates MedChem-similar molecules.

Among other things, I covered how we tweak the training set to encourage stereo to be retained, why training on aligned SMILES is important, and why we don't canonicalise the input.

I ended with some practical tips: 

• Changes tend to occur on the RHS of the generated SMILES
• It follows that by presenting the input SMILES string in different ways, the user can selectively focus ANNalog’s attention on different parts of the molecule
• If changes should be distributed evenly, an option “-e variants” generates random SMILES variants and pools the results
• ANNalog allows the user to supply a fixed SMILES prefix to ensure that a particular portion of the molecule is unchanged
• Sampling can be used to create more diverse output; 
• Feeding the result through more than once can increase diversity further (see option “-e recursive”)
• Not everything is gold
• chembl_gen_check results are appended to the output to flag unusual structures

The wrong type of SMILES

When I wrote the blog post on ANNalog, there was an instant flurry of interest. ANNalog is a generative model that suggests medchem analogs given an input SMILES string. I was confident in the work that Wei Dai had done, but you can't help feeling nervous that there might be something you've missed. And sure enough, an issue was filed on GitHub the next day and instantly my heart sank:

I had a shot at ANNalog by submitting losartan as query:

annalog-generate -i losartan.smi -n 50 -m beam -o losartan_50_beam.tsv --device cuda

Frankly I find the results unexpected to say the least (see below). The top-ranked generated molecule looks very odd IMO. I did a substructure search in ChEMBL on the 4-membered ring with 3 nitrogens (C1=NN=N1) and (as expected) get zero hits. 

This was filed by Evert Homan on a Saturday morning. I was away from my computer, but pasting a few SMILES strings into CDK Depict instantly showed that, if anything, he was being generous:

 

But how could this be? I had never seen any results this bad, especially in beam search which typically stays very close to the input string. Then I realised something...

The SMILES string provided was "CCCCC1=NC(=C(N1CC2=CC=C(C=C2)C3=CC=CC=C3C4=NNN=N4)CO)Cl". Do you notice anything? There are no aromatic atom symbols (e.g. lowercase 'c') despite losartan having multiple aromatic rings - this is a Kekulé SMILES string. And therein lay the problem. ANNalog was trained on aromatic SMILES strings and had never seen aromatic rings written as conjugated bonds. When presented with this, there were no high probability predictions, only bad choices.

The potential fix was simple. Before making the prediction, we just needed to convert to an aromatic SMILES string. So I ran it through 'obabel' to see if this worked:

$ obabel -:"CCCCC1=NC(=C(N1CC2=CC=C(C=C2)C3=CC=CC=C3C4=NNN=N4)CO)Cl" -osmi
CCCCc1nc(c(n1Cc1ccc(cc1)c1ccccc1c1n[nH]nn1)CO)Cl

Using this string as input brought the results back into the realm of reality: 

But something was still a bit off. Almost all of the top ten results have fewer rings than the input. Also, the probabilities of these structures (given the input) is still in the same range (the highest score is -7.2 this time versus -7.4 last time).

This is perhaps a good moment to explain that ANNalog is trained on RDKit SMILES, after randomly shuffling the atoms (as we want to handle SMILES strings written in different ways). Here's the corresponding RDKit SMILES string:

CCCCc1nc(Cl)c(CO)n1Cc1ccc(-c2ccccc2-c2nn[nH]n2)cc1

Comparing to the Open Babel SMILES above there are various atom order shuffles, but if you look beyond that, you will notice two hyphens where OB had none. RDKit adds hyphens to indicate a single bond between aromatic rings that is not itself aromatic, whereas Open Babel only does this if the bond is itself in a ring. This doesn't matter in a semantic sense (*), but ANNalog has learnt to expect to see a hyphen separating aromatic rings, and so gives poor results when the provided SMILES string does not match this pattern. Using the RDKit one, the results are finally as boring as expected from beam search (and the highest score is now -3.8):

 

There's a lesson here around how it's important to normalise input to the same form used in training. Also, how subtle differences in SMILES strings can have an outsized effect. And a useful reminder not to release on a Friday. But ultimately, users of ANNalog can now provide SMILES in whatever format they like, and it should just work as expected.

 

* When reading a SMILES string, an implicit bond joining two aromatic atoms is considered a single bond except where it is in a ring itself, in which case it is considered an aromatic bond. Hence why you *have* to use a single bond symbol when it is in a ring and joins two aromatic atoms, but it is optional otherwise.

ANNalog, a generative model for MedChem analogs

I'm delighted to announce that Wei Dai's work on ANNalog, has just been published in Journal of Cheminformatics (currently early access). This is a Python application that takes a molecule represented by a SMILES string and generates MedChem analogs using a deep neural network trained on pairs of molecules from the same ChEMBL assay. This work comes from Wei's PhD with Arianna Fornili at QMUL and Nxera as industry partner (Jon Tyzack), where I continue to act as co-supervisor.

The code is available on GitHub. I won't recapitulate the README in the repo, but I'll mention a few points which are not covered there fully.

To begin with, in my case I need to use 'uv' instead of 'conda' to install (due to conda's licensing conditions). Here's how I do it:

$ uv venv annalog_env --python=3.12 
$ source annalog_env/bin/activate
(annalog_env) $ uv pip install numpy==2.4.3 pandas==3.0.1 tqdm==4.67.3 torch==2.10.0 torchvision==0.25.0 rdkit==2025.9.6 scikit-learn==1.8.0 annalog

Once installed, whether with conda or uv, here's a basic example of use, that generates the 10 most probable analogs given a single SMILES as input (a SMILES file is also accepted):

$ annalog-generate -i "CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12" -n 10
input_smiles	rank	generated_smiles	score
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	1	CCCc1nn(C)c2c(=O)[nH]c(-c3ccccc3OCC)nc12	-4.036181999828045
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	2	CCCc1nn(C)c2c(=O)[nH]c(-c3cc(CN4CCN(C)CC4)ccc3OCC)nc12	-5.025642307042602
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	3	CCCc1nn(C)c2c(=O)[nH]c(-c3cc(N4CCN(C)CC4)ccc3OCC)nc12	-5.148663511712925
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	4	CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N)ccc3OCC)nc12	-5.308277614323515
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	5	CCCc1nn(C)c2c(=O)[nH]c(-c3cc(NS(=O)(=O)C)ccc3OCC)nc12	-5.468033235034966
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	6	CCCc1nn(C)c2c(=O)[nH]c(-c3cc(CN4CCOCC4)ccc3OCC)nc12	-5.6292664217390325
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	7	CCCc1nn(C)c2c(=O)[nH]c(-c3c(OC)cccc3)nc12	-5.676750207183716
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	8	CCCc1nn(C)c2c(=O)[nH]c(-c3ccccc3)nc12	-5.692219721810034
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	9	CCCc1nn(C)c2c(=O)[nH]c(-c3cc(Cl)ccc3OCC)nc12	-5.694125330995178
CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12	10	CCCc1nn(C)c2c(=O)[nH]c(-c3ccc(N4CCN(C)CC4)cc3OCC)nc12	-5.8864854959632 

The example above uses (classic) beam search. In the course of implementing the library, we realised that the term "beam search" seems to mean different things to different people, and typically you need to look at the code to see what they actually meant. ANNalog implements two variants of beam search. With classic beam search (--method beam), at each token position, all beam_width candidates are expanded simultaneously as a batch, then pruned back to beam_width. This makes it fast but greedy, and it can miss the globally optimal sequence by pruning it too early. With best-first beam search (--method BF-beam), a priority queue is used to always expand the single highest-probability partial sequence first, regardless of length. This is slower and more memory-intensive, but more likely to find the globally best sequences.

As a general rule, modifications tend to occur on the right-hand side of the SMILES string. By rewriting the SMILES string in a particular way, this can be used to direct modifications to a certain part of the molecule; this can be enforced more formally by required the start of the string, the prefix to be fixed (see "--prefix"). Conversely, if you want to spread modifications evenly across the molecule, you may wish to pass in multiple SMILES variants or have the script do this for you (see "--exploration-method variants" and "--variant-number"). This is shown by the following example that samples from the distribution:

annalog-generate -i "CCCc1nn(C)c2c(=O)[nH]c(-c3cc(S(=O)(=O)N4CCN(C)CC4)ccc3OCC)nc12" -n 10 --method sample --exploration-method variants

When sampling, the default temperature is 1.2. Increasing this too far may not be a good idea as it increases the chance of unlikely tokens being sampled. Instead, if you wish to explore the search space further, it might make sense to reduce the temperature a bit (e.g. 1.1) and feed the output of one run back in again. The "--exploration-method recursive" does this automatically, but if you want to want to combine this with variants the easiest way is to write the output to a file, pull out the generated SMILES, and feed them back in as a file.

By default, invalid SMILES are filtered. This uses the partialsmiles library (which I've previously described here) to avoid selecting tokens during the generation process that would result in semantically or syntactically invalid SMILES. Obviously, if the model had the ability to perfectly understand SMILES, this would not be necessary but at least this approach is more efficient than filtering after-the-fact. You can turn off this filter if you wish ("--keep-invalid") to see how much the results change.

Similarly, not every molecule generated is gold - there's a certain (I would say low) percentage of dubious structures generated. Structures are checked automatically by default using Eloy Félix's chembl_gen_check to run a set of tests against the generated structures (including code adapted from Wim Dehaen's LACAN). Depending on your starting structure and use case, you wish to use these as hard filters or to rank prior to visual inspection.

We have tried to make sure that this tool is of practical use in a drug discovery setting. Let us know how you get on.

Improvement in reasoning performance of LLMs over time

If you tried using ChatGPT when it first came out and concluded that it wasn't much use for a scientific reasoning task, it might be time to try it again. For the task I semi-described in the previous post, it worked exceptionally well. Here I've tried to gather some information on what's changed between the original release and today, and how you can keep track of improvements in performance.

Background

The LLM developed by OpenAI is called GPT ("Generative Pre-trained Transformer"). It was released as a chatbot publicly in Nov 2022 and it changed the world (for better or worse). This release was already version 3.5. Since then, LLMs have been developed by several companies/groups, often with new versions multiple times a year.

For scientific problems, a particular landmark was the development of LLMs with 'reasoning' capability; the ability to perform multi-step logical problem-solving. OpenAI's o1 model (Sept 2024) was the poster child for this. Currently OpenAI's GPT-5.2 represents the state-of-the-art, though the upcoming Gemini 3 Pro (currently in preview) may have the edge on it. But how do we know this?

The benchmark

LLM benchmarks exist for coding tasks, solving maths problems, medical exams and so forth, but for reasoning, a single benchmark - the GPQA benchmark - has established itself as the gold standard, in particular the "Diamond" subset thereof. GPQA stands for "Google-Proof Q&A". In Nov 2023, Rein et al (New York Uni) published a preprint describing how they employed incentivised PhD-level contractors to put together a set of scientific multiple-choice questions (4 choices) covering biology, chemistry and physics that could not be answered by Google searches and required domain knowledge. Random answering scores 25%; PhD-level non-scientists (note: my interpretation) with access to the internet scored 34%; PhD-level scientists scored 69.7% (note: I don't see this figure in the preprint, but it is widely quoted).

Just a note about the 69.7%. The exact figure doesn't matter, but there's some idea that beating this figure means superhuman results. Um, no. An expert in biology is not expected to answer the physics questions very well. Are coding LLMs superhuman if they are better than a Python expert at answering Java questions? If they beat an expert at their own area of expertise, now that would be impressive. Given that it is expected that three experts sitting together would get 100%, it follows that the benchmark is not capable of measuring superhuman performance (this is actually mentioned at the end of the preprint). A more interesting number perhaps is 55% which roughly gives a measure of how well a scientist does when answering questions outside of their domain (I estimated this value by subtracting the 33.3% of correct answers inside the expert's domain from 69.7% and scaling it up, i.e. (69.7 - 100/3) / 2 * 3).

The results

The key point is that this is a set of difficult scientific questions that require reasoning, and that we can track the performance of LLMs over time against this benchmark. Here are the results for OpenAI models taken from an analysis by Epoch AI. Each individual model was tested 16 times, with the typical stderr around 2%:

• GPT-3.5 (Nov 2022): 28% (taken from here) 
• GPT-4 (Mar 2023): 36%
• GPT-4 Turbo (Nov 2023): 42%
• GPT-4o (Aug 2024): 49%
• o1-mini (Sep 2024): 62% (at "high" reasoning setting)
• o1 (Dec 2024): 76% (medium)
• o3-mini (Jan 2025): 77% (high)
• GPT-4.1 (Apr 2025): 67% (note: the 4.1s are a non-reasoning model)
• GPT-4.1 mini: 66% (this is a smaller/faster/cheaper version of GPT-4.1)
• GPT-4.1 nano: 49%
• GPT-5 (Aug 2025): 85% (medium)
• GPT-5-mini: 72%
• GPT-5-nano: 67%
• GPT-5.1 (Nov 2025): 85% (medium) 88% (high)
• GPT-5.2 (Dec 2025): 88% (medium) 88% (high) 91% (xhigh)

Conclusion

In summary, there have been major improvements over the last three years in terms of reasoning. Nor has it reached a peak. Gemini 3 Pro Preview (Nov 2025) is reported as 93%. In August, a new kid on the block, Autopoiesis, reported 92.4% for their Aristotle X1 Verify model (but this is not publicly available nor independently verified). It feels like we will soon get to the point where the field will need a more difficult benchmark.

Finally, it should be noted that pure performance is not the only relevant criteria. Price, energy usage, and ability to run locally will also be important considerations. Let's see what 2026 brings.

Image credit: Epoch AI, CC BY, https://epoch.ai/benchmarks/gpqa-diamond.

Classifying PubMed Abstracts with LLMs

I've just spent the Christmas playing with the OpenAI API, and I am impressed at what is now possible. There are 38M PubMed Abstracts (*); through judicious use of keyword filters (bringing it down to 250K) followed by a pass through a cheap model (bringing it down to 3.5K), I have essentially categorised all of the abstracts of interest with GPT-5.2 for less than $20 (**). Even leaving out the cheaper model and going directly to GPT-5.2 would still only be less than $200.

(If you read this before Jan 11 2026, and think that this topic sounds interesting, I encourage you to apply for openings in my group !)

This is what's known as a reasoning task. Here is the general prompt that I used, except that I have hidden some of the details (the definition section was suggested by ChatGPT as an improvement but I don't know that it helped).

Based on this article title and abstract, does the associated paper describe X?

Definitions:
- Y
- Z

Instructions:
- Reply with exactly one of: Yes, No, or Maybe.
- If the answer is No, reply with only: No
- If the answer is Yes or Maybe, briefly justify your answer by quoting specific words or phrases from the title or abstract as evidence.
- Do not infer beyond what is explicitly stated in the text.

When running it through GPT-5.2, I used slightly more complicated instructions in order to identify primary research data::

Instructions:
- Reply with exactly one of: Yes (primary), Yes (secondary), No, or Maybe.
  - Yes (primary) indicates that the paper is a primary reference for X, rather than being a review article or mentioning it as background
- If the answer is No, reply with only: No
- Otherwise, briefly justify your answer by quoting specific words or phrases from the title or abstract as evidence.
- Do not infer beyond what is explicitly stated in the text.

The cheaper/faster/less accurate model I used was GPT5-Nano. As a validation, I put all of the 2023 and 2024 abstracts through GPT-5.2 and only found a single additional "Yes" (on top of the 53 found by GPT5-Nano), and I'm not convinced about that particular case. GPT-5.2 was definitely required though as a final pass to reduce false positives; I passed all of the Yeses and Maybes from the Nano through to GPT-5.2 and it mostly trimmed down the Yeses and didn't promote the Maybes.

Run-time is somewhat dependant on how many batches you need to use as there is an overhead associated with each batch. When you sign up for API access, you start on Tier 1 and need to use lots of batches which must be run serially to avoid queue limits. In one case, I ran 75 batches of 3000 with GPT5-Nano; each one took 20 mins on average which was fine, but the API SLA only guarantees an answer within 24h would would be bit of a problem. After 7 days, I was bumped up to Tier 3 and it's possible that I could do the same thing in a single batch now. For example, on Tier 3, I ran 4 batches of 10000 with GPT5.2, each one taking 1h.

Given that the prompt above lends itself to many questions, you can use the details I've provided to get an idea of what might be possible for your questions of interest and for what budget. The trickiest part right now is dealing with the batch API, queue limits, and JSON, which can all be a bit tedious until you have something in place. Indeed, it feels like the sort of thing where there's a business opportunity to provide a biologist-friendly interface to abstract (***) this away. In any case, it's easy to predict that in 2026 we will be seeing less of BioBERT and similar models, and more adoption of LLMs.

Footnotes:

* in the 2025 baseline file - soon to be updated.

** $17.07 for GPT5-Nano over 250K abstracts and $1.82 for GPT-5.2 over 3.5K abstracts.

*** for once no pun intended!

Openings for Scientific Developers in the Chemical Biology Services team

As part of a collaboration with Open Targets on the development of a Side Effect Resource to guide target selection, we have two positions across NLP engineer/scientific developer/data engineer (i.e. your expertise might be in one or the other). Cheminformatics and bioinformatics would both be within scope.

We also have a more senior role, the position of Technical Lead in our team. This could be someone from either a scientific or developer background, but with appropriate experience dealing with technical infrastructure (Kubernetes, etc.). From previous experience, we know that this is a difficult role to fill so I would very much encourage suitable candidates to apply - if you are unsure, then reach out (oboyle@ebi.ac.uk) and we can discuss.

And as I always point out, due to EMBL-EBI's special status, the quoted salaries are tax-free and so are equivalent to the net salary from another job. Benefits, especially for non-residents, are generous also. Both positions have 11th Jan as the closing date, so apply now.

A new job, a postdoc opportunity, an open biological curator role, and a user group meeting

Almost exactly 6 months ago, I took over the leadership of the Chemical Biology Services team at EMBL-EBI, Hinxton, UK. This is the team that looks after ChEMBL, ChEBI, SureChEMBL, UniChem and OPSIN (check out our webpage, blog and follow us on LinkedIn, Bluesky). I could say a lot more about the honour of being appointed to this role, the responsibility I feel to the community, how this is the dream job for a cheminformatician, but let's keep it short!

And you can join this team too! If you want to be part of improving these services for the community, then please get in touch regarding the ARISE2 postdoc call described on our blog. There are a lot of cool things that we could do.

And that's not the only opportunity. We have just opened a role for a biological data curator. This is a key role in terms of maintaining the quality of our data. Currently this work is not very visible to the community, and we want to change that, which is one of the reasons we want to build closer links with everyone via....

...the ChEMBL UGM (!), a user group meeting around our services. This will take place on Jun 10-11 2026. You will be hearing more about this on our blog, likely in September. Apart from everything else you might expect from a UGM, this particular meeting at the start of my 9-year tenure will play a key role in influencing future development of our services. If you want to be informed when more information is available and registration opens, please email chembl-ugm@ebi.ac.uk.

A 16-year time loop

Next week I wlil be attending the Computer Aided Drug Design Gordon Research Conference (CADD GRC) in the US (Maine).

The last (and only) time I attended this meeting was 16 years ago. At that time I was a postdoc at the CCDC working on the GOLD docking software, and presented a poster "Why multiple scoring functions can improve docking performance". Right after my 5-minute flash presentation (see program), some random English guy presented his poster about a database with a funny name that was about to be released.

This was July 2009; three months later ChEMBL 01 would be released.

Even at the time, I recognised that this was a big deal. In academia, research was difficult if not impossible because of the lack of data. Meanwhile industry published research that used their internal data and couldn't be challenged, compared, or built upon. To quote a recent conversation with an industry figure, showing results on ChEMBL keeps everyone honest.

But as well as learning about ChEMBL, I got to know John Overington, who subsequently invited me to present on protein ligand docking and cheminformatics at various training courses he organised at the EBI, even after I moved back to Ireland. For me, this was a really great connection to have as it increased my profile, and I got to meet many of the leading figures in the field. In return, John got the occassional bug report emailed directly to his inbox which I'm sure was exactly what he wanted. :-) We've kept in touch over the years, and he has been a great help when I've needed it.

And so, I am returning to where it all started. Only this time I am the random English Irish guy with the poster about ChEMBL. And who knows, maybe I, in my turn, will meet a future custodian of ChEMBL...?

MMPA: The case of the missing group Part II

I was reading over my recent blogpost on MMPA when I started wondering whether the approaches I described actually were the best ways of finding the correspondences with the 'missing groups'. Just looking at the diagrams, in each case there's another approach that jumps out...

In the original post, I described taking the third case and manufacturing matches to the first two. What if, instead, we take the first two and look for matches to the third? That is, chop off the asterisks (where attached to a ring), adjust the hydrogens and then see whether there exists a match to a full molecule:

This approach may well be more efficient as the number of implicit hydrogens in a typical molecule is significant.

In the original blogpost, I again described taking the third case and manufacturing a match to the first two. An alternative would be to take the constant parts (i.e. the singly-attached groups) of all potential series, join the two together and then look for a match to a full molecule.

In this case, it's not so clear-cut that this would be more efficient than the original approach described. But if I get around to implementing these alternative approaches, I'll report back...

SmiZip 2.0 released

SmiZip is a compressor for SMILES strings that uses the unused characters from the extended ASCII set. So, for example, while the character "C" is already used for Carbon, part of Chlorine, etc.,  the character 'J' could be used to indicate "C(=O)".

While most of the v2.0 changes (listed below) are tidy-ups or minor improvements, there is one in particular that simplifies a particular use case, to encode only into ASCII printable characters. This means that the compressed string can be stored in places where the full 256 characters are not supported, e.g. a JavaScript file (you can encode such characters with Unicode but that would increase the file size).

As an example, I trained on a dataset derived from ChEMBL and was able to compress to 30.8% on a hold-out test set using all available characters (202 multigrams). But even sticking to the printable ASCII character set, I was able to compress to 43.7% despite only having 43 characters available for multigrams. We benefit here from the fact that there are diminishing returns with each additional multigram - the majority of the benefit is derived from the initial multigrams.

As a specific example, here's CHEMBL505931 which is a polysaccharide attached to a sterol and compresses from 322 characters to 138 (42.9%):

[C@]12(O[C@@H](C[C@H]1C)C(=O)CC)[C@]1(CCC3=C(CC[C@H]4[C@@]([C@H](CC[C@]34C)O[C@@H]3O[C@@H]([C@H]([C@@H]([C@H]3O)O)O)CO[C@@H]3OC[C@@H]([C@@H]([C@H]3O[C@@H]3O[C@@H]([C@H]([C@@H]([C@H]3O[C@@H]3O[C@H]([C@@H]([C@H]([C@H]3O)O)O)C)O[C@@H]3OC[C@@H]([C@@H]([C@H]3O[C@@H]3OC[C@H]([C@@H]([C@H]3O)O)O)O)O)O)CO)O)O)(C)CO)[C@@]1(CC2)C)C
compresses to
]12(O>,D1C$ *):]1,*3=~*D4:@](D,C:]34CV>3O>(D(>(D3OVV$O>3OC>(>(D3O>3O>(D(>(D3O>3OD(>(D(D3OVV$V>3OC>(>(D3O>3OCD(>(D3OVVVVV$OVV),$O):@]1,C2$$

Changes: v2.0 (2025-01)

• find_best_ngrams.py: new option --non-printable to facilitate encoding into printable ASCII charactersr
• find_best_ngrams.py: --chars is now required (help text provides a reasonable starting point) to force the user to consider the list
• find_best_ngrams.py: if the end of the training .SMI file is reached, the script wraps around to the start
• find_best_ngrams.py: --cr corrected to --lf.
• compress.py: A better error message is generated if an attempt is made to encode a character not present in the JSON file
• compress.py: support added for .SMI files without titles.

Thanks to Adriano Rutz (@adafede) and Charles Tapley Hoyt (@cthoyt) for feedback.

MMPA: The case of the missing group

An efficient approach to Matched Molecular Pair Analysis (MMPA) consists of fragmenting molecules at a subset of single bonds and then collating the results to find cases where there's the same scaffold (or linker) but a different R group (or groups). If my memory serves me correctly this algorithm was first described by Wagener and Lommerse, but it was the publication by Hussain and Rea that really popularised it.

What I want to discuss here is the importance of handling the 'missing group'. To simplify things, let's fragment only at non-ring single bonds that are attached to rings.

Single cut

Consider the fragmentation of the following three molecules:

Molecule A: 

 Molecule B:

Molecule C:

For molecules A and B, the algorithm will break the bond to the F and Cl, replacing it by a dummy atom; the two scaffolds are identical and so A and B will be identified as a matched pair.

We also want C to be identified as part of the same matched series, but because the hydrogens are implicit, there's nothing to break off that would yield the matching scaffold. You could convert all implicit Hs to explicit which, well, would work; however you would have to hand in your cheminformatics badge at the desk due to the performance hit :-). Instead, we handle this empty group separately by iterating over the ring atoms that have an implicit H and generating the same scaffold fragments that would have been generated in the explicit case, each time associating it with the *[H] R group:

This is shown here for Molecule C, but the same procedure would apply to A and B. In fact, a significant proportion of the run time will just be dealing with this case as hydrogens on rings are not uncommon.

Double cut

The typical use of double cuts is to find molecules that are the same except for replacement of a linker between two rings. After a first cut, a second cut is done on the results if possible.

Molecule D

Molecule E

Molecule F

Molecules D and E share the same scaffolds, and so their linkers are identified as a matched pair. This isn't the case for Molecule F as it only gets as far as a single cut after which there is nothing left to cut for the next step.

The solution is to treat all single cut results as an additional outcome of the double cut process if (and only if) both ends are attached to a ring. After renumbering one of the locants, we associate the resulting scaffold pair with the empty linker: *-*. This procedure doesn't change anything for D and E, but gives the necessary result for F:

Are you finished?

Yes. Fortunately, these corner cases don't affect higher numbers of cuts. The smallest group that can have three or more cuts is a single atom, and this will be handled automatically.

ReplaceR - common R group and linker replacements

Introducing ReplaceR, a webapp that suggests R group or linker replacements that a medicinal chemist might consider making after the query, as well as those typically made before the query. This builds upon the ideas I described in a previous post.

In the course of a medicinal chemistry project, molecules progress through a series of changes that often involve replacement of R groups (or linkers). A key insight is that there is a certain order to these changes; in general, R groups that are simpler (in size, or in ease of synthesis, or just cheaper) will be made before those are more complex. It would be unusual to make the bromide without having already tried chloride; if you see a tetrazole, then a carboxylic acid was likely already made.

Given true time series data from projects, we could work out this information. However, even without this, we can infer similar information from co-occurrence data in matched molecular series across publications and patents in ChEMBL (see this 2019 talk and blogpost of mine for more details). To give a sense of the method, consider that butyl co-occurs with propyl in matched series in 801 papers, but given that propyl occurs in 2765 series versus butyl's 1905 we infer that butyl is typically made after propyl (if at all). Groups in the 'after' panel are sorted by co-occurence with the query. Groups in the 'before' panel are instead sorted by co-occurrence/frequency - if we didn't do this, hydrogen or methyl would always be top of the list!

Note that no biological activity data or physicochemical property data (e.g. MWs) is used in this analysis. The suggestions are not based on bioisosteres, SAR transfer or QSAR models, but are simply a reflection of historical choices made by medicinal chemists across a large number of publications. I note that the 2021 papers by Mahendra Awale et al and Kosuke Takeuchi et al touch on many of the same issues.

This webapp may be useful to medicinal chemists. But in particular, computational chemists and cheminformaticians will find the app and the data behind it useful to support automatic enumeration and virtual screening of likely next steps. Note that the webapp just lists the 'nearest' 12 groups before and after. The full dataset is available on Zenodo. You can find me at baoilleach@gmail.com or file an issue on GitHub.

Cheminformatics Credits: EMBL-EBI Chemical Biology Services Team for ChEMBL data, John Mayfield et al for Chemistry Development Kit and CDKDepict, Geoff Hutchison et al for Open Babel, Greg Landrum et al for RDKit, Chen Jiang for openbabel.js, Michel Moreau for rdkit.js, Peter Ertl and Bruno Bienfait for JSME, Noel O'Boyle (me) for ReplaceR.

Release of partialsmiles v2.0 - A validating parser for partial SMILES

A while back, I wrote a library called partialsmiles to parse and validate prefixes of SMILES strings. I've just updated it with a minor enhancement (it now returns the full state of the parser, not just the molecule) and thought I would show an example of typical use in the context of generative models.

For a change, let's not use neural nets. Using molecules from ChEMBL converted to random SMILES, I looked at all runs of 8 tokens and counted up the frequency of occurrence of the following token. For example, following Cc1ccccc, '1' was observed 340389 times, '2' 582 times (etc.); following COc1ccc(, 'c' was observed 33288 times, then 'C', 'N', 'O', 'Cl', etc. To handle the start, I padded with 8 start tokens.

Once these frequencies are in hand, we can use this to power a generative model where we sample from the distribution of frequencies at each step. Just like with a neural net, we can apply a temperature factor; this determines to what extent we exaggerate or downplay the difference between frequencies. A value of 1.0 means use the frequencies/probabilities as provided, while an extreme value of 10 would mean treat all (non-zero probability) tokens as equally likely. The other extreme of 0.1 would mean only ever pick the most likely token.

The rationale for partialsmiles is that we can perform checks for validity while generating the tokens, rather than waiting until we complete the SMILES string (i.e. when we generate an END token). If, at some point, the partial SMILES string turns out to be invalid we could just discard the string; this would speed things up but not change the result. An alternative, which I show below, is to avoid sampling tokens that will lead to an invalid SMILES; this will increase the number of valid SMILES thus making the process more efficient. Here I do this by setting the frequencies of those tokens to 0.

If we use the code below to generate 1000 SMILES strings and use a temperature factor of 1.0, the number of valid strings increases from 230 to 327 when partialsmiles is used. A smaller increase is observed for a temperature factor of 1.2 (205 to 247) and larger for 0.8 (263 to 434).

I note in passing that it's possible to force the model to avoid adding an END token unless it results in a valid full SMILES string (see PREVENT_FULL_INVALID_STRING below). I'm afraid that doesn't work out very well. Some really long SMILES string are generated that just keep going until the monkeys typing it eventually close the ring or parenthesis or whatever is the problem. Some things you just gotta let go, or patch up afterwards.

Oh, you want to see some of the generated molecules? No, you don't; you really don't :-) There's a reason neural nets are popular. :-) You can find the associated code on GitHub.

import pickle
import tqdm
import numpy as np
import partialsmiles as ps

np.random.seed(1)

PREVENT_FULL_INVALID_STRING = False

SOS, EOS = range(2)

TOKENS = [
  '^', '$', # i.e. SOS, EOS
  'c', 'C', '(', ')', 'O', '1', '2', '=', 'N', 'n', '3', '[C@H]',
  '[C@@H]', '4', 'F', '-', 'S', '/', 'Cl', '[nH]', 's', 'o', '5', '[C@]', '#',
  '[C@@]', '\\', '[O-]', '[N+]', 'Br', '6', 'P', '7', '8', '9']

TOKEN_IDXS = list(range(len(TOKENS)))

def find_disallowed_tokens(seq, freqs):
    disallowed = set()
    smi = "".join(TOKENS[x] for x in seq)
    for i, x in enumerate(TOKENS[2:], 2):
        if freqs[i] > 0:
            try:
                ps.ParseSmiles(smi + x, partial=True)
            except ps.Error:
                disallowed.add(i)

    if PREVENT_FULL_INVALID_STRING:
        if freqs[EOS] > 0:
            try:
                ps.ParseSmiles(smi, partial=False)
            except:
                disallowed.add(EOS)
    return disallowed

def generate(all_freqs, prefix_length, temperature, avoid_invalid):
    seq = [SOS] * prefix_length
    i = 0
    while True:
        idx = (seq[i]<<6) + (seq[i+1]<<12) + (seq[i+2]<<18) + (seq[i+3]<<24) + (seq[i+4]<<30) + (seq[i+5]<<36) + (seq[i+6]<<42) + (seq[i+7]<<48)
        freqs = [all_freqs.get(idx+j, 0) for j in TOKEN_IDXS]
        if avoid_invalid:
            disallowed_tokens = find_disallowed_tokens(seq[prefix_length:], freqs)
            freqs = [all_freqs.get(idx+j, 0) if j not in disallowed_tokens else 0 for j in TOKEN_IDXS]
            if all(freq==0 for freq in freqs):
                return ""
        probs = freqs/np.sum(freqs)
        p = np.power(probs, 1.0/temperature)
        q = p / np.sum(p)
        chosen = np.random.choice(TOKEN_IDXS, 1, p=q)
        seq.append(chosen[0])
        if chosen == 1: # EOS
            break
        i += 1
    return seq[PREFIX_LENGTH:-1]

if __name__ == "__main__":
    PREFIX_LENGTH = 8
    TEMPERATURE = 0.8
    AVOID_INVALID = True

    with open(f"../nbu/probabilities.{PREFIX_LENGTH}.2.pkl", "rb") as inp:
        all_freqs = pickle.load(inp)
        print("Loaded...")

    ofile = f"../nbu/out.{PREFIX_LENGTH}_{TEMPERATURE}_{AVOID_INVALID}.smi"
    with open(ofile, "w") as out:
        for i in tqdm.tqdm(range(1000)):
            tokens = generate(all_freqs, PREFIX_LENGTH, TEMPERATURE, AVOID_INVALID)
            smi = "".join(TOKENS[x] for x in tokens)
            out.write(f"{smi}\n")

push and pop directories

Here's a tool that I've been using daily since 2005 at least. I had a write-up on my old website, but with its recent disappearance it seems like a good time to update it and publish it here.

If you add the code below to your .bashrc, then you can type "push here" in one terminal window to record the current directory under the name 'here', and then "pop here" in another terminal window to change to the same directory.

In other words, it's a simple bookmark system for directories backed by a persistent on-disk file-based database (i.e. a text file :-) ). You may find this useful to support sshing into multiple machines that have a shared home folder, or to synchronise windows in a screen session or tabs in a terminal emulator.

See you all next year...🎉

popnamelist="$HOME/.popnames.txt"

# Supporting code for 'pop'
function popname () {

  if [ -z "$1" ]; then
      echo "  Syntax: popname TAG"
      return 1
  fi

  if [ -f $popnamelist ]; then
      grep "^$1" $popnamelist > $popnamelist.tmp
      if [ ! -s $popnamelist.tmp ]; then
          echo "No such tag"
      else
          awk '{print $2}' $popnamelist.tmp
      fi
  else
      echo "  There isn't anything to pop!"
      return 1
  fi
}

# Pushes the current directory into a 'memory slot' indexed by a tag
# See also 'pop'
function push() {

  if [ -z "$1" ]; then
      echo "  Syntax: push TAG"
      return 1
  fi

  touch $popnamelist
  grep -v "^$1" $popnamelist > $popnamelist.tmp

  echo "$1 `pwd`" >> $popnamelist.tmp
  sort $popnamelist.tmp > $popnamelist
  rm $popnamelist.tmp
}

# Pops the directory associated with 'tag' and makes it the current
# directory.
# Then you can use: pop mytag
# or echo `popname mytag`, or cat `popname mytag`/readme.txt

function pop() {
  if [ -z "$1" ]; then
    echo "  Syntax: pop TAG"
    echo
    cat $popnamelist
  else
    cd `popname $1`
  fi
}

The effect of ties on evaluation of virtual screens

For a virtual screening method where the associated score is prone to ties (e.g. a fingerprint comparison), care must be taken to handle these appropriately or the results will be biased towards better performance of either inactives or actives.

Let's work through a specific example, starting with some setup code assigning scores to a set of actives and separately to a set of inactives:

import random
random.seed(42)

INACTIVE, ACTIVE = range(2)

def mean_rank_of_actives(data):
    ranks = [i for (i, x) in enumerate(data) if x[1]==ACTIVE]
    return sum(ranks)/len(ranks)

if __name__ == "__main__":
    # Generate scores for actives and inactives
    actives = [1.0]*10 + [0.5]*10
    random.shuffle(actives)
    inactives = [0.5]*10 + [0.0]*70
    random.shuffle(inactives)

I've baked in a large number of ties. Half of the actives have a score of 0.5, a value that is shared by an equal number of inactives. We will use mean rank of the actives as a proxy here for ROC AUC - if I remember correctly, one is proportional to the other.

Our first attempt at evaluating performance is as follows:

def first(actives, inactives):
    everything = [(x, ACTIVE) for x in actives] + [(y, INACTIVE) for y in inactives]
    everything.sort(reverse=True) # rank the entries - highest score first
    return mean_rank_of_actives(everything)

It turns out that all of the actives are sorted ahead of the inactives, giving an overoptimistic mean rank of 9.5. Why? Because the ties are resolved based on the second item in the tuples, and ACTIVE (i.e. 1) is greater than INACTIVE (i.e. 0). In fact, swapping the values of these flags changes the results to the overpessimistic value of 14.5. Using a piece of text, e.g. "active" or "inactive", has the same problem.

No worries - when we sort, let's make sure we sort using the score values only:

def second(actives, inactives):
    everything = [(x, ACTIVE) for x in actives] + [(y, INACTIVE) for y in inactives]
    everything.sort(reverse=True, key=lambda x:x[0])
    return mean_rank_of_actives(everything)

Unfortunately, the result is still 9.5. Python's sort is a stable sort, which means that items retain their original order if there are ties. Since all of the actives come first in the original list, the tied actives still come before the tied inactives after sorting. To get rid of the influence of the original order, we need to shuffle the list into a random order:

def third(actives, inactives):
    everything = [(x, ACTIVE) for x in actives] + [(y, INACTIVE) for y in inactives]
    random.shuffle(everything)
    everything.sort(reverse=True, key=lambda x:x[0])
    return mean_rank_of_actives(everything)

This gives 12.0 +/- 0.67 (from 1000 repetitions), a more accurate assessment of the performance.

I think that this is an interesting little problem because it's rather subtle; even when you think you've solved it (as in the second solution above), some other issue rears its head. The only way to be sure is to test with tied data.

Note: If we are interested specifically in handling ties correctly in ranks, a common approach is to assign the same rank to all members of the tie. That's not really the point I want to make here. In practice, the evaluation metric might something quite different such as enrichment at 1%, which would not be amenable to such a solution.

Engagement on X vs Mastodon vs BlueSky

Egon posted recently on the situation with scientists leaving X and moving to another microblogging platform. He pointed out some pros and cons of BlueSky vs Mastodon and how it may be possible to bridge between the two.

This post doesn't really speak to the points that Egon raises, but rather gives a n=1 measure of the current level of community engagement on different platforms. I posted a link to my previous blog post on the three platforms mentioned in the title, and here's what I found:

• X - Since Sept 2011 - I have 1190 followers - 270 views (34 link clicks), 4 comments, 11 additional likes (including one from the poteen brewery I often get confused with :-) )
• Mastodon - Since Dec 2022 - I have 147 followers - 2 comments, 1 additional like
• BlueSky - Since 15 Dec 2024, i.e. one day before the blog post - I have 28 followers - 2 comments, 11 additional likes/reposts (including one from "Low Quality Facts" - make of that what you will)
• My blog - Since April 2007 - I no longer have visibility of the number of followers, but there are 57 people who use the follow.it link to follow by email - 428 views, 3 comments

First of all, the numbers confirm that X is no longer where the community is; despite being only 5 days on BlueSky and a handful of followers, the engagement just about matches X. Mastodon is an interesting platform, and great to have as an alternative should BlueSky not pan out, but the numbers and engagement have not taken off; it feels quiet over there, in a way that BlueSky doesn't. In each case, to actually find out the news, the user had to click through to the blog, so the figure of 428 there should give the overall total across everything (including LinkedIn where I also posted).

From the little time I have spent on BlueSky, it has the feel of the old Twitter. The Android app needs a bit of work (scrolling is very jerky on the main feed) but I'm sure it will be sorted. I moved early to Mastodon as I would have liked it to take off, but that hasn't happened. I note that BlueSky has plans to monetise somehow and time will tell whether this turns it into something closer to the current X. In the meanwhile, I guess I'll see you over there.

@baoilleach.bsky.social/@baoilleach@mstdn.science/@baoilleach

Moving to pastures new, but still in the same field Part IV

Following on from my previous post on the topic, I will shortly be leaving Nxera Pharma (Sosei Heptares as was).

The last five years at Nxera have passed quickly. I've gone from "what is this GPCR you speak of?" to providing tools and resources to support structure-based GPCR drug discovery, and indeed doing it myself. I'll miss working with my colleagues in CompChem and throughout Nxera, but I'm leaving in the knowledge that the Cheminformatics line is in very capable hands.

Which brings me to my next move...

I am honoured and super excited to take up the role of Chemical Biology Resources Team Leader at the EBI, leading the team responsible for ChEMBL, SureChEMBL, UniChem and ChEBI. There are some big shoes to fill, but I can't wait to start working with the fantastic team behind these resources. The start date is early February, so wish me luck!

Do NOT use SDFdb

It's not even a real database. All it does is:

import sdfdb
db = sdfdb.SDFdb("mysdf.sdf")
molfile = db.get_mol("mytitle")

Ok, so it's quick to index even quite a large SD file. But if Andrew Dalke ever sees the corners it cuts...! I mean, it doesn't even support Windows line endings.

And the code? Well, it speaks for itself:

import re

class SDFdb:
    def __init__(self, fname):
        self.fname = fname
        self.file = open(self.fname, "rb")
        self._create_index()

    def _create_index(self):
        patt = re.compile(b"\$\$\$\$\n(.*)\n")
        chunksize = 100000
        self.lookup = {}
        self.start = []
        idx = 0
        position = 0
        title = self.file.readline().rstrip().decode("ascii")
        self.lookup[title] = idx
        idx += 1
        self.start.append(position)
        position = self.file.tell()
        while text := self.file.read(chunksize):
            if text[-1] != '\n':
                text += self.file.readline()
            # Invariant: text ends with "\n"
            if text.endswith(b"\n$$$$\n"):
                text += self.file.readline()
            # Invariant: text never ends with "\n$$$$\n"
            for m in patt.finditer(text):
                title = m.groups()[0].decode("ascii")
                if title in self.lookup:
                    print(f"WARNING: Duplicate title {title}")
                self.lookup[title] = idx
                offset = m.start() + 5
                idx += 1
                self.start.append(position + offset)
            position = self.file.tell()
        self.start.append(self.file.tell()) # should be EOF

    def get_mol(self, title):
        idx = self.lookup.get(title, None)
        if idx is None:
            raise KeyError(f"Title '{title}' not found")
        self.file.seek(self.start[idx])
        text = self.file.read(self.start[idx+1]-self.start[idx])
        return text.decode("ascii")

    def close(self):
        self.file.close()

...though maybe I can see it would be useful for a large number of random accesses into a large (well-behaved) SD file.

That is, if they weren't already in a database. (Or *cough* they were in a database but you thought this was simpler than writing the code to batch up queries.)

Every ChEMBL everywhere, all at once

The team over at ChEMBL has just announced a symposium to celebrate 15 years of ChEMBLdb and 10 of SureChEMBL. With this in mind (or perhaps for other reasons), I have attempted to extract data from all 34 releases of ChEMBLdb. Let's see how I got on...

One of the things I love about ChEMBL is that they still provide all of the download files for all previous releases. But before w'get wget, let's see if there's an easier way, namely Charles Tapley Hoyt's ChEMBL Downloader - if it can do the work for me, then that'll save a whole bunch of work. Looking into it, however, the ability to run SQL queries relies on an SQLite database being available, which is only the case for ChEMBL 21 or later.

So old skool it is, which leads to the question of which database format to use. Well, it turns out that there's only one database format that is available for every single release of ChEMBL, and that's MySQL. One bash loop later, and I have all 34 chembl_nn_mysql.tar.gz files (well, two extra wgets for 22_1 and 24_1). A bit of unzipping later and my home folder is 300G heavier and we are ready to roll.

That is, if we can navigate the slightly differently named/arranged files in every case. Is the unzipped folder chembl_nn_mysql, or is it chembl_nn with the chembl_nn_mysql folder inside that? Are the setup instructions in INSTALL or INSTALL_mysql? Is the SQL data in a .sql file or a .dmp file? It is of course nitpicky to even mention this given the nature of what I'm getting for free from ChEMBL, but have you ever tried to install 34 ChEMBLs? :-) In the end, a carefully crafted find/grep was able to pull out the correct create command:

for d in `find . | grep INSTALL | sort`; do grep "create database" $d |sed "s/mysql>//g" ; done

This is simply 'create database chembl_nn' from version 1 until version 22, and then from 23 onwards it changes to specify the character set as UTF-8. For the actual import, I used a Python script to work out the right commands, wrote them to a shell script and then ran them. Just another 700G later, and it's all available in MySQL.

What I wanted to do was simply export a set of SMILES, InChIs and CHEMBL Ids from each of the versions. What I thought would be a problem turned out not to be at all; while the schema as a whole has grown and changed quite a bit over the years, the only effect for me was that 'compounds' became 'compound_structures' in ChEMBL 9.

More of an issue was that the data stored has changed over the years. I had forgotten that ChEMBL 1 did not have ChEMBL Ids, nor version 2; in fact, it wasn't until version 8 that they appeared on the scene. Before that, there was an attempt to use ChEBI ids and at the very start there was just the molregno. Furthermore, version 1 doesn't have Standard InChI - it just has a non-standard one (why, or what the settings are, I don't know ...Added 30/06/2024: Standard InChI was not available until shortly before ChEMBL 1 - that would explain it!). This at least we can pull in from version 2 by matching on molregno and SMILES, and then calculate any missing Standard InChIs off-line.

There was a point to all this, but that'll wait for another day. What I've shown here is that despite 15 years covering 34 releases, it's possible to recreate every single release and relive those glory days (ah, ChEMBL 16, will I ever forget you?).