![](\"http://datacite.org/wp-content/uploads/2021/11/direct_member.jpg\"/)
As a membership organization, members are at the core of everything we do. Understanding our members and ensuring good 2-way communication is therefore extremely important. Over the last year, we have been working on clarifying and aligning DataCite's membership model. We believe this will enable us to work with all organizations in the same way and give us more insight into the organizations using our services.
\nTypes of Membership
\nPreviously, we distinguished between Members, Providers, Clients, and Consortia. This often led to confusion among organizations about whether or not they are a DataCite Member. Therefore, going forward, we distinguish three clear membership categories: Member-only, Direct Member, and Consortium Member. The Member-only category remains unchanged and consists of Members that support DataCite's data sharing mission, want to collaborate with DataCite and/or be part of DataCite's governance. These Members do not register DOIs.
\nWhen it comes to Members that do consume DOI services, we distinguish two Member types: the Direct Member and the Consortium Member.
\nDirect Member
\nThis type of Member is a single organization that joins DataCite as a Member in its own right and consumes DOI services. This single organization may have one or more repositories under their umbrella, but those repositories are under the same administrative structure as the organization. A Direct Member also takes on membership responsibilities such as participating in DataCite's governance.
\n
Consortium Member
\nA consortium is a group of like-minded organizations that have come together to collectively participate in DataCite's community and governance activities and to use DataCite's DOI services. A consortium is composed of two or more non-profit organizations that are under different administrative structures. Consortia are generally located in a single country or subject-based. Organizations within a consortium can work with one or more repositories that are under the same administrative structure as that organization.
\nThe entire consortium is a Member collectively, meaning that the consortium has a single representative in DataCite's governance structure, regardless of the number of organizations that make up the consortium. One organization will take on the role of Consortium Lead, with the others participating as Consortium Organizations.
\n![](\"http://datacite.org/wp-content/uploads/2021/11/consortium_member.jpg\"/)
You can find more information about the different types of membership on our membership page.
\nRepositories
\nRepositories play a key role in DataCite services and are therefore a key component of the DataCite membership model. We define a repository as a service operated by research organizations, where research materials are stored, managed, and made accessible. A repository is a single unit and DataCite links the repository to information in re3data, where additional repository metadata are available.
\nHowever, we realize that not all content is hosted in repositories and therefore DataCite also includes periodicals as a designation. The term periodicals include journals, proceedings, books, blogs, and working paper series.
\nThe repository information provided will be used and displayed throughout DataCite services including search, data usage statistics, Member DOI statistics, and third-party indexing services.
\nDOI Fabrica
\nGoing forward, we'll be making some changes to DOI Fabrica to align with the membership categories described above.
\nRepositories
\nWe'll be replacing the old \"Clients\" terminology with \"Repositories\" and making a stronger connection to the repository concept by adding metadata elements from re3data to the repository settings information displayed in Fabrica. Periodicals, as their own sub-category of research material service, will be designated by a periodicals badge for easy identification.
\nConsortia
\nWe'll implement the Consortium Member structure so that what Consortium Organizations can see and what they have permissions to administer in Fabrica aligns with their unique status according to the model. Implementing this properly will mean a brief transfer period (on the order of hours) for the DOIs belonging to these organizations. We'll be reaching out to our current Consortium Members and their constituent organizations over the next few weeks to confirm their repository details and schedule their transfer periods. For existing Members that are planning on forming Consortia, we will work with you to ensure that your current Clients are successfully migrated to Consortium Organizations with Repositories in Fabrica.
\nTimeline
\nWe're planning on rolling the consortium and repository changes into the test system during the last week of August, and we plan to go live in production by mid-September. In the interim, there will be a brief period when the test system will reflect the new changes while the production system will not, so what you see in the two systems will not match. During this period we'd welcome your feedback in addition to our more formal testing.
\nWhat will I see?
\nMember Information
\nWith all these changes happening, it's critical for us that we have the correct contact information for each organization. This information is used on our Member page, to administer Fabrica, to communicate with you, and for billing purposes. Therefore, we want to encourage all our Members, both our long-time Members and those going through changes, to update your contact information in your Member settings. We've recently added the ability to add multiple contacts per organization, based on relevant themes of interest. Please see DataCite's privacy policy to better understand how we use personal information. Expect reminders from us because at the end of the year we'll be replacing our current mailing lists with the information obtained through this form!
\nWe realize this is a lot of information and we want to support you throughout these changes in any way we can. We created an FAQ section on our support site and are happy to talk to you at any time to find the best model for your organization. Please don't hesitate to contact support@datacite.org at any point during this process.
\nThe post The DataCite Membership model: consortia, repositories, and more appeared first on DataCite.
","doi":"https://doi.org/10.5438/gk09-ba24","guid":"http://datacite.org/the-datacite-membership-model-consortia-repositories-and-more/","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1566777600,"rid":"405s4-00m63","summary":"As a membership organization, members are at the core of everything we do. Understanding our members and ensuring good 2-way communication is therefore extremely important. Over the last year, we [\u2026] The post The DataCite Membership model: consortia, repositories, and more appeared first on DataCite.","title":"The DataCite Membership model: consortia, repositories, and more","updated_at":1788766955,"url":"https://datacite.org/blog/the-datacite-membership-model-consortia-repositories-and-more/","version":"v1"}},{"document":{"authors":[{"contributor_roles":[],"family":"Eden","given":"Terence","url":"https://orcid.org/0000-0002-9265-9069"}],"blog":{"authors":null,"community_id":"61ce553a-bafd-4aba-a952-d3bab5e85bcc","created":1788652800,"current_feed_url":null,"description":"Regular nonsense about tech and its effects \ud83d\ude43","doi":"https://doi.org/10.59350/shkspr","favicon":"https://rogue-scholar.org/api/communities/61ce553a-bafd-4aba-a952-d3bab5e85bcc/logo","feed_format":"application/atom+xml","feed_url":"https://shkspr.mobi/blog/feed/atom/","filter":null,"generator":"WordPress","home_page_url":"https://shkspr.mobi/blog","issn":"2753-1570","language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"shkspr","status":"active","subfield":"1712","title":"Terence Eden's Blog","updated":1788694460,"use_api":false},"blog_name":"Terence Eden's Blog","blog_slug":"shkspr","content_html":"I imagine everyone here has received an unsolicited message telling them that their tax is overdue and that they urgently need to visit Genuine-Tax-Payment-Website.fart or that a parcel is delayed at customs and you can pay a small sum for its release at Almost-The-Right-Acronym.ak
\nYou know it is a scam. Most people just mark as spam and move on with their day. But a significant number of people don't. They hastily visit the site, tap in their credit card details, give it their mother's maiden name, confirm address, upload a nude selfie, and only then realise that they've been had.
\nThe Internet works at pretty close to the speed of light. You can register a .uk domain and a minute later it's accessible from the other side of the planet. Brilliant for users who want to quickly launch a website. Also brilliant for abusers who want to launch a spam campaign.
\nBy the time enough people have reported the scammers' domain as suspicious, it is too late. In the time it takes for a registrar to disable the domain, or for its name to make its way to the Safe Browsing List, a million messages have already been sent and enough people have handed over their details.
\nWe're told that \"the purpose of a system is what it does\". At the moment, the Domain Name System's purpose seems to be a vector for criminals to run scams on people at a terrifyingly high rate.
\nBIG!
\nThere's a great blog post by Andrew Campling which reports on this startling claim:
\n\nThe study found that at least 10% of all new gTLD domain names registered during the year had subsequently appeared on security blocklists by the time of analysis. It estimated that, taking account of subsequent blocklisting and associated domains not themselves blocklisted, the share of names registered by malicious actors may be closer to 20%.
\nDNS Abuse and Criminal Infrastructure: Beyond Definitions and Blocklists (emphasis added)
That links to a presentation by Interisle which contains some rather shocking statistics (albeit with disputed methodology). It looks at generic Top Level Domains (gTLD) - those are things like .com and .fun rather than country code TLDs (ccTLD) like .uk and .de.
\nIt says 85 million new registrations of gTLDs were made in 2025. Of those 8.5 million were added to blocklists by May 2025. It reckons that a 10% abuse rate is the likely floor for these numbers and it's probably closer to 20%. One in five newly registered domains with a gTLD are scams. That's a bloody crisis.
\n13 TLDs had more than 50% of their registrations blocklisted.
\n![\"Table](\"https://shkspr.mobi/blog/wp-content/uploads/2026/08/gTLDs.webp\")
\nI can understand why .bid and .loan are popular with scammers. But why .mobi?! What did I ever do to you, eh?
\nWho are the scammers registering these through?
\n![\"List](\"https://shkspr.mobi/blog/wp-content/uploads/2026/08/registrars.webp\")
\nAh, our old friends at NameCheap. See Why do scammers love NameCheap?
\nIf those five registrars had more effective policies, it might significantly dent the scammers' ability to ply their devious wares. Or they might just move on to other registrars.
\nAs the report points out:
\n\nsuspension rates for blocklisted domains were 7.4% to 16.3%.
The full report is on the Interisle website.
\nI don't know.
\nIn the first instance, it might make sense for registrars to do strong Know Your Customer (KYC) checks on anyone buying a domain. But that stops anyone who wants to anonymously register I-Hate-Nintendo.whatever without risking the wrath of Intellectual Property lawyers.
Also, criminals have access to stolen money and stolen cards. They can convince a hapless mule to register a domain on the criminals' behalf.
\nRegistrars could ask for an escrow payment. Pay \u20ac9 for the domain name put \u20ac900 in escrow. If your domain appears on a blocklist within the year, you forfeit the money. Criminals with stolen funds are unlikely to care but it would probably put off lots of people from getting a new domain.
\nThere are various banned words and phrases depending on the TLD. For example, South Sudan has a list of political words which they don't want associated with their .ss ccTLD.
\nBut if one gTLD bans a word, a different one might not. A scammer doesn't care if the gTLD is .arse or .elbow - they just want the start of the domain to look legitimate.
\nSome registrars have strings that they don't allow. In fairness to NameCheap, when I tried to register dwp-payments-gov-uk.pizza it told me that domain was banned. It wouldn't let me get any gTLD with that name.
But all it takes is one registrar to be slightly lax and the scammers get through. Increasing the complexity of the rules is also a hell of a burden on smaller registrars.
\nBesides, it's pretty easy to get a generic enough looking domain and stick the confusing bit on a subdomain. Here are a clutch mentioned in the report:
\nhttps://gov.uk-dwpaph.bond/uk/https://gov.uk-dwpcjh.bond/uk/https://gov.uk-dwpclc.bond/ukhttps://gov.uk-dwpclw.bond/ukhttps://gov.uk-dwpclj.bond/uk/Perhaps there ought to be a delay before a new domain goes live to allow people to object to it? That would give governments, banks, delivery companies, and a dozen more \"important\" organisations a right to veto any \"dodgy\" looking domain.
\nBut suppose someone wants to register gov-uk-stole-my-horse.horse to protest the government's cruel policy of stealing horses - is that a legitimate use of a domain? What if the Darwin Pensioner Divas - a group of elderly singers - want to take payments for their new album of goth/punk covers, can the DPD delivery company veto dpd-payments.music?
Do we want a domain name system where powerful companies control exactly which domains we can register? If I have an idea for a domain on a Friday night do I have to wait until Monday before it can be launched? Are those companies realistically able to parse millions of domains per year and have a low false-positive rate?
\nAll of these things are possible - but all of them come with an impact on legitimate users. To be clear, I don't know what the right answer is.
\nLots! It has been a few years since I've been to an ICANN meeting, but even back then the topic of abuse was high on the agenda. They appear to be looking at ways to coordinate abuse reports between various entities, along with some other policies which should hopefully work.
\nThere are two salient points from one of the discussions held at the recent meeting
\n\nIf anybody thinks that in our current age of AI and as we move into different kinds of computing, DNS abuse is going to numerically stay steady and we will have a downward effect on that baseline 2027 number. I'm not sure that that's an accurate assumption. I think it's going to be the other thing, which is [\u2026] it's going to be easier to abuse the DNS.
And
\n\nAbusers are going to abuse because it's just too lucrative, because no matter what we do, they will find the way to make profit off of that, and will try to circumvent everything that we do. That is not a reason not to do it, though.
Quite!
\nAs I said, I don't know the answer to this. What I do know is, much like Android's app ecosystem being a haven for scammers, DNS is facing a crisis. When trust in a system goes, only chaos follows.
\nI don't want to live in a world where I have to show my passport and pay thousands of pounds to register a domain which is only available after being vetted by private interests. But I also don't want to live in a world where scammers have effectively no deterrent from abusing millions of people.
\nThe purpose of a system is what it does. I hope DNS's purpose can become less dangerous while still remaining open.
\n![\"\"](\"https://shkspr.mobi/blog/wp-content/themes/edent-wordpress-theme/info/okgo.php?ID=74588&HTTP_REFERER=Atom\")
","doi":"https://doi.org/10.59350/395ha-fss97","guid":"https://shkspr.mobi/blog/?p=74588","image":"https://shkspr.mobi/blog/wp-content/uploads/2026/08/gTLDs.webp","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1788652800,"rid":"8hyr0-9p774","summary":"I imagine everyone here has received an unsolicited message telling them that their tax is overdue and that they urgently need to visit Genuine-Tax-Payment-Website.fart or that a parcel is delayed at customs and you can pay a small sum for its release at Almost-The-Right-Acronym.ak You know it is a scam. Most people just mark as spam and move on with their day.","tags":["/etc/","ICANN","Internet","Scam","Spam"],"title":"The purpose of DNS is to spread scams","updated_at":1788766212,"url":"https://shkspr.mobi/blog/2026/09/the-purpose-of-dns-is-to-spread-scams/","version":"v1"}},{"document":{"authors":[{"contributor_roles":[],"family":"Swain","given":"Jonathan","url":"https://orcid.org/0000-0003-4457-1481"}],"blog":{"authors":[{"name":"Jon Swain","url":"https://orcid.org/0000-0003-4457-1481"}],"community_id":"13f55986-f209-443c-ae0d-2f9f3f521e5a","created":1788652800,"current_feed_url":null,"description":"I am a cheminformatician and data scientist, originally from the UK, but often found in Aotearoa (New Zealand). I'm interested in using data science and machine learning to solve problems in drug discovery.","doi":"https://doi.org/10.59350/jonswain","favicon":"https://rogue-scholar.org/api/communities/13f55986-f209-443c-ae0d-2f9f3f521e5a/logo","feed_format":"application/atom+xml","feed_url":"https://jonswain.github.io/feed.xml","filter":null,"generator":"Jekyll","home_page_url":"https://jonswain.github.io/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"jonswain","status":"active","subfield":"3002","title":"Jon Swain","updated":1778680800,"use_api":null},"blog_name":"Jon Swain","blog_slug":"jonswain","content_html":"![\"A](\"https://jonswain.github.io/images/career_change/chemistry_to_data_science.png\"/)
In 2020 I was working as a synthetic chemist in a laboratory in Aotearoa (New Zealand). I was getting disillusioned with the work (and academia in general) and started looking for a new career direction and new challenges. After chatting with a few friends who were working as data scientists I realised that my current interests in scientific computing, statistics, and drug discovery would transfer nicely to a job as a data scientist or cheminformatician.
\nIt felt like an emerging field within biomedical research, with huge potential. Alphafold 2 had just been released and was causing a lot of excitement, but the LLM hype was still yet to fully start (ChatGPT was still a few years away). It seemed like a perfect time to make the jump to a new, exciting field that would grow in the future.
\nWhen I started searching for information about making the switch I somehow felt both completely lost and unguided, and overwhelmed by the resources out there. Part of the thinking behind this post was: What do I wish I knew before setting out on this journey?
\nI've been asked a couple of times about making the change, mostly by old colleagues looking to move from chemistry research to data science, but I hope this advice is useful whatever your background is. I'm always happy to answer questions but thought it would be useful to have something to point people to that contained my main ideas and could be used as a starting point for further questions.
\nThere's a huge amount resources out there, and no right or wrong way to make the change. This just lists my experience and what worked for me, this will be different for everyone. One of my favourite things about data science is the low bar for entry, there's so much open-source software and data, YouTube videos, and Kaggle competitions that you only really need a basic computer to get started. On top of this I realise I was also lucky to have friends and colleagues that I could ask for advice, and a supportive family who were willing to house me whilst I re-trained. Through 2020 and 2021 we were occasionally confined at home due to COVID lockdowns (though thankfully in NZ these were fairly limited compared to the rest of the world), and with not much else going on this was an idea time for some re-training.
\nAs I asked around my network for advice on making this change, the first suggestion was always the same, learn to code in Python. So much data science work is done using Python libraries (pandas, numpy, scikit-learn, pytorch, and rdkit for chemistry). I now spend a significant part of each day writing and reading Python code and documentation.
\nI did what I often do (probably a bit too much) and scrolled reddit. I kept seeing two online Python courses recommended. Automate the Boring Stuff with Python Programming by Al Sweigart, and 100 Days of Code: The Complete Python Pro Bootcamp by Dr. Angela Yu. Both were available on Udemy, which seems to regularly have special offers and I completed both courses for around $10 (US) each. There are definitely similar courses available for free so it's a personal preference what to go for. I chose the paid courses for two reasons:
\nThe first course I completed was Automate the Boring Stuff with Python Programming by Al Sweigart. It's a fairly short course but very well taught, designed for anyone who works regularly with a computer and wants to automate repetitive tasks, so they never have to do them again. It doesn't go too deep into theory or style and gets stuck in with practical applications quickly. Within a couple of days, I was building small programs and could start to see how I could use it for to automate parts of my job. Before starting I had no idea if I would even enjoy programming in Python, and this was an ideal introduction (turns out I did enjoy it) and gave me enough understanding to look at other courses that might be useful. It seemed perfect to make my current job less boring, but not enough for the new job I wanted, so I started looking for something more in-depth.
\nThe next course I completed was 100 Days of Code: The Complete Python Pro Bootcamp by Dr. Angela Yu. This is a much more in-depth course that covers a huge range of applications for Python. Every day has a couple of hours of lectures that cover theory or usage, and a few problem sets to apply what you've learnt. Every 10 days there's a small project that uses everything you've learnt, and at the end there's about 20 days of projects to construct a personal Python portfolio. I found some sections a little confusing (the first time OOP is introduced), and there's a large section borrowed from a course on web development that I felt dragged on a bit (it's useful to understand HTML, CSS, and JS for things like web-scraping, but I didn't need quite that much detail on building websites). There's a significant section towards the end on data science that was particularly useful. Overall, this course probably made the biggest difference in my re-training.
\n![\"I](\"https://jonswain.github.io/images/career_change/cant_get_job.png\"/)
At this stage I started thinking about jobs, but everyone I spoke to gave the same response: \"What you're doing is the right idea, but we're just looking for someone with more experience.\" This seems to be a classic problem when re-training, how do you get your first experience when everyone wants someone with experience? There are a few ways to go about this, but the main two seem to be either an internship or a bootcamp. Both have pros and cons, and I decided to go down the bootcamp route. To me the main advantages were the ease of organisation (I was relocating to the UK at the time and having something waiting for me was nice), and the shorter duration (most internships seemed to be multiple months, whereas the bootcamp I did was 5 weeks).
\nAt the recommendation of a friend, I applied to Science to Data Science (S2DS), a 5-week, intensive, project-based bootcamp for research scientists with a PhD or MSc looking to move into data science. The application involved a short technical project to test your exploratory data analysis skills and a behavioural interview (that seemed mostly to filter out psychopaths).
\nFor the 5 weeks you're part of a small group that gets paired with a company, working on a data science project. There's a couple of talks in the first week about best practices, but you're working on your project straight away. My team was paired with Deutsche Welle (DW), a German broadcaster. They were interested in the gender breakdown of people mentioned and quoted in their articles, so the project involved Natural Language Processing (NLP) and Named Entity Recognition (NER), as well as using techniques like web-scraping for data gathering and machine learning classifiers for predicting gender from the name.
\nWorking in a team made the experience a lot of fun. We all came from a range of backgrounds and had different skillsets, but all started the bootcamp from similar positions so were learning together. It also gave a first experience using git and github collaboratively. In addition to the team, we had a mentor from S2DS who helped with problems, and contacts with DW who helped with project direction.
\nWhen I applied S2DS was \u00a3800 for the entire course, not cheap, but much cheaper than many bootcamps I saw advertised online. I made the calculation that getting a job and starting to get paid more quickly would pay off in the long run. I've heard bad things about some bootcamps found online, so always do your research before signing up. If you've got the connections or are willing to put the work in organising it yourself, an internship is another great way to gain experience without the up-front cost.
\nMy experience with S2DS was positive, and that seems to invariably be the feeling among other alumni. Having this experience on my CV and a project to talk about in interviews was invaluable during my job hunt. If you do decide to sign up, feel free to put me as your reference, I think I get some Amazon vouchers or similar. One unique advantage from completing S2DS is the ongoing career advice and the community of alumni, which will be useful going forwards in my career.
\nTwo other online courses I completed were CS50x and CS50 SQL from Harvard.
\nCS50x is Harvard's introductory computer science course. It's an 11-week course with lectures, short videos on specific subjects, and problem sets each week. Whilst not essential information for a data scientist, I feel having a much better understanding of how a computer works has made me a better data scientist. The first few weeks working through problems using C made me appreciate how easy Python is!
\nCS50 SQL is a 7-week SQL course with the same structure that teaches you all the basics you'll need for a job in data science. Most of the 7 weeks are done using SQLite, but it later moves onto ProgreSQL and MySQL.
\nWith a background in chemistry research, cheminformatics seemed like a natural field for me. Giving me the opportunity to combine by experience in chemistry and drug discovery with my new skills in Python programming and data science. I looked for short online courses and textbooks to get experience with cheminformatics problems and gain an understanding of the basics.
\nI started with the Cheminformatics OLCC. This is an introductory course with 8 sections that cover the basics of cheminformatics. It starts with help setting up a Python environment, and covers topics such as representing molecules, chemical databases, QSAR modelling, and simple machine learning for chemistry.
\nI then worked through Deep Learning for the Life Sciences. As datasets in the life sciences get larger, deep learning becomes more powerful. This covers how deep learning is used on molecules, proteins, and nucleic acids, including code so you can follow along.
\nI also attended the AI4SD Machine Learning Summer School at the University of Southampton. This was a week of lectures on topics such as machine learning, github, and LaTeX, with a focus on applications in chemistry. The week ended with a hackathon where we worked as a team to solve a chemical property prediction problem.
\nTeachOpenCADD looks like another great resource. I've not yet had the chance to go through it in detail, but it has a wide range of Jupyter Notebooks on all sorts of cheminformatics topics.
\nI found attending meetings and conferences on AI and cheminformatics was a great way to learn about current areas of research and meet interesting people doing similar research. This list is mostly limited to my local area (Cambridgeshire, UK).
\nThere's a load of useful blogs and newsletters for keeping up to date with cheminformatics:
\nAnd a few useful journals to add to your RSS feed:
\nWhen I started writing this post, I didn't expect it to be quite so long! Looking back at the process of re-training I realise how much work it took, and how lucky I was to have family to help support me through it all. I'm also incredibly glad I did it! My current job suits me so much more than the lab work I was doing before.
\nI think the general process I went through worked well, but there are still a few changes I would make if I could do it all again. I've heard that getting your first data science job has three equally important parts:
\nIn the search for my first job, I over-prioritised improving my skills, spending a lot of time doing courses. Whilst this is important for succeeding at a job, getting the job requires a wider focus. If I were to do it again, I would spend more time on personal projects to put on my github, and networking with other data scientists, asking about what problems their companies have and how I might be able to help fix them. When I finally did get my first data science job, it's no surprise it was on the recommendation of a previous colleague.
","doi":"https://doi.org/10.59350/qfdrh-57d64","guid":"https://jonswain.github.io/I-want-to-become-a-data-scientist-but-I-have-no-idea-where-to-start","image":"https://jonswain.github.io/images/career_change/chemistry_to_data_science.png","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1714867200,"rid":"s8tvz-zy751","summary":":seedling: Time for a change In 2020 I was working as a synthetic chemist in a laboratory in Aotearoa (New Zealand). I was getting disillusioned with the work (and academia in general) and started looking for a new career direction and new challenges.","tags":["Data-science","Cheminformatics","Career-development"],"title":"I Want to Become a Data Scientist, but I Have No Idea Where to Start\u2026","updated_at":1788763540,"url":"https://jonswain.github.io/become-a-data-scientist-starting-guide/","version":"v1"}},{"document":{"authors":[{"contributor_roles":[],"family":"Swain","given":"Jonathan","url":"https://orcid.org/0000-0003-4457-1481"}],"blog":{"authors":[{"name":"Jon Swain","url":"https://orcid.org/0000-0003-4457-1481"}],"community_id":"13f55986-f209-443c-ae0d-2f9f3f521e5a","created":1788652800,"current_feed_url":null,"description":"I am a cheminformatician and data scientist, originally from the UK, but often found in Aotearoa (New Zealand). I'm interested in using data science and machine learning to solve problems in drug discovery.","doi":"https://doi.org/10.59350/jonswain","favicon":"https://rogue-scholar.org/api/communities/13f55986-f209-443c-ae0d-2f9f3f521e5a/logo","feed_format":"application/atom+xml","feed_url":"https://jonswain.github.io/feed.xml","filter":null,"generator":"Jekyll","home_page_url":"https://jonswain.github.io/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"jonswain","status":"active","subfield":"3002","title":"Jon Swain","updated":1778680800,"use_api":null},"blog_name":"Jon Swain","blog_slug":"jonswain","content_html":"This is part 1 of a planned three post series on working with large chemical libraries.\nThe notebook used to create this post and all the files can be found in this github repo.
\nEstimates of the size of chemical space (the number of different molecules that could exist) range from 1050 - 1080 possible compounds, with the upper limits approaching estimates of the total number of atoms in the universe (around 1080 atoms). It's never going to be possible to make and test every possible chemical compound, even if we limit it to only \"drug-like\" compounds. Even with fast computational methods it's not always possible to exhaustively screen virtual libraries. We need methods to find the best scoring compounds from within virtual libraries without scoring every single compound.
\nSmall companies may have hundreds of thousands to millions of physical compounds in their virtual library, and larger companies will have libraries of a few million. Virtual make-on-demand libraries, especially combinatorial libraries, are expanding into the billions of compounds (often called ultra-large libraries). Enamine REAL has 48 billion possible compounds. If your scoring function takes one second per compound, you could maybe screen one million compounds (but it will take 11.6 days), but a billion seconds is nearly 32 years. As the size of virtual libraries expands, it becomes impossible to screen them exhaustively.
\nThere are a range of methods for filtering the library for screening, such as similarity searching or Thompson sampling. In this post I'll discuss active learning.
\nActive learning is a machine learning method for searching large libraries when the scoring function that is too computationally expensive to label the full library of compounds. A selection of data is labelled with the expensive scoring function, and a machine learning model is trained on these labels and used to score all compounds from within the library. The compounds with the best scores from the machine learning model are labelled using the more expensive function, and the data from this pooled and used to train a new machine learning model. This cycle is repeated until a finish criterion is met.
\nThe initial steps are:
\nFollowed then by a repeating cycle:
\nFirst we need to import the libraries we will be using.
\nimport gzip\nimport math\nimport re\nimport shutil\nimport subprocess\nfrom functools import partial\nfrom itertools import product\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem, DataStructs, Descriptors\nfrom sklearn.ensemble import RandomForestRegressor\nfrom tqdm import tqdm\nThe first function we will define is the expensive scoring function, this will take a list of SMILES and return a list of scores. For this example, I'm going to try find the compound from within the library with the lowest calculated Log P. This is actually a very fast calculation and can be done exhaustively, which means we can confirm if active learning is finding the lowest value.
\ndef calc_logp(smiles: list[str]) -> list[float]:\n \"\"\"Calculates the logP value for a list of compounds.\n\n Args:\n smiles (pd.Series): The input molecules.\n\n Returns:\n np.arary: The scores of the molecules.\n \"\"\"\n return [Descriptors.MolLogP(Chem.MolFromSmiles(s)) for s in smiles]\nNext, we need to define some useful functions for the active learning pipeline.
\nThe first function creates a virtual library by combining three sets of building blocks. The smi files used here were borrowed from Pat Walters repository on Thompson sampling. The second creates the Morgan fingerprints for all compounds within the library for training the machine learning model, which is trained using the third function. The final function scores the entire library using the machine learning model.
\ndef build_virtual_library() -> pd.DataFrame:\n \"\"\"Builds a virtual library by coupling building blocks from the input smi files.\n\n Returns:\n pd.DataFrame: A DataFrame containing the molecular objects and SMILES strings of\n the products.\n \"\"\"\n try:\n library = pd.read_csv(\"data/library.csv\", index_col=\"smiles\")\n library[\"mol\"] = [Chem.MolFromSmiles(s) for s in tqdm(library.index.to_list())]\n except FileNotFoundError:\n reaction_smarts = \"N[c:4][c:3]C(O)=O.[#6:1][NH2].[#6:2]C(=O)[OH]>>[C:2]c1n[c:4][c:3]c(=O)n1[C:1]\"\n bb_types = [\"aminobenzoic\", \"carboxylic_acids\", \"primary_amines\"]\n rxn = AllChem.ReactionFromSmarts(reaction_smarts)\n\n building_blocks = []\n for bb in bb_types:\n smil = []\n with open(Path(f\"data/{bb}_100.smi\"), \"r\") as f:\n for line in f.readlines():\n smiles, _ = line.split()\n smil.append(smiles)\n building_blocks.append(smil)\n\n total_prods = math.prod([len(x) for x in building_blocks])\n\n product_list = []\n for reagents in tqdm(product(*building_blocks), total=total_prods):\n reagent_mol_list = [Chem.MolFromSmiles(x) for x in reagents]\n products = rxn.RunReactants(reagent_mol_list)\n if products:\n Chem.SanitizeMol(products[0][0])\n product_list.append(products[0][0])\n\n library = pd.DataFrame(\n product_list,\n index=[Chem.MolToSmiles(m) for m in product_list],\n columns=[\"mol\"],\n )\n library.index.name = \"smiles\"\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n library.to_csv(\"data/library.csv\")\n return library\n\n\ndef create_morgan_fingerprints(library: pd.DataFrame) -> pd.DataFrame:\n \"\"\"Creates Morgan fingerprints for the input library.\n\n Args:\n library (pd.DataFrame): The input library.\n\n Returns:\n pd.DataFrame: The Morgan fingerprints of the input library.\n \"\"\"\n try:\n fps_df = pd.read_csv(\"data/fingerprints.csv\").set_index(\"smiles\")\n except FileNotFoundError:\n fps = [\n list(AllChem.GetMorganFingerprintAsBitVect(mol, radius=2))\n for mol in tqdm(library.mol)\n ]\n fps_df = pd.DataFrame(\n fps, columns=[f\"fp_{x}\" for x in range(len(fps[0]))], index=library.index\n )\n fps_df.to_csv(\"data/fingerprints.csv\")\n return fps_df\n\n\ndef train_ml_model(\n library: pd.DataFrame, fingerprints: pd.DataFrame\n) -> RandomForestRegressor:\n \"\"\"Trains a random forest regressor model on slow scores for the input library.\n\n Args:\n library (pd.DataFrame): The input library.\n fingerprints (pd.DataFrame): The Morgan fingerprints of the input library.\n\n Returns:\n RandomForestRegressor: The trained random forest regressor model.\n \"\"\"\n scored = library[~library[\"slow_scores\"].isna()]\n X = fingerprints.loc[scored.index]\n y = scored[\"slow_scores\"]\n regressor = RandomForestRegressor(max_depth=10)\n regressor.fit(X, y)\n return regressor\n\n\ndef score_library(\n library: pd.DataFrame,\n regressor: RandomForestRegressor,\n fingerprints: pd.DataFrame,\n) -> pd.DataFrame:\n \"\"\"Scores the entire library with the trained model.\n\n Args:\n library (pd.DataFrame): The input library.\n regressor (RandomForestRegressor): The trained random forest regressor model.\n fingerprints (pd.DataFrame): The Morgan fingerprints of the input library.\n\n Returns:\n pd.DataFrame: The input library with the model scores.\n \"\"\"\n library[\"model_scores\"] = regressor.predict(fingerprints)\n return library\nWe then need to set up the active learning pipeline. This follows the steps from above:
\nFollowed then by a repeating cycle:
\nBecause we sometimes know the best value for a scoring function, I've added an early stopping condition if the best compound is found to save time.
\ndef run_active_learning(\n library: pd.DataFrame,\n fingerprints: pd.DataFrame,\n compounds_per_round: int,\n number_of_rounds: int,\n scoring_function: callable,\n minimize: bool,\n early_stopping_value: float | None = None,\n) -> pd.DataFrame:\n \"\"\"Runs active learning on the virtual library.\n\n Args:\n library (pd.DataFrame): The input virtual library.\n fingerprints (pd.DataFrame): The Morgan fingerprints of the input library.\n compounds_per_round (int): Number of compounds to select per round.\n number_of_rounds (int): Number of active learning rounds.\n scoring_function (callable): The slow scoring function.\n minimize (bool): Whether to minimize or maximize the slow scoring function.\n early_stopping_value (float | None): If the slow scoring function finds a\n compound with this value, stops the active learning.\n\n Returns:\n pd.DataFrame: The input library with the slow scores added.\n \"\"\"\n # Select initial random sample\n initial_sample = library.sample(compounds_per_round)\n library[\"slow_scores\"] = np.NaN\n\n # Score the initial sample\n initial_scores = scoring_function(initial_sample.index.to_list())\n\n # Save the slow scores\n library.loc[initial_sample.index, \"slow_scores\"] = initial_scores\n\n # Run active learning\n al_round = 0\n while al_round < number_of_rounds:\n # Train the ML model\n model = train_ml_model(library, fingerprints)\n\n # Use the model to score the entire virtual library\n library = score_library(library, model, fingerprints)\n\n # Select the top scoring molecules with no slow scores\n top_compounds = (\n library[library[\"slow_scores\"].isna()]\n .sort_values(\"slow_scores\", ascending=minimize)\n .head(compounds_per_round)\n )\n\n # Score the top molecules with the slow function\n slow_scores = scoring_function(top_compounds.index.to_list())\n # Save the slow scores\n library.loc[top_compounds.index, \"slow_scores\"] = slow_scores\n library.loc[top_compounds.index, \"scored_round\"] = al_round\n\n al_round += 1\n if early_stopping_value is not None:\n if minimize:\n if library[\"slow_scores\"].min() == early_stopping_value:\n break\n else:\n if library[\"slow_scores\"].max() == early_stopping_value:\n break\n\n return library\nCalculated LogP (cLogP) is an important property in drug discovery. Whilst we wouldn't usually be aiming to find the minimum or maximum value in a library, it's a useful example as it can be calculated quickly enough to exhaustively screen the library to confirm we're finding the best compound. In the example below, two compounds are randomly selected from the library and their Morgan fingerprints and cLogP values are used to train a machine learning model. This model scores the entire library and the cLogP for two best scoring compounds are pooled with the earlier compounds and the process repeated. After only six rounds (12 compounds scored), the active learning algorithm can find the lowest cLogP in the library, much quicker than exhaustively screening all 132,500 compounds! The experiment was repeated ten times and the number of compounded needed is constant across the repeats. cLogP is a simple property, but it's still impressive how little data the machine learning model needed to find the best scores.
\n# Active learning parameters\nCOMPOUNDS_PER_ROUND = 2\nNUMBER_OF_ROUNDS = 50\nMINIMIZE = True\nNUMBER_OF_REPEATS = 10\n\n# Create the virtual library\nprint(\"Creating virtual library\")\nlibrary = build_virtual_library()\n\n# Create morgan fingerprints for the library\nprint(\"Creating Morgan fingerprints\")\nfingerprints = create_morgan_fingerprints(library)\n\n# Find the minimum logP in the library\nall_clogp_values = calc_logp(library.index.to_list())\nprint(f\"Minimum logP in the library: {min(all_clogp_values):.2f}\")\n\n# Run active learning\nprint(\"Running active learning\")\nscored_compounds_per_run = []\nfor i in tqdm(range(NUMBER_OF_REPEATS)):\n # Reset library between rounds\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n\n # Active learning\n library = run_active_learning(\n library=library,\n fingerprints=fingerprints,\n compounds_per_round=COMPOUNDS_PER_ROUND,\n number_of_rounds=NUMBER_OF_ROUNDS,\n scoring_function=calc_logp,\n minimize=MINIMIZE,\n early_stopping_value=min(all_clogp_values),\n )\n\n # Log results\n scored_compounds_per_run.append(len(library.loc[library.slow_scores.notna()]))\n\nave_scored_compounds = np.mean(scored_compounds_per_run) - COMPOUNDS_PER_ROUND\nprint(\n f\"Average number of scored compounds before finding lowest value: {ave_scored_compounds:.2f}\"\n)\nCreating virtual library\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 1000000/1000000 [01:42<00:00, 9747.08it/s]\n\nCreating Morgan fingerprints\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 132500/132500 [04:11<00:00, 526.70it/s]\n\nMinimum logP in the library: -5.00\nRunning active learning\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:21<00:00, 2.11s/it]\n\nAverage number of scored compounds before finding lowest value: 12.00\nMolecular similarity is often defined as the Tanimoto similarity between Morgan fingerprints. In the next experiment we'll try to maximise the Tanimoto similarity to a reference molecule. The reference molecule will be randomly selected from the library, so the maximum similarity will be 1. This turned out to be a significantly more difficult problem than the previous example minimizing cLogP. The number of compounds per round and number of rounded were increased (5000 compounds screened), and even still the reference molecule was only recovered once out of the ten repeats. This was a slightly surprising result considering the machine learning models were trained using the Morgan fingerprints that were also used in the Tanimoto similarity.
\ndef calc_similarity(comparison_smiles: list[str], ref_smiles: str) -> list[float]:\n \"\"\"Calculates the Tanimoto similarity of a reference compound to a list of\n comparison compounds.\n\n Args:\n comparison_smiles (list[str]): List of SMILES to compare to\n ref_smiles (str): SMILES of the reference compound\n\n Returns:\n list[float]: List of similarity scores\n \"\"\"\n ref_mol = Chem.MolFromSmiles(ref_smiles)\n comparison_mols = [Chem.MolFromSmiles(x) for x in comparison_smiles]\n fpgen = AllChem.GetMorganGenerator()\n ref_fp = fpgen.GetFingerprint(ref_mol)\n comparison_fps = [fpgen.GetFingerprint(x) for x in comparison_mols]\n return [DataStructs.FingerprintSimilarity(ref_fp, x) for x in comparison_fps]\n\n\n# Active learning parameters\nCOMPOUNDS_PER_ROUND = 100\nNUMBER_OF_ROUNDS = 50\nMINIMIZE = False\nNUMBER_OF_REPEATS = 10\n\n# Run active learning\nprint(\"Running active learning\")\nscored_compounds_per_run = []\nfor i in tqdm(range(NUMBER_OF_REPEATS)):\n # Choose a random reference molecule\n ref_smiles = library.sample(1).index[0]\n\n # Reset library between rounds\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n\n # Active learning\n library = run_active_learning(\n library=library,\n fingerprints=fingerprints,\n compounds_per_round=COMPOUNDS_PER_ROUND,\n number_of_rounds=NUMBER_OF_ROUNDS,\n scoring_function=partial(calc_similarity, ref_smiles=ref_smiles),\n minimize=MINIMIZE,\n early_stopping_value=1,\n )\n\n # Log results\n top = library[~library.slow_scores.isna()].sort_values(\n \"slow_scores\", ascending=MINIMIZE\n )\n scored_compounds_per_run.append(\n {\n \"round\": i,\n \"reference_molecule\": ref_smiles,\n \"compounds_scored\": len(library[~library.slow_scores.isna()]),\n \"top_scorer\": top.index.to_list()[0],\n \"top_score\": top.slow_scores.values[0],\n }\n )\n\npd.DataFrame(scored_compounds_per_run).set_index(\"round\").style.format(\n {\"top_score\": \"{:.3f}\"}\n)\nRunning active learning\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [22:30<00:00, 135.02s/it]\n| Round | \nReference compound | \nCompounds scored | \nTop score | \n
|---|---|---|---|
| 0 | \nCOc1ccc2c(=O)n(C(=N)NCC(=O)O)c([C@@H](CN)CO)nc2c1OC | \n5100 | \n0.466 | \n
| 1 | \nCC(=O)OC[C@H](C(=O)O)n1c([C@H](N)CN(C)C)nc2[nH][nH]c(=O)c2c1=O | \n3300 | \n1.000 | \n
| 2 | \nCc1cc2nc([C@H](C)NC(=O)CN)n(CC[C@H](N)C(=O)O)c(=O)c2cc1O | \n5100 | \n0.506 | \n
| 3 | \nNc1nncn1Cc1nc2c(cnc3ccnn32)c(=O)n1C(=O)[C@H](N)CC(=O)O | \n5100 | \n0.495 | \n
| 4 | \nCn1nnc2nc(C3(O)CC(N)C3)n([C@H](CO)C(=O)O)c(=O)c21 | \n5100 | \n0.488 | \n
| 5 | \nCn1cc2c(=O)n(CC[C@@H](N)C(=O)O)c([C@H](CO)NC(N)=O)nc2n1 | \n5100 | \n0.518 | \n
| 6 | \nCOC(=O)c1ccc2c(=O)n(CC[C@@H](N)C(=O)O)c([C@@H](N)CCO)nc2c1 | \n5100 | \n0.449 | \n
| 7 | \nC[C@H](O)[C@H](N)c1nc2cnc(Cl)nc2c(=O)n1[C@H](C(=O)O)[C@@H](C)O | \n5100 | \n0.434 | \n
| 8 | \nCn1cnc2c(=O)n(C(=O)C[C@H](N)C(=O)O)c([C@@H](N)CO)nc21 | \n5100 | \n0.519 | \n
| 9 | \nC[C@@H](O)[C@H](N)c1nc2nc(N)ncc2c(=O)n1C(=O)NCC(=O)O | \n5100 | \n0.458 | \n
Looking at the final run, we can see that in general the maximum similarity is not increasing significantly with further rounds, apart from the occasional spike as the model finds a similar compound. This suggests that increasing the number of rounds may not increase the performance.
\nsns.lineplot(\n library[[\"slow_scores\", \"scored_round\"]].groupby(\"scored_round\").max().values\n).set(xlabel=\"Round\", ylabel=\"Max similarity to reference compound\")\nplt.show()\n\n![\"A](\"https://jonswain.github.io/images/active_learning/active_learning_13_0.png\"/)
\n
AL was very successful at finding the lowest cLogP value in a library but struggled to recover a randomly chosen reference molecule by trying to maximize the Tanimoto similarity. The reason behind this may be due to the underlying distribution of data within the library. By plotting the cLogP values of all compounds within the library, we can see that the data is close to being normally distributed. The machine learning model is able to learn to predict cLogP well from the data it is given. For the Tanimoto similarity, the scores for the whole library of similarity to the reference compound from the final run are plotted, which have a large positive skew. There is lots of data for predicting Tanimoto similarities between 0.1 and 0.4, but very little data for scores above 0.4. So even though the machine learning model is using the same representation as is used to calculate the Tanimoto similarity (Morgan fingerprint), it doesn't have enough information to learn what makes a high score.
\nall_similarities = calc_similarity(library.index.to_list(), ref_smiles)\nfig, ax = plt.subplots(1, 2, figsize=(12, 5))\nsns.histplot(all_clogp_values, bins=50, alpha=0.5, label=\"True scores\", ax=ax[0])\nax[0].set(title=\"cLogP values for library\", xlabel=\"cLogP\", ylabel=\"Frequency\")\nsns.histplot(all_similarities, bins=50, alpha=0.5, label=\"Predicted scores\", ax=ax[1])\nax[1].set(\n title=\"Tanimoto similarity to reference compound\",\n xlabel=\"Similarity\",\n ylabel=\"Frequency\",\n)\nplt.show()\n\n![\"A](\"https://jonswain.github.io/images/active_learning/active_learning_15_0.png\"/)
\n
Molecular docking is a common technique used in drug discovery to predict how a small molecule ligand will bind to a protein target. This can be combined with a scoring function to estimate the binding affinity, and these scores can be used to rank virtual compounds for experimental validation. I am far from an expert on molecular docking, but here we'll use it as an example of how active learning could be used in a real drug discovery project.
\nDocking is slow. In this example we'll use SMINA, which seem to take around one minute for each compound. We're docking compounds to the SARS MPro protein and using active learning to try find the highest binding affinity by minimizing the predicted binding energy. At one minute per docking calculation, it would take 13 weeks to exhaustively screen the virtual library used here.
\nThe experiment was repeated three times, twice it improved on the docking scores of the initial random sample, and on one run it found a very good scoring compound in the initial random sample and wasn't able to improve on that. Due to running this on my laptop, I kept the number of compounds per round (five) and the number of rounds (10) low. With more computing power, it would be interesting to see how the performance of the active learning improves with increasing these values.
\ndef dock_mols(smiles: list[str]) -> list[float]:\n \"\"\"Dock the molecules to the SARS MPro protein and return the affinity of the best\n pose for each.\n\n Args:\n smiles: list[str]: The molecules to dock\n\n Returns:\n list[float]: The affinities of the best pose for each input molecule\n \"\"\"\n mols = [Chem.MolFromSmiles(s) for s in smiles]\n scores = []\n\n for mol in mols:\n with Chem.SDWriter(\"data/docking/tmp_conf.sdf\") as w:\n m = Chem.AddHs(mol)\n _ = AllChem.EmbedMultipleConfs(m, numConfs=5, numThreads=0)\n confs = m.GetConformers()\n for c in confs:\n w.write(m, confId=c.GetId())\n\n subprocess.run(\n [\n \"./smina.osx.12\",\n \"--exhaustiveness\",\n \"10\",\n \"--cpu\",\n \"10\",\n \"--seed\",\n \"0\",\n \"--autobox_ligand\",\n \"data/docking/ligand_only.pdb\",\n \"-r\",\n \"data/docking/protein_minus_ligand.pdb\",\n \"-l\",\n \"data/docking/tmp_conf.sdf\",\n \"-o\",\n \"data/docking/tmp_conf_docked.sdf.gz\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n )\n\n with gzip.open(\"data/docking/tmp_conf_docked.sdf.gz\", \"rb\") as f_in:\n with open(\"data/docking/tmp_conf_docked.sdf\", \"wb\") as f_out:\n shutil.copyfileobj(f_in, f_out)\n\n # The output sdf wasn't loading so this is a workaround\n with open(\"data/docking/tmp_conf_docked.sdf\", \"r\") as f:\n text = f.read()\n affinities = re.findall(r\"\\n(-\\d.\\d+)\" , text)\n affinities = [float(x) for x in affinities]\n\n scores.append(min(affinities))\n\n return scores\n\n\n# Active learning parameters\nCOMPOUNDS_PER_ROUND = 5\nNUMBER_OF_ROUNDS = 10\nMINIMIZE = True\nNUMBER_OF_REPEATS = 3\n\n# Run active learning\nprint(\"Running active learning\")\nscored_compounds_per_run = []\ncumulative_scores = {}\nfor i in tqdm(range(NUMBER_OF_REPEATS)):\n # Reset library between rounds\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n\n # Active learning\n library = run_active_learning(\n library=library,\n fingerprints=fingerprints,\n compounds_per_round=COMPOUNDS_PER_ROUND,\n number_of_rounds=NUMBER_OF_ROUNDS,\n scoring_function=dock_mols,\n minimize=MINIMIZE,\n )\n\n # Log results\n top = library[~library.slow_scores.isna()].sort_values(\n \"slow_scores\", ascending=MINIMIZE\n )\n scored_compounds_per_run.append(library.slow_scores.min())\n cumulative_scores[f\"run_{i}\"] = [\n library[library.scored_round <= j].slow_scores.min()\n for j in range(NUMBER_OF_ROUNDS)\n ]\n\nprint(\n f\"Average lowest docking score after screening {COMPOUNDS_PER_ROUND * NUMBER_OF_ROUNDS} compounds: {sum(scored_compounds_per_run) / len(scored_compounds_per_run):.2f}\"\n)\n\n# Plot the cumulative best score found in each round for the final run\nsns.lineplot(cumulative_scores).set(\n xlabel=\"Round\", ylabel=\"Minimum docking score\", xticks=range(NUMBER_OF_ROUNDS)\n)\nplt.show()\nRunning active learning\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 3/3 [4:54:58<00:00, 5899.39s/it] \n\nAverage lowest docking score after screening 50 compounds: -9.14\n\n![\"A](\"https://jonswain.github.io/images/active_learning/active_learning_17_4.png\"/)
\n
In this post we've discussed active learning and seen three examples of how it could be used when our virtual library has grown too large to exhaustively screen every compound.
\nI've tried to visualise and compare distributions using violin plots for reports and presentations in the past, and the feedback I've got was generally\u2026 not great. When searching for better methods I came across this excellent blog post by Alex Belengeanu on raincloud plots and I'm now a big fan.
\nThere's lots of methods for visualising distributions to compare different populations, and each has its own pros and cons. Histograms (often with a kernel density estimation) are good method for showing complex distributions, but it's not possible to easily compare statistics such as the mean or quantiles from the raw graph, and there's also some loss of information from binning or smoothing, which could misrepresent the data, such as making discrete values appear continuous. Box and whisker plots are great for being able to compare statistics such as means and quartiles, but they don't represent complex distributions such as bimodal data well. A jitter plot is good for showing the structure of the underlying data, but difficult to compare statistics or absolute densities of points.
\nA violin plot with box plot combines the strengths of the two methods to visualise complex distributions with easy to read off statistics. The main downside is they're ugly (IMO). Every time I've used them in a report or presentation the discussion has veered off to whether they look more like sea creatures or medieval weapons than violins, with very little discussion of the actual data!
\nWhilst searching for better methods I came across this blog post by Alex Belengeanu on raincloud plots, which seem to be the best alternative.
\nRaincloud plots are pretty simple, they're just a combination of half a violin plot (sometimes called a ridgeline plot), a box plot, and a jitter plot. They're intuitive to read and show distributions clearly, as well as being much more aesthetically pleasing than many alternatives. I've customised Alex Belengeanu's code to suit my needs, and it's all in a function in this GitHub repo.
\n![\"A](\"https://jonswain.github.io/images/raincloudplot/petals.png\"/)
This is part 2 of a a planned three post series on working with large chemical libraries.\nThe notebook used to create this post, and all the files can be found in this github repo.
\nCombinatorial libraries can grow quickly, combining three sets of 1,000 building blocks allows you to access one billion possibilities. For example, 1,000 amines, 1,000 halide functionalised carboxylic acids, and 1,000 boronic acids can be combined to make a billion-member library. That would be a lot of synthesis and a lot of testing if you wanted to screen every member of the library! Even if a computational method took 1 second to score each virtual compound, it would take 32 years to score the entire library. Make-on-demand chemical suppliers such as Enamine have libraries that contain multiple billions of compounds, often created by combining sets of building blocks using reliable chemistry. We need methods for effectively sampling this chemical space.
\nIn this post I'll be looking into using a genetic algorithm to search a small combinatorial library. These are quick and incredibly simple methods that can be used to search combinatorial libraries for the combination of building blocks that maximises or minimises a scoring function.
\nWhilst this method uses a combinatorial library that contains all possible configurations of all three building blocks, it's possible to extend it to more complex combinatorial libraries by storing the building blocks and reactions as a graph.
\nGenetic algorithms are biologically inspired, based on biological natural selection. They use a very simple algorithm, with no machine learning or complex statistics. One of the big advantages is that you don't need to fully enumerate of score the entire library, only certain combinations of building blocks need to be enumerated and scores, reducing memory and computational costs. Working in building block space rather than with enumerated structures means the complexity will scale with number of building blocks, which will increase much more slowly than the total library size. The downside of the genetic algorithm is that it requires combinatorial libraries and won't work with large collections of diverse molecules not made up from the same building blocks.
\nInitially a random population of the unlabelled data is selected by randomly choosing building blocks and enumerating the reaction products from these building blocks.
\nThe genetic algorithm has a cycle that is repeated, and each round is called a generation.
\nFirst, we need to import the libraries we will be using.
\nimport math\nimport time\nfrom functools import partial\nfrom itertools import product\nfrom pathlib import Path\nfrom typing import Callable\n\nimport numpy as np\nimport pandas as pd\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem, DataStructs, Descriptors\nfrom tqdm import tqdm\nThe first function we will define is the expensive scoring function, this will take a list of RDKit molecules and return a list of scores. As in the previous example, I'm going to try find the compound from within the library with the lowest calculated Log P (cLogP). This is actually a very fast calculation and can be done exhaustively, which means we can confirm if the genetic algorithm is finding the lowest value and triggering the early stopping.
\ndef calc_logp(mols: list[Chem.Mol]) -> list[float]:\n \"\"\"Calculate the logP value for a list of compounds.\n\n Args:\n mols (list[Chem.Mol]): The molecules.\n\n Returns:\n list[float]: The scores of the molecules.\n \"\"\"\n return [Descriptors.MolLogP(mol) for mol in mols]\nNext, we need to define some useful functions for the active learning pipeline.
\nThe first function takes the indices of three building blocks and enumerates the product of their reaction. Since not all combinations are possible, it will return None if the reaction fails. The second function takes the surviving population from the selection pressure and creates a new population by shuffling the building blocks. The third takes the new shuffled population and adds in random mutations, this prevents the algorithm getting stuck in a local minimum by ensuring building block diversity.
\ndef make_molecule(\n r1: int,\n r2: int,\n r3: int,\n building_blocks: list[list[str]],\n rxn: Chem.rdChemReactions.ChemicalReaction,\n) -> Chem.Mol | None:\n \"\"\"React building blocks to make a molecule.\n\n Args:\n r1 (int): Index of the first building block.\n r2 (int): Index of the second building block.\n r3 (int): Index of the third building block.\n building_blocks (list[list[str]]): The building blocks.\n rxn (Chem.rdChemReactions.ChemicalReaction): The reaction to combine the\n building blocks.\n\n\n Returns:\n Chem.Mol | None: The molecule or None if the reaction fails.\n \"\"\"\n bbs = [building_blocks[0][r1], building_blocks[1][r2], building_blocks[2][r3]]\n reagent_mol_list = [Chem.MolFromSmiles(x) for x in bbs]\n products = rxn.RunReactants(reagent_mol_list)\n if products:\n Chem.SanitizeMol(products[0][0])\n return products[0][0]\n return None\n\n\ndef shuffle_population(data: pd.DataFrame, population_size: int) -> pd.DataFrame:\n \"\"\"Shuffle the population building blocks.\n\n Args:\n data (pd.DataFrame): The surviving population data.\n population_size (int): The size of the population.\n\n Returns:\n pd.DataFrame: The shuffled population.\n \"\"\"\n r1s = data[\"r1\"].to_list()\n r2s = data[\"r2\"].to_list()\n r3s = data[\"r3\"].to_list()\n new_population = pd.DataFrame(\n {\n \"r1\": [np.random.choice(r1s) for _ in range(population_size)],\n \"r2\": [np.random.choice(r2s) for _ in range(population_size)],\n \"r3\": [np.random.choice(r3s) for _ in range(population_size)],\n }\n )\n return new_population\n\n\ndef mutate_population(\n data: pd.DataFrame, mutation_rate: float, building_blocks: list[list[str]]\n) -> pd.DataFrame:\n \"\"\"Mutate the population building blocks.\n\n Args:\n data (pd.DataFrame): The surviving population data.\n mutation_rate (float): The fraction of building blocks to mutate.\n building_blocks (list[list[str]]): The building blocks.\n\n Returns:\n pd.DataFrame: The mutated population.\n \"\"\"\n for count, column in enumerate(data.columns):\n selection = data.sample(frac=mutation_rate)\n data.loc[selection.index, column] = np.random.randint(\n 0, len(building_blocks[count]), len(selection)\n )\n return data\nThe genetic algorithm needs the lists of building blocks, a reaction to couple them to enumerate products, and a function to score the molecules. It also has a number of hyperparameters than can be tuned:
\npopulation_size: The number of molecules in each generation to enumerate. A larger population will have a greater diversity of building blocks but will take longer to enumerate and score.num_generations: The number of cycles of the genetic algorithm to run. More generations will give a better chance of finding the optimum value but will take longer to run.selection_pressure: The fraction of the population to remove each cycle. A stronger selection pressure will speed up selection but increase the chance of being caught in a local minimum.mutation_rate: The proportion of building blocks in the new population to replace with a random building block. A too high value will stop the algorithm finding the best combination, but a too low value will also increase the chance of being caught in a local minimum.def run_genetic_algorithm(\n building_blocks: list[list[str]],\n rnx: Chem.rdChemReactions.ChemicalReaction,\n population_size: int,\n num_generations: int,\n selection_pressure: float,\n mutation_rate: float,\n scoring_function: Callable,\n minimize: bool,\n early_stopping_value: float | None = None,\n) -> pd.DataFrame:\n \"\"\"Run the genetic algorithm.\n\n Args:\n building_blocks (list[list[str]]): The building blocks.\n rnx (Chem.rdChemReactions.ChemicalReaction): The reaction to combine the\n building blocks.\n population_size (int): The size of the population in each generation.\n num_generations (int): The number of generations to run.\n selection_pressure (float): The fraction of the population to discard.\n mutation_rate (float): The fraction of the population to mutate.\n scoring_function (Callable): The function to score the molecules.\n minimize (bool): Whether to minimize the scoring function.\n early_stopping_value (float | None, optional): The value to stop early at.\n Defaults to None.\n\n Returns:\n pd.DataFrame: The history of the population.\n \"\"\"\n # To keep track of the population each generation\n history = pd.DataFrame()\n\n # Choose initial population\n population = pd.DataFrame(\n {\n \"r1\": np.random.randint(0, len(building_blocks[0]), population_size),\n \"r2\": np.random.randint(0, len(building_blocks[1]), population_size),\n \"r3\": np.random.randint(0, len(building_blocks[2]), population_size),\n }\n )\n\n for generation in range(num_generations):\n # Generate molecules\n population[\"ROMol\"] = population.apply(\n lambda x: make_molecule(x[\"r1\"], x[\"r2\"], x[\"r3\"], building_blocks, rnx),\n axis=1,\n )\n\n # Kill off fatal mutants\n population = population.dropna().reset_index(drop=True)\n\n # Score molecules\n population[\"score\"] = scoring_function(population[\"ROMol\"].to_list())\n\n # Save population for analysis\n population[\"generation\"] = generation\n history = pd.concat([history, population])\n\n # Early stopping\n if early_stopping_value:\n if minimize:\n if round(population[\"score\"].min(), 5) <= round(\n early_stopping_value, 5\n ):\n break\n else:\n if round(population[\"score\"].max(), 5) >= round(\n early_stopping_value, 5\n ):\n break\n\n # Select top performing molecules\n population = (\n population.sort_values(\"score\", ascending=minimize)\n .head(int(population_size * (1 - selection_pressure)))\n .reset_index(drop=True)\n )\n\n # Shuffle and mutate\n population = shuffle_population(population, population_size)\n population = mutate_population(population, mutation_rate, building_blocks)\n population = population.drop_duplicates().reset_index(drop=True)\n\n return history\nThis function fully enumerates a virtual library by combining three sets of building blocks. The smi files used here were borrowed from Pat Walters repository on Thompson sampling. This is not necessary for the genetic algorithm but will allow us to see if the genetic algorithm if finding the best values.
\ndef build_virtual_library() -> pd.DataFrame:\n \"\"\"Build a virtual library by coupling building blocks from the input smi files.\n\n Returns:\n pd.DataFrame: The virtual library.\n \"\"\"\n try:\n library = pd.read_csv(\"data/library.csv\", index_col=\"smiles\")\n library[\"mol\"] = [Chem.MolFromSmiles(s) for s in tqdm(library.index.to_list())]\n except FileNotFoundError:\n reaction_smarts = \"N[c:4][c:3]C(O)=O.[#6:1][NH2].[#6:2]C(=O)[OH]>>[C:2]c1n[c:4][c:3]c(=O)n1[C:1]\"\n bb_types = [\"aminobenzoic\", \"carboxylic_acids\", \"primary_amines\"]\n rxn = AllChem.ReactionFromSmarts(reaction_smarts)\n\n building_blocks = []\n for bb in bb_types:\n smil = []\n with open(Path(f\"data/{bb}_100.smi\"), \"r\") as f:\n for line in f.readlines():\n smiles, _ = line.split()\n smil.append(smiles)\n building_blocks.append(smil)\n\n total_prods = math.prod([len(x) for x in building_blocks])\n\n product_list = []\n for reagents in tqdm(product(*building_blocks), total=total_prods):\n reagent_mol_list = [Chem.MolFromSmiles(x) for x in reagents]\n products = rxn.RunReactants(reagent_mol_list)\n if products:\n Chem.SanitizeMol(products[0][0])\n product_list.append(products[0][0])\n\n library = pd.DataFrame(\n product_list,\n index=[Chem.MolToSmiles(m) for m in product_list],\n columns=[\"mol\"],\n )\n library.to_csv(\"data/library.csv\")\n return library\n\n\n# Chemistry parameters\nRXN = AllChem.ReactionFromSmarts(\n \"N[c:4][c:3]C(O)=O.[#6:1][NH2].[#6:2]C(=O)[OH]>>[C:2]c1n[c:4][c:3]c(=O)n1[C:1]\"\n)\nBUILDING_BLOCKS = []\nfor bb in [\"aminobenzoic\", \"carboxylic_acids\", \"primary_amines\"]:\n smil = []\n with open(Path(f\"data/{bb}_100.smi\"), \"r\") as f:\n for line in f.readlines():\n smiles, _ = line.split()\n smil.append(smiles)\n BUILDING_BLOCKS.append(smil)\nThe genetic algorithm was used to find the lowest cLogP from within the library. This was repeated 10 times with no tuning of the hyperparameters which may improve performance. The genetic algorithm reliably finds the combination of building blocks with the lowest cLogP in the combinatorial library, nearly 100 times faster than enumerating the entire library. On average the genetic algorithm only had to enumerate and score 600 combinations before it found the best scoring combination.
\n# GA parameters\nPOPULATION_SIZE = 500\nNUM_GENERATIONS = 20\nSELECTION_PRESSURE = 0.5\nMUTATION_RATE = 0.1\nMINIMIZE = True\nNUMBER_OF_REPEATS = 10\n\n# Create the virtual library\nvirtual_library_start = time.time()\nprint(\"Creating virtual library\")\nlibrary = build_virtual_library()\nvirtual_library_end = time.time()\nprint(\n f\"Virtual library created in {virtual_library_end - virtual_library_start:.2f} seconds\",\n)\n\n# Find the minimum logP in the library\nlogp_start = time.time()\nprint(\"Calculating logP values for the library\")\nall_clogp_values = [Descriptors.MolLogP(mol) for mol in tqdm(library.mol.to_list())]\nlogp_end = time.time()\nprint(f\"Minimum logP in the library: {min(all_clogp_values):.2f}\")\nprint(f\"LogP calculations took {logp_end - logp_start:.2f} seconds\")\n\n# Run the genetic algorithm\nprint(f\"Running genetic algorithm {NUMBER_OF_REPEATS} times\")\nga_times = []\nmax_generations = []\nfor _ in tqdm(range(NUMBER_OF_REPEATS)):\n ga_start = time.time()\n logp_ga_run = run_genetic_algorithm(\n building_blocks=BUILDING_BLOCKS,\n rnx=RXN,\n population_size=POPULATION_SIZE,\n num_generations=NUM_GENERATIONS,\n selection_pressure=SELECTION_PRESSURE,\n mutation_rate=MUTATION_RATE,\n scoring_function=calc_logp,\n minimize=MINIMIZE,\n early_stopping_value=min(all_clogp_values),\n )\n ga_end = time.time()\n ga_times.append(ga_end - ga_start)\n max_generations.append(logp_ga_run.generation.max())\nprint(\n f\"Genetic algorithm run took {np.average(ga_times):.2f} seconds on average, with a maximum of {np.max(max_generations)} generations\"\n)\n\n# Calculate improvement\nga_time = ga_end - ga_start\nvirtual_library_time = virtual_library_end - virtual_library_start\nlogp_time = logp_end - logp_start\nimprovement = (virtual_library_time + logp_time) / ga_time\nprint(\n f\"GA was {improvement:.2f} times faster than building and scoring the virtual library\"\n)\nCreating virtual library\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 1000000/1000000 [01:48<00:00, 9183.87it/s]\nVirtual library created in 120.60 seconds\n\nCalculating logP values for the library\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 132500/132500 [00:27<00:00, 4803.96it/s]\nMinimum logP in the library: -5.00\nLogP calculations took 27.59 seconds\n\nRunning genetic algorithm 10 times\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:14<00:00, 1.41s/it]\n\nGenetic algorithm run took 1.41 seconds on average, with a maximum of 12 generations\nGA was 81.90 times faster than building and scoring the virtual library\nA reference molecule was randomly selected from the library, and the genetic algorithm was used to recover it from the library by maximising the Tanimoto similarity. In the active learning experiments, this was failing, possibly due to the low number of compounds from within the library with high Tanimoto similarities. The genetic algorithm was able to recover the reference molecule every time, taking less than a second to do so each time.
\ndef calc_similarity(comparison_mols: list[Chem.Mol], ref_mol: Chem.Mol) -> list[float]:\n \"\"\"Calculate the Tanimoto similarity to a reference compound.\n\n Args:\n comparison_mols (list[Chem.Mol]): List of molecules to compare.\n ref_mol (Chem.Mol): The reference molecule.\n\n Returns:\n list[float]: List of similarity scores\n \"\"\"\n fpgen = AllChem.GetMorganGenerator()\n ref_fp = fpgen.GetFingerprint(ref_mol)\n comparison_fps = [fpgen.GetFingerprint(x) for x in comparison_mols]\n return [DataStructs.FingerprintSimilarity(ref_fp, x) for x in comparison_fps]\n\n\n# Score the virtual library\nprint(\"Scoring the virtual library\")\ntanimoto_start = time.time()\nall_tanimoto_values = calc_similarity(\n library.mol.to_list(), library.mol.sample().values[0]\n)\ntanimoto_end = time.time()\nprint(f\"Tanimoto calculations took {tanimoto_end - tanimoto_start:.2f} seconds\")\n\n# GA parameters\nPOPULATION_SIZE = 500\nNUM_GENERATIONS = 10\nSELECTION_PRESSURE = 0.5\nMUTATION_RATE = 0.1\nMINIMIZE = False\nNUMBER_OF_REPEATS = 10\n\n# Run the genetic algorithm\nprint(f\"Running genetic algorithm {NUMBER_OF_REPEATS} times\")\nga_times = []\nmax_generations = []\nfor _ in tqdm(range(NUMBER_OF_REPEATS)):\n # Choose random reference molecule\n ref_mol = library.mol.sample().values[0]\n\n # Run the genetic algorithm\n ga_start = time.time()\n tanimoto_ga_run = run_genetic_algorithm(\n building_blocks=BUILDING_BLOCKS,\n rnx=RXN,\n population_size=POPULATION_SIZE,\n num_generations=NUM_GENERATIONS,\n selection_pressure=SELECTION_PRESSURE,\n mutation_rate=MUTATION_RATE,\n scoring_function=partial(calc_similarity, ref_mol=ref_mol),\n minimize=MINIMIZE,\n early_stopping_value=1,\n )\n ga_end = time.time()\n ga_times.append(ga_end - ga_start)\n max_generations.append(tanimoto_ga_run.generation.max())\nprint(\n f\"Genetic algorithm run took {np.average(ga_times):.2f} seconds on average, with a maximum of {np.max(max_generations)} generations\"\n)\n\n# Calculate improvement\ntanimoto_time = tanimoto_end - tanimoto_start\nimprovement = tanimoto_time / np.average(ga_times)\nprint(\n f\"GA was {improvement:.2f} times faster than scoring the virtual library with Tanimoto similarity\"\n)\nScoring the virtual library\nTanimoto calculations took 7.41 seconds\nRunning genetic algorithm 10 times\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:07<00:00, 1.29it/s]\n\nGenetic algorithm run took 0.77 seconds on average, with a maximum of 9 generations\nGA was 9.60 times faster than scoring the virtual library with Tanimoto similarity\nDeep Learning models have traditionally performed well on unstructured data such as text and images, but poorly on structured tabular data, and are usually outperformed by Gradient Boosted Decision Trees (GBDTs) on tabular chemical data. TabPFN (Tabular Prior-data Fitted Network) is a transformer-based foundation model for tabular data, pre-trained on millions of synthetic datasets to solve supervised learning tasks, with state-of-the-art performance on benchmarks. But does it work for cheminformatics?
\nA recent paper published in Nature, Accurate predictions on small data with a tabular foundation model by Hollmann et al., describes a new Deep Learning (DL) algorithm for making predictions on tabular data, TabPFN (Tabular Prior-data Fitted Network). TabPFN is a transformer-based foundation model for tabular data that operates via in-context learning (ICL), enabling it to train and predict on an entire dataset in a single forward pass.
\nThe authors generated over 100 million synthetic datasets using structural causal models (SCMs). These models simulate causal relationships and mimic real-world tabular data challenges, such as non-linear relationships, missing values, outliers, and diverse feature types (e.g., categorical, ordinal, numerical). These datasets sampled high-level hyperparameters (such as dataset size and number of features), both classification and regression tasks, and Gaussian noise was added to mimic real-world complexities.
\nThese synthetic datasets were used to pre-train a transformer model. During the training, parts of the datasets were masked, and the model is trained to predict masked target values in synthetic datasets, given features and labelled samples as context. The parameters of the neural network were updated until the predictions matched the masked values. Through this training, the model learns to fill in missing data from a dataset. The model learns a generic prediction algorithm that approximates Bayesian inference for the synthetic data prior, enabling robust handling of unseen datasets.
\nWhen making predictions TabPFN uses ICL and processes an entire dataset (both labelled and test samples) in one pass, performing training and inference simultaneously.
\nThe Therapeutic data commons (TDC) aims to help development of AI/ML tools for therapeutic science by providing datasets and curated benchmarks to assess the performance of new methods. The TDC ADMET benchmark group contains 22 datasets for molecular property prediction, ranging from 475 to 13,130 entries, with both classification and regression tasks. The datasets contain SMILES strings for each chemical compound, and a target variable to be predicted. These can all be easily downloaded using their Python library.
\nTabPFN is designed for small tabular datasets, and the default parameters limit to training on 10,000 entries with 500 features. If a dataset had more than 10,000 entries in the training and validation sets, a random selection of 10,000 entries was used. For each dataset, the training and validation data were combined, and the test set put aside. For all entries, the 210 RDKit descriptors were calculated and used as the features for training. Predictions were made on the test set and compared to the true values to evaluate the performance of TabPFN. The training was repeated using molecular fingerprints (MACCS keys and ECPF folded to 500 bits), but the performance was found to be lower than using RDKit calculated properties.
\nNOTE: Due to memory limits on my computer, I had to limit to datasets with fewer than 1,800 entries in the training data. Training models for the larger datasets is ongoing.
\nUsing the RDKit calculated descriptors as features, TabPFN come in the top 10 models for all TDC datasets apart from \"Clearance_Hepatocyte_Az\". It comes 3rd for the \"Vdss_Lombardo\" dataset, 2nd for \"Caco2_Wang\", \"Pgp_Broccatelli\", and \"Bbb_Martins\" datasets, and is the highest performing model for the \"Clearance_Microsome_Az\" dataset. The high performing datasets include both classification and regression tasks and doesn't seem to be a clear link between classification performance and imbalanced datasets. \"Cyp2C9_Substrate_Carbonmangels\" is somewhat imbalanced (19.3% positive - TDC rank: 10th) but so are \"HIA_Hou\" (11.1% negative - TDC rank: 5th), and \"Bioavailability_Ma\", (22.9% negative - TDC rank: 5th). The strongest link seems to be between dataset size and performance, with all the highest performing models having >900 entries.
\nThe training time is an average of five repeats, using WSL on a Windows machine with 8 GB RAM allocated to WSL.
\n| Dataset | \nSize | \nTask | \nMetric | \nTraining time (min) | \nTabFPN performance | \nCurrent TDC best performance | \nTabPFN TDC leaderboard rank | \n
|---|---|---|---|---|---|---|---|
| Caco2_Wang | \n906 | \nRegression | \nMAE | \n8.25 | \n0.282 \u00b1 0.005 | \n0.276 \u00b1 0.005 | \n2nd | \n
| HIA_Hou | \n578 | \nClassification | \nAUROC | \n3.17 | \n0.987 \u00b1 0.001 | \n0.990 \u00b1 0.002 | \n5th | \n
| Pgp_Broccatelli | \n1218 | \nClassification | \nAUROC | \n12.77 | \n0.936 \u00b1 0.004 | \n0.938 \u00b1 0.002 | \n2th | \n
| Bioavailability_Ma | \n640 | \nClassification | \nAUROC | \n4.32 | \n0.735 \u00b1 0.016 | \n0.753 \u00b1 0.000 | \n5th | \n
| Bbb_Martins | \n2030 | \nClassification | \nAUROC | \n35.51 | \n0.917 \u00b1 0.003 | \n0.920 \u00b1 0.006 | \n2nd | \n
| Vdss_Lombardo | \n1130 | \nRegression | \nSpearman | \n12.91 | \n0.693 \u00b1 0.004 | \n0.713 \u00b1 0.007 | \n3rd | \n
| Cyp2D6_Substrate_Carbonmangels | \n667 | \nClassification | \nAUPRC | \n4.82 | \n0.714 \u00b1 0.009 | \n0.736 | \n6th | \n
| Cyp3A4_Substrate_Carbonmangels | \n670 | \nClassification | \nAUROC | \n3.97 | \n0.641 \u00b1 0.004 | \n0.667 \u00b1 0.019 | \n7th | \n
| Cyp2C9_Substrate_Carbonmangels | \n669 | \nClassification | \nAUPRC | \n4.28 | \n0.400 \u00b1 0.013 | \n0.441 \u00b1 0.033 | \n10th | \n
| Half_Life_Obach | \n667 | \nRegression | \nSpearman | \n4.20 | \n0.546 \u00b1 0.013 | \n0.576 \u00b1 0.025 | \n6th | \n
| Clearance_Microsome_Az | \n1102 | \nRegression | \nSpearman | \n12.71 | \n0.632 \u00b1 0.006 | \n0.630 \u00b1 0.010 | \n1st | \n
| Clearance_Hepatocyte_Az | \n1213 | \nRegression | \nSpearman | \n11.27 | \n0.391 \u00b1 0.004 | \n0.536 \u00b1 0.02 | \n>10th | \n
| Herg | \n655 | \nClassification | \nAUROC | \n3.54 | \n0.850 \u00b1 0.002 | \n0.880 \u00b1 0.002 | \n6th | \n
| Dili | \n475 | \nClassification | \nAUROC | \n1.92 | \n0.910 \u00b1 0.005 | \n0.925 \u00b1 0.005 | \n6th | \n
On some of the TDC datasets, TabPFN exhibits state-of-the-art performance out of the box. Regularly outperforming other DL methods such as graph neural networks (GNNs) designed specifically for cheminformatics. It can learn complex relationships between features and the target on small datasets due to the pre-training on synthetic data. TabPFN uses a scikit-learn style API, making it very easy to use and integrate into existing workflows. Since TabPFN approximates Bayesian inference, it can include important features such as uncertainty estimation possible.
\nTabPFN is limited to 10,000 entries with 500 features as a default, as this is the limit of the synthetic data it was trained on. It is possible to use datasets larger than this using the ignore_pretraining_limits=True parameter, but this may lead to very long training times. The computational requirements for TabPFN scale quadratically with the number of samples (n) and the number of features (m), i.e. O(n2 + m2) so training and predicting on larger datasets gets increasing longer. The training time can be long, though not significantly worse than other DL methods such as Chemprop in my experience. The single pass training and prediction architecture is not great for situations when you need to make multiple predictions, but the authors have included an option to cache the effect of the training data (using the fit_mode=\"fit_with_cache\" parameter) when training the model, making the training and prediction API much more like a standard scikit-learn ML model.
TabPFN is definitely a method to consider when building QSAR models, especially when you have a small dataset. With more time (and computing power), I'd be interested to see if combining the RDKit calculated descriptors with a molecular fingerprint (i.e. adding MACCS keys would make 376 feature columns) improves performance. In my experience sometimes combining the local description of molecules using molecular fingerprints with global descriptions of molecules using calculated descriptors creates a very good embedding of your molecules.
\nIt's possible to fine tune TabPFN with specialised datasets. Perhaps fine tuning with a large number of chemical datasets, or synthetic datasets of calculated molecular properties would create a chemical tabular foundation model (ChemTabPFN?), capable of accurately predicting chemical property-structure relationships.
\n![\"An](\"https://jonswain.github.io/images/traffic_widget/traffic_notification.png\"/)
No data science or cheminformatics today!
\nI usually cycle or get public transport to work, but occasionally I do have to drive. The traffic in T\u0101maki Makaurau (Auckland) is very variable, especially on rainy or windy days, so I often find myself checking Google Maps every 5 minutes after 4 pm to work out when I need to leave to get home. As with anything that I have to do repeatedly, I decided to automate it.
\nMy work computer runs Windows so I use WSL, which adds a few extra complications. All the code needed can be found in this repository.
\n(I've also been experimenting with Google Gemini as a coding assistant. This wasn't completely vibe-coding, but as this was a quick personal project, I was much less vigorous in checking the code it generated compared to more important work!)
\nAll the Python code is kept in a file called traffic-widget.py which is stored on my WSL disk, I first needed to import the necessary libraries:
import datetime\nimport json\nimport os\nimport subprocess\n\nimport requests\nfrom dotenv import load_dotenv\nGoogle maps doesn't seem to have a completely free API. They do offer some free usage, but you still have to sign up and give credit card details, which is always a worry in case you accidentally go over the free limits (or accidentally leak your API key to the internet). TomTom on the other hand does have a free API, you have to sign up and get an API key, and it comes with plenty of free requests. After signing up and getting an API key, I first created a TomTomAPI Python class, this uses the API key to make calls to the route calculation endpoint, providing the GPS coordinates for the start and end of your route, and returns the travel time.
\nclass TomTomAPI:\n \"\"\"Encapsulates interactions with the TomTom Routing API.\"\"\"\n\n BASE_URL = \"https://api.tomtom.com/routing/1/calculateRoute\"\n\n def __init__(self, api_key: str):\n \"\"\"Initialize the TomTomAPI with the API key.\"\"\"\n self.api_key = api_key\n\n def get_travel_time(\n self, start_lat: float, start_lon: float, end_lat: float, end_lon: float\n ) -> int | None:\n \"\"\"Calculate the travel time between two points.\n\n Args:\n start_lat (float): Latitude of the starting point.\n start_lon (float): Longitude of the starting point.\n end_lat (float): Latitude of the destination point.\n end_lon (float): Longitude of the destination point.\n\n Returns:\n int: The travel time in seconds, or None if an error occurs.\n \"\"\"\n start_point = f\"{start_lat},{start_lon}\"\n end_point = f\"{end_lat},{end_lon}\"\n url = f\"{self.BASE_URL}/{start_point}:{end_point}/json?key={self.api_key}\"\n\n try:\n response = requests.get(url)\n response.raise_for_status()\n data = response.json()\n travel_time_seconds = data[\"routes\"][0][\"summary\"][\"travelTimeInSeconds\"]\n return travel_time_seconds\n\n except requests.exceptions.RequestException as e:\n print(f\"TomTom API Error: {e}\")\n return None\n except (json.JSONDecodeError, KeyError) as e:\n print(f\"TomTom API Error: Invalid response format or missing data: {e}\")\n return None\nTo display the notifications I used BurntToast, a Windows PowerShell module for displaying Toast Notifications. To install, I needed to run PowerShell as admin and enter:
\nInstall-Module -Name BurntToast\nBurntToast can be called from a PowerShell file (.ps1). The PowerShell file to display notifications is called show_notification.ps1 and is stored on my Windows disk. The -Sound Alarm5 adds a sound to the notification and makes it last longer and -AppLogo gives an icon beside the notification.
param(\n [string]$Title,\n [string]$Message\n)\n\n$ImagePath = 'C:\\Path\\to\\icon.jpg'\n\nNew-BurntToastNotification -Text $Title, $Message -AppLogo $ImagePath -Sound Alarm5\nTo allow the script to run, it may be necessary to run this command in PowerShell as admin:
\nSet-ExecutionPolicy -ExecutionPolicy RemoteSigned -Scope CurrentUser\nBack in the Python file (traffic-widget.py), I next needed a class to run the PowerShell file to display the notification.
class WindowsNotifier:\n \"\"\"Handles displaying Windows toast notifications.\"\"\"\n\n def __init__(self, powershell_path: str, powershell_script: str):\n \"\"\"Initialize the WindowsNotifier with the paths to PowerShell.\"\"\"\n self.powershell_path = powershell_path\n self.powershell_script = powershell_script\n\n def show_notification(self, title: str, message: str):\n \"\"\"Show a Windows toast notification.\n\n Args:\n title (str): The title of the notification.\n message (str): The body of the notification.\n \"\"\"\n try:\n subprocess.run(\n [\n self.powershell_path,\n \"-ExecutionPolicy\",\n \"Bypass\",\n \"-File\",\n self.powershell_script,\n \"-Title\",\n title,\n \"-Message\",\n message,\n ],\n check=True,\n capture_output=True,\n text=True,\n )\n except subprocess.CalledProcessError as e:\n print(f\"Error showing notification: {e}\")\n print(f\"PowerShell Output:\\n{e.stderr}\")\n except FileNotFoundError:\n print(f\"Error: PowerShell executable or script not found. Check the paths.\")\nThis function formats the data from the TomTomAPI to be clearer.
\ndef format_travel_time(label: str, travel_time: int | None) -> str:\n \"\"\"Format the travel time into a readable string.\n\n Args:\n label (str): The label for the travel time (e.g., \"Home\").\n travel_time (int | None): The travel time in seconds, or None if an error occurred.\n\n Returns:\n str: The formatted travel time string.\n \"\"\"\n now = datetime.datetime.now()\n arrival_time = (\n now + datetime.timedelta(seconds=travel_time) if travel_time else None\n )\n if travel_time is not None:\n return f\"{label}: {travel_time / 60:.1f} minutes (arrive at {arrival_time.strftime('%H:%M')})\"\n else:\n return f\"{label}: ERROR minutes.\"\nTo prevent sharing sensitive data such as my API key and home address, I stored these in a .env file. This also stores the paths to Windows PowerShell and the PowerShell script. Since this is working between two operating systems, the paths are slightly more complicated than usual. The PowerShell path is the path to the powershell.exe file on your Windows disk, from your Linux environment. This is usually something like: /mnt/c/Windows/System32/WindowsPowerShell/v1.0/powershell.exe. The PowerShell script is in Windows format with escaped backslashes.
API_KEY = \"your_api_key\"\nWORK_LATITUDE = \"work_latitude\"\nWORK_LONGITUDE = \"work_longitude\"\nHOME_LATITUDE = \"home_latitude\"\nHOME_LONGITUDE = \"home_longitude\"\n\nPOWERSHELL_PATH = \"/mnt/c/path/to/powershell.exe\"\nPOWERSHELL_SCRIPT = \"C:\\\\windows\\\\path\\\\to\\\\show_notification.ps1\"\nTo run the Python script I used a conda environment. It can be created and activated with:
\nconda env create -f environment.yml\nconda activate traffic-widget\nFinally I needed some Python code to run the whole process:
\nif __name__ == \"__main__\":\n load_dotenv()\n api_key = os.environ[\"API_KEY\"]\n work_lat = float(os.environ[\"WORK_LATITUDE\"])\n work_lon = float(os.environ[\"WORK_LONGITUDE\"])\n home_lat = float(os.environ[\"HOME_LATITUDE\"])\n home_lon = float(os.environ[\"HOME_LONGITUDE\"])\n powershell_path = os.environ.get(\"POWERSHELL_PATH\")\n powershell_script = os.environ.get(\"POWERSHELL_SCRIPT\")\n if not powershell_script:\n print(\"Error: POWERSHELL_SCRIPT environment variable not set.\")\n exit(1)\n\n tomtom_api = TomTomAPI(api_key)\n notifier = WindowsNotifier(powershell_path, powershell_script)\n\n home_travel_time = tomtom_api.get_travel_time(\n work_lat, work_lon, home_lat, home_lon\n )\n\n message_lines = [\n format_travel_time(\"Home\", home_travel_time),\n ]\n notification_message = \"\\n\".join(message_lines)\n\n notifier.show_notification(\"Driving times:\", notification_message)\nManually running the script with python traffic-widget.py should cause the following pop-up:
To make the script run automatically, back in WSL, I ran:
\ncrontab -e\nAnd added the details for the cron job. I wanted mine to run every 5 minutes from 4-5 pm on weekdays.
\n0-55/5 16 * * 1-5 /home//miniconda3/envs/traffic-widget/bin/python /home//path/to/traffic-widget/traffic-widget.py\n An unwelcome cyclone may have cancelled my hiking plans, but it gave me the opportunity to delve into the fascinating world of drug repurposing. This blog post explores how identifying new uses for existing drugs can dramatically accelerate and reduce the cost of bringing vital treatments to patients, and how recent advances in machine learning can further streamline this process.
\nOver the Easter long weekend, T\u0101maki Makaurau (Auckland) was visited by Cyclone Tam, which meant my hiking trip was cancelled. At least it gave me some time to work on a few projects and finish this blog post I'd been meaning to do for a while!
\nTraditional drug discovery is a marathon, often taking over a decade and requiring an investment of more than a billion dollars for each successful drug. This lengthy and costly process highlights the urgent need for more efficient strategies to combat disease. One such promising strategy is drug repurposing: identifying new therapeutic uses for existing drugs or those already undergoing clinical trials. By bypassing much of the early-stage development, and sometimes significant parts of clinical trials, drug repurposing offers the potential to significantly reduce both the time and financial cost associated with bringing new treatments to patients. Drug repurposing is ideal for addressing therapeutic needs in disease areas where financial incentives are limited (such as neglected tropical diseases and antibiotic development), or for rapid deployment in situations demanding swift treatment interventions for imminent public health emergencies, such as during the COVID-19 pandemic.
\nOne well-known example of drug repurposing is sildenafil, originally developed to treat high blood pressure and angina (chest pain), and famously repurposed for erectile dysfunction after clinical trials revealed this unexpected side effect. Another is thalidomide, developed to treat morning sickness and found to cause birth defects, which is now used to treat cancers and skin disorders associated with leprosy. The crucial role of drug repurposing was further highlighted during the recent COVID-19 pandemic, where the corticosteroid dexamethasone emerged as a vital intervention in reducing mortality among severely ill patients.
\nMachine learning offers a novel method for drug repurposing. Using high-quality activity data for a target of interest, a machine learning model can be trained to predict active compounds based on their molecular structure. This trained model can then be used to screen libraries of compounds for potential repurposing.
\nThe Drug Repurposing Hub is a curated and annotated collection of drugs that have already received FDA approval or are currently undergoing clinical trials. Maintained by the Broad Institute, this resource serves as both a virtual library, providing comprehensive information on each compound, and a physical library, offering access to compound plates for experimental screening.
\nThis post was inspired by the work of the Collins Lab at MIT, who used deep learning techniques to identify potential compounds for repurposing as treatments for antimicrobial-resistant bacteria. Their research addresses the escalating crisis of antibiotic resistance, a looming threat projected to cause 10 million deaths annually by 2050.
\nTo tackle this challenge, the researchers assembled a primary training set of 2,335 molecules, experimentally screening them against E. coli BW25113 and identifying 120 compounds with antimicrobial activity. With the goal of creating a robust and generalizable model, they ensured the training set was structurally diverse. This data was then used to train a directed message-passing neural network (dMPNN) to create a classification model capable of predicting antimicrobial activity based on a molecule's chemical structure. The trained model was subsequently used to make predictions on virtual libraries of compounds, including the Drug Repurposing Hub.
\nNinety-nine molecules from the Drug Repurposing Hub that were predicted to be active were experimentally tested, with 51 displaying growth inhibition against E. coli. This process led to the identification of a particularly promising compound, the c-Jun N-terminal kinase inhibitor SU3327, which the researchers named halicin. Notably, halicin possesses a structure distinctly different from conventional antibiotics and demonstrated potent inhibitory activity against E. coli growth.
\nCentral to this research is the open-source software package Chemprop. Developed primarily by researchers at MIT using the PyTorch framework, Chemprop harnesses the capabilities of message-passing neural networks (MPNNs) for molecular property prediction. Chemprop can be used both as a command-line interface (CLI) and a Python API. It implements a d-MPNN architecture, a type of graph neural network particularly effective for extracting relationships from molecular structures represented as graphs. An advantage of Chemprop is its ability to learn directly from molecular representations, such as SMILES strings converted into molecular graphs, eliminating the need for manual feature engineering, which can also introduce biases.
\nLacking access to a physical lab for experimental data, I needed an alternative approach for data gathering. For this, I used ChEMBL, a freely accessible and meticulously curated chemical database housing bioactivity data for drug-like molecules. It is maintained by the European Bioinformatics Institute (EBI), part of the European Molecular Biology Laboratory (EMBL) in the UK. Using the ChEMBL Webservice, I was able to specifically download data relevant to the target of interest. This retrieved data then served as the foundational dataset for training a machine learning model, which was subsequently employed to virtually screen the compounds within The Drug Repurposing Hub.
\nimport subprocess\nfrom pathlib import Path\n\nimport pandas as pd\nfrom chembl_webresource_client.new_client import new_client\nfrom rdkit import Chem\nfrom rdkit.Chem import Draw\nfrom tqdm import tqdm\n\nPath(\"data\").mkdir(exist_ok=True)\nPath(\"models\").mkdir(exist_ok=True)\nUsing the ChEMBL Webservice, I downloaded the data associated with a target for training a machine learning model using the three functions below. The first searches ChEMBL with a target name, and returns the target ChEMBL ID of the most similar target. The second downloads the activity data associated that the target ChEMBL ID, along with the molecular ChEMBL IDs. The third returns the canonical SMILES for the molecular ChEMBL IDs.
\ntarget = new_client.target\nactivity = new_client.activity\nmolecule = new_client.molecule\n\n\ndef get_target_chembl_id(target_name: str) -> str:\n \"\"\"Get the ChEMBL ID of a target given its name.\"\"\"\n print(f'Searching for target: \"{target_name}\"')\n target_query = target.search(target_name)\n name, chembl_id = target_query[0][\"pref_name\"], target_query[0][\"target_chembl_id\"]\n print(f\"Target name: {name}, ChEMBL ID: {chembl_id}\")\n return chembl_id\n\n\ndef get_target_activity(chembl_id: str) -> pd.DataFrame:\n \"\"\"Get the activity data for a target given its ChEMBL ID.\"\"\"\n print(f\"Getting activity data for {chembl_id}\")\n activities = activity.filter(\n target_chembl_id=chembl_id,\n standard_type=\"IC50\",\n relation=\"=\",\n standard_units=\"nM\",\n ).only(\n \"molecule_chembl_id\",\n \"standard_value\",\n )\n print(f\"Found {len(activities)} activities\")\n activities_df = pd.DataFrame(\n list(tqdm(activities[:100], desc=\"Processing Activities\"))\n )[[\"molecule_chembl_id\", \"standard_value\"]]\n activities_df[\"standard_value\"] = pd.to_numeric(activities_df[\"standard_value\"])\n return activities_df\n\n\ndef get_molecule_data(molecule_chembl_ids: list[str]) -> pd.DataFrame:\n \"\"\"Get the molecule data for a list of ChEMBL IDs.\"\"\"\n print(f\"Getting molecule data for {len(molecule_chembl_ids)} ChEMBL IDs\")\n compounds_provider = molecule.filter(\n molecule_chembl_id__in=molecule_chembl_ids\n ).only(\"molecule_chembl_id\", \"molecule_structures\")\n compounds_df = pd.json_normalize(\n list(tqdm(compounds_provider, desc=\"Processing Compounds\"))\n )[[\"molecule_chembl_id\", \"molecule_structures.canonical_smiles\"]].rename(\n columns={\"molecule_structures.canonical_smiles\": \"canonical_smiles\"}\n )\n return compounds_df\nChemprop takes the data in the form of a table containing one column of SMILES strings, and one columns of activities. For this example I am creating a classification model, so the activity column is a binary classification (1 or 0). The first function takes the activity data and molecule data and merges them on the molecule ChEMBL ID, before cleaning the data up and creating the binary variable based on a cutoff. This cutoff can be changed depending on the target and dataset. I've always found the CLI usage of Chemprop easier, so the data is saved as a CSV file for training.
\ndef create_training_data(\n activity_data: pd.DataFrame, molecule_data: pd.DataFrame, nm_cutoff: float\n) -> pd.DataFrame:\n \"\"\"Create the training data for the model.\"\"\"\n print(\"Creating training data\")\n training_data = (\n (\n pd.merge(\n left=activity_data,\n right=molecule_data,\n on=\"molecule_chembl_id\",\n how=\"left\",\n )\n .dropna(how=\"any\")\n .drop(columns=\"molecule_chembl_id\")\n )\n .groupby(\"canonical_smiles\", as_index=False)\n .agg(\"mean\")\n )\n training_data[\"activity\"] = (training_data[\"standard_value\"] < nm_cutoff).astype(\n int\n )\n actives = training_data[\"activity\"].sum()\n inactives = len(training_data) - actives\n percentage = 100 * actives / len(training_data)\n print(f\"Actives: {actives} ({percentage:.2f}%), Inactives: {inactives}\")\n return training_data.drop(columns=\"standard_value\")\n\n\ndef save_training_data(\n training_data: pd.DataFrame, target_name: str, output_dir: Path\n) -> None:\n \"\"\"Save the training data to a CSV file.\"\"\"\n print(f\"Saving training data to {output_dir}\")\n output_dir.mkdir(parents=True, exist_ok=True)\n training_data.to_csv(output_dir / f\"{target_name}_training_data.csv\", index=False)\n print(f\"Training data saved to {output_dir / f'{target_name}_training_data.csv'}\")\nThis function downloaded the Drug Repurposing Hub dataset, tidies the data up and saves it for making predictions using Chemprop.
\ndef canon_smiles(smiles: str) -> str | None:\n \"\"\"Convert SMILES string to canonical SMILES.\"\"\"\n mol = Chem.MolFromSmiles(smiles)\n if mol is None:\n return None\n return Chem.MolToSmiles(mol, canonical=True)\n\n\ndef download_drug_repurposing_hub() -> None:\n \"\"\"Download the Drug Repurposing Hub dataset.\"\"\"\n if Path(\"data/drug_repurposing_hub.csv\").exists():\n print(\"Drug Repurposing Hub dataset already downloaded\")\n return None\n url = \"https://storage.googleapis.com/cdot-general-storage/repurposing_samples_20240610.txt\"\n drug_repurposing_df = pd.read_csv(\n url,\n sep=\"\\t\",\n skiprows=9,\n )\n drug_repurposing_df[\"smiles\"] = drug_repurposing_df[\"smiles\"].apply(canon_smiles)\n drug_repurposing_df.dropna(subset=[\"smiles\"], inplace=True)\n drug_repurposing_df.drop_duplicates(subset=[\"smiles\"], keep=\"first\", inplace=True)\n drug_repurposing_df.to_csv(\"data/drug_repurposing_hub.csv\", index=False)\n print(\"Drug Repurposing Hub dataset downloaded and processed\")\n print(f\"Number of unique SMILES: {len(drug_repurposing_df['smiles'].unique())}\")\nAs I mentioned above, I find the CLI usage of Chemprop easier (especially since the move from v1 to v2), so here I'm using subprocess to run Chemprop. These two functions train a Chemprop model, and make predictions on the Drug Repurposing Hub data.
\ndef train_chemprop_model(target_name: str, data_dir: str, output_dir: str) -> None:\n \"\"\"Train a ChemProp model using the training data.\"\"\"\n print(f\"Training ChemProp model for {target_name}\")\n chemprop_train_args = [\n \"chemprop\",\n \"train\",\n \"--data-path\",\n f\"{data_dir}/{target_name.replace(\" \", \"_\")}_training_data.csv\",\n \"--task-type\",\n \"classification\",\n \"--save-dir\",\n f\"{output_dir}/{target_name.replace(\" \", \"_\")}\",\n \"--split-type\",\n \"scaffold_balanced\",\n ]\n subprocess.run(chemprop_train_args, check=True)\n\n\ndef predict_with_chemprop(target_name: str, preds_dir: str, model_dir: str) -> None:\n \"\"\"Predict using the ChemProp model.\"\"\"\n print(f\"Predicting with ChemProp model for {target_name}\")\n chemprop_predict_args = [\n \"chemprop\",\n \"predict\",\n \"--test-path\",\n \"data/drug_repurposing_hub.csv\",\n \"--model-paths\",\n f\"{model_dir}/{target_name.replace(\" \", \"_\")}\",\n \"--smiles-columns\",\n \"smiles\",\n \"--preds-path\",\n f\"{preds_dir}/{target_name.replace(\" \", \"_\")}_predictions.csv\",\n ]\n subprocess.run(chemprop_predict_args, check=True)\nOnce we've scored the compounds, we want to take a look at the top scoring compounds.
\ndef visualize_top_predictions(search: str) -> None:\n \"\"\"Visualize the top predictions.\"\"\"\n predictions_df = (\n pd.read_csv(f\"data/{search.replace(\" \", \"_\")}_predictions.csv\")\n .sort_values(\"activity\", ascending=False)\n .head(9)\n )\n return Draw.MolsToGridImage(\n [Chem.MolFromSmiles(s) for s in predictions_df[\"smiles\"]],\n legends=[f\"{a:.2f}\" for a in predictions_df[\"activity\"]],\n molsPerRow=3,\n )\nThis function runs the whole process for a specified target search.
\ndef run_repurposing_pipeline(search: str, nm_cutoff=25) -> None:\n \"\"\"Run the drug repurposing pipeline.\"\"\"\n download_drug_repurposing_hub()\n target_chembl_id = get_target_chembl_id(search)\n activity_data = get_target_activity(target_chembl_id)\n molecule_data = get_molecule_data(activity_data[\"molecule_chembl_id\"].tolist())\n training_data = create_training_data(\n activity_data, molecule_data, nm_cutoff=nm_cutoff\n )\n save_training_data(training_data, search.replace(\" \", \"_\"), Path(\"data\"))\n train_chemprop_model(target_name=search, data_dir=\"data\", output_dir=\"models/\")\n predict_with_chemprop(target_name=search, preds_dir=\"data\", model_dir=\"models/\")\nMalaria kills around half a million people per year, with 76% of global malaria deaths in children under 5 years old. It is considered a neglected tropical disease (NTD) due to its significant impact on global health, particularly in low-income countries. It disproportionately affects some of the world's poorest people, and is often overlooked in global health efforts. Drug repurposing could offer alternative therapies significantly cheaper than developing an entirely new drug. Here I look to target the most deadly species of parasite that carries malaria - Plasmodium falciparum. After training the Chemprop model and screening the Drug Repurposing Hub, some compounds are flagged as potential hits and could be investigated experimentally.
\nOf the top 10 predictions, pyronaridine and KAF-156 are both malaria treatments, and were both in the training data, but interestingly pyronaridine was not classified as a hit in the training data.\nCarfilzomib, TPPS4, nanchangmycin, oprozomib, and lindane have all been investigated for treating malaria, so the model appears to be picking out the right compounds. Abemaciclib appears to have been recently investigated for repurposing with malaria. GS-967 appears to have no link to malaria yet, but maybe should be investigated experimentally!
\nsearch = \"Plasmodium falciparum\"\n\nrun_repurposing_pipeline(search)\nvisualize_top_predictions(search)\nDrug Repurposing Hub dataset downloaded and processed\nNumber of unique SMILES: 6734\nSearching for target: \"Plasmodium falciparum\"\nTarget name: Plasmodium falciparum, ChEMBL ID: CHEMBL364\nGetting activity data for CHEMBL364\nFound 45314 activities\nProcessing Activities: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 45314/45314 [00:02<00:00, 15928.83it/s]\nGetting molecule data for 45314 ChEMBL IDs\nProcessing Molecules: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 21821/21821 [1:32:46<00:00, 3.92it/s]\nCreating training data\nActives: 2677 (12.30%), Inactives: 19082\nSaving training data to data\nTraining data saved to data/Plasmodium_falciparum_training_data.csv\nTraining ChemProp model for Plasmodium falciparum\nPredicting with ChemProp model for Plasmodium falciparum\n![\"Grid](\"https://jonswain.github.io/images/drug_repurposing/malaria-hits.png\"/)
The COVID-19 pandemic has caused the deaths of over 7 million people so far, and in 2020 more than half of the world's population was under some form of lockdown. With the speed that COVID spread around the globe, there was no time to develop new drug from scratch, and drug repurposing was explored in clinics to find treatments. In this example, I didn't find any promising compounds in the Drug Repurposing Hub, but with a much smaller dataset than for malaria, the machine learning model is potentially much less powerfull.
\nsearch = \"SARS-CoV-2\"\n\nrun_repurposing_pipeline(search, nm_cutoff=500)\nvisualize_top_predictions(search)\nDrug Repurposing Hub dataset already downloaded\nSearching for target: \"SARS-CoV-2\"\nTarget name: SARS-CoV-2, ChEMBL ID: CHEMBL4303835\nGetting activity data for CHEMBL4303835\nFound 802 activities\nProcessing Activities: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 802/802 [00:00<00:00, 12406.99it/s]\nGetting molecule data for 802 ChEMBL IDs\nProcessing Compounds: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 490/490 [00:00<00:00, 6081.22it/s]\nCreating training data\nActives: 71 (14.55%), Inactives: 417\nSaving training data to data\nTraining data saved to data/SARS-CoV-2_training_data.csv\nTraining ChemProp model for SARS-CoV-2\nPredicting with ChemProp model for SARS-CoV-2\n![\"Grid](\"https://jonswain.github.io/images/drug_repurposing/covid-hits.png\"/)
Whilst neural networks (NNs) have done amazing things with unstructured data such as text and images, they've traditionally been outperformed on tabular data by Gradient-Boosted Decision Trees (GBDTs), although recent advances such as TabPFN and TabICL suggest that the performance gap may have closed. I've written a bit about TabPFN for chemical datasets here.
\nOne of the most popular deep learning architectures for chemical property prediction is Chemprop, developed by researchers at MIT. A Chemprop model is made up of two NNs. The first is a directed Message Passing Neural Network (d-MPNN), a type of Graph Convolutional Neural Network (GCNN), that takes a graphical representation of a molecule and converts it to a molecular embedding, a vector that describes the original molecule. This molecular embedding is then put through a Feed-Forward Neural Network (FFN), a fully connected NN used to make the final prediction. The model is trained end-to-end, with both NNs being updated at the same time to minimise a loss function.
\nOnce the model is trained, the d-MPNN can be used to calculate a molecular embedding for a collection of molecules. This is a learned representation optimised for the specific chemical prediction problem. This learned representation is a vector with a length equal to the number of input nodes in the FFN. This is very similar to a traditional molecular fingerprint, such as a Morgan fingerprint.
\nA collection of these learned fingerprints forms a 2-dimensional array, which is tabular data. Since GBDTs regularly outperform NNs on tabular data, can we improve the performance of Chemprop by replacing the FFN with a GBDT after training the d-MPNN?
\nWhilst writing this post, a paper was published in The Journal of Chemical Information and Modeling by Pat Walters et al., Practically Significant Method Comparison Protocols for Machine Learning in Small Molecule Drug Discovery. This is probably my favourite cheminformatics paper so far this year, with comprehensive guidelines on how to compare machine learning models in a statistically robust way. I've tried to follow these recommendations when comparing models here.
\nBelow I compare three models: the first is a Random Forest (this is not gradient-boosted, but provides a good baseline without hyperparameter tuning) trained on Morgan fingerprints; the second is a standard Chemprop model; and the third generates learned fingerprints using the Chemprop d-MPNN and then trains a Random Forest model on these. I used two datasets, one for a classification problem and one for a regression problem. The classification data is the Bile Salt Export Pump dataset and the regression data is the Polaris LogS dataset, both ideas 'borrowed' from Practical Cheminformatics posts.
\nimport subprocess\nfrom collections import defaultdict\nfrom pathlib import Path\nfrom typing import Callable\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom numpy.typing import NDArray\nfrom rdkit import Chem\nfrom rdkit.Chem import rdFingerprintGenerator\nfrom scipy.stats import f_oneway\nfrom sklearn.base import BaseEstimator\nfrom sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\nfrom sklearn.metrics import (\n matthews_corrcoef,\n mean_absolute_error,\n r2_score,\n root_mean_squared_error,\n)\nfrom sklearn.model_selection import ShuffleSplit\nfrom statsmodels.stats.multicomp import pairwise_tukeyhsd\nfrom tqdm.notebook import tqdm\n\nfrom visualisation import confusion_matricies, scatter_plots\nThe first function prepare_data takes the path to a CSV file, creates RDKit molecular objects, removes any that are causing issues, and calculates the Morgan fingerprints for each model.
The second, save_data_for_chemprop, formats the data and saves it for training Chemprop models.
The next three functions: train_chemprop_model, make_chemprop_predictions, and calculate_chemprop_fingerprints, all use subprocess to call the Chemprop CLI. I've generally found it easier to use than the Chemprop Python module.
The next two: train_morgan_sklearn_model and train_chemprop_sklearn_model, use Scikit-learn to train Random Forest models, the first using the previously generated Morgan fingerprints, and the second uses the Chemprop learned fingerprints.
Finally, gather_predictions uses all of the trained models to make predictions on a test dataset.
def prepare_data(\n data_path: Path, smiles_col: str, target_col: str\n) -> tuple[pd.DataFrame, pd.DataFrame]:\n \"\"\"Load and clean data, generate Morgan fingerprints.\n\n Args:\n data_path (Path): Path to CSV file containing data.\n smiles_col (str): Name of the column containing SMILES strings.\n target_col (str): Name of the column containing target variable.\n\n Returns:\n tuple[pd.DataFrame, pd.DataFrame]: Cleaned dataframe and dataframe of Morgan\n fingerprints.\n \"\"\"\n df = pd.read_csv(data_path)\n df[\"ROMol\"] = df[smiles_col].apply(lambda x: Chem.MolFromSmiles(x))\n df.dropna(subset=[\"ROMol\"], inplace=True)\n mfpgen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)\n fps = [mfpgen.GetFingerprint(mol) for mol in df[\"ROMol\"]]\n fps_df = pd.DataFrame(np.array(fps))\n df = df[[smiles_col, target_col]]\n return df, fps_df\n\n\ndef save_data_for_chemprop(\n data: pd.DataFrame,\n project_name: str,\n train_idx: NDArray[np.int64],\n test_idx: NDArray[np.int64],\n seed: int,\n fold_number: int,\n) -> None:\n \"\"\"Save training and test data for Chemprop CLI.\n\n Args:\n data (pd.DataFrame): Dataframe containing data.\n project_name (str): Name of the project.\n train_idx (NDArray[np.int64]): List of indices for training data.\n test_idx (NDArray[np.int64]): List of indices for test data.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n \"\"\"\n Path(f\"data/{project_name}/seed_{seed}/fold_{fold_number}/\").mkdir(\n parents=True, exist_ok=True\n )\n data.iloc[train_idx].to_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/train.csv\", index=False\n )\n data.iloc[test_idx].to_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/test.csv\", index=False\n )\n\n\ndef train_chemprop_model(\n project_name: str, seed: int, fold_number: int, task_type: str\n) -> None:\n \"\"\"Train Chemprop model using CLI.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n task_type (str): Type of task ('classification' or 'regression').\n \"\"\"\n subprocess.run(\n [\n \"chemprop\",\n \"train\",\n \"--data-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/train.csv\",\n \"--task-type\",\n task_type,\n \"--output-dir\",\n f\"models/{project_name}/seed_{seed}/fold_{fold_number}\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n stderr=subprocess.DEVNULL,\n )\n\n\ndef make_chemprop_predictions(project_name: str, seed: int, fold_number: int) -> None:\n \"\"\"Make predictions using trained Chemprop model.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n \"\"\"\n subprocess.run(\n [\n \"chemprop\",\n \"predict\",\n \"--test-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/test.csv\",\n \"--model-paths\",\n f\"models/{project_name}/seed_{seed}/fold_{fold_number}/\",\n \"--preds-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/preds.csv\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n stderr=subprocess.DEVNULL,\n )\n\n\ndef calculate_chemprop_fingerprints(\n project_name: str, seed: int, fold_number: int, dataset_type: str\n) -> None:\n \"\"\"Calculate Chemprop fingerprints using trained Chemprop model.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n dataset_type (str): Type of dataset ('train' or 'test').\n \"\"\"\n subprocess.run(\n [\n \"chemprop\",\n \"fingerprint\",\n \"--test-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/{dataset_type}.csv\",\n \"--output\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/{dataset_type}_fps.csv\",\n \"--model-path\",\n f\"models/{project_name}/seed_{seed}/fold_{fold_number}/\",\n \"--ffn-block-index\",\n \"0\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n stderr=subprocess.DEVNULL,\n )\n\n\ndef train_morgan_sklearn_model(\n fps_df: pd.DataFrame, y: pd.Series, train_idx: NDArray[np.int64], task_type: str\n) -> BaseEstimator:\n \"\"\"Train Scikit-learn model using Morgan fingerprints.\n\n Args:\n fps_df (pd.DataFrame): Dataframe of Morgan fingerprints.\n y (pd.Series): Series of target variable.\n train_idx (NDArray[np.int64]): List of indices for training data.\n task_type (str): Type of task ('classification' or 'regression').\n\n Returns:\n BaseEstimator: Trained Scikit-learn model.\n \"\"\"\n if task_type == \"regression\":\n rf = RandomForestRegressor(n_estimators=100, random_state=0)\n else:\n rf = RandomForestClassifier(n_estimators=100, random_state=0)\n X_train = fps_df.iloc[train_idx]\n y_train = y.iloc[train_idx]\n return rf.fit(X_train, y_train)\n\n\ndef train_chemprop_sklearn_model(\n project_name: str,\n seed: int,\n fold_number: int,\n y: pd.Series,\n train_idx: NDArray[np.int64],\n task_type: str,\n) -> BaseEstimator:\n \"\"\"Train Scikit-learn model using Chemprop fingerprints.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n y (pd.Series): Series of target variable.\n train_idx (NDArray[np.int64]): List of indices for training data.\n task_type (str): Type of task ('classification' or 'regression').\n\n Returns:\n BaseEstimator: Trained Scikit-learn model.\n \"\"\"\n X_train = pd.read_csv(\n Path(f\"data/{project_name}/seed_{seed}/fold_{fold_number}/train_fps_0.csv\")\n )\n if task_type == \"regression\":\n rf = RandomForestRegressor(n_estimators=100, random_state=0)\n else:\n rf = RandomForestClassifier(n_estimators=100, random_state=0)\n y_train = y.iloc[train_idx]\n return rf.fit(X_train, y_train)\n\n\ndef gather_predictions(\n morgan_rf: BaseEstimator,\n chemprop_rf: BaseEstimator,\n project_name: str,\n seed: int,\n fold_number: int,\n df: pd.DataFrame,\n fps_df: pd.DataFrame,\n test_idx: NDArray[np.int64],\n target_col: str,\n task_type: str,\n) -> dict:\n \"\"\"Gather predictions from all models.\n\n Args:\n morgan_rf (BaseEstimator): Trained Scikit-learn model on Morgan fingerprints.\n chemprop_rf (BaseEstimator): Trained Scikit-learn model on Chemprop fingerprints.\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n df (pd.DataFrame): Dataframe containing data.\n fps_df (pd.DataFrame): Dataframe of Morgan fingerprints.\n test_idx (NDArray[np.int64]): List of indices for test data.\n target_col (str): Name of the column containing target variable.\n task_type (str): Type of task ('classification' or 'regression').\n\n Returns:\n dict: Updated dictionary with predictions.\n \"\"\"\n predictions = {}\n rf_preds = morgan_rf.predict(fps_df.iloc[test_idx])\n chemprop_rf_preds = chemprop_rf.predict(\n pd.read_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/test_fps_0.csv\"\n )\n )\n chemprop_probas = pd.read_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/preds.csv\"\n )[target_col]\n if task_type == \"regression\":\n chemprop_preds = chemprop_probas\n else:\n chemprop_preds = (chemprop_probas >= 0.5).astype(int)\n y_true = df.iloc[test_idx][target_col]\n predictions[f\"seed_{seed}_fold_{fold_number}\"] = {\n \"rf\": rf_preds,\n \"chemprop_rf\": chemprop_rf_preds,\n \"chemprop\": chemprop_preds,\n \"true\": y_true,\n }\n return predictions\nTo compare the models, I used 5x5 cross-validation. This approach involves using five different random seeds to generate data splits for 5-fold cross-validation, resulting in 25 unique train-test splits. For each split, the models are trained on the training data and then used to make predictions on the test data.
\ndef run_model_comparison(\n data_path: Path,\n project_name: str,\n smiles_col: str,\n target_col: str,\n task_type: str,\n n_seeds: int = 5,\n n_folds: int = 5,\n test_frac: float = 0.2,\n) -> dict:\n \"\"\"Run model comparison pipeline.\n\n Args:\n data_path (Path): Path to CSV file containing data.\n project_name (str): Name of the project.\n smiles_col (str): Name of the column containing SMILES strings.\n target_col (str): Name of the column containing target variable.\n task_type (str): Type of task ('classification' or 'regression').\n n_seeds (int, optional): Number of random seeds for cross-validation. Defaults to 5.\n n_folds (int, optional): Number of folds for cross-validation. Defaults to 5.\n test_frac (float, optional): Fraction of data to use as test set. Defaults to 0.2.\n\n Returns:\n dict: Dictionary containing predictions from all models.\n \"\"\"\n df, fps_df = prepare_data(data_path, smiles_col, target_col)\n\n all_predictions = {}\n for seed in tqdm(range(n_seeds), total=n_seeds, desc=\"Seeds\", unit=\"seed\"):\n cv = ShuffleSplit(n_splits=n_folds, test_size=test_frac, random_state=seed)\n for fold_number, (train_idx, test_idx) in tqdm(\n enumerate(cv.split(df)),\n total=n_folds,\n desc=f\"Seed {seed}: Folds\",\n unit=\"fold\",\n ):\n # Save data for training models using Chemprop CLI\n save_data_for_chemprop(\n df, project_name, train_idx, test_idx, seed, fold_number\n )\n\n # Train Chemprop model using CLI\n train_chemprop_model(project_name, seed, fold_number, task_type)\n\n # Make predictions on test set using trained Chemprop model\n make_chemprop_predictions(project_name, seed, fold_number)\n\n # Use trained model to create fingerprints\n calculate_chemprop_fingerprints(project_name, seed, fold_number, \"train\")\n calculate_chemprop_fingerprints(project_name, seed, fold_number, \"test\")\n\n # Train Scikit-learn models on Morgan fingerprints\n morgan_rf = train_morgan_sklearn_model(\n fps_df, df[target_col], train_idx, task_type\n )\n\n # Train Scikit-learn model on Chemprop fingerprints\n chemprop_rf = train_chemprop_sklearn_model(\n project_name,\n seed,\n fold_number,\n df[target_col],\n train_idx,\n task_type,\n )\n\n # Gather predictions\n predictions = gather_predictions(\n morgan_rf,\n chemprop_rf,\n project_name,\n seed,\n fold_number,\n df,\n fps_df,\n test_idx,\n target_col,\n task_type,\n )\n all_predictions.update(predictions)\n return all_predictions\nThe next two cells run the training and prediction process for the classification and regression datasets.
\n# Set up workspace\nPath(\"data\").mkdir(exist_ok=True)\nPath(\"models\").mkdir(exist_ok=True)\nnp.random.seed(42)\n\nn_seeds = 5\nn_folds = 5\ntest_frac = 0.2\ndata_path = Path(\"data/raw/bsep.csv\")\nproject_name = \"bsep\"\nsmiles_col = \"smiles\"\ntarget_col = \"bsep\"\ntask_type = \"classification\"\n\nclassification_predictions = run_model_comparison(\n data_path,\n project_name,\n smiles_col,\n target_col,\n task_type,\n n_seeds=n_seeds,\n n_folds=n_folds,\n test_frac=test_frac,\n)\nSeeds: 0%| | 0/5 [00:00n_seeds = 5\nn_folds = 5\ntest_frac = 0.2\ndata_path = Path(\"data/raw/biogen_logS.csv\")\nproject_name = \"logS\"\nsmiles_col = \"SMILES\"\ntarget_col = \"logS\"\ntask_type = \"regression\"\n\nregression_predictions = run_model_comparison(\n data_path,\n project_name,\n smiles_col,\n target_col,\n task_type,\n n_seeds=n_seeds,\n n_folds=n_folds,\n test_frac=test_frac,\n)\nSeeds: 0%| | 0/5 [00:00The next function takes the predicted values and true values from the test set and calculates a relevant performance metric.
\nFor the moderately imbalanced BSEP classification data (17% positive class), I used the Matthews Correlation Coefficient (MCC). MCC is a value between -1 and 1 and is a robust metric for imbalanced datasets. For the LogS regression data, I used the coefficient of determination (R2). An R2 value of 1 indicates perfect agreement between predicted and true values, with no lower limit for poor predictions.
\ndef calculate_performance(preds_dict: dict, metric: Callable) -> dict:\n \"\"\"Analyse performance of models.\n\n Args:\n preds_dict (dict): Dictionary containing predictions from all models.\n metric (Callable): Metric function to evaluate performance.\n\n Returns:\n dict: Dictionary containing performance metrics for each model.\n \"\"\"\n metrics = defaultdict(list)\n for predictions in preds_dict.values():\n metrics[f\"rf_{metric.__name__}\"].append(\n metric(predictions[\"true\"], predictions[\"rf\"])\n )\n metrics[f\"chemprop_rf_{metric.__name__}\"].append(\n metric(predictions[\"true\"], predictions[\"chemprop_rf\"])\n )\n metrics[f\"chemprop_{metric.__name__}\"].append(\n metric(predictions[\"true\"], predictions[\"chemprop\"])\n )\n return metrics\n\n\nclassification_metrics = calculate_performance(\n classification_predictions, matthews_corrcoef\n)\n\nregression_metrics = calculate_performance(regression_predictions, r2_score)\nTo visualise and inspect the classification results, I've used confusion matricies. All three models have a significant number of True Negatives (TN) in the top left square. The Random Forest model and the Chemprop model have a significant number of False Negatives (FN), whereas the Chemprop-RF model has more True Positives (TP) and more False Positives (FP).
\nconfusion_matricies(classification_predictions)\n![\"Confusion](\"https://jonswain.github.io/images/chemprop-rf/chemprop-rf-1.png\"/)
To visualise and inspect the regression results, I've used x-y scatter plots of the predicted values against the true values. The Random Forest regressor has poor predictions, with more values being overestimated. The Chemprop model and the Chemprop-RF models perform better, especially on these lower LogS values. Some performance metrics are annotated on the plots in the top left: R2, RMSE (root mean squared error), and MAE (mean absolute error).
\nscatter_plots(regression_predictions)\n![\"Predicted](\"https://jonswain.github.io/images/chemprop-rf/chemprop-rf-2.png\"/)
As recommended in the paper I've linked, I used one-way ANOVA to compare the performance distributions from the 5x5 cross-validation and determine if there was a statistically significant difference between the means. This was followed by Tukey's Honest Significant Differences as a post-hoc pairwise test. The results are plotted using the Simultaneous Confidence Interval plot from the statsmodels Python library.
\nFor both the classification and regression datasets, the ANOVA results showed a statistically significant difference in model performance (p<0.05).
\nFor the classification dataset, post-hoc tests revealed a statistically significant difference in performance between all pairs of models, with the Chemprop-RF model performing the best.
\nFor the regression dataset, the Random Forest model performed statistically significantly worse than the other two. However, there was no statistically significant difference between the performance of the Chemprop and Chemprop-RF models.
\nThe classification dataset (BSEP, 807 entries) is much smaller than the regression dataset (LogS, 2,173 entries). This might suggest that Chemprop-RF models are particularly effective for low-data problems, where there isn't enough data to train a high-performing Feed-Forward Neural Network (FFN) alone.
\n# Adapted from https://github.com/PatWalters/practical_cheminformatics_posts/blob/main/adme_comparison/\ndef run_anova(df_in: pd.DataFrame, col: str) -> float:\n \"\"\"Run one-way ANOVA on model performance.\"\"\"\n res_list = []\n for _, v in df_in.groupby(\"method\"):\n res_list.append(v[col].values)\n return f_oneway(*res_list)[1]\n\n\ndef plot_metric_comparison(metric_dict: dict, metric: Callable) -> None:\n \"\"\"Plot comparison of model performance using a specified metric.\n\n Args:\n metric_dict (dict): Dictionary containing performance metrics for each model.\n metric (Callable): Metric function used to evaluate performance.\n \"\"\"\n _, ax = plt.subplots()\n melt_df = pd.DataFrame(metric_dict).melt()\n melt_df[\"method\"] = melt_df.variable.map(\n {\n f\"rf_{metric.__name__}\": \"Random Forest\",\n f\"chemprop_rf_{metric.__name__}\": \"Chemprop-RF\",\n f\"chemprop_{metric.__name__}\": \"Chemprop\",\n }\n )\n best_model = (\n melt_df.groupby(\"method\")[\"value\"]\n .mean()\n .reset_index()\n .sort_values(\"value\", ascending=False)\n .method.values[0]\n )\n tukey = pairwise_tukeyhsd(\n endog=melt_df[\"value\"], groups=melt_df[\"method\"], alpha=0.05\n )\n tukey.plot_simultaneous(comparison_name=best_model, ax=ax, figsize=(8, 5))\n anova_p_value = run_anova(melt_df, \"value\")\n ax.set_title(\n f\"{metric.__name__.replace('_', ' ').title()} \" + f\"ANOVA p={anova_p_value:.3f}\"\n )\n ax.set_xlabel(f\"{metric.__name__.replace('_', ' ').title()}\")\n\n plt.tight_layout()\n plt.show()\n\n\nplot_metric_comparison(classification_metrics, matthews_corrcoef)\nplot_metric_comparison(regression_metrics, r2_score)\n![\"Tukey](\"https://jonswain.github.io/images/chemprop-rf/chemprop-rf-3.png\"/)
![\"Tukey](\"https://jonswain.github.io/images/chemprop-rf/chemprop-rf-4.png\"/)
As a membership organization, members are at the core of everything we do. Understanding our members and ensuring good 2-way communication is therefore extremely important. Over the last year, we have been working on clarifying and aligning DataCite's membership model. We believe this will enable us to work with all organizations in the same way and give us more insight into the organizations using our services.
\nTypes of Membership
\nPreviously, we distinguished between Members, Providers, Clients, and Consortia. This often led to confusion among organizations about whether or not they are a DataCite Member. Therefore, going forward, we distinguish three clear membership categories: Member-only, Direct Member, and Consortium Member. The Member-only category remains unchanged and consists of Members that support DataCite's data sharing mission, want to collaborate with DataCite and/or be part of DataCite's governance. These Members do not register DOIs.
\nWhen it comes to Members that do consume DOI services, we distinguish two Member types: the Direct Member and the Consortium Member.
\nDirect Member
\nThis type of Member is a single organization that joins DataCite as a Member in its own right and consumes DOI services. This single organization may have one or more repositories under their umbrella, but those repositories are under the same administrative structure as the organization. A Direct Member also takes on membership responsibilities such as participating in DataCite's governance.
\n
Consortium Member
\nA consortium is a group of like-minded organizations that have come together to collectively participate in DataCite's community and governance activities and to use DataCite's DOI services. A consortium is composed of two or more non-profit organizations that are under different administrative structures. Consortia are generally located in a single country or subject-based. Organizations within a consortium can work with one or more repositories that are under the same administrative structure as that organization.
\nThe entire consortium is a Member collectively, meaning that the consortium has a single representative in DataCite's governance structure, regardless of the number of organizations that make up the consortium. One organization will take on the role of Consortium Lead, with the others participating as Consortium Organizations.
\n
You can find more information about the different types of membership on our membership page.
\nRepositories
\nRepositories play a key role in DataCite services and are therefore a key component of the DataCite membership model. We define a repository as a service operated by research organizations, where research materials are stored, managed, and made accessible. A repository is a single unit and DataCite links the repository to information in re3data, where additional repository metadata are available.
\nHowever, we realize that not all content is hosted in repositories and therefore DataCite also includes periodicals as a designation. The term periodicals include journals, proceedings, books, blogs, and working paper series.
\nThe repository information provided will be used and displayed throughout DataCite services including search, data usage statistics, Member DOI statistics, and third-party indexing services.
\nDOI Fabrica
\nGoing forward, we'll be making some changes to DOI Fabrica to align with the membership categories described above.
\nRepositories
\nWe'll be replacing the old \"Clients\" terminology with \"Repositories\" and making a stronger connection to the repository concept by adding metadata elements from re3data to the repository settings information displayed in Fabrica. Periodicals, as their own sub-category of research material service, will be designated by a periodicals badge for easy identification.
\nConsortia
\nWe'll implement the Consortium Member structure so that what Consortium Organizations can see and what they have permissions to administer in Fabrica aligns with their unique status according to the model. Implementing this properly will mean a brief transfer period (on the order of hours) for the DOIs belonging to these organizations. We'll be reaching out to our current Consortium Members and their constituent organizations over the next few weeks to confirm their repository details and schedule their transfer periods. For existing Members that are planning on forming Consortia, we will work with you to ensure that your current Clients are successfully migrated to Consortium Organizations with Repositories in Fabrica.
\nTimeline
\nWe're planning on rolling the consortium and repository changes into the test system during the last week of August, and we plan to go live in production by mid-September. In the interim, there will be a brief period when the test system will reflect the new changes while the production system will not, so what you see in the two systems will not match. During this period we'd welcome your feedback in addition to our more formal testing.
\nWhat will I see?
\nMember Information
\nWith all these changes happening, it's critical for us that we have the correct contact information for each organization. This information is used on our Member page, to administer Fabrica, to communicate with you, and for billing purposes. Therefore, we want to encourage all our Members, both our long-time Members and those going through changes, to update your contact information in your Member settings. We've recently added the ability to add multiple contacts per organization, based on relevant themes of interest. Please see DataCite's privacy policy to better understand how we use personal information. Expect reminders from us because at the end of the year we'll be replacing our current mailing lists with the information obtained through this form!
\nWe realize this is a lot of information and we want to support you throughout these changes in any way we can. We created an FAQ section on our support site and are happy to talk to you at any time to find the best model for your organization. Please don't hesitate to contact support@datacite.org at any point during this process.
\nThe post The DataCite Membership model: consortia, repositories, and more appeared first on DataCite.
","doi":"https://doi.org/10.5438/gk09-ba24","guid":"http://datacite.org/the-datacite-membership-model-consortia-repositories-and-more/","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1566777600,"rid":"405s4-00m63","summary":"As a membership organization, members are at the core of everything we do. Understanding our members and ensuring good 2-way communication is therefore extremely important. Over the last year, we [\u2026] The post The DataCite Membership model: consortia, repositories, and more appeared first on DataCite.","title":"The DataCite Membership model: consortia, repositories, and more","updated_at":1788766955,"url":"https://datacite.org/blog/the-datacite-membership-model-consortia-repositories-and-more/","version":"v1"},{"authors":[{"contributor_roles":[],"family":"Eden","given":"Terence","url":"https://orcid.org/0000-0002-9265-9069"}],"blog":{"authors":null,"community_id":"61ce553a-bafd-4aba-a952-d3bab5e85bcc","created":1788652800,"current_feed_url":null,"description":"Regular nonsense about tech and its effects \ud83d\ude43","doi":"https://doi.org/10.59350/shkspr","favicon":"https://rogue-scholar.org/api/communities/61ce553a-bafd-4aba-a952-d3bab5e85bcc/logo","feed_format":"application/atom+xml","feed_url":"https://shkspr.mobi/blog/feed/atom/","filter":null,"generator":"WordPress","home_page_url":"https://shkspr.mobi/blog","issn":"2753-1570","language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"shkspr","status":"active","subfield":"1712","title":"Terence Eden's Blog","updated":1788694460,"use_api":false},"blog_name":"Terence Eden's Blog","blog_slug":"shkspr","content_html":"I imagine everyone here has received an unsolicited message telling them that their tax is overdue and that they urgently need to visit Genuine-Tax-Payment-Website.fart or that a parcel is delayed at customs and you can pay a small sum for its release at Almost-The-Right-Acronym.ak
\nYou know it is a scam. Most people just mark as spam and move on with their day. But a significant number of people don't. They hastily visit the site, tap in their credit card details, give it their mother's maiden name, confirm address, upload a nude selfie, and only then realise that they've been had.
\nThe Internet works at pretty close to the speed of light. You can register a .uk domain and a minute later it's accessible from the other side of the planet. Brilliant for users who want to quickly launch a website. Also brilliant for abusers who want to launch a spam campaign.
\nBy the time enough people have reported the scammers' domain as suspicious, it is too late. In the time it takes for a registrar to disable the domain, or for its name to make its way to the Safe Browsing List, a million messages have already been sent and enough people have handed over their details.
\nWe're told that \"the purpose of a system is what it does\". At the moment, the Domain Name System's purpose seems to be a vector for criminals to run scams on people at a terrifyingly high rate.
\nBIG!
\nThere's a great blog post by Andrew Campling which reports on this startling claim:
\n\nThe study found that at least 10% of all new gTLD domain names registered during the year had subsequently appeared on security blocklists by the time of analysis. It estimated that, taking account of subsequent blocklisting and associated domains not themselves blocklisted, the share of names registered by malicious actors may be closer to 20%.
\nDNS Abuse and Criminal Infrastructure: Beyond Definitions and Blocklists (emphasis added)
That links to a presentation by Interisle which contains some rather shocking statistics (albeit with disputed methodology). It looks at generic Top Level Domains (gTLD) - those are things like .com and .fun rather than country code TLDs (ccTLD) like .uk and .de.
\nIt says 85 million new registrations of gTLDs were made in 2025. Of those 8.5 million were added to blocklists by May 2025. It reckons that a 10% abuse rate is the likely floor for these numbers and it's probably closer to 20%. One in five newly registered domains with a gTLD are scams. That's a bloody crisis.
\n13 TLDs had more than 50% of their registrations blocklisted.
\n\nI can understand why .bid and .loan are popular with scammers. But why .mobi?! What did I ever do to you, eh?
\nWho are the scammers registering these through?
\n\nAh, our old friends at NameCheap. See Why do scammers love NameCheap?
\nIf those five registrars had more effective policies, it might significantly dent the scammers' ability to ply their devious wares. Or they might just move on to other registrars.
\nAs the report points out:
\n\nsuspension rates for blocklisted domains were 7.4% to 16.3%.
The full report is on the Interisle website.
\nI don't know.
\nIn the first instance, it might make sense for registrars to do strong Know Your Customer (KYC) checks on anyone buying a domain. But that stops anyone who wants to anonymously register I-Hate-Nintendo.whatever without risking the wrath of Intellectual Property lawyers.
Also, criminals have access to stolen money and stolen cards. They can convince a hapless mule to register a domain on the criminals' behalf.
\nRegistrars could ask for an escrow payment. Pay \u20ac9 for the domain name put \u20ac900 in escrow. If your domain appears on a blocklist within the year, you forfeit the money. Criminals with stolen funds are unlikely to care but it would probably put off lots of people from getting a new domain.
\nThere are various banned words and phrases depending on the TLD. For example, South Sudan has a list of political words which they don't want associated with their .ss ccTLD.
\nBut if one gTLD bans a word, a different one might not. A scammer doesn't care if the gTLD is .arse or .elbow - they just want the start of the domain to look legitimate.
\nSome registrars have strings that they don't allow. In fairness to NameCheap, when I tried to register dwp-payments-gov-uk.pizza it told me that domain was banned. It wouldn't let me get any gTLD with that name.
But all it takes is one registrar to be slightly lax and the scammers get through. Increasing the complexity of the rules is also a hell of a burden on smaller registrars.
\nBesides, it's pretty easy to get a generic enough looking domain and stick the confusing bit on a subdomain. Here are a clutch mentioned in the report:
\nhttps://gov.uk-dwpaph.bond/uk/https://gov.uk-dwpcjh.bond/uk/https://gov.uk-dwpclc.bond/ukhttps://gov.uk-dwpclw.bond/ukhttps://gov.uk-dwpclj.bond/uk/Perhaps there ought to be a delay before a new domain goes live to allow people to object to it? That would give governments, banks, delivery companies, and a dozen more \"important\" organisations a right to veto any \"dodgy\" looking domain.
\nBut suppose someone wants to register gov-uk-stole-my-horse.horse to protest the government's cruel policy of stealing horses - is that a legitimate use of a domain? What if the Darwin Pensioner Divas - a group of elderly singers - want to take payments for their new album of goth/punk covers, can the DPD delivery company veto dpd-payments.music?
Do we want a domain name system where powerful companies control exactly which domains we can register? If I have an idea for a domain on a Friday night do I have to wait until Monday before it can be launched? Are those companies realistically able to parse millions of domains per year and have a low false-positive rate?
\nAll of these things are possible - but all of them come with an impact on legitimate users. To be clear, I don't know what the right answer is.
\nLots! It has been a few years since I've been to an ICANN meeting, but even back then the topic of abuse was high on the agenda. They appear to be looking at ways to coordinate abuse reports between various entities, along with some other policies which should hopefully work.
\nThere are two salient points from one of the discussions held at the recent meeting
\n\nIf anybody thinks that in our current age of AI and as we move into different kinds of computing, DNS abuse is going to numerically stay steady and we will have a downward effect on that baseline 2027 number. I'm not sure that that's an accurate assumption. I think it's going to be the other thing, which is [\u2026] it's going to be easier to abuse the DNS.
And
\n\nAbusers are going to abuse because it's just too lucrative, because no matter what we do, they will find the way to make profit off of that, and will try to circumvent everything that we do. That is not a reason not to do it, though.
Quite!
\nAs I said, I don't know the answer to this. What I do know is, much like Android's app ecosystem being a haven for scammers, DNS is facing a crisis. When trust in a system goes, only chaos follows.
\nI don't want to live in a world where I have to show my passport and pay thousands of pounds to register a domain which is only available after being vetted by private interests. But I also don't want to live in a world where scammers have effectively no deterrent from abusing millions of people.
\nThe purpose of a system is what it does. I hope DNS's purpose can become less dangerous while still remaining open.
\n","doi":"https://doi.org/10.59350/395ha-fss97","guid":"https://shkspr.mobi/blog/?p=74588","image":"https://shkspr.mobi/blog/wp-content/uploads/2026/08/gTLDs.webp","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1788652800,"rid":"8hyr0-9p774","summary":"I imagine everyone here has received an unsolicited message telling them that their tax is overdue and that they urgently need to visit Genuine-Tax-Payment-Website.fart or that a parcel is delayed at customs and you can pay a small sum for its release at Almost-The-Right-Acronym.ak You know it is a scam. Most people just mark as spam and move on with their day.","tags":["/etc/","ICANN","Internet","Scam","Spam"],"title":"The purpose of DNS is to spread scams","updated_at":1788766212,"url":"https://shkspr.mobi/blog/2026/09/the-purpose-of-dns-is-to-spread-scams/","version":"v1"},{"authors":[{"contributor_roles":[],"family":"Swain","given":"Jonathan","url":"https://orcid.org/0000-0003-4457-1481"}],"blog":{"authors":[{"name":"Jon Swain","url":"https://orcid.org/0000-0003-4457-1481"}],"community_id":"13f55986-f209-443c-ae0d-2f9f3f521e5a","created":1788652800,"current_feed_url":null,"description":"I am a cheminformatician and data scientist, originally from the UK, but often found in Aotearoa (New Zealand). I'm interested in using data science and machine learning to solve problems in drug discovery.","doi":"https://doi.org/10.59350/jonswain","favicon":"https://rogue-scholar.org/api/communities/13f55986-f209-443c-ae0d-2f9f3f521e5a/logo","feed_format":"application/atom+xml","feed_url":"https://jonswain.github.io/feed.xml","filter":null,"generator":"Jekyll","home_page_url":"https://jonswain.github.io/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"jonswain","status":"active","subfield":"3002","title":"Jon Swain","updated":1778680800,"use_api":null},"blog_name":"Jon Swain","blog_slug":"jonswain","content_html":"
In 2020 I was working as a synthetic chemist in a laboratory in Aotearoa (New Zealand). I was getting disillusioned with the work (and academia in general) and started looking for a new career direction and new challenges. After chatting with a few friends who were working as data scientists I realised that my current interests in scientific computing, statistics, and drug discovery would transfer nicely to a job as a data scientist or cheminformatician.
\nIt felt like an emerging field within biomedical research, with huge potential. Alphafold 2 had just been released and was causing a lot of excitement, but the LLM hype was still yet to fully start (ChatGPT was still a few years away). It seemed like a perfect time to make the jump to a new, exciting field that would grow in the future.
\nWhen I started searching for information about making the switch I somehow felt both completely lost and unguided, and overwhelmed by the resources out there. Part of the thinking behind this post was: What do I wish I knew before setting out on this journey?
\nI've been asked a couple of times about making the change, mostly by old colleagues looking to move from chemistry research to data science, but I hope this advice is useful whatever your background is. I'm always happy to answer questions but thought it would be useful to have something to point people to that contained my main ideas and could be used as a starting point for further questions.
\nThere's a huge amount resources out there, and no right or wrong way to make the change. This just lists my experience and what worked for me, this will be different for everyone. One of my favourite things about data science is the low bar for entry, there's so much open-source software and data, YouTube videos, and Kaggle competitions that you only really need a basic computer to get started. On top of this I realise I was also lucky to have friends and colleagues that I could ask for advice, and a supportive family who were willing to house me whilst I re-trained. Through 2020 and 2021 we were occasionally confined at home due to COVID lockdowns (though thankfully in NZ these were fairly limited compared to the rest of the world), and with not much else going on this was an idea time for some re-training.
\nAs I asked around my network for advice on making this change, the first suggestion was always the same, learn to code in Python. So much data science work is done using Python libraries (pandas, numpy, scikit-learn, pytorch, and rdkit for chemistry). I now spend a significant part of each day writing and reading Python code and documentation.
\nI did what I often do (probably a bit too much) and scrolled reddit. I kept seeing two online Python courses recommended. Automate the Boring Stuff with Python Programming by Al Sweigart, and 100 Days of Code: The Complete Python Pro Bootcamp by Dr. Angela Yu. Both were available on Udemy, which seems to regularly have special offers and I completed both courses for around $10 (US) each. There are definitely similar courses available for free so it's a personal preference what to go for. I chose the paid courses for two reasons:
\nThe first course I completed was Automate the Boring Stuff with Python Programming by Al Sweigart. It's a fairly short course but very well taught, designed for anyone who works regularly with a computer and wants to automate repetitive tasks, so they never have to do them again. It doesn't go too deep into theory or style and gets stuck in with practical applications quickly. Within a couple of days, I was building small programs and could start to see how I could use it for to automate parts of my job. Before starting I had no idea if I would even enjoy programming in Python, and this was an ideal introduction (turns out I did enjoy it) and gave me enough understanding to look at other courses that might be useful. It seemed perfect to make my current job less boring, but not enough for the new job I wanted, so I started looking for something more in-depth.
\nThe next course I completed was 100 Days of Code: The Complete Python Pro Bootcamp by Dr. Angela Yu. This is a much more in-depth course that covers a huge range of applications for Python. Every day has a couple of hours of lectures that cover theory or usage, and a few problem sets to apply what you've learnt. Every 10 days there's a small project that uses everything you've learnt, and at the end there's about 20 days of projects to construct a personal Python portfolio. I found some sections a little confusing (the first time OOP is introduced), and there's a large section borrowed from a course on web development that I felt dragged on a bit (it's useful to understand HTML, CSS, and JS for things like web-scraping, but I didn't need quite that much detail on building websites). There's a significant section towards the end on data science that was particularly useful. Overall, this course probably made the biggest difference in my re-training.
\n
At this stage I started thinking about jobs, but everyone I spoke to gave the same response: \"What you're doing is the right idea, but we're just looking for someone with more experience.\" This seems to be a classic problem when re-training, how do you get your first experience when everyone wants someone with experience? There are a few ways to go about this, but the main two seem to be either an internship or a bootcamp. Both have pros and cons, and I decided to go down the bootcamp route. To me the main advantages were the ease of organisation (I was relocating to the UK at the time and having something waiting for me was nice), and the shorter duration (most internships seemed to be multiple months, whereas the bootcamp I did was 5 weeks).
\nAt the recommendation of a friend, I applied to Science to Data Science (S2DS), a 5-week, intensive, project-based bootcamp for research scientists with a PhD or MSc looking to move into data science. The application involved a short technical project to test your exploratory data analysis skills and a behavioural interview (that seemed mostly to filter out psychopaths).
\nFor the 5 weeks you're part of a small group that gets paired with a company, working on a data science project. There's a couple of talks in the first week about best practices, but you're working on your project straight away. My team was paired with Deutsche Welle (DW), a German broadcaster. They were interested in the gender breakdown of people mentioned and quoted in their articles, so the project involved Natural Language Processing (NLP) and Named Entity Recognition (NER), as well as using techniques like web-scraping for data gathering and machine learning classifiers for predicting gender from the name.
\nWorking in a team made the experience a lot of fun. We all came from a range of backgrounds and had different skillsets, but all started the bootcamp from similar positions so were learning together. It also gave a first experience using git and github collaboratively. In addition to the team, we had a mentor from S2DS who helped with problems, and contacts with DW who helped with project direction.
\nWhen I applied S2DS was \u00a3800 for the entire course, not cheap, but much cheaper than many bootcamps I saw advertised online. I made the calculation that getting a job and starting to get paid more quickly would pay off in the long run. I've heard bad things about some bootcamps found online, so always do your research before signing up. If you've got the connections or are willing to put the work in organising it yourself, an internship is another great way to gain experience without the up-front cost.
\nMy experience with S2DS was positive, and that seems to invariably be the feeling among other alumni. Having this experience on my CV and a project to talk about in interviews was invaluable during my job hunt. If you do decide to sign up, feel free to put me as your reference, I think I get some Amazon vouchers or similar. One unique advantage from completing S2DS is the ongoing career advice and the community of alumni, which will be useful going forwards in my career.
\nTwo other online courses I completed were CS50x and CS50 SQL from Harvard.
\nCS50x is Harvard's introductory computer science course. It's an 11-week course with lectures, short videos on specific subjects, and problem sets each week. Whilst not essential information for a data scientist, I feel having a much better understanding of how a computer works has made me a better data scientist. The first few weeks working through problems using C made me appreciate how easy Python is!
\nCS50 SQL is a 7-week SQL course with the same structure that teaches you all the basics you'll need for a job in data science. Most of the 7 weeks are done using SQLite, but it later moves onto ProgreSQL and MySQL.
\nWith a background in chemistry research, cheminformatics seemed like a natural field for me. Giving me the opportunity to combine by experience in chemistry and drug discovery with my new skills in Python programming and data science. I looked for short online courses and textbooks to get experience with cheminformatics problems and gain an understanding of the basics.
\nI started with the Cheminformatics OLCC. This is an introductory course with 8 sections that cover the basics of cheminformatics. It starts with help setting up a Python environment, and covers topics such as representing molecules, chemical databases, QSAR modelling, and simple machine learning for chemistry.
\nI then worked through Deep Learning for the Life Sciences. As datasets in the life sciences get larger, deep learning becomes more powerful. This covers how deep learning is used on molecules, proteins, and nucleic acids, including code so you can follow along.
\nI also attended the AI4SD Machine Learning Summer School at the University of Southampton. This was a week of lectures on topics such as machine learning, github, and LaTeX, with a focus on applications in chemistry. The week ended with a hackathon where we worked as a team to solve a chemical property prediction problem.
\nTeachOpenCADD looks like another great resource. I've not yet had the chance to go through it in detail, but it has a wide range of Jupyter Notebooks on all sorts of cheminformatics topics.
\nI found attending meetings and conferences on AI and cheminformatics was a great way to learn about current areas of research and meet interesting people doing similar research. This list is mostly limited to my local area (Cambridgeshire, UK).
\nThere's a load of useful blogs and newsletters for keeping up to date with cheminformatics:
\nAnd a few useful journals to add to your RSS feed:
\nWhen I started writing this post, I didn't expect it to be quite so long! Looking back at the process of re-training I realise how much work it took, and how lucky I was to have family to help support me through it all. I'm also incredibly glad I did it! My current job suits me so much more than the lab work I was doing before.
\nI think the general process I went through worked well, but there are still a few changes I would make if I could do it all again. I've heard that getting your first data science job has three equally important parts:
\nIn the search for my first job, I over-prioritised improving my skills, spending a lot of time doing courses. Whilst this is important for succeeding at a job, getting the job requires a wider focus. If I were to do it again, I would spend more time on personal projects to put on my github, and networking with other data scientists, asking about what problems their companies have and how I might be able to help fix them. When I finally did get my first data science job, it's no surprise it was on the recommendation of a previous colleague.
","doi":"https://doi.org/10.59350/qfdrh-57d64","guid":"https://jonswain.github.io/I-want-to-become-a-data-scientist-but-I-have-no-idea-where-to-start","image":"https://jonswain.github.io/images/career_change/chemistry_to_data_science.png","language":"en","license":"https://creativecommons.org/licenses/by/4.0/legalcode","published_at":1714867200,"rid":"s8tvz-zy751","summary":":seedling: Time for a change In 2020 I was working as a synthetic chemist in a laboratory in Aotearoa (New Zealand). I was getting disillusioned with the work (and academia in general) and started looking for a new career direction and new challenges.","tags":["Data-science","Cheminformatics","Career-development"],"title":"I Want to Become a Data Scientist, but I Have No Idea Where to Start\u2026","updated_at":1788763540,"url":"https://jonswain.github.io/become-a-data-scientist-starting-guide/","version":"v1"},{"authors":[{"contributor_roles":[],"family":"Swain","given":"Jonathan","url":"https://orcid.org/0000-0003-4457-1481"}],"blog":{"authors":[{"name":"Jon Swain","url":"https://orcid.org/0000-0003-4457-1481"}],"community_id":"13f55986-f209-443c-ae0d-2f9f3f521e5a","created":1788652800,"current_feed_url":null,"description":"I am a cheminformatician and data scientist, originally from the UK, but often found in Aotearoa (New Zealand). I'm interested in using data science and machine learning to solve problems in drug discovery.","doi":"https://doi.org/10.59350/jonswain","favicon":"https://rogue-scholar.org/api/communities/13f55986-f209-443c-ae0d-2f9f3f521e5a/logo","feed_format":"application/atom+xml","feed_url":"https://jonswain.github.io/feed.xml","filter":null,"generator":"Jekyll","home_page_url":"https://jonswain.github.io/","issn":null,"language":"eng","license":"https://creativecommons.org/licenses/by/4.0/legalcode","prefix":"10.59350","relative_url":null,"secure":true,"slug":"jonswain","status":"active","subfield":"3002","title":"Jon Swain","updated":1778680800,"use_api":null},"blog_name":"Jon Swain","blog_slug":"jonswain","content_html":"This is part 1 of a planned three post series on working with large chemical libraries.\nThe notebook used to create this post and all the files can be found in this github repo.
\nEstimates of the size of chemical space (the number of different molecules that could exist) range from 1050 - 1080 possible compounds, with the upper limits approaching estimates of the total number of atoms in the universe (around 1080 atoms). It's never going to be possible to make and test every possible chemical compound, even if we limit it to only \"drug-like\" compounds. Even with fast computational methods it's not always possible to exhaustively screen virtual libraries. We need methods to find the best scoring compounds from within virtual libraries without scoring every single compound.
\nSmall companies may have hundreds of thousands to millions of physical compounds in their virtual library, and larger companies will have libraries of a few million. Virtual make-on-demand libraries, especially combinatorial libraries, are expanding into the billions of compounds (often called ultra-large libraries). Enamine REAL has 48 billion possible compounds. If your scoring function takes one second per compound, you could maybe screen one million compounds (but it will take 11.6 days), but a billion seconds is nearly 32 years. As the size of virtual libraries expands, it becomes impossible to screen them exhaustively.
\nThere are a range of methods for filtering the library for screening, such as similarity searching or Thompson sampling. In this post I'll discuss active learning.
\nActive learning is a machine learning method for searching large libraries when the scoring function that is too computationally expensive to label the full library of compounds. A selection of data is labelled with the expensive scoring function, and a machine learning model is trained on these labels and used to score all compounds from within the library. The compounds with the best scores from the machine learning model are labelled using the more expensive function, and the data from this pooled and used to train a new machine learning model. This cycle is repeated until a finish criterion is met.
\nThe initial steps are:
\nFollowed then by a repeating cycle:
\nFirst we need to import the libraries we will be using.
\nimport gzip\nimport math\nimport re\nimport shutil\nimport subprocess\nfrom functools import partial\nfrom itertools import product\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem, DataStructs, Descriptors\nfrom sklearn.ensemble import RandomForestRegressor\nfrom tqdm import tqdm\nThe first function we will define is the expensive scoring function, this will take a list of SMILES and return a list of scores. For this example, I'm going to try find the compound from within the library with the lowest calculated Log P. This is actually a very fast calculation and can be done exhaustively, which means we can confirm if active learning is finding the lowest value.
\ndef calc_logp(smiles: list[str]) -> list[float]:\n \"\"\"Calculates the logP value for a list of compounds.\n\n Args:\n smiles (pd.Series): The input molecules.\n\n Returns:\n np.arary: The scores of the molecules.\n \"\"\"\n return [Descriptors.MolLogP(Chem.MolFromSmiles(s)) for s in smiles]\nNext, we need to define some useful functions for the active learning pipeline.
\nThe first function creates a virtual library by combining three sets of building blocks. The smi files used here were borrowed from Pat Walters repository on Thompson sampling. The second creates the Morgan fingerprints for all compounds within the library for training the machine learning model, which is trained using the third function. The final function scores the entire library using the machine learning model.
\ndef build_virtual_library() -> pd.DataFrame:\n \"\"\"Builds a virtual library by coupling building blocks from the input smi files.\n\n Returns:\n pd.DataFrame: A DataFrame containing the molecular objects and SMILES strings of\n the products.\n \"\"\"\n try:\n library = pd.read_csv(\"data/library.csv\", index_col=\"smiles\")\n library[\"mol\"] = [Chem.MolFromSmiles(s) for s in tqdm(library.index.to_list())]\n except FileNotFoundError:\n reaction_smarts = \"N[c:4][c:3]C(O)=O.[#6:1][NH2].[#6:2]C(=O)[OH]>>[C:2]c1n[c:4][c:3]c(=O)n1[C:1]\"\n bb_types = [\"aminobenzoic\", \"carboxylic_acids\", \"primary_amines\"]\n rxn = AllChem.ReactionFromSmarts(reaction_smarts)\n\n building_blocks = []\n for bb in bb_types:\n smil = []\n with open(Path(f\"data/{bb}_100.smi\"), \"r\") as f:\n for line in f.readlines():\n smiles, _ = line.split()\n smil.append(smiles)\n building_blocks.append(smil)\n\n total_prods = math.prod([len(x) for x in building_blocks])\n\n product_list = []\n for reagents in tqdm(product(*building_blocks), total=total_prods):\n reagent_mol_list = [Chem.MolFromSmiles(x) for x in reagents]\n products = rxn.RunReactants(reagent_mol_list)\n if products:\n Chem.SanitizeMol(products[0][0])\n product_list.append(products[0][0])\n\n library = pd.DataFrame(\n product_list,\n index=[Chem.MolToSmiles(m) for m in product_list],\n columns=[\"mol\"],\n )\n library.index.name = \"smiles\"\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n library.to_csv(\"data/library.csv\")\n return library\n\n\ndef create_morgan_fingerprints(library: pd.DataFrame) -> pd.DataFrame:\n \"\"\"Creates Morgan fingerprints for the input library.\n\n Args:\n library (pd.DataFrame): The input library.\n\n Returns:\n pd.DataFrame: The Morgan fingerprints of the input library.\n \"\"\"\n try:\n fps_df = pd.read_csv(\"data/fingerprints.csv\").set_index(\"smiles\")\n except FileNotFoundError:\n fps = [\n list(AllChem.GetMorganFingerprintAsBitVect(mol, radius=2))\n for mol in tqdm(library.mol)\n ]\n fps_df = pd.DataFrame(\n fps, columns=[f\"fp_{x}\" for x in range(len(fps[0]))], index=library.index\n )\n fps_df.to_csv(\"data/fingerprints.csv\")\n return fps_df\n\n\ndef train_ml_model(\n library: pd.DataFrame, fingerprints: pd.DataFrame\n) -> RandomForestRegressor:\n \"\"\"Trains a random forest regressor model on slow scores for the input library.\n\n Args:\n library (pd.DataFrame): The input library.\n fingerprints (pd.DataFrame): The Morgan fingerprints of the input library.\n\n Returns:\n RandomForestRegressor: The trained random forest regressor model.\n \"\"\"\n scored = library[~library[\"slow_scores\"].isna()]\n X = fingerprints.loc[scored.index]\n y = scored[\"slow_scores\"]\n regressor = RandomForestRegressor(max_depth=10)\n regressor.fit(X, y)\n return regressor\n\n\ndef score_library(\n library: pd.DataFrame,\n regressor: RandomForestRegressor,\n fingerprints: pd.DataFrame,\n) -> pd.DataFrame:\n \"\"\"Scores the entire library with the trained model.\n\n Args:\n library (pd.DataFrame): The input library.\n regressor (RandomForestRegressor): The trained random forest regressor model.\n fingerprints (pd.DataFrame): The Morgan fingerprints of the input library.\n\n Returns:\n pd.DataFrame: The input library with the model scores.\n \"\"\"\n library[\"model_scores\"] = regressor.predict(fingerprints)\n return library\nWe then need to set up the active learning pipeline. This follows the steps from above:
\nFollowed then by a repeating cycle:
\nBecause we sometimes know the best value for a scoring function, I've added an early stopping condition if the best compound is found to save time.
\ndef run_active_learning(\n library: pd.DataFrame,\n fingerprints: pd.DataFrame,\n compounds_per_round: int,\n number_of_rounds: int,\n scoring_function: callable,\n minimize: bool,\n early_stopping_value: float | None = None,\n) -> pd.DataFrame:\n \"\"\"Runs active learning on the virtual library.\n\n Args:\n library (pd.DataFrame): The input virtual library.\n fingerprints (pd.DataFrame): The Morgan fingerprints of the input library.\n compounds_per_round (int): Number of compounds to select per round.\n number_of_rounds (int): Number of active learning rounds.\n scoring_function (callable): The slow scoring function.\n minimize (bool): Whether to minimize or maximize the slow scoring function.\n early_stopping_value (float | None): If the slow scoring function finds a\n compound with this value, stops the active learning.\n\n Returns:\n pd.DataFrame: The input library with the slow scores added.\n \"\"\"\n # Select initial random sample\n initial_sample = library.sample(compounds_per_round)\n library[\"slow_scores\"] = np.NaN\n\n # Score the initial sample\n initial_scores = scoring_function(initial_sample.index.to_list())\n\n # Save the slow scores\n library.loc[initial_sample.index, \"slow_scores\"] = initial_scores\n\n # Run active learning\n al_round = 0\n while al_round < number_of_rounds:\n # Train the ML model\n model = train_ml_model(library, fingerprints)\n\n # Use the model to score the entire virtual library\n library = score_library(library, model, fingerprints)\n\n # Select the top scoring molecules with no slow scores\n top_compounds = (\n library[library[\"slow_scores\"].isna()]\n .sort_values(\"slow_scores\", ascending=minimize)\n .head(compounds_per_round)\n )\n\n # Score the top molecules with the slow function\n slow_scores = scoring_function(top_compounds.index.to_list())\n # Save the slow scores\n library.loc[top_compounds.index, \"slow_scores\"] = slow_scores\n library.loc[top_compounds.index, \"scored_round\"] = al_round\n\n al_round += 1\n if early_stopping_value is not None:\n if minimize:\n if library[\"slow_scores\"].min() == early_stopping_value:\n break\n else:\n if library[\"slow_scores\"].max() == early_stopping_value:\n break\n\n return library\nCalculated LogP (cLogP) is an important property in drug discovery. Whilst we wouldn't usually be aiming to find the minimum or maximum value in a library, it's a useful example as it can be calculated quickly enough to exhaustively screen the library to confirm we're finding the best compound. In the example below, two compounds are randomly selected from the library and their Morgan fingerprints and cLogP values are used to train a machine learning model. This model scores the entire library and the cLogP for two best scoring compounds are pooled with the earlier compounds and the process repeated. After only six rounds (12 compounds scored), the active learning algorithm can find the lowest cLogP in the library, much quicker than exhaustively screening all 132,500 compounds! The experiment was repeated ten times and the number of compounded needed is constant across the repeats. cLogP is a simple property, but it's still impressive how little data the machine learning model needed to find the best scores.
\n# Active learning parameters\nCOMPOUNDS_PER_ROUND = 2\nNUMBER_OF_ROUNDS = 50\nMINIMIZE = True\nNUMBER_OF_REPEATS = 10\n\n# Create the virtual library\nprint(\"Creating virtual library\")\nlibrary = build_virtual_library()\n\n# Create morgan fingerprints for the library\nprint(\"Creating Morgan fingerprints\")\nfingerprints = create_morgan_fingerprints(library)\n\n# Find the minimum logP in the library\nall_clogp_values = calc_logp(library.index.to_list())\nprint(f\"Minimum logP in the library: {min(all_clogp_values):.2f}\")\n\n# Run active learning\nprint(\"Running active learning\")\nscored_compounds_per_run = []\nfor i in tqdm(range(NUMBER_OF_REPEATS)):\n # Reset library between rounds\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n\n # Active learning\n library = run_active_learning(\n library=library,\n fingerprints=fingerprints,\n compounds_per_round=COMPOUNDS_PER_ROUND,\n number_of_rounds=NUMBER_OF_ROUNDS,\n scoring_function=calc_logp,\n minimize=MINIMIZE,\n early_stopping_value=min(all_clogp_values),\n )\n\n # Log results\n scored_compounds_per_run.append(len(library.loc[library.slow_scores.notna()]))\n\nave_scored_compounds = np.mean(scored_compounds_per_run) - COMPOUNDS_PER_ROUND\nprint(\n f\"Average number of scored compounds before finding lowest value: {ave_scored_compounds:.2f}\"\n)\nCreating virtual library\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 1000000/1000000 [01:42<00:00, 9747.08it/s]\n\nCreating Morgan fingerprints\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 132500/132500 [04:11<00:00, 526.70it/s]\n\nMinimum logP in the library: -5.00\nRunning active learning\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:21<00:00, 2.11s/it]\n\nAverage number of scored compounds before finding lowest value: 12.00\nMolecular similarity is often defined as the Tanimoto similarity between Morgan fingerprints. In the next experiment we'll try to maximise the Tanimoto similarity to a reference molecule. The reference molecule will be randomly selected from the library, so the maximum similarity will be 1. This turned out to be a significantly more difficult problem than the previous example minimizing cLogP. The number of compounds per round and number of rounded were increased (5000 compounds screened), and even still the reference molecule was only recovered once out of the ten repeats. This was a slightly surprising result considering the machine learning models were trained using the Morgan fingerprints that were also used in the Tanimoto similarity.
\ndef calc_similarity(comparison_smiles: list[str], ref_smiles: str) -> list[float]:\n \"\"\"Calculates the Tanimoto similarity of a reference compound to a list of\n comparison compounds.\n\n Args:\n comparison_smiles (list[str]): List of SMILES to compare to\n ref_smiles (str): SMILES of the reference compound\n\n Returns:\n list[float]: List of similarity scores\n \"\"\"\n ref_mol = Chem.MolFromSmiles(ref_smiles)\n comparison_mols = [Chem.MolFromSmiles(x) for x in comparison_smiles]\n fpgen = AllChem.GetMorganGenerator()\n ref_fp = fpgen.GetFingerprint(ref_mol)\n comparison_fps = [fpgen.GetFingerprint(x) for x in comparison_mols]\n return [DataStructs.FingerprintSimilarity(ref_fp, x) for x in comparison_fps]\n\n\n# Active learning parameters\nCOMPOUNDS_PER_ROUND = 100\nNUMBER_OF_ROUNDS = 50\nMINIMIZE = False\nNUMBER_OF_REPEATS = 10\n\n# Run active learning\nprint(\"Running active learning\")\nscored_compounds_per_run = []\nfor i in tqdm(range(NUMBER_OF_REPEATS)):\n # Choose a random reference molecule\n ref_smiles = library.sample(1).index[0]\n\n # Reset library between rounds\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n\n # Active learning\n library = run_active_learning(\n library=library,\n fingerprints=fingerprints,\n compounds_per_round=COMPOUNDS_PER_ROUND,\n number_of_rounds=NUMBER_OF_ROUNDS,\n scoring_function=partial(calc_similarity, ref_smiles=ref_smiles),\n minimize=MINIMIZE,\n early_stopping_value=1,\n )\n\n # Log results\n top = library[~library.slow_scores.isna()].sort_values(\n \"slow_scores\", ascending=MINIMIZE\n )\n scored_compounds_per_run.append(\n {\n \"round\": i,\n \"reference_molecule\": ref_smiles,\n \"compounds_scored\": len(library[~library.slow_scores.isna()]),\n \"top_scorer\": top.index.to_list()[0],\n \"top_score\": top.slow_scores.values[0],\n }\n )\n\npd.DataFrame(scored_compounds_per_run).set_index(\"round\").style.format(\n {\"top_score\": \"{:.3f}\"}\n)\nRunning active learning\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [22:30<00:00, 135.02s/it]\n| Round | \nReference compound | \nCompounds scored | \nTop score | \n
|---|---|---|---|
| 0 | \nCOc1ccc2c(=O)n(C(=N)NCC(=O)O)c([C@@H](CN)CO)nc2c1OC | \n5100 | \n0.466 | \n
| 1 | \nCC(=O)OC[C@H](C(=O)O)n1c([C@H](N)CN(C)C)nc2[nH][nH]c(=O)c2c1=O | \n3300 | \n1.000 | \n
| 2 | \nCc1cc2nc([C@H](C)NC(=O)CN)n(CC[C@H](N)C(=O)O)c(=O)c2cc1O | \n5100 | \n0.506 | \n
| 3 | \nNc1nncn1Cc1nc2c(cnc3ccnn32)c(=O)n1C(=O)[C@H](N)CC(=O)O | \n5100 | \n0.495 | \n
| 4 | \nCn1nnc2nc(C3(O)CC(N)C3)n([C@H](CO)C(=O)O)c(=O)c21 | \n5100 | \n0.488 | \n
| 5 | \nCn1cc2c(=O)n(CC[C@@H](N)C(=O)O)c([C@H](CO)NC(N)=O)nc2n1 | \n5100 | \n0.518 | \n
| 6 | \nCOC(=O)c1ccc2c(=O)n(CC[C@@H](N)C(=O)O)c([C@@H](N)CCO)nc2c1 | \n5100 | \n0.449 | \n
| 7 | \nC[C@H](O)[C@H](N)c1nc2cnc(Cl)nc2c(=O)n1[C@H](C(=O)O)[C@@H](C)O | \n5100 | \n0.434 | \n
| 8 | \nCn1cnc2c(=O)n(C(=O)C[C@H](N)C(=O)O)c([C@@H](N)CO)nc21 | \n5100 | \n0.519 | \n
| 9 | \nC[C@@H](O)[C@H](N)c1nc2nc(N)ncc2c(=O)n1C(=O)NCC(=O)O | \n5100 | \n0.458 | \n
Looking at the final run, we can see that in general the maximum similarity is not increasing significantly with further rounds, apart from the occasional spike as the model finds a similar compound. This suggests that increasing the number of rounds may not increase the performance.
\nsns.lineplot(\n library[[\"slow_scores\", \"scored_round\"]].groupby(\"scored_round\").max().values\n).set(xlabel=\"Round\", ylabel=\"Max similarity to reference compound\")\nplt.show()\n\n\n
AL was very successful at finding the lowest cLogP value in a library but struggled to recover a randomly chosen reference molecule by trying to maximize the Tanimoto similarity. The reason behind this may be due to the underlying distribution of data within the library. By plotting the cLogP values of all compounds within the library, we can see that the data is close to being normally distributed. The machine learning model is able to learn to predict cLogP well from the data it is given. For the Tanimoto similarity, the scores for the whole library of similarity to the reference compound from the final run are plotted, which have a large positive skew. There is lots of data for predicting Tanimoto similarities between 0.1 and 0.4, but very little data for scores above 0.4. So even though the machine learning model is using the same representation as is used to calculate the Tanimoto similarity (Morgan fingerprint), it doesn't have enough information to learn what makes a high score.
\nall_similarities = calc_similarity(library.index.to_list(), ref_smiles)\nfig, ax = plt.subplots(1, 2, figsize=(12, 5))\nsns.histplot(all_clogp_values, bins=50, alpha=0.5, label=\"True scores\", ax=ax[0])\nax[0].set(title=\"cLogP values for library\", xlabel=\"cLogP\", ylabel=\"Frequency\")\nsns.histplot(all_similarities, bins=50, alpha=0.5, label=\"Predicted scores\", ax=ax[1])\nax[1].set(\n title=\"Tanimoto similarity to reference compound\",\n xlabel=\"Similarity\",\n ylabel=\"Frequency\",\n)\nplt.show()\n\n\n
Molecular docking is a common technique used in drug discovery to predict how a small molecule ligand will bind to a protein target. This can be combined with a scoring function to estimate the binding affinity, and these scores can be used to rank virtual compounds for experimental validation. I am far from an expert on molecular docking, but here we'll use it as an example of how active learning could be used in a real drug discovery project.
\nDocking is slow. In this example we'll use SMINA, which seem to take around one minute for each compound. We're docking compounds to the SARS MPro protein and using active learning to try find the highest binding affinity by minimizing the predicted binding energy. At one minute per docking calculation, it would take 13 weeks to exhaustively screen the virtual library used here.
\nThe experiment was repeated three times, twice it improved on the docking scores of the initial random sample, and on one run it found a very good scoring compound in the initial random sample and wasn't able to improve on that. Due to running this on my laptop, I kept the number of compounds per round (five) and the number of rounds (10) low. With more computing power, it would be interesting to see how the performance of the active learning improves with increasing these values.
\ndef dock_mols(smiles: list[str]) -> list[float]:\n \"\"\"Dock the molecules to the SARS MPro protein and return the affinity of the best\n pose for each.\n\n Args:\n smiles: list[str]: The molecules to dock\n\n Returns:\n list[float]: The affinities of the best pose for each input molecule\n \"\"\"\n mols = [Chem.MolFromSmiles(s) for s in smiles]\n scores = []\n\n for mol in mols:\n with Chem.SDWriter(\"data/docking/tmp_conf.sdf\") as w:\n m = Chem.AddHs(mol)\n _ = AllChem.EmbedMultipleConfs(m, numConfs=5, numThreads=0)\n confs = m.GetConformers()\n for c in confs:\n w.write(m, confId=c.GetId())\n\n subprocess.run(\n [\n \"./smina.osx.12\",\n \"--exhaustiveness\",\n \"10\",\n \"--cpu\",\n \"10\",\n \"--seed\",\n \"0\",\n \"--autobox_ligand\",\n \"data/docking/ligand_only.pdb\",\n \"-r\",\n \"data/docking/protein_minus_ligand.pdb\",\n \"-l\",\n \"data/docking/tmp_conf.sdf\",\n \"-o\",\n \"data/docking/tmp_conf_docked.sdf.gz\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n )\n\n with gzip.open(\"data/docking/tmp_conf_docked.sdf.gz\", \"rb\") as f_in:\n with open(\"data/docking/tmp_conf_docked.sdf\", \"wb\") as f_out:\n shutil.copyfileobj(f_in, f_out)\n\n # The output sdf wasn't loading so this is a workaround\n with open(\"data/docking/tmp_conf_docked.sdf\", \"r\") as f:\n text = f.read()\n affinities = re.findall(r\"\\n(-\\d.\\d+)\" , text)\n affinities = [float(x) for x in affinities]\n\n scores.append(min(affinities))\n\n return scores\n\n\n# Active learning parameters\nCOMPOUNDS_PER_ROUND = 5\nNUMBER_OF_ROUNDS = 10\nMINIMIZE = True\nNUMBER_OF_REPEATS = 3\n\n# Run active learning\nprint(\"Running active learning\")\nscored_compounds_per_run = []\ncumulative_scores = {}\nfor i in tqdm(range(NUMBER_OF_REPEATS)):\n # Reset library between rounds\n library[\"slow_scores\"] = np.NaN\n library[\"model_scores\"] = np.NaN\n library[\"scored_round\"] = np.NaN\n\n # Active learning\n library = run_active_learning(\n library=library,\n fingerprints=fingerprints,\n compounds_per_round=COMPOUNDS_PER_ROUND,\n number_of_rounds=NUMBER_OF_ROUNDS,\n scoring_function=dock_mols,\n minimize=MINIMIZE,\n )\n\n # Log results\n top = library[~library.slow_scores.isna()].sort_values(\n \"slow_scores\", ascending=MINIMIZE\n )\n scored_compounds_per_run.append(library.slow_scores.min())\n cumulative_scores[f\"run_{i}\"] = [\n library[library.scored_round <= j].slow_scores.min()\n for j in range(NUMBER_OF_ROUNDS)\n ]\n\nprint(\n f\"Average lowest docking score after screening {COMPOUNDS_PER_ROUND * NUMBER_OF_ROUNDS} compounds: {sum(scored_compounds_per_run) / len(scored_compounds_per_run):.2f}\"\n)\n\n# Plot the cumulative best score found in each round for the final run\nsns.lineplot(cumulative_scores).set(\n xlabel=\"Round\", ylabel=\"Minimum docking score\", xticks=range(NUMBER_OF_ROUNDS)\n)\nplt.show()\nRunning active learning\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 3/3 [4:54:58<00:00, 5899.39s/it] \n\nAverage lowest docking score after screening 50 compounds: -9.14\n\n\n
In this post we've discussed active learning and seen three examples of how it could be used when our virtual library has grown too large to exhaustively screen every compound.
\nI've tried to visualise and compare distributions using violin plots for reports and presentations in the past, and the feedback I've got was generally\u2026 not great. When searching for better methods I came across this excellent blog post by Alex Belengeanu on raincloud plots and I'm now a big fan.
\nThere's lots of methods for visualising distributions to compare different populations, and each has its own pros and cons. Histograms (often with a kernel density estimation) are good method for showing complex distributions, but it's not possible to easily compare statistics such as the mean or quantiles from the raw graph, and there's also some loss of information from binning or smoothing, which could misrepresent the data, such as making discrete values appear continuous. Box and whisker plots are great for being able to compare statistics such as means and quartiles, but they don't represent complex distributions such as bimodal data well. A jitter plot is good for showing the structure of the underlying data, but difficult to compare statistics or absolute densities of points.
\nA violin plot with box plot combines the strengths of the two methods to visualise complex distributions with easy to read off statistics. The main downside is they're ugly (IMO). Every time I've used them in a report or presentation the discussion has veered off to whether they look more like sea creatures or medieval weapons than violins, with very little discussion of the actual data!
\nWhilst searching for better methods I came across this blog post by Alex Belengeanu on raincloud plots, which seem to be the best alternative.
\nRaincloud plots are pretty simple, they're just a combination of half a violin plot (sometimes called a ridgeline plot), a box plot, and a jitter plot. They're intuitive to read and show distributions clearly, as well as being much more aesthetically pleasing than many alternatives. I've customised Alex Belengeanu's code to suit my needs, and it's all in a function in this GitHub repo.
\n
This is part 2 of a a planned three post series on working with large chemical libraries.\nThe notebook used to create this post, and all the files can be found in this github repo.
\nCombinatorial libraries can grow quickly, combining three sets of 1,000 building blocks allows you to access one billion possibilities. For example, 1,000 amines, 1,000 halide functionalised carboxylic acids, and 1,000 boronic acids can be combined to make a billion-member library. That would be a lot of synthesis and a lot of testing if you wanted to screen every member of the library! Even if a computational method took 1 second to score each virtual compound, it would take 32 years to score the entire library. Make-on-demand chemical suppliers such as Enamine have libraries that contain multiple billions of compounds, often created by combining sets of building blocks using reliable chemistry. We need methods for effectively sampling this chemical space.
\nIn this post I'll be looking into using a genetic algorithm to search a small combinatorial library. These are quick and incredibly simple methods that can be used to search combinatorial libraries for the combination of building blocks that maximises or minimises a scoring function.
\nWhilst this method uses a combinatorial library that contains all possible configurations of all three building blocks, it's possible to extend it to more complex combinatorial libraries by storing the building blocks and reactions as a graph.
\nGenetic algorithms are biologically inspired, based on biological natural selection. They use a very simple algorithm, with no machine learning or complex statistics. One of the big advantages is that you don't need to fully enumerate of score the entire library, only certain combinations of building blocks need to be enumerated and scores, reducing memory and computational costs. Working in building block space rather than with enumerated structures means the complexity will scale with number of building blocks, which will increase much more slowly than the total library size. The downside of the genetic algorithm is that it requires combinatorial libraries and won't work with large collections of diverse molecules not made up from the same building blocks.
\nInitially a random population of the unlabelled data is selected by randomly choosing building blocks and enumerating the reaction products from these building blocks.
\nThe genetic algorithm has a cycle that is repeated, and each round is called a generation.
\nFirst, we need to import the libraries we will be using.
\nimport math\nimport time\nfrom functools import partial\nfrom itertools import product\nfrom pathlib import Path\nfrom typing import Callable\n\nimport numpy as np\nimport pandas as pd\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem, DataStructs, Descriptors\nfrom tqdm import tqdm\nThe first function we will define is the expensive scoring function, this will take a list of RDKit molecules and return a list of scores. As in the previous example, I'm going to try find the compound from within the library with the lowest calculated Log P (cLogP). This is actually a very fast calculation and can be done exhaustively, which means we can confirm if the genetic algorithm is finding the lowest value and triggering the early stopping.
\ndef calc_logp(mols: list[Chem.Mol]) -> list[float]:\n \"\"\"Calculate the logP value for a list of compounds.\n\n Args:\n mols (list[Chem.Mol]): The molecules.\n\n Returns:\n list[float]: The scores of the molecules.\n \"\"\"\n return [Descriptors.MolLogP(mol) for mol in mols]\nNext, we need to define some useful functions for the active learning pipeline.
\nThe first function takes the indices of three building blocks and enumerates the product of their reaction. Since not all combinations are possible, it will return None if the reaction fails. The second function takes the surviving population from the selection pressure and creates a new population by shuffling the building blocks. The third takes the new shuffled population and adds in random mutations, this prevents the algorithm getting stuck in a local minimum by ensuring building block diversity.
\ndef make_molecule(\n r1: int,\n r2: int,\n r3: int,\n building_blocks: list[list[str]],\n rxn: Chem.rdChemReactions.ChemicalReaction,\n) -> Chem.Mol | None:\n \"\"\"React building blocks to make a molecule.\n\n Args:\n r1 (int): Index of the first building block.\n r2 (int): Index of the second building block.\n r3 (int): Index of the third building block.\n building_blocks (list[list[str]]): The building blocks.\n rxn (Chem.rdChemReactions.ChemicalReaction): The reaction to combine the\n building blocks.\n\n\n Returns:\n Chem.Mol | None: The molecule or None if the reaction fails.\n \"\"\"\n bbs = [building_blocks[0][r1], building_blocks[1][r2], building_blocks[2][r3]]\n reagent_mol_list = [Chem.MolFromSmiles(x) for x in bbs]\n products = rxn.RunReactants(reagent_mol_list)\n if products:\n Chem.SanitizeMol(products[0][0])\n return products[0][0]\n return None\n\n\ndef shuffle_population(data: pd.DataFrame, population_size: int) -> pd.DataFrame:\n \"\"\"Shuffle the population building blocks.\n\n Args:\n data (pd.DataFrame): The surviving population data.\n population_size (int): The size of the population.\n\n Returns:\n pd.DataFrame: The shuffled population.\n \"\"\"\n r1s = data[\"r1\"].to_list()\n r2s = data[\"r2\"].to_list()\n r3s = data[\"r3\"].to_list()\n new_population = pd.DataFrame(\n {\n \"r1\": [np.random.choice(r1s) for _ in range(population_size)],\n \"r2\": [np.random.choice(r2s) for _ in range(population_size)],\n \"r3\": [np.random.choice(r3s) for _ in range(population_size)],\n }\n )\n return new_population\n\n\ndef mutate_population(\n data: pd.DataFrame, mutation_rate: float, building_blocks: list[list[str]]\n) -> pd.DataFrame:\n \"\"\"Mutate the population building blocks.\n\n Args:\n data (pd.DataFrame): The surviving population data.\n mutation_rate (float): The fraction of building blocks to mutate.\n building_blocks (list[list[str]]): The building blocks.\n\n Returns:\n pd.DataFrame: The mutated population.\n \"\"\"\n for count, column in enumerate(data.columns):\n selection = data.sample(frac=mutation_rate)\n data.loc[selection.index, column] = np.random.randint(\n 0, len(building_blocks[count]), len(selection)\n )\n return data\nThe genetic algorithm needs the lists of building blocks, a reaction to couple them to enumerate products, and a function to score the molecules. It also has a number of hyperparameters than can be tuned:
\npopulation_size: The number of molecules in each generation to enumerate. A larger population will have a greater diversity of building blocks but will take longer to enumerate and score.num_generations: The number of cycles of the genetic algorithm to run. More generations will give a better chance of finding the optimum value but will take longer to run.selection_pressure: The fraction of the population to remove each cycle. A stronger selection pressure will speed up selection but increase the chance of being caught in a local minimum.mutation_rate: The proportion of building blocks in the new population to replace with a random building block. A too high value will stop the algorithm finding the best combination, but a too low value will also increase the chance of being caught in a local minimum.def run_genetic_algorithm(\n building_blocks: list[list[str]],\n rnx: Chem.rdChemReactions.ChemicalReaction,\n population_size: int,\n num_generations: int,\n selection_pressure: float,\n mutation_rate: float,\n scoring_function: Callable,\n minimize: bool,\n early_stopping_value: float | None = None,\n) -> pd.DataFrame:\n \"\"\"Run the genetic algorithm.\n\n Args:\n building_blocks (list[list[str]]): The building blocks.\n rnx (Chem.rdChemReactions.ChemicalReaction): The reaction to combine the\n building blocks.\n population_size (int): The size of the population in each generation.\n num_generations (int): The number of generations to run.\n selection_pressure (float): The fraction of the population to discard.\n mutation_rate (float): The fraction of the population to mutate.\n scoring_function (Callable): The function to score the molecules.\n minimize (bool): Whether to minimize the scoring function.\n early_stopping_value (float | None, optional): The value to stop early at.\n Defaults to None.\n\n Returns:\n pd.DataFrame: The history of the population.\n \"\"\"\n # To keep track of the population each generation\n history = pd.DataFrame()\n\n # Choose initial population\n population = pd.DataFrame(\n {\n \"r1\": np.random.randint(0, len(building_blocks[0]), population_size),\n \"r2\": np.random.randint(0, len(building_blocks[1]), population_size),\n \"r3\": np.random.randint(0, len(building_blocks[2]), population_size),\n }\n )\n\n for generation in range(num_generations):\n # Generate molecules\n population[\"ROMol\"] = population.apply(\n lambda x: make_molecule(x[\"r1\"], x[\"r2\"], x[\"r3\"], building_blocks, rnx),\n axis=1,\n )\n\n # Kill off fatal mutants\n population = population.dropna().reset_index(drop=True)\n\n # Score molecules\n population[\"score\"] = scoring_function(population[\"ROMol\"].to_list())\n\n # Save population for analysis\n population[\"generation\"] = generation\n history = pd.concat([history, population])\n\n # Early stopping\n if early_stopping_value:\n if minimize:\n if round(population[\"score\"].min(), 5) <= round(\n early_stopping_value, 5\n ):\n break\n else:\n if round(population[\"score\"].max(), 5) >= round(\n early_stopping_value, 5\n ):\n break\n\n # Select top performing molecules\n population = (\n population.sort_values(\"score\", ascending=minimize)\n .head(int(population_size * (1 - selection_pressure)))\n .reset_index(drop=True)\n )\n\n # Shuffle and mutate\n population = shuffle_population(population, population_size)\n population = mutate_population(population, mutation_rate, building_blocks)\n population = population.drop_duplicates().reset_index(drop=True)\n\n return history\nThis function fully enumerates a virtual library by combining three sets of building blocks. The smi files used here were borrowed from Pat Walters repository on Thompson sampling. This is not necessary for the genetic algorithm but will allow us to see if the genetic algorithm if finding the best values.
\ndef build_virtual_library() -> pd.DataFrame:\n \"\"\"Build a virtual library by coupling building blocks from the input smi files.\n\n Returns:\n pd.DataFrame: The virtual library.\n \"\"\"\n try:\n library = pd.read_csv(\"data/library.csv\", index_col=\"smiles\")\n library[\"mol\"] = [Chem.MolFromSmiles(s) for s in tqdm(library.index.to_list())]\n except FileNotFoundError:\n reaction_smarts = \"N[c:4][c:3]C(O)=O.[#6:1][NH2].[#6:2]C(=O)[OH]>>[C:2]c1n[c:4][c:3]c(=O)n1[C:1]\"\n bb_types = [\"aminobenzoic\", \"carboxylic_acids\", \"primary_amines\"]\n rxn = AllChem.ReactionFromSmarts(reaction_smarts)\n\n building_blocks = []\n for bb in bb_types:\n smil = []\n with open(Path(f\"data/{bb}_100.smi\"), \"r\") as f:\n for line in f.readlines():\n smiles, _ = line.split()\n smil.append(smiles)\n building_blocks.append(smil)\n\n total_prods = math.prod([len(x) for x in building_blocks])\n\n product_list = []\n for reagents in tqdm(product(*building_blocks), total=total_prods):\n reagent_mol_list = [Chem.MolFromSmiles(x) for x in reagents]\n products = rxn.RunReactants(reagent_mol_list)\n if products:\n Chem.SanitizeMol(products[0][0])\n product_list.append(products[0][0])\n\n library = pd.DataFrame(\n product_list,\n index=[Chem.MolToSmiles(m) for m in product_list],\n columns=[\"mol\"],\n )\n library.to_csv(\"data/library.csv\")\n return library\n\n\n# Chemistry parameters\nRXN = AllChem.ReactionFromSmarts(\n \"N[c:4][c:3]C(O)=O.[#6:1][NH2].[#6:2]C(=O)[OH]>>[C:2]c1n[c:4][c:3]c(=O)n1[C:1]\"\n)\nBUILDING_BLOCKS = []\nfor bb in [\"aminobenzoic\", \"carboxylic_acids\", \"primary_amines\"]:\n smil = []\n with open(Path(f\"data/{bb}_100.smi\"), \"r\") as f:\n for line in f.readlines():\n smiles, _ = line.split()\n smil.append(smiles)\n BUILDING_BLOCKS.append(smil)\nThe genetic algorithm was used to find the lowest cLogP from within the library. This was repeated 10 times with no tuning of the hyperparameters which may improve performance. The genetic algorithm reliably finds the combination of building blocks with the lowest cLogP in the combinatorial library, nearly 100 times faster than enumerating the entire library. On average the genetic algorithm only had to enumerate and score 600 combinations before it found the best scoring combination.
\n# GA parameters\nPOPULATION_SIZE = 500\nNUM_GENERATIONS = 20\nSELECTION_PRESSURE = 0.5\nMUTATION_RATE = 0.1\nMINIMIZE = True\nNUMBER_OF_REPEATS = 10\n\n# Create the virtual library\nvirtual_library_start = time.time()\nprint(\"Creating virtual library\")\nlibrary = build_virtual_library()\nvirtual_library_end = time.time()\nprint(\n f\"Virtual library created in {virtual_library_end - virtual_library_start:.2f} seconds\",\n)\n\n# Find the minimum logP in the library\nlogp_start = time.time()\nprint(\"Calculating logP values for the library\")\nall_clogp_values = [Descriptors.MolLogP(mol) for mol in tqdm(library.mol.to_list())]\nlogp_end = time.time()\nprint(f\"Minimum logP in the library: {min(all_clogp_values):.2f}\")\nprint(f\"LogP calculations took {logp_end - logp_start:.2f} seconds\")\n\n# Run the genetic algorithm\nprint(f\"Running genetic algorithm {NUMBER_OF_REPEATS} times\")\nga_times = []\nmax_generations = []\nfor _ in tqdm(range(NUMBER_OF_REPEATS)):\n ga_start = time.time()\n logp_ga_run = run_genetic_algorithm(\n building_blocks=BUILDING_BLOCKS,\n rnx=RXN,\n population_size=POPULATION_SIZE,\n num_generations=NUM_GENERATIONS,\n selection_pressure=SELECTION_PRESSURE,\n mutation_rate=MUTATION_RATE,\n scoring_function=calc_logp,\n minimize=MINIMIZE,\n early_stopping_value=min(all_clogp_values),\n )\n ga_end = time.time()\n ga_times.append(ga_end - ga_start)\n max_generations.append(logp_ga_run.generation.max())\nprint(\n f\"Genetic algorithm run took {np.average(ga_times):.2f} seconds on average, with a maximum of {np.max(max_generations)} generations\"\n)\n\n# Calculate improvement\nga_time = ga_end - ga_start\nvirtual_library_time = virtual_library_end - virtual_library_start\nlogp_time = logp_end - logp_start\nimprovement = (virtual_library_time + logp_time) / ga_time\nprint(\n f\"GA was {improvement:.2f} times faster than building and scoring the virtual library\"\n)\nCreating virtual library\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 1000000/1000000 [01:48<00:00, 9183.87it/s]\nVirtual library created in 120.60 seconds\n\nCalculating logP values for the library\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 132500/132500 [00:27<00:00, 4803.96it/s]\nMinimum logP in the library: -5.00\nLogP calculations took 27.59 seconds\n\nRunning genetic algorithm 10 times\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:14<00:00, 1.41s/it]\n\nGenetic algorithm run took 1.41 seconds on average, with a maximum of 12 generations\nGA was 81.90 times faster than building and scoring the virtual library\nA reference molecule was randomly selected from the library, and the genetic algorithm was used to recover it from the library by maximising the Tanimoto similarity. In the active learning experiments, this was failing, possibly due to the low number of compounds from within the library with high Tanimoto similarities. The genetic algorithm was able to recover the reference molecule every time, taking less than a second to do so each time.
\ndef calc_similarity(comparison_mols: list[Chem.Mol], ref_mol: Chem.Mol) -> list[float]:\n \"\"\"Calculate the Tanimoto similarity to a reference compound.\n\n Args:\n comparison_mols (list[Chem.Mol]): List of molecules to compare.\n ref_mol (Chem.Mol): The reference molecule.\n\n Returns:\n list[float]: List of similarity scores\n \"\"\"\n fpgen = AllChem.GetMorganGenerator()\n ref_fp = fpgen.GetFingerprint(ref_mol)\n comparison_fps = [fpgen.GetFingerprint(x) for x in comparison_mols]\n return [DataStructs.FingerprintSimilarity(ref_fp, x) for x in comparison_fps]\n\n\n# Score the virtual library\nprint(\"Scoring the virtual library\")\ntanimoto_start = time.time()\nall_tanimoto_values = calc_similarity(\n library.mol.to_list(), library.mol.sample().values[0]\n)\ntanimoto_end = time.time()\nprint(f\"Tanimoto calculations took {tanimoto_end - tanimoto_start:.2f} seconds\")\n\n# GA parameters\nPOPULATION_SIZE = 500\nNUM_GENERATIONS = 10\nSELECTION_PRESSURE = 0.5\nMUTATION_RATE = 0.1\nMINIMIZE = False\nNUMBER_OF_REPEATS = 10\n\n# Run the genetic algorithm\nprint(f\"Running genetic algorithm {NUMBER_OF_REPEATS} times\")\nga_times = []\nmax_generations = []\nfor _ in tqdm(range(NUMBER_OF_REPEATS)):\n # Choose random reference molecule\n ref_mol = library.mol.sample().values[0]\n\n # Run the genetic algorithm\n ga_start = time.time()\n tanimoto_ga_run = run_genetic_algorithm(\n building_blocks=BUILDING_BLOCKS,\n rnx=RXN,\n population_size=POPULATION_SIZE,\n num_generations=NUM_GENERATIONS,\n selection_pressure=SELECTION_PRESSURE,\n mutation_rate=MUTATION_RATE,\n scoring_function=partial(calc_similarity, ref_mol=ref_mol),\n minimize=MINIMIZE,\n early_stopping_value=1,\n )\n ga_end = time.time()\n ga_times.append(ga_end - ga_start)\n max_generations.append(tanimoto_ga_run.generation.max())\nprint(\n f\"Genetic algorithm run took {np.average(ga_times):.2f} seconds on average, with a maximum of {np.max(max_generations)} generations\"\n)\n\n# Calculate improvement\ntanimoto_time = tanimoto_end - tanimoto_start\nimprovement = tanimoto_time / np.average(ga_times)\nprint(\n f\"GA was {improvement:.2f} times faster than scoring the virtual library with Tanimoto similarity\"\n)\nScoring the virtual library\nTanimoto calculations took 7.41 seconds\nRunning genetic algorithm 10 times\n100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 10/10 [00:07<00:00, 1.29it/s]\n\nGenetic algorithm run took 0.77 seconds on average, with a maximum of 9 generations\nGA was 9.60 times faster than scoring the virtual library with Tanimoto similarity\nDeep Learning models have traditionally performed well on unstructured data such as text and images, but poorly on structured tabular data, and are usually outperformed by Gradient Boosted Decision Trees (GBDTs) on tabular chemical data. TabPFN (Tabular Prior-data Fitted Network) is a transformer-based foundation model for tabular data, pre-trained on millions of synthetic datasets to solve supervised learning tasks, with state-of-the-art performance on benchmarks. But does it work for cheminformatics?
\nA recent paper published in Nature, Accurate predictions on small data with a tabular foundation model by Hollmann et al., describes a new Deep Learning (DL) algorithm for making predictions on tabular data, TabPFN (Tabular Prior-data Fitted Network). TabPFN is a transformer-based foundation model for tabular data that operates via in-context learning (ICL), enabling it to train and predict on an entire dataset in a single forward pass.
\nThe authors generated over 100 million synthetic datasets using structural causal models (SCMs). These models simulate causal relationships and mimic real-world tabular data challenges, such as non-linear relationships, missing values, outliers, and diverse feature types (e.g., categorical, ordinal, numerical). These datasets sampled high-level hyperparameters (such as dataset size and number of features), both classification and regression tasks, and Gaussian noise was added to mimic real-world complexities.
\nThese synthetic datasets were used to pre-train a transformer model. During the training, parts of the datasets were masked, and the model is trained to predict masked target values in synthetic datasets, given features and labelled samples as context. The parameters of the neural network were updated until the predictions matched the masked values. Through this training, the model learns to fill in missing data from a dataset. The model learns a generic prediction algorithm that approximates Bayesian inference for the synthetic data prior, enabling robust handling of unseen datasets.
\nWhen making predictions TabPFN uses ICL and processes an entire dataset (both labelled and test samples) in one pass, performing training and inference simultaneously.
\nThe Therapeutic data commons (TDC) aims to help development of AI/ML tools for therapeutic science by providing datasets and curated benchmarks to assess the performance of new methods. The TDC ADMET benchmark group contains 22 datasets for molecular property prediction, ranging from 475 to 13,130 entries, with both classification and regression tasks. The datasets contain SMILES strings for each chemical compound, and a target variable to be predicted. These can all be easily downloaded using their Python library.
\nTabPFN is designed for small tabular datasets, and the default parameters limit to training on 10,000 entries with 500 features. If a dataset had more than 10,000 entries in the training and validation sets, a random selection of 10,000 entries was used. For each dataset, the training and validation data were combined, and the test set put aside. For all entries, the 210 RDKit descriptors were calculated and used as the features for training. Predictions were made on the test set and compared to the true values to evaluate the performance of TabPFN. The training was repeated using molecular fingerprints (MACCS keys and ECPF folded to 500 bits), but the performance was found to be lower than using RDKit calculated properties.
\nNOTE: Due to memory limits on my computer, I had to limit to datasets with fewer than 1,800 entries in the training data. Training models for the larger datasets is ongoing.
\nUsing the RDKit calculated descriptors as features, TabPFN come in the top 10 models for all TDC datasets apart from \"Clearance_Hepatocyte_Az\". It comes 3rd for the \"Vdss_Lombardo\" dataset, 2nd for \"Caco2_Wang\", \"Pgp_Broccatelli\", and \"Bbb_Martins\" datasets, and is the highest performing model for the \"Clearance_Microsome_Az\" dataset. The high performing datasets include both classification and regression tasks and doesn't seem to be a clear link between classification performance and imbalanced datasets. \"Cyp2C9_Substrate_Carbonmangels\" is somewhat imbalanced (19.3% positive - TDC rank: 10th) but so are \"HIA_Hou\" (11.1% negative - TDC rank: 5th), and \"Bioavailability_Ma\", (22.9% negative - TDC rank: 5th). The strongest link seems to be between dataset size and performance, with all the highest performing models having >900 entries.
\nThe training time is an average of five repeats, using WSL on a Windows machine with 8 GB RAM allocated to WSL.
\n| Dataset | \nSize | \nTask | \nMetric | \nTraining time (min) | \nTabFPN performance | \nCurrent TDC best performance | \nTabPFN TDC leaderboard rank | \n
|---|---|---|---|---|---|---|---|
| Caco2_Wang | \n906 | \nRegression | \nMAE | \n8.25 | \n0.282 \u00b1 0.005 | \n0.276 \u00b1 0.005 | \n2nd | \n
| HIA_Hou | \n578 | \nClassification | \nAUROC | \n3.17 | \n0.987 \u00b1 0.001 | \n0.990 \u00b1 0.002 | \n5th | \n
| Pgp_Broccatelli | \n1218 | \nClassification | \nAUROC | \n12.77 | \n0.936 \u00b1 0.004 | \n0.938 \u00b1 0.002 | \n2th | \n
| Bioavailability_Ma | \n640 | \nClassification | \nAUROC | \n4.32 | \n0.735 \u00b1 0.016 | \n0.753 \u00b1 0.000 | \n5th | \n
| Bbb_Martins | \n2030 | \nClassification | \nAUROC | \n35.51 | \n0.917 \u00b1 0.003 | \n0.920 \u00b1 0.006 | \n2nd | \n
| Vdss_Lombardo | \n1130 | \nRegression | \nSpearman | \n12.91 | \n0.693 \u00b1 0.004 | \n0.713 \u00b1 0.007 | \n3rd | \n
| Cyp2D6_Substrate_Carbonmangels | \n667 | \nClassification | \nAUPRC | \n4.82 | \n0.714 \u00b1 0.009 | \n0.736 | \n6th | \n
| Cyp3A4_Substrate_Carbonmangels | \n670 | \nClassification | \nAUROC | \n3.97 | \n0.641 \u00b1 0.004 | \n0.667 \u00b1 0.019 | \n7th | \n
| Cyp2C9_Substrate_Carbonmangels | \n669 | \nClassification | \nAUPRC | \n4.28 | \n0.400 \u00b1 0.013 | \n0.441 \u00b1 0.033 | \n10th | \n
| Half_Life_Obach | \n667 | \nRegression | \nSpearman | \n4.20 | \n0.546 \u00b1 0.013 | \n0.576 \u00b1 0.025 | \n6th | \n
| Clearance_Microsome_Az | \n1102 | \nRegression | \nSpearman | \n12.71 | \n0.632 \u00b1 0.006 | \n0.630 \u00b1 0.010 | \n1st | \n
| Clearance_Hepatocyte_Az | \n1213 | \nRegression | \nSpearman | \n11.27 | \n0.391 \u00b1 0.004 | \n0.536 \u00b1 0.02 | \n>10th | \n
| Herg | \n655 | \nClassification | \nAUROC | \n3.54 | \n0.850 \u00b1 0.002 | \n0.880 \u00b1 0.002 | \n6th | \n
| Dili | \n475 | \nClassification | \nAUROC | \n1.92 | \n0.910 \u00b1 0.005 | \n0.925 \u00b1 0.005 | \n6th | \n
On some of the TDC datasets, TabPFN exhibits state-of-the-art performance out of the box. Regularly outperforming other DL methods such as graph neural networks (GNNs) designed specifically for cheminformatics. It can learn complex relationships between features and the target on small datasets due to the pre-training on synthetic data. TabPFN uses a scikit-learn style API, making it very easy to use and integrate into existing workflows. Since TabPFN approximates Bayesian inference, it can include important features such as uncertainty estimation possible.
\nTabPFN is limited to 10,000 entries with 500 features as a default, as this is the limit of the synthetic data it was trained on. It is possible to use datasets larger than this using the ignore_pretraining_limits=True parameter, but this may lead to very long training times. The computational requirements for TabPFN scale quadratically with the number of samples (n) and the number of features (m), i.e. O(n2 + m2) so training and predicting on larger datasets gets increasing longer. The training time can be long, though not significantly worse than other DL methods such as Chemprop in my experience. The single pass training and prediction architecture is not great for situations when you need to make multiple predictions, but the authors have included an option to cache the effect of the training data (using the fit_mode=\"fit_with_cache\" parameter) when training the model, making the training and prediction API much more like a standard scikit-learn ML model.
TabPFN is definitely a method to consider when building QSAR models, especially when you have a small dataset. With more time (and computing power), I'd be interested to see if combining the RDKit calculated descriptors with a molecular fingerprint (i.e. adding MACCS keys would make 376 feature columns) improves performance. In my experience sometimes combining the local description of molecules using molecular fingerprints with global descriptions of molecules using calculated descriptors creates a very good embedding of your molecules.
\nIt's possible to fine tune TabPFN with specialised datasets. Perhaps fine tuning with a large number of chemical datasets, or synthetic datasets of calculated molecular properties would create a chemical tabular foundation model (ChemTabPFN?), capable of accurately predicting chemical property-structure relationships.
\n
No data science or cheminformatics today!
\nI usually cycle or get public transport to work, but occasionally I do have to drive. The traffic in T\u0101maki Makaurau (Auckland) is very variable, especially on rainy or windy days, so I often find myself checking Google Maps every 5 minutes after 4 pm to work out when I need to leave to get home. As with anything that I have to do repeatedly, I decided to automate it.
\nMy work computer runs Windows so I use WSL, which adds a few extra complications. All the code needed can be found in this repository.
\n(I've also been experimenting with Google Gemini as a coding assistant. This wasn't completely vibe-coding, but as this was a quick personal project, I was much less vigorous in checking the code it generated compared to more important work!)
\nAll the Python code is kept in a file called traffic-widget.py which is stored on my WSL disk, I first needed to import the necessary libraries:
import datetime\nimport json\nimport os\nimport subprocess\n\nimport requests\nfrom dotenv import load_dotenv\nGoogle maps doesn't seem to have a completely free API. They do offer some free usage, but you still have to sign up and give credit card details, which is always a worry in case you accidentally go over the free limits (or accidentally leak your API key to the internet). TomTom on the other hand does have a free API, you have to sign up and get an API key, and it comes with plenty of free requests. After signing up and getting an API key, I first created a TomTomAPI Python class, this uses the API key to make calls to the route calculation endpoint, providing the GPS coordinates for the start and end of your route, and returns the travel time.
\nclass TomTomAPI:\n \"\"\"Encapsulates interactions with the TomTom Routing API.\"\"\"\n\n BASE_URL = \"https://api.tomtom.com/routing/1/calculateRoute\"\n\n def __init__(self, api_key: str):\n \"\"\"Initialize the TomTomAPI with the API key.\"\"\"\n self.api_key = api_key\n\n def get_travel_time(\n self, start_lat: float, start_lon: float, end_lat: float, end_lon: float\n ) -> int | None:\n \"\"\"Calculate the travel time between two points.\n\n Args:\n start_lat (float): Latitude of the starting point.\n start_lon (float): Longitude of the starting point.\n end_lat (float): Latitude of the destination point.\n end_lon (float): Longitude of the destination point.\n\n Returns:\n int: The travel time in seconds, or None if an error occurs.\n \"\"\"\n start_point = f\"{start_lat},{start_lon}\"\n end_point = f\"{end_lat},{end_lon}\"\n url = f\"{self.BASE_URL}/{start_point}:{end_point}/json?key={self.api_key}\"\n\n try:\n response = requests.get(url)\n response.raise_for_status()\n data = response.json()\n travel_time_seconds = data[\"routes\"][0][\"summary\"][\"travelTimeInSeconds\"]\n return travel_time_seconds\n\n except requests.exceptions.RequestException as e:\n print(f\"TomTom API Error: {e}\")\n return None\n except (json.JSONDecodeError, KeyError) as e:\n print(f\"TomTom API Error: Invalid response format or missing data: {e}\")\n return None\nTo display the notifications I used BurntToast, a Windows PowerShell module for displaying Toast Notifications. To install, I needed to run PowerShell as admin and enter:
\nInstall-Module -Name BurntToast\nBurntToast can be called from a PowerShell file (.ps1). The PowerShell file to display notifications is called show_notification.ps1 and is stored on my Windows disk. The -Sound Alarm5 adds a sound to the notification and makes it last longer and -AppLogo gives an icon beside the notification.
param(\n [string]$Title,\n [string]$Message\n)\n\n$ImagePath = 'C:\\Path\\to\\icon.jpg'\n\nNew-BurntToastNotification -Text $Title, $Message -AppLogo $ImagePath -Sound Alarm5\nTo allow the script to run, it may be necessary to run this command in PowerShell as admin:
\nSet-ExecutionPolicy -ExecutionPolicy RemoteSigned -Scope CurrentUser\nBack in the Python file (traffic-widget.py), I next needed a class to run the PowerShell file to display the notification.
class WindowsNotifier:\n \"\"\"Handles displaying Windows toast notifications.\"\"\"\n\n def __init__(self, powershell_path: str, powershell_script: str):\n \"\"\"Initialize the WindowsNotifier with the paths to PowerShell.\"\"\"\n self.powershell_path = powershell_path\n self.powershell_script = powershell_script\n\n def show_notification(self, title: str, message: str):\n \"\"\"Show a Windows toast notification.\n\n Args:\n title (str): The title of the notification.\n message (str): The body of the notification.\n \"\"\"\n try:\n subprocess.run(\n [\n self.powershell_path,\n \"-ExecutionPolicy\",\n \"Bypass\",\n \"-File\",\n self.powershell_script,\n \"-Title\",\n title,\n \"-Message\",\n message,\n ],\n check=True,\n capture_output=True,\n text=True,\n )\n except subprocess.CalledProcessError as e:\n print(f\"Error showing notification: {e}\")\n print(f\"PowerShell Output:\\n{e.stderr}\")\n except FileNotFoundError:\n print(f\"Error: PowerShell executable or script not found. Check the paths.\")\nThis function formats the data from the TomTomAPI to be clearer.
\ndef format_travel_time(label: str, travel_time: int | None) -> str:\n \"\"\"Format the travel time into a readable string.\n\n Args:\n label (str): The label for the travel time (e.g., \"Home\").\n travel_time (int | None): The travel time in seconds, or None if an error occurred.\n\n Returns:\n str: The formatted travel time string.\n \"\"\"\n now = datetime.datetime.now()\n arrival_time = (\n now + datetime.timedelta(seconds=travel_time) if travel_time else None\n )\n if travel_time is not None:\n return f\"{label}: {travel_time / 60:.1f} minutes (arrive at {arrival_time.strftime('%H:%M')})\"\n else:\n return f\"{label}: ERROR minutes.\"\nTo prevent sharing sensitive data such as my API key and home address, I stored these in a .env file. This also stores the paths to Windows PowerShell and the PowerShell script. Since this is working between two operating systems, the paths are slightly more complicated than usual. The PowerShell path is the path to the powershell.exe file on your Windows disk, from your Linux environment. This is usually something like: /mnt/c/Windows/System32/WindowsPowerShell/v1.0/powershell.exe. The PowerShell script is in Windows format with escaped backslashes.
API_KEY = \"your_api_key\"\nWORK_LATITUDE = \"work_latitude\"\nWORK_LONGITUDE = \"work_longitude\"\nHOME_LATITUDE = \"home_latitude\"\nHOME_LONGITUDE = \"home_longitude\"\n\nPOWERSHELL_PATH = \"/mnt/c/path/to/powershell.exe\"\nPOWERSHELL_SCRIPT = \"C:\\\\windows\\\\path\\\\to\\\\show_notification.ps1\"\nTo run the Python script I used a conda environment. It can be created and activated with:
\nconda env create -f environment.yml\nconda activate traffic-widget\nFinally I needed some Python code to run the whole process:
\nif __name__ == \"__main__\":\n load_dotenv()\n api_key = os.environ[\"API_KEY\"]\n work_lat = float(os.environ[\"WORK_LATITUDE\"])\n work_lon = float(os.environ[\"WORK_LONGITUDE\"])\n home_lat = float(os.environ[\"HOME_LATITUDE\"])\n home_lon = float(os.environ[\"HOME_LONGITUDE\"])\n powershell_path = os.environ.get(\"POWERSHELL_PATH\")\n powershell_script = os.environ.get(\"POWERSHELL_SCRIPT\")\n if not powershell_script:\n print(\"Error: POWERSHELL_SCRIPT environment variable not set.\")\n exit(1)\n\n tomtom_api = TomTomAPI(api_key)\n notifier = WindowsNotifier(powershell_path, powershell_script)\n\n home_travel_time = tomtom_api.get_travel_time(\n work_lat, work_lon, home_lat, home_lon\n )\n\n message_lines = [\n format_travel_time(\"Home\", home_travel_time),\n ]\n notification_message = \"\\n\".join(message_lines)\n\n notifier.show_notification(\"Driving times:\", notification_message)\nManually running the script with python traffic-widget.py should cause the following pop-up:
To make the script run automatically, back in WSL, I ran:
\ncrontab -e\nAnd added the details for the cron job. I wanted mine to run every 5 minutes from 4-5 pm on weekdays.
\n0-55/5 16 * * 1-5 /home//miniconda3/envs/traffic-widget/bin/python /home//path/to/traffic-widget/traffic-widget.py\n An unwelcome cyclone may have cancelled my hiking plans, but it gave me the opportunity to delve into the fascinating world of drug repurposing. This blog post explores how identifying new uses for existing drugs can dramatically accelerate and reduce the cost of bringing vital treatments to patients, and how recent advances in machine learning can further streamline this process.
\nOver the Easter long weekend, T\u0101maki Makaurau (Auckland) was visited by Cyclone Tam, which meant my hiking trip was cancelled. At least it gave me some time to work on a few projects and finish this blog post I'd been meaning to do for a while!
\nTraditional drug discovery is a marathon, often taking over a decade and requiring an investment of more than a billion dollars for each successful drug. This lengthy and costly process highlights the urgent need for more efficient strategies to combat disease. One such promising strategy is drug repurposing: identifying new therapeutic uses for existing drugs or those already undergoing clinical trials. By bypassing much of the early-stage development, and sometimes significant parts of clinical trials, drug repurposing offers the potential to significantly reduce both the time and financial cost associated with bringing new treatments to patients. Drug repurposing is ideal for addressing therapeutic needs in disease areas where financial incentives are limited (such as neglected tropical diseases and antibiotic development), or for rapid deployment in situations demanding swift treatment interventions for imminent public health emergencies, such as during the COVID-19 pandemic.
\nOne well-known example of drug repurposing is sildenafil, originally developed to treat high blood pressure and angina (chest pain), and famously repurposed for erectile dysfunction after clinical trials revealed this unexpected side effect. Another is thalidomide, developed to treat morning sickness and found to cause birth defects, which is now used to treat cancers and skin disorders associated with leprosy. The crucial role of drug repurposing was further highlighted during the recent COVID-19 pandemic, where the corticosteroid dexamethasone emerged as a vital intervention in reducing mortality among severely ill patients.
\nMachine learning offers a novel method for drug repurposing. Using high-quality activity data for a target of interest, a machine learning model can be trained to predict active compounds based on their molecular structure. This trained model can then be used to screen libraries of compounds for potential repurposing.
\nThe Drug Repurposing Hub is a curated and annotated collection of drugs that have already received FDA approval or are currently undergoing clinical trials. Maintained by the Broad Institute, this resource serves as both a virtual library, providing comprehensive information on each compound, and a physical library, offering access to compound plates for experimental screening.
\nThis post was inspired by the work of the Collins Lab at MIT, who used deep learning techniques to identify potential compounds for repurposing as treatments for antimicrobial-resistant bacteria. Their research addresses the escalating crisis of antibiotic resistance, a looming threat projected to cause 10 million deaths annually by 2050.
\nTo tackle this challenge, the researchers assembled a primary training set of 2,335 molecules, experimentally screening them against E. coli BW25113 and identifying 120 compounds with antimicrobial activity. With the goal of creating a robust and generalizable model, they ensured the training set was structurally diverse. This data was then used to train a directed message-passing neural network (dMPNN) to create a classification model capable of predicting antimicrobial activity based on a molecule's chemical structure. The trained model was subsequently used to make predictions on virtual libraries of compounds, including the Drug Repurposing Hub.
\nNinety-nine molecules from the Drug Repurposing Hub that were predicted to be active were experimentally tested, with 51 displaying growth inhibition against E. coli. This process led to the identification of a particularly promising compound, the c-Jun N-terminal kinase inhibitor SU3327, which the researchers named halicin. Notably, halicin possesses a structure distinctly different from conventional antibiotics and demonstrated potent inhibitory activity against E. coli growth.
\nCentral to this research is the open-source software package Chemprop. Developed primarily by researchers at MIT using the PyTorch framework, Chemprop harnesses the capabilities of message-passing neural networks (MPNNs) for molecular property prediction. Chemprop can be used both as a command-line interface (CLI) and a Python API. It implements a d-MPNN architecture, a type of graph neural network particularly effective for extracting relationships from molecular structures represented as graphs. An advantage of Chemprop is its ability to learn directly from molecular representations, such as SMILES strings converted into molecular graphs, eliminating the need for manual feature engineering, which can also introduce biases.
\nLacking access to a physical lab for experimental data, I needed an alternative approach for data gathering. For this, I used ChEMBL, a freely accessible and meticulously curated chemical database housing bioactivity data for drug-like molecules. It is maintained by the European Bioinformatics Institute (EBI), part of the European Molecular Biology Laboratory (EMBL) in the UK. Using the ChEMBL Webservice, I was able to specifically download data relevant to the target of interest. This retrieved data then served as the foundational dataset for training a machine learning model, which was subsequently employed to virtually screen the compounds within The Drug Repurposing Hub.
\nimport subprocess\nfrom pathlib import Path\n\nimport pandas as pd\nfrom chembl_webresource_client.new_client import new_client\nfrom rdkit import Chem\nfrom rdkit.Chem import Draw\nfrom tqdm import tqdm\n\nPath(\"data\").mkdir(exist_ok=True)\nPath(\"models\").mkdir(exist_ok=True)\nUsing the ChEMBL Webservice, I downloaded the data associated with a target for training a machine learning model using the three functions below. The first searches ChEMBL with a target name, and returns the target ChEMBL ID of the most similar target. The second downloads the activity data associated that the target ChEMBL ID, along with the molecular ChEMBL IDs. The third returns the canonical SMILES for the molecular ChEMBL IDs.
\ntarget = new_client.target\nactivity = new_client.activity\nmolecule = new_client.molecule\n\n\ndef get_target_chembl_id(target_name: str) -> str:\n \"\"\"Get the ChEMBL ID of a target given its name.\"\"\"\n print(f'Searching for target: \"{target_name}\"')\n target_query = target.search(target_name)\n name, chembl_id = target_query[0][\"pref_name\"], target_query[0][\"target_chembl_id\"]\n print(f\"Target name: {name}, ChEMBL ID: {chembl_id}\")\n return chembl_id\n\n\ndef get_target_activity(chembl_id: str) -> pd.DataFrame:\n \"\"\"Get the activity data for a target given its ChEMBL ID.\"\"\"\n print(f\"Getting activity data for {chembl_id}\")\n activities = activity.filter(\n target_chembl_id=chembl_id,\n standard_type=\"IC50\",\n relation=\"=\",\n standard_units=\"nM\",\n ).only(\n \"molecule_chembl_id\",\n \"standard_value\",\n )\n print(f\"Found {len(activities)} activities\")\n activities_df = pd.DataFrame(\n list(tqdm(activities[:100], desc=\"Processing Activities\"))\n )[[\"molecule_chembl_id\", \"standard_value\"]]\n activities_df[\"standard_value\"] = pd.to_numeric(activities_df[\"standard_value\"])\n return activities_df\n\n\ndef get_molecule_data(molecule_chembl_ids: list[str]) -> pd.DataFrame:\n \"\"\"Get the molecule data for a list of ChEMBL IDs.\"\"\"\n print(f\"Getting molecule data for {len(molecule_chembl_ids)} ChEMBL IDs\")\n compounds_provider = molecule.filter(\n molecule_chembl_id__in=molecule_chembl_ids\n ).only(\"molecule_chembl_id\", \"molecule_structures\")\n compounds_df = pd.json_normalize(\n list(tqdm(compounds_provider, desc=\"Processing Compounds\"))\n )[[\"molecule_chembl_id\", \"molecule_structures.canonical_smiles\"]].rename(\n columns={\"molecule_structures.canonical_smiles\": \"canonical_smiles\"}\n )\n return compounds_df\nChemprop takes the data in the form of a table containing one column of SMILES strings, and one columns of activities. For this example I am creating a classification model, so the activity column is a binary classification (1 or 0). The first function takes the activity data and molecule data and merges them on the molecule ChEMBL ID, before cleaning the data up and creating the binary variable based on a cutoff. This cutoff can be changed depending on the target and dataset. I've always found the CLI usage of Chemprop easier, so the data is saved as a CSV file for training.
\ndef create_training_data(\n activity_data: pd.DataFrame, molecule_data: pd.DataFrame, nm_cutoff: float\n) -> pd.DataFrame:\n \"\"\"Create the training data for the model.\"\"\"\n print(\"Creating training data\")\n training_data = (\n (\n pd.merge(\n left=activity_data,\n right=molecule_data,\n on=\"molecule_chembl_id\",\n how=\"left\",\n )\n .dropna(how=\"any\")\n .drop(columns=\"molecule_chembl_id\")\n )\n .groupby(\"canonical_smiles\", as_index=False)\n .agg(\"mean\")\n )\n training_data[\"activity\"] = (training_data[\"standard_value\"] < nm_cutoff).astype(\n int\n )\n actives = training_data[\"activity\"].sum()\n inactives = len(training_data) - actives\n percentage = 100 * actives / len(training_data)\n print(f\"Actives: {actives} ({percentage:.2f}%), Inactives: {inactives}\")\n return training_data.drop(columns=\"standard_value\")\n\n\ndef save_training_data(\n training_data: pd.DataFrame, target_name: str, output_dir: Path\n) -> None:\n \"\"\"Save the training data to a CSV file.\"\"\"\n print(f\"Saving training data to {output_dir}\")\n output_dir.mkdir(parents=True, exist_ok=True)\n training_data.to_csv(output_dir / f\"{target_name}_training_data.csv\", index=False)\n print(f\"Training data saved to {output_dir / f'{target_name}_training_data.csv'}\")\nThis function downloaded the Drug Repurposing Hub dataset, tidies the data up and saves it for making predictions using Chemprop.
\ndef canon_smiles(smiles: str) -> str | None:\n \"\"\"Convert SMILES string to canonical SMILES.\"\"\"\n mol = Chem.MolFromSmiles(smiles)\n if mol is None:\n return None\n return Chem.MolToSmiles(mol, canonical=True)\n\n\ndef download_drug_repurposing_hub() -> None:\n \"\"\"Download the Drug Repurposing Hub dataset.\"\"\"\n if Path(\"data/drug_repurposing_hub.csv\").exists():\n print(\"Drug Repurposing Hub dataset already downloaded\")\n return None\n url = \"https://storage.googleapis.com/cdot-general-storage/repurposing_samples_20240610.txt\"\n drug_repurposing_df = pd.read_csv(\n url,\n sep=\"\\t\",\n skiprows=9,\n )\n drug_repurposing_df[\"smiles\"] = drug_repurposing_df[\"smiles\"].apply(canon_smiles)\n drug_repurposing_df.dropna(subset=[\"smiles\"], inplace=True)\n drug_repurposing_df.drop_duplicates(subset=[\"smiles\"], keep=\"first\", inplace=True)\n drug_repurposing_df.to_csv(\"data/drug_repurposing_hub.csv\", index=False)\n print(\"Drug Repurposing Hub dataset downloaded and processed\")\n print(f\"Number of unique SMILES: {len(drug_repurposing_df['smiles'].unique())}\")\nAs I mentioned above, I find the CLI usage of Chemprop easier (especially since the move from v1 to v2), so here I'm using subprocess to run Chemprop. These two functions train a Chemprop model, and make predictions on the Drug Repurposing Hub data.
\ndef train_chemprop_model(target_name: str, data_dir: str, output_dir: str) -> None:\n \"\"\"Train a ChemProp model using the training data.\"\"\"\n print(f\"Training ChemProp model for {target_name}\")\n chemprop_train_args = [\n \"chemprop\",\n \"train\",\n \"--data-path\",\n f\"{data_dir}/{target_name.replace(\" \", \"_\")}_training_data.csv\",\n \"--task-type\",\n \"classification\",\n \"--save-dir\",\n f\"{output_dir}/{target_name.replace(\" \", \"_\")}\",\n \"--split-type\",\n \"scaffold_balanced\",\n ]\n subprocess.run(chemprop_train_args, check=True)\n\n\ndef predict_with_chemprop(target_name: str, preds_dir: str, model_dir: str) -> None:\n \"\"\"Predict using the ChemProp model.\"\"\"\n print(f\"Predicting with ChemProp model for {target_name}\")\n chemprop_predict_args = [\n \"chemprop\",\n \"predict\",\n \"--test-path\",\n \"data/drug_repurposing_hub.csv\",\n \"--model-paths\",\n f\"{model_dir}/{target_name.replace(\" \", \"_\")}\",\n \"--smiles-columns\",\n \"smiles\",\n \"--preds-path\",\n f\"{preds_dir}/{target_name.replace(\" \", \"_\")}_predictions.csv\",\n ]\n subprocess.run(chemprop_predict_args, check=True)\nOnce we've scored the compounds, we want to take a look at the top scoring compounds.
\ndef visualize_top_predictions(search: str) -> None:\n \"\"\"Visualize the top predictions.\"\"\"\n predictions_df = (\n pd.read_csv(f\"data/{search.replace(\" \", \"_\")}_predictions.csv\")\n .sort_values(\"activity\", ascending=False)\n .head(9)\n )\n return Draw.MolsToGridImage(\n [Chem.MolFromSmiles(s) for s in predictions_df[\"smiles\"]],\n legends=[f\"{a:.2f}\" for a in predictions_df[\"activity\"]],\n molsPerRow=3,\n )\nThis function runs the whole process for a specified target search.
\ndef run_repurposing_pipeline(search: str, nm_cutoff=25) -> None:\n \"\"\"Run the drug repurposing pipeline.\"\"\"\n download_drug_repurposing_hub()\n target_chembl_id = get_target_chembl_id(search)\n activity_data = get_target_activity(target_chembl_id)\n molecule_data = get_molecule_data(activity_data[\"molecule_chembl_id\"].tolist())\n training_data = create_training_data(\n activity_data, molecule_data, nm_cutoff=nm_cutoff\n )\n save_training_data(training_data, search.replace(\" \", \"_\"), Path(\"data\"))\n train_chemprop_model(target_name=search, data_dir=\"data\", output_dir=\"models/\")\n predict_with_chemprop(target_name=search, preds_dir=\"data\", model_dir=\"models/\")\nMalaria kills around half a million people per year, with 76% of global malaria deaths in children under 5 years old. It is considered a neglected tropical disease (NTD) due to its significant impact on global health, particularly in low-income countries. It disproportionately affects some of the world's poorest people, and is often overlooked in global health efforts. Drug repurposing could offer alternative therapies significantly cheaper than developing an entirely new drug. Here I look to target the most deadly species of parasite that carries malaria - Plasmodium falciparum. After training the Chemprop model and screening the Drug Repurposing Hub, some compounds are flagged as potential hits and could be investigated experimentally.
\nOf the top 10 predictions, pyronaridine and KAF-156 are both malaria treatments, and were both in the training data, but interestingly pyronaridine was not classified as a hit in the training data.\nCarfilzomib, TPPS4, nanchangmycin, oprozomib, and lindane have all been investigated for treating malaria, so the model appears to be picking out the right compounds. Abemaciclib appears to have been recently investigated for repurposing with malaria. GS-967 appears to have no link to malaria yet, but maybe should be investigated experimentally!
\nsearch = \"Plasmodium falciparum\"\n\nrun_repurposing_pipeline(search)\nvisualize_top_predictions(search)\nDrug Repurposing Hub dataset downloaded and processed\nNumber of unique SMILES: 6734\nSearching for target: \"Plasmodium falciparum\"\nTarget name: Plasmodium falciparum, ChEMBL ID: CHEMBL364\nGetting activity data for CHEMBL364\nFound 45314 activities\nProcessing Activities: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 45314/45314 [00:02<00:00, 15928.83it/s]\nGetting molecule data for 45314 ChEMBL IDs\nProcessing Molecules: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 21821/21821 [1:32:46<00:00, 3.92it/s]\nCreating training data\nActives: 2677 (12.30%), Inactives: 19082\nSaving training data to data\nTraining data saved to data/Plasmodium_falciparum_training_data.csv\nTraining ChemProp model for Plasmodium falciparum\nPredicting with ChemProp model for Plasmodium falciparum\n
The COVID-19 pandemic has caused the deaths of over 7 million people so far, and in 2020 more than half of the world's population was under some form of lockdown. With the speed that COVID spread around the globe, there was no time to develop new drug from scratch, and drug repurposing was explored in clinics to find treatments. In this example, I didn't find any promising compounds in the Drug Repurposing Hub, but with a much smaller dataset than for malaria, the machine learning model is potentially much less powerfull.
\nsearch = \"SARS-CoV-2\"\n\nrun_repurposing_pipeline(search, nm_cutoff=500)\nvisualize_top_predictions(search)\nDrug Repurposing Hub dataset already downloaded\nSearching for target: \"SARS-CoV-2\"\nTarget name: SARS-CoV-2, ChEMBL ID: CHEMBL4303835\nGetting activity data for CHEMBL4303835\nFound 802 activities\nProcessing Activities: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 802/802 [00:00<00:00, 12406.99it/s]\nGetting molecule data for 802 ChEMBL IDs\nProcessing Compounds: 100%|\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588\u2588| 490/490 [00:00<00:00, 6081.22it/s]\nCreating training data\nActives: 71 (14.55%), Inactives: 417\nSaving training data to data\nTraining data saved to data/SARS-CoV-2_training_data.csv\nTraining ChemProp model for SARS-CoV-2\nPredicting with ChemProp model for SARS-CoV-2\n
Whilst neural networks (NNs) have done amazing things with unstructured data such as text and images, they've traditionally been outperformed on tabular data by Gradient-Boosted Decision Trees (GBDTs), although recent advances such as TabPFN and TabICL suggest that the performance gap may have closed. I've written a bit about TabPFN for chemical datasets here.
\nOne of the most popular deep learning architectures for chemical property prediction is Chemprop, developed by researchers at MIT. A Chemprop model is made up of two NNs. The first is a directed Message Passing Neural Network (d-MPNN), a type of Graph Convolutional Neural Network (GCNN), that takes a graphical representation of a molecule and converts it to a molecular embedding, a vector that describes the original molecule. This molecular embedding is then put through a Feed-Forward Neural Network (FFN), a fully connected NN used to make the final prediction. The model is trained end-to-end, with both NNs being updated at the same time to minimise a loss function.
\nOnce the model is trained, the d-MPNN can be used to calculate a molecular embedding for a collection of molecules. This is a learned representation optimised for the specific chemical prediction problem. This learned representation is a vector with a length equal to the number of input nodes in the FFN. This is very similar to a traditional molecular fingerprint, such as a Morgan fingerprint.
\nA collection of these learned fingerprints forms a 2-dimensional array, which is tabular data. Since GBDTs regularly outperform NNs on tabular data, can we improve the performance of Chemprop by replacing the FFN with a GBDT after training the d-MPNN?
\nWhilst writing this post, a paper was published in The Journal of Chemical Information and Modeling by Pat Walters et al., Practically Significant Method Comparison Protocols for Machine Learning in Small Molecule Drug Discovery. This is probably my favourite cheminformatics paper so far this year, with comprehensive guidelines on how to compare machine learning models in a statistically robust way. I've tried to follow these recommendations when comparing models here.
\nBelow I compare three models: the first is a Random Forest (this is not gradient-boosted, but provides a good baseline without hyperparameter tuning) trained on Morgan fingerprints; the second is a standard Chemprop model; and the third generates learned fingerprints using the Chemprop d-MPNN and then trains a Random Forest model on these. I used two datasets, one for a classification problem and one for a regression problem. The classification data is the Bile Salt Export Pump dataset and the regression data is the Polaris LogS dataset, both ideas 'borrowed' from Practical Cheminformatics posts.
\nimport subprocess\nfrom collections import defaultdict\nfrom pathlib import Path\nfrom typing import Callable\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom numpy.typing import NDArray\nfrom rdkit import Chem\nfrom rdkit.Chem import rdFingerprintGenerator\nfrom scipy.stats import f_oneway\nfrom sklearn.base import BaseEstimator\nfrom sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\nfrom sklearn.metrics import (\n matthews_corrcoef,\n mean_absolute_error,\n r2_score,\n root_mean_squared_error,\n)\nfrom sklearn.model_selection import ShuffleSplit\nfrom statsmodels.stats.multicomp import pairwise_tukeyhsd\nfrom tqdm.notebook import tqdm\n\nfrom visualisation import confusion_matricies, scatter_plots\nThe first function prepare_data takes the path to a CSV file, creates RDKit molecular objects, removes any that are causing issues, and calculates the Morgan fingerprints for each model.
The second, save_data_for_chemprop, formats the data and saves it for training Chemprop models.
The next three functions: train_chemprop_model, make_chemprop_predictions, and calculate_chemprop_fingerprints, all use subprocess to call the Chemprop CLI. I've generally found it easier to use than the Chemprop Python module.
The next two: train_morgan_sklearn_model and train_chemprop_sklearn_model, use Scikit-learn to train Random Forest models, the first using the previously generated Morgan fingerprints, and the second uses the Chemprop learned fingerprints.
Finally, gather_predictions uses all of the trained models to make predictions on a test dataset.
def prepare_data(\n data_path: Path, smiles_col: str, target_col: str\n) -> tuple[pd.DataFrame, pd.DataFrame]:\n \"\"\"Load and clean data, generate Morgan fingerprints.\n\n Args:\n data_path (Path): Path to CSV file containing data.\n smiles_col (str): Name of the column containing SMILES strings.\n target_col (str): Name of the column containing target variable.\n\n Returns:\n tuple[pd.DataFrame, pd.DataFrame]: Cleaned dataframe and dataframe of Morgan\n fingerprints.\n \"\"\"\n df = pd.read_csv(data_path)\n df[\"ROMol\"] = df[smiles_col].apply(lambda x: Chem.MolFromSmiles(x))\n df.dropna(subset=[\"ROMol\"], inplace=True)\n mfpgen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)\n fps = [mfpgen.GetFingerprint(mol) for mol in df[\"ROMol\"]]\n fps_df = pd.DataFrame(np.array(fps))\n df = df[[smiles_col, target_col]]\n return df, fps_df\n\n\ndef save_data_for_chemprop(\n data: pd.DataFrame,\n project_name: str,\n train_idx: NDArray[np.int64],\n test_idx: NDArray[np.int64],\n seed: int,\n fold_number: int,\n) -> None:\n \"\"\"Save training and test data for Chemprop CLI.\n\n Args:\n data (pd.DataFrame): Dataframe containing data.\n project_name (str): Name of the project.\n train_idx (NDArray[np.int64]): List of indices for training data.\n test_idx (NDArray[np.int64]): List of indices for test data.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n \"\"\"\n Path(f\"data/{project_name}/seed_{seed}/fold_{fold_number}/\").mkdir(\n parents=True, exist_ok=True\n )\n data.iloc[train_idx].to_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/train.csv\", index=False\n )\n data.iloc[test_idx].to_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/test.csv\", index=False\n )\n\n\ndef train_chemprop_model(\n project_name: str, seed: int, fold_number: int, task_type: str\n) -> None:\n \"\"\"Train Chemprop model using CLI.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n task_type (str): Type of task ('classification' or 'regression').\n \"\"\"\n subprocess.run(\n [\n \"chemprop\",\n \"train\",\n \"--data-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/train.csv\",\n \"--task-type\",\n task_type,\n \"--output-dir\",\n f\"models/{project_name}/seed_{seed}/fold_{fold_number}\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n stderr=subprocess.DEVNULL,\n )\n\n\ndef make_chemprop_predictions(project_name: str, seed: int, fold_number: int) -> None:\n \"\"\"Make predictions using trained Chemprop model.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n \"\"\"\n subprocess.run(\n [\n \"chemprop\",\n \"predict\",\n \"--test-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/test.csv\",\n \"--model-paths\",\n f\"models/{project_name}/seed_{seed}/fold_{fold_number}/\",\n \"--preds-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/preds.csv\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n stderr=subprocess.DEVNULL,\n )\n\n\ndef calculate_chemprop_fingerprints(\n project_name: str, seed: int, fold_number: int, dataset_type: str\n) -> None:\n \"\"\"Calculate Chemprop fingerprints using trained Chemprop model.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n dataset_type (str): Type of dataset ('train' or 'test').\n \"\"\"\n subprocess.run(\n [\n \"chemprop\",\n \"fingerprint\",\n \"--test-path\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/{dataset_type}.csv\",\n \"--output\",\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/{dataset_type}_fps.csv\",\n \"--model-path\",\n f\"models/{project_name}/seed_{seed}/fold_{fold_number}/\",\n \"--ffn-block-index\",\n \"0\",\n ],\n check=True,\n stdout=subprocess.DEVNULL,\n stderr=subprocess.DEVNULL,\n )\n\n\ndef train_morgan_sklearn_model(\n fps_df: pd.DataFrame, y: pd.Series, train_idx: NDArray[np.int64], task_type: str\n) -> BaseEstimator:\n \"\"\"Train Scikit-learn model using Morgan fingerprints.\n\n Args:\n fps_df (pd.DataFrame): Dataframe of Morgan fingerprints.\n y (pd.Series): Series of target variable.\n train_idx (NDArray[np.int64]): List of indices for training data.\n task_type (str): Type of task ('classification' or 'regression').\n\n Returns:\n BaseEstimator: Trained Scikit-learn model.\n \"\"\"\n if task_type == \"regression\":\n rf = RandomForestRegressor(n_estimators=100, random_state=0)\n else:\n rf = RandomForestClassifier(n_estimators=100, random_state=0)\n X_train = fps_df.iloc[train_idx]\n y_train = y.iloc[train_idx]\n return rf.fit(X_train, y_train)\n\n\ndef train_chemprop_sklearn_model(\n project_name: str,\n seed: int,\n fold_number: int,\n y: pd.Series,\n train_idx: NDArray[np.int64],\n task_type: str,\n) -> BaseEstimator:\n \"\"\"Train Scikit-learn model using Chemprop fingerprints.\n\n Args:\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n y (pd.Series): Series of target variable.\n train_idx (NDArray[np.int64]): List of indices for training data.\n task_type (str): Type of task ('classification' or 'regression').\n\n Returns:\n BaseEstimator: Trained Scikit-learn model.\n \"\"\"\n X_train = pd.read_csv(\n Path(f\"data/{project_name}/seed_{seed}/fold_{fold_number}/train_fps_0.csv\")\n )\n if task_type == \"regression\":\n rf = RandomForestRegressor(n_estimators=100, random_state=0)\n else:\n rf = RandomForestClassifier(n_estimators=100, random_state=0)\n y_train = y.iloc[train_idx]\n return rf.fit(X_train, y_train)\n\n\ndef gather_predictions(\n morgan_rf: BaseEstimator,\n chemprop_rf: BaseEstimator,\n project_name: str,\n seed: int,\n fold_number: int,\n df: pd.DataFrame,\n fps_df: pd.DataFrame,\n test_idx: NDArray[np.int64],\n target_col: str,\n task_type: str,\n) -> dict:\n \"\"\"Gather predictions from all models.\n\n Args:\n morgan_rf (BaseEstimator): Trained Scikit-learn model on Morgan fingerprints.\n chemprop_rf (BaseEstimator): Trained Scikit-learn model on Chemprop fingerprints.\n project_name (str): Name of the project.\n seed (int): Random seed used for splitting data.\n fold_number (int): Fold number for cross-validation.\n df (pd.DataFrame): Dataframe containing data.\n fps_df (pd.DataFrame): Dataframe of Morgan fingerprints.\n test_idx (NDArray[np.int64]): List of indices for test data.\n target_col (str): Name of the column containing target variable.\n task_type (str): Type of task ('classification' or 'regression').\n\n Returns:\n dict: Updated dictionary with predictions.\n \"\"\"\n predictions = {}\n rf_preds = morgan_rf.predict(fps_df.iloc[test_idx])\n chemprop_rf_preds = chemprop_rf.predict(\n pd.read_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/test_fps_0.csv\"\n )\n )\n chemprop_probas = pd.read_csv(\n f\"data/{project_name}/seed_{seed}/fold_{fold_number}/preds.csv\"\n )[target_col]\n if task_type == \"regression\":\n chemprop_preds = chemprop_probas\n else:\n chemprop_preds = (chemprop_probas >= 0.5).astype(int)\n y_true = df.iloc[test_idx][target_col]\n predictions[f\"seed_{seed}_fold_{fold_number}\"] = {\n \"rf\": rf_preds,\n \"chemprop_rf\": chemprop_rf_preds,\n \"chemprop\": chemprop_preds,\n \"true\": y_true,\n }\n return predictions\nTo compare the models, I used 5x5 cross-validation. This approach involves using five different random seeds to generate data splits for 5-fold cross-validation, resulting in 25 unique train-test splits. For each split, the models are trained on the training data and then used to make predictions on the test data.
\ndef run_model_comparison(\n data_path: Path,\n project_name: str,\n smiles_col: str,\n target_col: str,\n task_type: str,\n n_seeds: int = 5,\n n_folds: int = 5,\n test_frac: float = 0.2,\n) -> dict:\n \"\"\"Run model comparison pipeline.\n\n Args:\n data_path (Path): Path to CSV file containing data.\n project_name (str): Name of the project.\n smiles_col (str): Name of the column containing SMILES strings.\n target_col (str): Name of the column containing target variable.\n task_type (str): Type of task ('classification' or 'regression').\n n_seeds (int, optional): Number of random seeds for cross-validation. Defaults to 5.\n n_folds (int, optional): Number of folds for cross-validation. Defaults to 5.\n test_frac (float, optional): Fraction of data to use as test set. Defaults to 0.2.\n\n Returns:\n dict: Dictionary containing predictions from all models.\n \"\"\"\n df, fps_df = prepare_data(data_path, smiles_col, target_col)\n\n all_predictions = {}\n for seed in tqdm(range(n_seeds), total=n_seeds, desc=\"Seeds\", unit=\"seed\"):\n cv = ShuffleSplit(n_splits=n_folds, test_size=test_frac, random_state=seed)\n for fold_number, (train_idx, test_idx) in tqdm(\n enumerate(cv.split(df)),\n total=n_folds,\n desc=f\"Seed {seed}: Folds\",\n unit=\"fold\",\n ):\n # Save data for training models using Chemprop CLI\n save_data_for_chemprop(\n df, project_name, train_idx, test_idx, seed, fold_number\n )\n\n # Train Chemprop model using CLI\n train_chemprop_model(project_name, seed, fold_number, task_type)\n\n # Make predictions on test set using trained Chemprop model\n make_chemprop_predictions(project_name, seed, fold_number)\n\n # Use trained model to create fingerprints\n calculate_chemprop_fingerprints(project_name, seed, fold_number, \"train\")\n calculate_chemprop_fingerprints(project_name, seed, fold_number, \"test\")\n\n # Train Scikit-learn models on Morgan fingerprints\n morgan_rf = train_morgan_sklearn_model(\n fps_df, df[target_col], train_idx, task_type\n )\n\n # Train Scikit-learn model on Chemprop fingerprints\n chemprop_rf = train_chemprop_sklearn_model(\n project_name,\n seed,\n fold_number,\n df[target_col],\n train_idx,\n task_type,\n )\n\n # Gather predictions\n predictions = gather_predictions(\n morgan_rf,\n chemprop_rf,\n project_name,\n seed,\n fold_number,\n df,\n fps_df,\n test_idx,\n target_col,\n task_type,\n )\n all_predictions.update(predictions)\n return all_predictions\nThe next two cells run the training and prediction process for the classification and regression datasets.
\n# Set up workspace\nPath(\"data\").mkdir(exist_ok=True)\nPath(\"models\").mkdir(exist_ok=True)\nnp.random.seed(42)\n\nn_seeds = 5\nn_folds = 5\ntest_frac = 0.2\ndata_path = Path(\"data/raw/bsep.csv\")\nproject_name = \"bsep\"\nsmiles_col = \"smiles\"\ntarget_col = \"bsep\"\ntask_type = \"classification\"\n\nclassification_predictions = run_model_comparison(\n data_path,\n project_name,\n smiles_col,\n target_col,\n task_type,\n n_seeds=n_seeds,\n n_folds=n_folds,\n test_frac=test_frac,\n)\nSeeds: 0%| | 0/5 [00:00n_seeds = 5\nn_folds = 5\ntest_frac = 0.2\ndata_path = Path(\"data/raw/biogen_logS.csv\")\nproject_name = \"logS\"\nsmiles_col = \"SMILES\"\ntarget_col = \"logS\"\ntask_type = \"regression\"\n\nregression_predictions = run_model_comparison(\n data_path,\n project_name,\n smiles_col,\n target_col,\n task_type,\n n_seeds=n_seeds,\n n_folds=n_folds,\n test_frac=test_frac,\n)\nSeeds: 0%| | 0/5 [00:00The next function takes the predicted values and true values from the test set and calculates a relevant performance metric.
\nFor the moderately imbalanced BSEP classification data (17% positive class), I used the Matthews Correlation Coefficient (MCC). MCC is a value between -1 and 1 and is a robust metric for imbalanced datasets. For the LogS regression data, I used the coefficient of determination (R2). An R2 value of 1 indicates perfect agreement between predicted and true values, with no lower limit for poor predictions.
\ndef calculate_performance(preds_dict: dict, metric: Callable) -> dict:\n \"\"\"Analyse performance of models.\n\n Args:\n preds_dict (dict): Dictionary containing predictions from all models.\n metric (Callable): Metric function to evaluate performance.\n\n Returns:\n dict: Dictionary containing performance metrics for each model.\n \"\"\"\n metrics = defaultdict(list)\n for predictions in preds_dict.values():\n metrics[f\"rf_{metric.__name__}\"].append(\n metric(predictions[\"true\"], predictions[\"rf\"])\n )\n metrics[f\"chemprop_rf_{metric.__name__}\"].append(\n metric(predictions[\"true\"], predictions[\"chemprop_rf\"])\n )\n metrics[f\"chemprop_{metric.__name__}\"].append(\n metric(predictions[\"true\"], predictions[\"chemprop\"])\n )\n return metrics\n\n\nclassification_metrics = calculate_performance(\n classification_predictions, matthews_corrcoef\n)\n\nregression_metrics = calculate_performance(regression_predictions, r2_score)\nTo visualise and inspect the classification results, I've used confusion matricies. All three models have a significant number of True Negatives (TN) in the top left square. The Random Forest model and the Chemprop model have a significant number of False Negatives (FN), whereas the Chemprop-RF model has more True Positives (TP) and more False Positives (FP).
\nconfusion_matricies(classification_predictions)\n
To visualise and inspect the regression results, I've used x-y scatter plots of the predicted values against the true values. The Random Forest regressor has poor predictions, with more values being overestimated. The Chemprop model and the Chemprop-RF models perform better, especially on these lower LogS values. Some performance metrics are annotated on the plots in the top left: R2, RMSE (root mean squared error), and MAE (mean absolute error).
\nscatter_plots(regression_predictions)\n
As recommended in the paper I've linked, I used one-way ANOVA to compare the performance distributions from the 5x5 cross-validation and determine if there was a statistically significant difference between the means. This was followed by Tukey's Honest Significant Differences as a post-hoc pairwise test. The results are plotted using the Simultaneous Confidence Interval plot from the statsmodels Python library.
\nFor both the classification and regression datasets, the ANOVA results showed a statistically significant difference in model performance (p<0.05).
\nFor the classification dataset, post-hoc tests revealed a statistically significant difference in performance between all pairs of models, with the Chemprop-RF model performing the best.
\nFor the regression dataset, the Random Forest model performed statistically significantly worse than the other two. However, there was no statistically significant difference between the performance of the Chemprop and Chemprop-RF models.
\nThe classification dataset (BSEP, 807 entries) is much smaller than the regression dataset (LogS, 2,173 entries). This might suggest that Chemprop-RF models are particularly effective for low-data problems, where there isn't enough data to train a high-performing Feed-Forward Neural Network (FFN) alone.
\n# Adapted from https://github.com/PatWalters/practical_cheminformatics_posts/blob/main/adme_comparison/\ndef run_anova(df_in: pd.DataFrame, col: str) -> float:\n \"\"\"Run one-way ANOVA on model performance.\"\"\"\n res_list = []\n for _, v in df_in.groupby(\"method\"):\n res_list.append(v[col].values)\n return f_oneway(*res_list)[1]\n\n\ndef plot_metric_comparison(metric_dict: dict, metric: Callable) -> None:\n \"\"\"Plot comparison of model performance using a specified metric.\n\n Args:\n metric_dict (dict): Dictionary containing performance metrics for each model.\n metric (Callable): Metric function used to evaluate performance.\n \"\"\"\n _, ax = plt.subplots()\n melt_df = pd.DataFrame(metric_dict).melt()\n melt_df[\"method\"] = melt_df.variable.map(\n {\n f\"rf_{metric.__name__}\": \"Random Forest\",\n f\"chemprop_rf_{metric.__name__}\": \"Chemprop-RF\",\n f\"chemprop_{metric.__name__}\": \"Chemprop\",\n }\n )\n best_model = (\n melt_df.groupby(\"method\")[\"value\"]\n .mean()\n .reset_index()\n .sort_values(\"value\", ascending=False)\n .method.values[0]\n )\n tukey = pairwise_tukeyhsd(\n endog=melt_df[\"value\"], groups=melt_df[\"method\"], alpha=0.05\n )\n tukey.plot_simultaneous(comparison_name=best_model, ax=ax, figsize=(8, 5))\n anova_p_value = run_anova(melt_df, \"value\")\n ax.set_title(\n f\"{metric.__name__.replace('_', ' ').title()} \" + f\"ANOVA p={anova_p_value:.3f}\"\n )\n ax.set_xlabel(f\"{metric.__name__.replace('_', ' ').title()}\")\n\n plt.tight_layout()\n plt.show()\n\n\nplot_metric_comparison(classification_metrics, matthews_corrcoef)\nplot_metric_comparison(regression_metrics, r2_score)\n