In this episode recorded 7 June 2011 at ICCS 2011, I interviewed Paolo Tosco and the Open3DQSAR project. It is a tool for high-throughput chemometric analysis of molecular interaction fields. As of 17 June 2011 this project became free software; before then it was restricted to non-Europeans because the ComFA patent for Europe did not expire until then.

I spent the last 4 hours transcribing the interview. Please let me know if you find it useful.

AD: Hello my name is Andrew Dalke. Welcome to episode two of "Molecular Coding." A few weeks ago I was at the International Conference on Chemical Structures in The Netherlands. During lunchtime I met Paolo Tosco. He was talking with others at the table about how the ComFA patent meant that his 3D-QSAR package was available at no cost to most people in the world but it couldn't be released as free software until the European patent expired on 17 June 2011.

AD: I asked him more about his project, and I was impressed by both what he's done and also with the amount of work he puts in to supporting the people who might use his software. He very kindly gave me a crash course in 3D-QSAR so I would know enough to be able to interview him about his project. The interview took place on the 7th of June, 2011.

AD: If you want to try out his software, Open3DQSAR.org is for the tool that does high-throughput chemometric analysis of molecular interaction fields, its sister site Open3DAlign.org does unsupervised molecular alignment, and sdf2xyz2sdf.sourceforge.net does atom type assignment for the Merck molecular forcefield implementation in Tinker.

music

AD: Welcome again to Molecular Coding. I'm here with Paolo Tosco who has written a 3D-QSAR program. The software is Open3DQSAR and it uses OpenBabel and it's an open source project, or almost open source. And I thought I would talk about what is going on [and] more about the software components that go into it.

PT: Yes, officially OpenQSAR will be fully open source in a few days because the infamous ComFA patent will expire on the 18th of June so I will be able to freely distribute the software under the GPLv3 license in a few days. Until now it was restricted to no European countries because the patent expired in the USA so I could distribute everywhere but in like twelve European countries. It will be truly open source.

AD: Okay, and since I didn't know much about 3D-QSAR I asked Paolo to go through and walk me through the steps, and so I think we'll just do that now. How do we use your tools and do science?

PT: So basically 3D-QSAR is a way to derive a correlation between the activities of a set of molecules and their 3D properties in terms usually of steric and electostatic potential fields. You try to spot out the zones of each compound which drive better or worse affinity for your target. This technique is especially useful when you don't know very much about the target or maybe you don't even know what the target it. So the only experimental information you can use is the one which comes from ligands.

AD: So you start off with what, a SD files that contains "these are active" and "these other.."

PT: Yes. You can import an SD file which of course contains the 3D structures and along with the 3D structures there could be a field in which occurs an affinity or pIC50 or whatever. If you don't have an SD file which encode the affinity you can import an external text file [with] the affinities provided; of course in the correct order and that they match the molecules in the SD file but this is obvious.

AD: I know the next step in this was to do the conformation generation. I was looking through the mailing list and it looked like the conformation generation wasn't even added until last fall some time.

PT: The point of the conformation generation is a bug problem, that every knows because it's not very difficult as long as you take into consideration open structures but it becomes a pain when you have to sample for instance the ring systems. It's definitely not easy to cut a ring and sample the conformation and then rebuild the ring. It's not to beat a good Monte Carlo algorithm for instance. The way I decided to implement conformation searching was using quenched molecular dynamics so that is molecular dynamics at high temperature. The high temperature should be high enough to overcome the torsional energy barriers but too high for instance to enantiomerize chiral centers or isomerize double bonds or these kinds of things that you don't usually want to happen.

AD: Most people who might want to do conformation search might use say Corina or the OpenEye tools to do this. You decided go ahead and implement your own conformation generation tool. How did you go about doing that and how good of a replacement is it.. I guess, for the purposes of doing 3D-QSAR?

PT: The point is that I wanted the whole suite did not include any kind of either closed source program or any non free program so I had to rely on either free available tools or something I was coding on my own. The nice point is since I saw there are not too many; basically there is no Merck forcefield based tool on the market which is completely free. My only choices were OpenBabel and Tinker. While OpenBabel can accomplish basic energy conformation calculations using the Merck forcefield of course Tinker offers a much larger palette of possible calculations including also implicit solvent models which proved to be very handy when you're dealing with conformations in a solution of molecules.

My choice fell on Tinker. Merck forcefield was implemented in Tinker a few years ago but a big part was missing in this implementation because basically there was no tool to do the automatic atom typing of the input coordinate files. The first step was to code something which was able to do the atom typing. Since I didn't want to reinvent the wheel from scratch of course I relied on the code which was already existing and already GPL'ed and that's the OpenBabel code. I fixed a few bugs which were present in the Merck forcefield assignment code and I added the part which converts the 99 forcefield atom types in the 214 Tinker MMFF-like atom types.

AD: I've never actually implemented MMFF but that defines the atom types. We have SD file comes in, you fixed the OpenBabel code so it does the correct MMFF atom type assignment, and then you had to translate the MMFF type assignment into the Tinker types.

PT: Exactly. The Tinker atom types enrich the plain Merck forcefield atom type because they add atom types which kind of encode part of the charge inside the atom type itself. It's a richer description of the original Merck forcefield atom type. I talked with Professor Ponder about this fact and he was interested as well in a tool which was able to produce an XYZ file containing Merk forcefield atom types according to Tinker format. That's why I also released as a side dish a stand-alone tool which is able to convert an SD file either single or multi-molecule in either one or a set of XYZ files ready for accomplishing calculations inside of Tinker. You can find this other tool which has a cumbersome name. I'm still asking myself why I used such a ... laughs

AD: That's the sdf2xyz2sdf.

PT: Yes. You will win a prize for being able to spell it without ... laughs

AD: It is a mouthful.

PT: Yes. And this is just a stand-alone tool to accomplish atom typing.

AD: So that doesn't even use OpenBabel at all?

PT: It does. To make things easier I did not really directly link to the libraries but actually what the tool does calling the real binary so you bypass many of the problems you can face; I have an old library on my system, so the tool includes the binaries for any operating system so you can use it quite comfortably and you don't have to care about the fact that you may have an updated version or not on your own system.

AD: I'm also curious about how it was to integrate with Tinker because the one thing I know is it's written in Fortran and the other thing is I also know most of the times you use Tinker you have to get it from ... It's not freely available. I mean, it's available at not code but it's not a GPL or BSD license.

PT: That was also my concern because of course you can alway use an external tool but then all users will have the pain of having to get the tool on their own and to get the license and to get the binaries and to compile the binaries for their own architecture and so on and so on. This is something which in my opinion is already quite discouraging about using your software. So of course if you are able to put together a full-featured package which has all the pieces it needs to work of course I think it's better for anyone. That's why I asked Professor Ponder whether I could actually include in my distribution their binaries even if they are not strictly open source, and actually he was positive about that since given the open source nature of my project and the non-profit nature of the whole thing he was very eager to allow me to do this because he was really very collaborative and very interesting in the whole thing. I think this was a very good point in the whole story.

AD: From what you were saying earlier, you said that he also knew there's now an atom typer for the Merck forcefield.

PT: Actually it was him encouraging me to take the atom typing part of Open3DQSAR and encode it as a stand-alone tool which has the awkward name of sdf2xyz and so on. It was him really promoting the fact that I released this tools as a stand-alone tool.

AD: Did you ever read The Hobbit?

PT: Yes, I did.

AD: "There and Back Again"

PT: Well, I read it a long time ago so I don't have very clear memories about that.

AD: The subtitle that Bilbo Baggins was "There and Back Again", so "sdf2xyz2sdf" is "There and Back Again."

PT: That's true.

AD: So you've got the whole work of talking with Tinker during the conformations... doing the quenched MD to get all the different conformations and get them back. You store them in some result you use then for the next stage. What's the intermediate conformer storage format? And the reason I'm asking this is because I know that you can save as an SD file but then you're just saving a bunch of stuff over and over again. Back in my MD days we just saved the coordinates so I was wondering if, oh! people should start doing this for their conformers.

PT: The point of why the SDF format can actually store multiple conformation the XYZ format can't, because each XYZ file is for a single molecule. The way I tried to go around it in the simplest way for all the users is that you read the SDF file from standard input. In the SDF files you have a molecule name. In the case that the molecule name is actually defined you will obtain an XYZ file with proper name for each molecule which was in the original SD file and you will also obtain with the XYZ file a key file with the keywords which are needed to accomplish the computations. Of course as a default you get the basic set of keys then you can add for instance keys for implicit solvent models and whatever you want.

PT: Actually there was a limitation in the original Tinker implementation of the Merck forcefield since it didn't work with charged heteraromatic rings and I could overcome this limitation since you could apply an external set of Merck forcefield charges and apparently OpenBabel can perfectly assign Merck forcefield charges while the original implementation of the Merck forcefield inside Tinker couldn't, so actually it looks like mixing the Open3DQSAR stuff with the Tinker stuff did well to both of the software so this this is a good point.

AD: Okay, so you generated all of the conformers ...

PT: Once you have the XYZ files you can get back to SDF file just using the original SD as a template and just the coordinates are updated in the original SD. That's quite straight-forward.

AD: So the next step you do for your QSAR model is first you you generate pharmocophores. I saw that Silicos was saying that you use their "fay-ro?", "far-o?", "Pharao?"

PT: "Pharao", I would say "Pharao", I think. I don't know if the correct pronounciation.

AD: Using the Silicos Pharao tool; do you do most of the pharmacophore filtering with their tool?

PT: At the beginning I thought I could use Pharao not that much for the pharmacorphore part but essentially it's an alignment tool. After testing Pharao with several data sets, especially with the data sets from the Southerland benchmark suite I could see that some of the data set were actually very fit to be aligned with a pharmacorphore based tool but some of the others, especially those which have an ill-defined pharmacorphore, they are not fit for that because all of the compounds making part of the data set have a different collection of pharmacorphore points so they are not really amenable to obtain a common template out of that pharmacophores. That's why I decided to enhance the alignment part with some tools coded on my own which are not pharmacophore based but atom based. That is, first you find the pairs of atoms which are most amenable to be matched and then you do the extra matching. That seems to enlarge the pool of possible molecules you can align by a large extent.

AD: Alright, I'm lost. I thought this was part of the pharmacophore alignment you were doing as part of the precursor step for doing the grid alignment but I don't quite understand the 3D-QSAR well enough.

PT: Well, the first tool which was born was the 3D-QSAR part. That implies that you already have some way to make the alignment.

AD: So we haven't gotten to the point where you're doing the 3D-QSAR. We're still getting to the point where you're making the grids.

PT: Yeah, so the reason I started from the 3D-QSAR tool and not from the alignment tool which would be the most logical way to do was that actually I already had the alignments made and I needed a tool to make the computation on a high-throughput basis. That is, without having to input each and every model from a graphical user interface as I used to do when I was using tools like GULP to make the chemometric analysis.

PT: But then of course the tool was very fun for me because I had other ways to accomplish the alignment but after the 20th email where people were asking "ahh, that's a nice tool but how the hell am I going to align the compounds to make 3D-QSAR" I thought the time was ripe to make an alignment tool.

AD: You started off with the [3D]-QSAR tool and from that you needed the alignment tool and for that you needed the conformer generation tool.

PT: Exactly.

AD: That's a lot of work.

PT: Yep, yeah, it does.

AD: So the pharmacophore alignment, it's the precursor to your 3D-QSAR code. For that are you using Pharao?

PT: Yep.

AD: And of course what's interesting to me is that Silicos came out last year saying "We're going to support OpenBabel. We'll release our stuff as open source. We'll make this available to the world." So here's an example where you said "the source code's available, I can use it for my tool" and it probably saved you a lot of time too.

PT: Yep, and actually I think my implementation was actually one of the first using Pharao because I met Gert Thijs at a meeting in Goslar last year and actually they were very happy about seeing that someone was using their tool after just after I think a few weeks after they released it as open source. They also put the PDF of my poster on their web site as an example application using it so I was quite pleased with that.

AD: How did you find out about their tool? Was it because they announced it to the OpenBabel list?

PT: I saw the paper in Journal of Molecular Graphics and Modelling and the way I found it was because I was scraping the web to find a free, open source alignment tool because I wanted to try actually a kind of brute-force approach to 3D-QSAR model building and before spending a lot of time in coding my own alignment tool I wanted to have a hint whether this could make sense or not. I needed a tool to be able to test the whole procedure in a short time without losing a lot of time in something that probably didn't have a chance to work. As long as I saw that actually there was a chance that this was going to work then I coded my own alignment tool which was, not better than Pharao intrinsically but probably more fit for 3D-QSAR because actually the alignment that you need for 3D-QSAR are particularly consistent because every small inconsistency of the alignment of course increases the background noise in your models and hides the signal. That's why you need a quite special "pure posed" alignment tool, I would call, for 3D-QSAR.

AD: You use the pharmacophore alignment to do the 6-degree of freedom alignment of your molecule to the reference, then you take that, transform your original conformation you generated, and then that's were you finally get to the point where you do the 3D-QSAR.

PT: Yes.

AD: Okay, so how does that go? That's embedded inside of a grid, you're working with grid space, and doing more work with forcefields there. Is that inside of OpenBabel?

PT: Well that part of the calculation is extremely trivial because what you compute is just the non-bonded part of the forcefield because you're just computing the van der Waals, the Lennard-Jones potential and the Coulombic part of the potential. That's really an extremely easy part of the job to do and so you don't need any fancy molecular mechanics tool. You can code it in 10 minutes I would say.

AD: And you're only talking about a few dozen atoms so the n-squared easy algorithm is not a problem.

PT: No, definitely.

AD: Okay, so you generate the grid, you assign energy points, you then do an overlay of the grids. This is the part where I'm not sure what goes on next. You subtract the energy values, or the difference of the energy values and then you ...?

PT: No, you take the energy value grid for each and every molecule and then you build a PLS model of all the grid values so you have thousands, tens of thousands of energy values and the PLS is very [something] in picking out the ones which have a variance across the data set and just throw away all the points which are basically either 0 or they have the cutoff value because there is a long [something] steric clash or because your point is inside the molecule or many molecules of the data set are quite similar in certain parts and dissimilar only in small moieties. The PLS is good at extracting only the real different information across your data set.

AD: Now this is the point where you hit the patent problem.

PT: Yep. I admit that I was very naive, but I am alway naive in my life, because I was very much concern about this point because I thought "ahh, but I'm a good guy", I'm not doing that for - yeah, you're laughing - I was really convinced about that because I was doing it on a non-profit basis and this was an academic tool and of course I'm sure you could not even think of competing with a full-featured 3D modelling suite like Sybyl. I mean, that's a tool, a complete tool, with a graphical user interface which allows you to do every step from conformation search to alignments to 3D-QSAR model building. My tool at the moment was just doing 3D-QSAR model building so I don't think it was a competitor.

PT: I was very confident when I wrote to Tripos that they would just tell me "ahh, yeah, just go along; it's not a problem." Especially because the patent had already expired in the USA and it was going to expire in Europe in one year anyway so I didn't think it would make a difference. But I received, I admit a very polite, but it was a "no", not a "yes." That forced me to use a very unpleasant geographic redistribution policy and also forced me not to be able to distribute it on a real open source basis. Of course the source was available, has always been available, but you could not modify or redistribute because you don't have a license to do that. This is going to end starting from June 18th, 2011, that is in 10 days it will be over because the patent will be expired everywhere in the world. I will finally be able to release it on a GPLv3 without issues.

AD: You're releasing this under GPLv3. Why v3 instead of v2? I know OpenBabel has the v2 license for historical reasons.

PT: I admit that I'm not a big expert of licenses. I read both the GPLv2 and v3 and I admit that I could not catch the fine differences between them and I thought, as any software guy is going, that I would take the most recent one. I admit that this is again probably very naive but you know I don't care very much about these things.

AD: I've actually read and listened to various things about the GPLv3 and it sounds like a better license all around than the v2 is. There's some pretty good advantages to it that were considered in v3 with lots of input that in the early 1990s were not fully thought out when v2 came out.

AD: So one of the things that you talk about on your web page is that this is an automated tool. How do you automate it? Is there a scripting language or commands or how?

PT: I think that the main thrust behind the whole story is making a suite which was suited for high-throughput modelling. That is, completely unsupervised to the highest degree and scriptable. I thought quite a long time of it was convenient to make a pure command-line tool, but since the software has many many options it appeared to me making a command-line tool would just mean to add a plethora of options that no one would remember and in the end it would not have been very [unknown] as a command-line tool. I decided rather to go for an input script with a very simple syntax that was keyword equal to parameter; in a CHARMM-like fashion, I think could be a good example, or AMBER-like fashion. Small input scripts, text-based which are very amenable to building in an automated way with awk or sed or whatever stream editing tools. I think that was the best way to go around it.

AD: Because it's for high-throughput, I think I saw something there about how it's parallelized.

PT: It's parallelized using just pthreads. Of course pthreads on the Unix platforms and native Windows threads on Windows. Actually here I mentioned that the software is of course available as source code but it's also available as precompiled binary distributions for Linux, Windows, Mac OSX and Solaris and FreeBSD so you don't have to care about ...

AD: I was looking at your gallery pages and I was surprised to see how many operating systems that you actually tested it on.

PT: Actually that takes a long time. Each time you have to make just a small bug fix you have to rebuild all the versions. That's very time consuming, but I think it's quite essential, especially in the first stages of your software life that it's immediately available to users without their having to build it, especially when you have a number of different appendicies and you easily get lost in getting all of them, especially if you don't have big experience in building software from source.

AD: I released a package last week, a week and a half ago now, and in the first few comments I got back were "are there prebuilt packages for Windows; I'm not used to compiling." So I agree.

AD: In addition to all these operating systems I also saw that you had support for visualizating some of the results in, not only PyMol, but MOE, Sybyl, and Maestro. You do a lot of integration work there and that takes a good amount of time.

PT: Of course as long as you are using coordinate files then any software suite can read MOL or SDF file so that's not a problem. The problem was about importing the binary grids, that is, the isocontours that visually show where you should add, for instance, steric [something], negative electrostatic potential in your structures. Those are all proprietary formats you have to output all the different formats for all the different software platforms that you want to target. For PyMol the format that PyMol can read is the same as the InsightII format. Talking about MOE, talking about Sybyl, talking about Schrodinger Maestro then all of them have different formats. I had to support all of them and test all of them. But again I think that it's essential that people who are used to certain tools don't have to switch just because your software cannot target the platform that you are used to. I think that one should always have users in mind when you release software rather than just your own convenience. Testing different platforms and different formats is also a very good way of debugging nasty things that otherwise you can overlook. Testing a couple of compilers many times shows you some weaknesses of your code that otherwise you can overlook, so it was not lost time I think.

AD: So how do you do the testing. Do you have, for instance, automated tests, or do you mostly run manual tests?

PT: I made basically a script which touches more or less all the features of the software in order to be able to ... Still I admit that sometimes it happens that adding some module to the software I could experience regressions in other parts that I forgot to test each and every time because you would never imagine that touching a part that you thought was miles away from that one would screw up another one. That comes with experience. I don't have big developing experience. I'm building it slowly.

AD: How did you get into this? What you've told me now says that you have some development background but you're doing chemistry. How did you learn to do the software development?

PT: I must say that I was born as an organic synthesis chemist. My career had gone from 1998 to 2006 as an organic chemist and the only thing I was using computers was for writing papers and sending emails. I had no kind of programming experience and I barely knew what a computer was. I didn't know about a world called Linux and all of this. I really started from scratch. I know that if some professional programmer will look at my code, probably especially the oldest part of my code, well, he would have some comments about it, but it works.

AD: When programmers look at other programmers they also have comments, so that's not saying much.

PT: Still, it was a big effort to, for instance, to parallelize things efficiently because the first release I made of my software was really slow. Before making the public release I was able just making better parallelization to improve like 10-per the original performance. I really understood how much space you have for improvements and optimization of the code. I think right now I made benchmarks against GULP, which is the only chemometrics tool that I had access to and the performance on one-to-on core is the same, but the nice point is that GULP is not parallel while my software is, so of course if you use an 8 core machine then you get 8-per performance.

AD: You actually get 8-fold speedup?

PT: Yep, because you don't have very much IO so it really scales linearly. That's a very good point. I'm really planning to do GPU porting of some of the code because it's extremely amneable to high degree of parallelization, so I expect to get very good speedups when I will do that, but I really didn't have the time.

AD: Well thank you very much for your time. It's been enjoyable talking with you and learning more about how the 3D-QSAR program you've worked works and to hear more about the background of just what goes into writing such a program.

PT: Thank you. I hope that I have somehow raised curiosity in the audience and that someone can actually try my tools and especially give me feedback about bugs, new features, and whatever you can imagine of.

AD: So the name of the program and web site is?

PT: The name of the web site is http://open3dqsar.org/ and http://open3dalign.org/ for the 3D-QSAR and alignment programs respectively so it's quite easy to remember.

music

AD: Thank you for listening to Molecular Coding. This podcast and transcript are distributed under the
Creative Commons Attribution-Share-Alike 3.0 Unported license. The theme music was composed and performed by Andreas Steffen.