new paper: "A transcriptomics data-driven gene space accurately predicts liver cytopathology and drug-induced liver injury"

Figure from the article. CC-BY.
One of the projects I worked on at Karolinska Institutet with Prof. Grafström was the idea of combining transcriptomics data with dose-response data. Because we wanted to know if there was a relation between the structures of chemicals (drugs, toxicants, etc) and how biological systems react to that. Basically, testing the whole idea behind quantitative-structure activity relationship (QSAR) modeling.

Using data from the Connectivity Map (Cmap, doi:10.1126/science.1132939) and NCI60, we set out to do just that. My role in this work was to explore the actual structure-activity relationship. The Chemistry Development Kit (doi:10.1186/s13321-017-0220-4) was used to calculate molecular descriptor, and we used various machine learning approaches to explore possible regression models. Bottom line was, it is not possible to correlate the chemical structures with the biological activities. We explored the reason and ascribe this to the high diversity of the chemical structures in the Cmap data set. In fact, they selected the chemicals in that study based on chemical diversity. All the details can be found in this new paper.

It's important to note that these findings does not validate the QSAR concept, but just that they very unfortunately selected their compounds, making exploration of this idea impossible, by design.

However, using the transcriptomics data and a method developed by Juuso Parkkinen it is able to find multivariate patterns. In fact, what we saw is more than is presented in this paper, as we have not been able to support further findings with supporting evidence yet. This paper, however, presents experimental confirmation that predictions based on this component model, coined the Predictive Toxicogenocics Gene Space, actually makes sense. Biological interpretation is presented using a variety of bioinformatics analyses. But a full mechanistic description of the components is yet to be developed. My expectation is that we will be able to link these components to key events in biological responses to exposure to toxicants.

 Kohonen, P., Parkkinen, J. A., Willighagen, E. L., Ceder, R., Wennerberg, K., Kaski, S., Grafström, R. C., Jul. 2017. A transcriptomics data-driven gene space accurately predicts liver cytopathology and drug-induced liver injury. Nature Communications 8. 
https://doi.org/10.1038/ncomms15932

Postdoc Position in Computational Metabolomics

A postdoc position is available in my research group at Friedrich-Schiller-University Jena, Germany.

Open Positions:

Postdoc: We are looking for a talented cheminformatician, bioinformatician or someone with comparable skills to work on the development cloud-based methods for computational metabolomics. The successful candidate will work closely with the H2020 e-infrastructure project PhenoMeNal, a European consortium of 14 partners. This position requires excellent skills in at least one modern, object-oriented programming language. A strong interest in metabolomics and cloud computing as well as the ability to work in a distributed team will be advantageous. The postdoc will also have the opportunity to participate in the day-to-day management of the group as well as in the organisation of seminars and practical courses for our students

The position requires a strong interest in metabolomics, molecular informatics and current IT technologies, programming skills a modern object oriented programming language and the ability to work in geographically distributed teams.

Please send applications in PDF format by email to christoph.steinbeck@uni-jena.de. We will accept applications until the position is filled.

Background information:

The Friedrich Schiller University Jena (FSU Jena), founded in 1558, is one of the oldest universities in Europe and a member in the COIMBRA group, a network of prestigious, traditional European universities. The University of Jena has a distinguished record of innovations and resulting educational strengths in  major fields such as optics, photonics and optical technologies, innovative materials and related technologies, dynamics of complex biological systems and humans in changing social environments. It has more than 18,000 students. The university’s friendly and stimulating atmosphere and state-of-the-art facilities boost academic careers and enable excellence in learning, teaching and research. Assistance with proposing and inaugurating new research projects and with establishing public-private partnerships is considered a crucial point.

About Christoph Steinbeck

 

New paper: "The Chemistry Development Kit (CDK) v2.0: atom typing, depiction, molecular formulas, and substructure searching"

This paper was long overdue. But software papers are not easy to write, particularly not follow up papers. That actually seems a lot easier for databases. Moreover, we already publish too much. However, the scholarly community does not track software citations (data citations neither, but there seems to be a bit more momentum there; larger user group?). So, we need these kind of papers, and just a version, archived software release (e.g. on Zenodo) is not enough. But, there it is, the third CDK paper (doi:10.1186/s13321-017-0220-4). Fifth, if you include the two papers about ring finding, also describing CDK functionality.


Of course, I could detail what this paper has to offer, but let's not spoil the article. It's CC-BY so everyone can readily access it. You don't even need Sci-Hub (but are allowed for this paper; it's legal!).

A huge thanks to all co-authors, John's work as release manager and great performance improvements as well as code improvement, code clean up, etc, and all developers who are not co-author on this paper but contributed bigger or smaller patches over time (plz check the full AUTHOR list!). That list does not include the companies that have been supporting the project in kind, tho. Also huge thanks to all the users, particularly those who have used the CDK in downstream projects, many of which are listed in the introduction of the paper.

And, make sure to follow John Mayfield's blog with tons of good stuff.

Sharp Tools for Java Refactoring: Byte Code Analysis

I'm currently refactoring parts of the CDK core classes. As part of this I often need to find specific patterns/idioms that need to be changed. Whilst source code inspections and an IDE can make this task easy sometimes the tools aren't quite sharp enough.

I needed to find all occurrences of a reference (instead of value) comparison on a particular class. In Java there is no operator overload and so you can have situations like this:

Integer a = new Integer(25);
Integer b = new Integer(25);
if (a == b) {} // false
if (a.equals(b)) {} // true
I mentioned operating overloading but it's more subtle and is more about comparing reference vs. value comparison. In C/C++ we can have similar behaviour:

int aval = 25, bval = 25;
int *a = &aval;
int *b = &bval;
if (a == b) {} // false
if (*a == *b) {} // true
Most IDE's and code inspection programs will warn about common occurrences (for example Integer) but I wanted to find places where the CDK's classes were used like this. A simple text grep will find some but will have false positives and negatives requiring lots of manual checking. Daniel suggested the well known FindBugs might be able help.

Rather than analyze source code like PMD and Checkstyle, FindBugs analyses Java byte code with a set of defaults detectors to find often subtle but critical mistakes. FindBugs can be configured with custom detectors (see here), however the inspection I needed (RC: Suspicious reference comparison to constant) was almost there. After digging around in the source code I found you can provide a list of custom classes to detect. However it took a bit of trial and error to get what I needed.

First up we turn off all inspections except for the one we're looking for (we need to fix many others reported but I was looking for something specific). To do this we create an XML config that will only run the specific inspection (RC for Reference Comparison):
findbugs-include.xml

<?xml version="1.0" encoding="UTF-8"?>
<FindBugsFilter>
<Match>
<Bug code="RC"/>
</Match>
</FindBugsFilter>
We then run findbugs with this configuration and provide the frc.suspicious property.
Running findbugs

$> findbugs -textui \
-include findbugs-include.xml \
-property "frc.suspicious=org.openscience.cdk.interfaces.IAtom" \
base/*/target/cdk-*-2.0-SNAPSHOT.jar
This produces an accurate report of all the places the references are compared. Here's a sample:

H C RC: Suspicious comparison of org.openscience.cdk.interfaces.IAtom references in org.openscience.cdk.Bond.getOther(IAtom) At Bond.java:[line 253]
H C RC: Suspicious comparison of org.openscience.cdk.interfaces.IAtom references in org.openscience.cdk.Bond.getConnectedAtom(IAtom) At Bond.java:[line 265]
H C RC: Suspicious comparison of org.openscience.cdk.interfaces.IAtom references in org.openscience.cdk.Bond.getConnectedAtoms(IAtom) At Bond.java:[line 281]
H C RC: Suspicious comparison of org.openscience.cdk.interfaces.IAtom references in org.openscience.cdk.Bond.contains(IAtom) At Bond.java:[line 300]
...

Version 1.5.15 was released: https://github.com/cdk/cdk/releases/tag/cdk-1.5.15

Version 1.5.15 was released: https://github.com/cdk/cdk/releases/tag/cdk-1.5.15

cdk

The ACS Spring disclosures of 2017 #1

At the American Chemical Society meetings drug companies disclose recent new drugs to the world. Normally, the chemical structures are already out in the open, often as part of patents. But because these patents commonly discuss many compounds, the disclosures are a big thing.

Now, these disclosure meetings are weird. You will not get InChIKeys (see doi:10.1186/s13321-015-0068-4) or something similar. No, people sit down with paper, manually redraw the structure. Like Carmen Drahl has done in the past. And Bethany Halford has taken over that role at some point. Great work from both! The Chemical & Engineering News has aggregated the tweets into this overview.

Of course, a drug structure disclosure is not complete if it does not lead to deposition in databases. The first thing is to convert the drawings into something machine readable. And thanks to the great work from John May on the Chemistry Development Kit and the OpenSMILES team, I'm happy with this being SMILES. So, we (Chris Southan and me) started a Google Spreadsheet with CCZero data:



I drew the structures in Bioclipse 2.6.2 (which has CDK 1.5.13) and copy-pasted the SMILES and InChIKey into the spreadsheet. Of course, it is essential to get the stereochemistry right. The stereochemistry of the compounds was discussed on Twitter, and we think we got it right. But we cannot be 100% sure. For that, it would have been hugely helpful if the disclosures included the InChIKeys!

As I wrote before, I see Wikidata as a central resource in a web of linked chemical data. So, using the same code I used previously to add disclosures to Wikidata, I created Wikidata items for these compounds, except for one that was already in the database (see the right image). The code also fetches PubChem compound IDs, which are also listed in this spreadsheet.

The Wikidata IDs link to the SQID interface, giving a friendly GUI, one that I actually brought up before too. That said, until people add more information, it may be a bit sparsely populated:


But others are working on this series of disclosures too, and keep an eye on this blog post, as others may follow up with further information!

CDK AtomContainer's are Slow - Lets fix that

The core class for molecule representation in CDK is the AtomContainer. The AtomContainer uses an edge-list data structure for storing the underlying connection table (see The Right Representation for the Job).

Essentially this edge-list representation is efficient in space. Atoms can be shared between and belong to multiple AtomContainers. Therefore querying connectivity (is this atom connected to this other atom) is linear time in the number of bonds.

The inefficiency of the AtomContainer can really sting. If someone was to describe Morgan's relaxation algorithm you may implement it like Code 1. The algorithm looks reasonable however it will run much slower than you expected. You may expect the runtime of this algorithm to be ~N2 but it's actually ~N3. I've annotated with XXX where the extra effort creeps in.
Code 1 - Naive Morgan-like Relaxation (AtomContainer/AtomIter)
// Step 1. Algorithm body
int[] prev = new int[mol.getAtomCount()];
int[] next = new int[mol.getAtomCount()];
for (int i = 0; i < mol.getAtomCount(); i++) {
next[i] = prev[i] = mol.getAtom(i).getAtomicNumber();
}
for (int rep = 0; rep < mol.getAtomCount(); rep++) { // 0..numAtoms
for (int j = 0; j < mol.getAtomCount(); j++) { // 0..numAtoms
IAtom atom = mol.getAtom(j);
// XXX: linear traversal! 0..numBonds
for (IBond bond : mol.getConnectedBondsList(atom)) {
IAtom nbr = bond.getConnectedAtom(atom);
// XXX: linear traversal! 0..numAtoms avg=numAtoms/2
next[j] += prev[mol.getAtomNumber(nbr)];
}
}o
System.arraycopy(next, 0, prev, 0, next.length);
}

A New Start: API Rewrite?


Ultimately to fix this problem correctly, would involve changing the core AtomContainer representation, unfortunately this would require an API change, optimally I think adding the constraint that atoms/bonds can not be in multiple molecules would be needed**. This would be a monumental change and not one I can stomach right now.

Existing Trade Off: The GraphUtil class


In 2013 I added the GraphUtil class for converting an AtomContainer to a more optimal adjacency list (int[][]) that was subsequently used to speed up many algorithms including: ring finding, canonicalisation, and substructure searching. Each time one of these algorithm is invoked with an IAtomContainer the first step is to build the adjacency list 2D array.

Code 2 - GraphUtil usage
IAtomContainer mol = ...;
int[][] adj = GraphUtil.toAdjList(mol);

// optional with lookup map to bonds
EdgeToBondMap e2b = EdgeToBondMap.withSpaceFor(mol);
int[][] adj = GraphUtil.toAdjList(mol, e2b);

Although useful the usage of GraphUtil is somewhat clunky requiring passing around not just the adjacency list but the original molecule and the EdgeToBondMap if needed.
Code 3 - GraphUtil Depth First Traversal

void visit(IAtomContainer mol, int[][] adj, EdgeToBondMap bondmap, int beg, int prev) {
mol.getAtom(beg).setFlag(CDKConstants.VISITED, true);
for (int end : adjlist[beg]) {
if (end == prev)
continue;
if (!mol.getAtom(end).getFlag(CDKConstants.VISITED))
visit(mol, adj, bondmap, end, beg);
else
bondmap.get(beg, end).setIsInRing(true); // back edge
}
}

Using the GraphUtil approach has been successful but due to the clunky-ness I've not felt comfortable exposing the option of passing these through to public APIs. It was only ever meant as an internal optimisation to be hidden from the caller. Beyond causing unintentional poor performance (Code 1) what often happens in a workflow is GraphUtil is invoked multiple times. A typical use case would be matching multiple SMARTS against one AtomContainer.

A New Public API: Atom and Bond References


I wanted something nicer to work with and came up with the idea of using object composition to extend the existing Atom and Bond APIs with methods to improve performance and connectivity checks.

Essentially the idea is to provide two classes, and AtomRef and BondRef that reference a given atom or bond in a particular AtomContainer. An AtomRef knows about the original atom it's connected bonds and the index, the BondRef knows about the original bond, it's index and the AtomRef for the connected atoms. The majority of methods (e.g. setSymbol, setImplicitHydrogenCount, setBondOrder) are passed straight through to the original atom. Some methods (such as setAtom on IBond) are blocked as being unmodifiable.

Code 4 - AtomRef and BondRef structure
class AtomRef implements IAtom {
IAtom atm;
int idx;
List<BondRef> bnds;
}

class BondRef implements IBond {
IBond bnd;
int idx;
AtomRef beg, end;
}

We can now re-write the Morgan-like relaxation (Code 1) using AtomRef and BondRef. The scaling of this algorithm is now ~N2 as you would expect.
Code 5 - Morgan-like Relaxation (AtomRef/AtomIter)
// Step 1. Initial up front conversion cost
AtomRef[] arefs = AtomRef.getAtomRefs(mol);

// Step 2. Algorithm body
int[] prev = new int[mol.getAtomCount()];
int[] next = new int[mol.getAtomCount()];
for (int i = 0; i < mol.getAtomCount(); i++) {
next[i] = prev[i] = mol.getAtom(i).getAtomicNumber();
}
for (int rep = 0; rep < mol.getAtomCount(); rep++) {
for (AtomRef aref : arefs) {
int idx = aref.getIndex();
for (BondRef bond : aref.getBonds()) {
next[idx] += prev[bond.getConnectedAtom(aref).getIndex()];
}
}
System.arraycopy(next, 0, prev, 0, next.length);
}

The depth first implementation also improves in readability and only requires two arguments.
Code 6 - AromRef Depth First (AtomRef/AtomFlags)
// Step 1. Initial up front conversion cost
void visit(AtomRef beg, BondRef prev) {
beg.setFlag(CDKConstants.VISITED, true);
for (BondRef bond : beg.getBonds()) {
if (bond == prev)
continue;
AtomRef nbr = bond.getConnectedAtom(beg);
if (!nbr.getFlag(CDKConstants.VISITED))
visit(nbr, bond);
else
bond.setIsInRing(true); // back edge
}
}


Benchmark


I like the idea of exposing the AtomRef and BondRef to public APIs. I wanted to check that the trade-off in calculating and using the AtomRef/BondRef vs the current internal GraphUtil. To test this I wrote a benchmark that implements some variants of a Depth First Search and Morgan-like algorithms. I varied the algorithm implementations and whether I used, IAtomContainer, GraphUtil, or AtomRef.

The performance was measured over ChEMBL 22 and averaged the run time performance over 1/10th (167,839 records). You can find the code on GitHub (Benchmark.java). Each algorithm computes a checksum to verify the same work is being done. Here are the raw results: depthfirst.tsv, and relaxation.tsv.


Depth First Traversal


A Depth first traversal is a linear time algorithm. I tested eight implementations that varied the graph data structure and whether I used an external visit array or atom flags to mark visited atoms. When looking just at initialisation time the AtomRef creation is about the same as GraphUtil. There was some variability between the different variants but I couldn't isolate where the different came from (maybe GC/JIT related). The runtime of the AtomRef was marginally slower than GraphUtil. Both were significantly faster (18-20x) than the AtomContainer to do the traversal. When we look at the total run-time (initialisation+traversal) we see that even for a linear algorithm, the AtomRef (and GraphUtil) were ~3x faster. Including the EdgeToBondMap adds a significant penalty.




Graph Relaxation


A more interesting test is a Morgan-like relaxation, as a more expensive algorithm (N2) it should emphasise any difference between the AtomRef and GraphUtil. The variability in this algorithm is whether we relax over atoms (AtomIter - see Code 1/5) or bonds (BondIter). We see a huge variability in AtomContainer/AtomIter implementation. This is because the algorithm is more susceptible to difference in input (molecule) size.



Clearly the AtomContainer/AtomIter is really bad (~80x slower). Excluding this results shows that as expected the AtomRef/AtomIter is slower than the GraphUtil/AtomIter equivalent (~2x slower). However because the AtomRef has a richer syntax, we can do a trick with XOR number storage to improve performance or iterate over bonds (BondIter) giving like-for-like speeds.



Conclusions


The proposed AtomRef and BondRef provide a convenience API to use the CDK in a natural way with efficient connectivity access. The conversion to an AtomRef is efficient and provides a speedup even for linear algorithms. The encapsulation facilities the passing as a public API parameter, users will be able to compute it ahead of time and pass it along to multiple algorithms.

I'm somewhat tempted to provide an equivalent AtomContainerRef allowing a drop-in replacement for methods that take the IAtomContainer interface. It is technically possible to implement writes (e.g. delete bond) efficiently in which case it would no longer be a 'Ref'. Maybe I'll omit that functionality or use a better name?

Footnotes


  • ** My colleague Daniel Lowe notes that OPSIN allows atoms to be in multiple molecules and know about their neighbours but it's a bit of a fudge. It's certainly possible with some extra book keeping but prevents some other optimisations from being applied.

Postdoc and PhD positions in Cheminformatics at Jena University, Germany

One postdoc position and three phd positions are available in my newly founded research group at Jena University, Germany.

I am currently moving from my previous position as Head of Cheminformatics and Metabolism at the European Bioinformatics Institute (EBI) to the Institute for Analytical Chemistry as Professor for Analytical Chemistry, Cheminformatics and Chemometrics at Jena University. The successful candidates will help forming the nucleus of the new research group and work in an exciting network of local and international collaborations, such as the PhenoMeNal project funded by the European Commission in their Horizon2020 framework program.

Open Positions:

  1. Postdoc: We are looking for a talented cheminformatician, bioinformatician or someone with comparable skills to work on the development cloud-based methods for computational metabolomics. The successful candidate will work closely with the H2020 e-infrastructure project PhenoMeNal, a European consortium of 14 partners. This position requires excellent skills in at least one modern, object-oriented programming language. A strong interest in metabolomics and cloud computing as well as the ability to work in a distributed team will be advantageous. The postdoc will also have the opportunity to participate in the day-to-day management of the group as well as in the organisation of seminars and practical courses for our students.
  2. PhD student, biomedical information mining: In this phd project the candidate will combine methods of text mining, image mining and cheminformatics to extract information about metabolites and natural products from the published primary literature. This includes opportunities to work with the OpenMinTed consortium, where we have been leading the biomedical use case in the last 1.5 years, as well as with the ContentMine team.
  3. PhD student, cheminformatic prediction of natural product structures: Depending on skills and interests of the successful candidate, this project can target the problem of structure prediction of natural products and metabolite from either the side of spectroscopic information which one might have about an unknown natural product or starting from the genome of a natural product producing organism. Two positions are available in this area.

All PhD positions require a strong interest in molecular informatics and current IT technologies, programming skills a modern object oriented programming language and the ability to work in geographically distributed teams.

Please send applications in PDF format by email to christoph.steinbeck@uni-jena.de. We will accept applications until the position is filled.

Background information:

The Friedrich Schiller University Jena (FSU Jena), founded in 1558, is one of the oldest universities in Europe and a member in the COIMBRA group, a network of prestigious, traditional European universities. The University of Jena has a distinguished record of innovations and resulting educational strengths in  major fields such as optics, photonics and optical technologies, innovative materials and related technologies, dynamics of complex biological systems and humans in changing social environments. It has more than 18,000 students. The university’s friendly and stimulating atmosphere and state-of-the-art facilities boost academic careers and enable excellence in learning, teaching and research. Assistance with proposing and inaugurating new research projects and with establishing public-private partnerships is considered a crucial point.

About Christoph Steinbeck

 

Alzheimer’s disease, PaDEL-Descriptor, CDK versions, and QSAR models

A new paper in PeerJ (doi:10.7717/peerj.2322) caught my eye for two reasons. First, it's nice to see a paper using the CDK in PeerJ, one of the journals of an innovative, gold Open Access publishing group. Second, that's what I call a graphical abstract (shown on the right)!

The paper describes a collection of Alzheimer-related QSAR models. It primarily uses fingerprints and the PaDeL-Descriptor software (doi:10.1002/jcc.21707) for it particularly. I just checked the (new) PaDeL-Descriptor website and it still seems to use CDK 1.4. The page has the note "Hence, there are some compatibility issues which will only be resolved when PaDEL-Descriptor updates to CDK 1.5.x, which will only happen when CDK 1.5.x becomes the new stable release." and I hope Yap Chun Wei will soon find time to make this update. I had a look at the source code, but with no NetBeans experience and no install instructions, I was unable to compile the source code. AMBIT is now up to speed with CDK 1.5, so the migration should not be too difficult.

Mind you, PaDEL is used quite a bit, so the impact of such an upgrade would be substantial. The Wiley webpage for the article mentions 184 citations, Google Scholar counts 369.

But there is another thing. The authors of the Alzheimer paper compare various fingerprints and the predictive powers of models based on them. I am really looking forward to a paper where the authors compare the same fingerprint (or set of descriptors) but with different CDK versions, particularly CDK 1.4 against 1.5. My guess is that the models based on 1.5 will be better, but I am not entirely convinced yet that the increased stability of 1.5 is actually going to make a significant impact on the QSAR performance... what do you think?


Simeon, S., Anuwongcharoen, N., Shoombuatong, W., Malik, A. A., Prachayasittikul, V., Wikberg, J. E. S., Nantasenamat, C., Aug. 2016. Probing the origins of human acetylcholinesterase inhibition via QSAR modeling and molecular docking. PeerJ 4, e2322+. 10.7717/peerj.2322

Yap, C. W., May 2011. PaDEL-descriptor: An open source software to calculate molecular descriptors and fingerprints. Journal of Computational Chemistry 32 (7), 1466-1474. 10.1002/jcc.21707

The Groovy Cheminformatics scripts are now online

My Groovy Cheminformatics with the Chemistry Development Kit book sold more than 100 times via Lulu.com now. An older release can be downloaded as CC-BY from Figshare and was "bought" 39 times. That does not really make a living, but does allow me to financially support CiteULike, for example, where you can find all the references I use in the book.

The content of the book is not unique. The book exists for convenience, it explains things around the APIs, gives tips and tricks. In the first place for myself, to help me quickly answer questions on the cdk-user mailing list. This list is a powerful source of answers, and the archive covers 14 years of user support:


One of the design goals of the book was to have many editions allowing me to keep all scripts updated. In fact, all scripts in the book are run each time I make a new release of the book, and, therefore, which each release of the CDK that I make a book release for. That also explains why a new release of the book currently takes quite a bit of time, because there are so many API updates at the moment, as you can read about in the draft CDK 3 paper.

Now, I had for a long time also the plan to make the scripts freely available. However, I never got around to making the website to go with that. I have given up on the idea of a website and now use GitHub. So, you can now, finally, find the scripts for the two active book releases on GitHub. Of course, without the explanations and context. For that you need the book.

Happy CDK hacking!