SeqALigSeqB v2: predicting which of two proteins a compound prefers

1. What goes in

One row is one question: given this compound, does protein A or protein B bind it more tightly? The two proteins occupy fixed positions, and the position is the question. Sequence A, then the ligand, then sequence B.

indices 0:480 · 480:1518 · 1518:1998  ·  total width 1,998

No protein structure, no docked pose, no model of the binding site. Sequence and compound only.

2. What comes out

Two numbers that sum to 1, and they are never collapsed into one.

A fully committed prediction reads 0, 1 or 1, 0. Anything between is the model reporting how sure it is, and that spread is not decoration: section 4 shows accuracy climbing from 58% to 99% across it on the test set. The winner is whichever of the two values is larger; both are kept.

3. A worked example

One comparison drawn from the training set, given in full so it can be redrawn and re-derived independently. Both proteins are human wild-type kinases whose sequences match the canonical UniProt entries exactly, verified by md5, and both have deposited structures. The compound is an approved drug.

Cc1nc(Nc2ncc(s2)C(=O)Nc3c(C)cccc3Cl)cc(n1)N4CCN(CCO)CC4

Presented the other way round, as GSK3B, dasatinib, ABL1, the two outputs exchange places and the model names GSK3B the loser. The answer follows the order the question is asked in, which is the whole design.

Sequence A, ABL1, P00519, complete

MLEICLKLVGCKSKKGLSSSSSCYLEEALQRPVASDFEPQGLSEAARWNSKENLLAGPSENDPNLFVALYDFVASGDNTLSITKGEKLRVLGYNHNGEWCEAQTKNGQGWVPSNYITPVNSLEKHSWYHGPVSRNAAEYLLSSGINGSFLVRESESSPGQRSISLRYEGRVYHYRINTASDGKLYVSSESRFNTLAELVHHHSTVADGLITTLHYPAPKRNKPTVYGVSPNYDKWEMERTDITMKHKLGGGQYGEVYEGVWKKYSLTVAVKTLKEDTMEVEEFLKEAAVMKEIKHPNLVQLLGVCTREPPFYIITEFMTYGNLLDYLRECNRQEVNAVVLLYMATQISSAMEYLEKKNFIHRDLAARNCLVGENHLVKVADFGLSRLMTGDTYTAHAGAKFPIKWTAPESLAYNKFSIKSDVWAFGVLLWEIATYGMSPYPGIDLSQVYELLEKDYRMERPEGCPEKVYELMRACWQWNPSDRPSFAEIHQAFETMFQESSISDEVEKELGKQGVRGAVSTLLQAPELPTKTRTSRRAAEHRDTTDVPEMPHSKGQGESDPLDHEPAVSPLLPRKERGPPEGGLNEDERLLPKDKKTNLFSALIKKKKKTAPTPPKRSSSFREMDGQPERRGAGEEEGRDISNGALAFTPLDTADPAKSPKPSNGAGVPNGALRESGGSGFRSPHLWKKSSTLTSSRLATGEEEGGGSSSKRFLRSCSASCVPHGAKDTEWRSVTLPRDLQSTGRQFDSSTFGGHKSEKPALPRKRAGENRSDQVTRGTVTPPPRLVKKNEEAADEVFKDIMESSPGSSPPNLTPKPLRRQVTVAPASGLPHKEEAGKGSALGTPAAAEPVTPTSKAGSGAPGGTSKGPAEESRVRRHKHSSESPGRDKGKLSRLKPAPPPPPAASAGKAGGKPSQSPSQEAAGEAVLGAKTKATSLVDAVNSDAAKPSQPGEGLKKPVLPATPKPQSAKPSGTPISPAPVPSTLPSASSALAGDQPSSTAFIPLISTRVSLRKTRQPPERIASGAITKGVVLDSTEALCLAISRNSEQMASHSAVLEAGKNLYTFCVSYVDSIQQMRNKFAFREAINKLENNLRELQICPATAGSGPAATQDFSKLLSSVKEISDIVQR

Sequence B, GSK3B, P49841, complete

MSGRPRTTSFAESCKPVQQPSAFGSMKVSRDKDGSKVTTVVATPGQGPDRPQEVSYTDTKVIGNGSFGVVYQAKLCDSGELVAIKKVLQDKRFKNRELQIMRKLDHCNIVRLRYFFYSSGEKKDEVYLNLVLDYVPETVYRVARHYSRAKQTLPVIYVKLYMYQLFRSLAYIHSFGICHRDIKPQNLLLDPDTAVLKLCDFGSAKQLVRGEPNVSYICSRYYRAPELIFGATDYTSSIDVWSAGCVLAELLLGQPIFPGDSGVDQLVEIIKVLGTPTREQIREMNPNYTEFKFPQIKAHPWTKVFRPRTPPEAIALCSRLLEYTPTARLTPLEACAHSFFDELRDPNVKLPNGRDTPALFNFTTQELSSNPPLATILIPPHARIQAAASTPTNATAASDANTGDRGQTNNAASASASNST

4. The data

A comparison exists wherever the same compound has been measured against two different protein sequences and the two readings can be ordered. Greater-than and less-than readings are used, not discarded: a compound that is weaker than 10 µM against one protein and 1 nM against another has a decidable winner. Pairs whose readings overlap, so that no winner can be named, are not used.

Every comparison is presented twice, once as A, ligand, B with its answer and once as B, ligand, A with the answer reversed. This is what forces the model to read the order: a model that ignored which slot a protein occupied would give the same answer to both versions and be wrong on exactly one of them every time, scoring 50% regardless of how well it understood the chemistry. Protein A is the winner in 49.0% of comparisons before the swap and exactly 50.0% after, so no model can profit from guessing a slot.

The exclusion operates on the whole triplet, protein-ligand-protein, because that is the unit the model is asked about. A ChEMBL comparison is removed only when that exact three-way combination is already in training. Nothing is removed for sharing a protein, and nothing for sharing a protein-and-compound pair: if ChEMBL pairs a familiar compound against a different second protein, that is a question the model has never been asked and it is kept. Potency is computed from the reported value in nanomolar as 9 - log10(nM) rather than from ChEMBL's pchembl_value, which is left empty on every greater-than and less-than row and would have silently discarded 311,000 qualified readings. Ki and Kd measurements are used, and both readings in a comparison must be of the same measurement type, so a Ki is never weighed against an IC50.

5. Train and test performance

The model reproduces its training data at 94.3% and generalises to unseen ChEMBL comparisons at 75.3%. Both are measured on the same two quantities, accuracy against the known winner, over the full sets rather than a sample.

Both curves run the same direction. The further apart the two proteins actually are in potency, the more reliably the model calls it, reaching 90.2% on the test set where the separation exceeds two logs. Where the two proteins are within half a log of each other, the test accuracy of 57.7% is close to a coin flip.

This is the practical result. The two outputs are not merely a label in disguise; the distance between them tracks how often the model is right. On unseen ChEMBL data, predictions where the two values are far apart are correct 99.3% of the time, and predictions near a coin flip are correct 57.9% of the time. The model flags its own uncertain calls. On the test set, 505,499 comparisons carry a spread above 0.6 and are correct 97.2% of the time taken together.

6. Order reversal

Exchanging the two proteins should reverse the answer. Measured on the output values rather than on the winner label, the mean of |p(1 | A,B) - p(0 | B,A)| would be zero for a perfectly antisymmetric model.

The swap augmentation did what it was there to do. Presenting the same comparison the other way round moves the two outputs almost exactly to their mirror positions and changes the declared winner in 95.1% of unseen cases.

7. Prediction strength and accuracy

Prediction strength is the larger of the two output probabilities, so it runs from 0.5, a coin flip, to 1.0, a fully committed call. This is the same definition used in the LigASeqLigB report, so the two models can be read side by side. Accuracy rises with strength on both sets, and it never falls at any step.

Acting only on predictions at or above a threshold. The point of this second table is the trade: accuracy rises as the threshold rises, but the number of comparisons you can still answer falls away. The right cutoff is where the accuracy is worth the predictions you give up.

On the test set the accuracy gained per prediction surrendered falls off after 0.70. Going from 0.50 to 0.70 buys 17 accuracy points and costs 63% of the answers. Going from 0.70 to 0.90 buys only 7 more points and costs 88% of what remained. On that reading 0.70 is the useful operating point: more than a third of all comparisons answered at 92.3%, against 75.3% if you answer everything.

Where the shortlist must not contain mistakes, 0.80 answers 505,499 comparisons at 97.2%. Beyond 0.90 the gain is fractions of a point for most of the remaining data, so there is little reason to go further.

Note that the strength value does not mean the same thing on the two sets. On training comparisons a 0.70 prediction is right 99% of the time; on ChEMBL the same 0.70 is right 92% of the time. The ordering is reliable, the absolute number is not transferable, which is why the cutoff must be set from the test column and not the training one.

8. Methods

Building a comparison

Each measurement is converted to a range on pIC50, where higher is more potent. An exact value gives a point. A less-than reading against a concentration establishes that the compound is more potent than that value and gives a lower bound. A greater-than reading gives an upper bound. A single-concentration screen is treated the same way, potent if it registered as active and weak if it did not. Repeat measurements of the same compound against the same sequence are combined into one range, with exact values taking precedence over bounds and the intersection taken otherwise.

Two ranges form a usable comparison when they cannot trade places, that is when the lower bound of one exceeds the upper bound of the other. Everything else is set aside, because no winner can be named. The recorded separation is the guaranteed minimum, the distance between the two bounds rather than between point estimates.

Featurisation

Protein sequences are keyed by the md5 of the sequence, so a mutant is a target in its own right rather than being folded into its wild type. Vectors were verified before use to be the plain full-sequence mean pool rather than a pocket-pooled or domain-pooled variant, reproducing a freshly computed mean pool to a maximum absolute difference of 0.0 across all 733 available. Nine sequences have no embedding and were excluded along with the 29 measurements that reference them; this is the only inclusion criterion applied.

Order augmentation

Every comparison is written into the design matrix twice, as A, ligand, B with its answer and as B, ligand, A with the answer inverted. The two copies are adjacent halves of one matrix and are fitted together in a single call.

Split policy

Nothing is withheld from the KKB. The test set is external, drawn from ChEMBL, and is excluded on the triplet: a ChEMBL comparison is discarded only when that same protein-ligand-protein combination is already present in training. Sequence order is ignored when matching, because a training comparison is presented in both directions, so encountering it either way counts as having seen it.

Fitting

Estimator

sklearn.ensemble.RandomForestClassifier

n_estimators

200

min_samples_leaf

8

max_features

sqrt

random_state

1

bootstrap

True (sklearn default)

class_weight

None; the swap makes the classes exactly balanced

Design matrix

8,680,234 × 1,998 float32, 64.6 GB, built in memory

Fit time

1,502 s on 28 cores

Peak memory

about 85 GB

Cap applied

none

Environment

Python 3.10.15, scikit-learn 1.7.2, numpy 2.2.6, rdkit 2025.09.5, macOS arm64

Reproducing or extending the fit

The pipeline is deterministic. With the same inputs and random_state=1 it produces the same forest, so a rerun is a check rather than a new result. The stages run in order: 300 builds the KKB comparisons, 301 builds the ChEMBL test set, 302 marks everything for training, 303 fits and scores the test set, 304 scores the training set, 305 measures output strength, 306 completes the export.

Adding trees to this forest without refitting it. The trees are independent, so an existing forest can be enlarged. Set warm_start=True on the loaded model, raise n_estimators, and call fit with the identical design matrix and labels. Only the new trees are grown; the existing 200 are left untouched. The matrix must be identical, since sklearn does not check and a changed matrix would give a forest whose trees disagree about what they were trained on. Rebuilding that matrix takes about 70 seconds, and each additional 100 trees costs roughly 750 seconds at this size.

Adding new measurements. A random forest has no incremental update for new rows; partial_fit does not exist for this estimator. New KKB or ChEMBL data therefore requires a complete refit, which is the 1,502 second path above. Rerun 300 so the comparison set is rebuilt from the enlarged source, then 302 and 303. Because the comparison count grows roughly with the square of the number of sequences a compound has been measured against, adding a broadly profiled compound adds far more comparisons than adding a singleton, and the design matrix should be re-sized before the run rather than discovered mid-fit.

9. The exported model

The release is named SeqALigSeqB_v2_selectivity and sits in export/, which is excluded from Dropbox so it does not sync. It is self-contained and packaged so a third party can test it and stand it up on another machine.

Standing it up elsewhere

Copy the directory, install the pinned versions, then run python selftest.py and require it to print ALL CHECKS PASSED before trusting any number from the model on that machine. The pinning is not a formality: a random forest is a pickled object graph, and a different scikit-learn will either refuse to load it or load it and score differently. The pins are scikit-learn 1.7.2, numpy 2.2.6, joblib 1.5.3 and rdkit 2025.09.5 on Python 3.10.15.

The self-test replays 50 reference comparisons chosen to span the confidence range, from near coin flips to fully committed calls, and requires each to match to within 1e-6. It then exchanges the two sequences on every case and checks that the outputs mirror. Run here against the finished model it reports a maximum absolute difference of 3.3e-16, which is floating-point noise, and passes.

Loading takes about a minute and holds roughly 12 GB of RAM resident, which is the practical constraint on where this can be deployed.

sequence_vectors.npz covers only the 728 sequences the model was trained on. A sequence outside that set raises an error rather than being scored against a guessed vector.

Every KKB comparison was used for training, including all 5,846 in which a wild type is compared against one of its own mutants. Because nothing was withheld and ChEMBL is keyed by gene and therefore contains no mutants, this release has no measured wild-type-versus-mutant accuracy and none can be quoted for it. The model is trained on those comparisons and will predict on them; that capability is unmeasured rather than absent. This is recorded in MANIFEST.json and in split_counts.json so it travels with the files.