Comparing compound sets

A compound table is a tab-separated file with an inchikey column, a smiles column, or both, and whatever else you keep about your compounds:

#ligase representation: inchikey; salt-stripped, charge-neutralised
inchikey	smiles	name	activity_nm
BSYNRYMUTXBXSQ-UHFFFAOYSA-N	CC(=O)Oc1ccccc1C(=O)O	aspirin	120
HEFNNWSXXWATRW-JTQLQIEISA-N	C[C@@H](c1ccc(CC(C)C)cc1)C(=O)O	ibuprofen	85

The first line is the part that matters. It says what a name in this file is — an InChIKey — and what was done before the keys were computed. Every other column is carried without being read.

Loading one

import ligase as lg
 
library = lg.molecule_set("vendor/catalogue.tsv", release="ChEMBL 34")
 
library.dimension("representation").value
# 'InChIKey, salt-stripped and charge-neutralised'
library.dimension("representation").source
# 'from the representation line the standardiser wrote: vendor/catalogue.tsv'

Four dimensions apply: representation, organism, subject — what the set is of, which for a screen is usually the target — and release. Asking for anything else gets None, because a table of substances has no coordinate system, no contigs and no samples to have.

There is deliberately no identifier_namespace. An InChIKey is an identifier space, and representation already says which one this file is in; two axes carrying one fact is one too many.

Letting it refuse

Two sets prepared differently do not join. This is the whole point, and it happens before anything runs:

library = lg.molecule_set("vendor/catalogue.tsv")
hits = lg.molecule_set("screen/hits.tsv")
 
lg.check_step("overlap", {"library": library, "hits": hits})

Both states in full, the operations that differ named on their own line, and a question to act on. Where the two differ by more than one operation that line is the only part worth reading:

Where a file says nothing

A downloaded table of InChIKeys usually carries no representation line, because nothing but Ligase writes one. Such a table leaves the dimension unknown, and because representation is demanded wherever two compound tables are joined, two of them are refused rather than allowed to agree about nothing:

The engine does not assume as given for a file that says nothing. A ChEMBL download is standardised and a vendor list is not, so assuming the weaker state would be inventing the most dangerous of the possible answers.

You have two ways out, and they are not equivalent.

Declare it, where you know

library = lg.molecule_set(
    "chembl/actives.tsv",
    representation="inchikey; salt-stripped, charge-neutralised",
)

Accepted as given and not checked against the molecules, because the engine has no chemistry to check it with. This is the same bargain a declared identifier space gets: a declaration is you telling the engine something it cannot read.

Standardise it, where you do not

library = lg.molecule_set("vendor/catalogue.tsv")  # smiles column
clean = lg.chem.standardize(
    library, to="inchikey; salt-stripped, charge-neutralised"
)

lg.chem.standardize needs the chemistry extra:

pip install 'ligase-bio[chem]'

The step declares RDKit as its tool, so two things follow without your doing anything: the toolkit’s version joins the cache key, and so an upgrade reruns the standardisation it could have changed; and ligase methods credits RDKit in its reference list, because a dependency an analysis rested on belongs in the write-up.

Its output carries its own representation line, so everything downstream is checked against a statement the step can stand behind rather than one you typed. The step is a registered conversion across representation, which means it is the only way across the axis and the error above sends you to it by name.

The trap worth knowing before you choose a state

Salt stripping alone does not match a neutral set. RDKit’s FragmentParent leaves the deprotonated carboxylate, so:

and plain aspirin is the second. Two files that both say “salts removed” can still fail to overlap if only one of them was also neutralised.

So salt-stripped on its own is not a state to standardise to. ChEMBL, PubChem and ZINC all hand out neutral parents, so the reference set you actually have is almost certainly keyed -N, whether or not whoever exported it said so — while lg.chem.standardize(to="inchikey; salt-stripped") emits -M for every acid in your own set. The two declarations are honest, they compare equal, the join passes, and the overlap is wrong. A consumer measured the same generated set scoring 0% novel against one such reference and 50% against another, with ligase check reporting no problems for both.

If you are choosing a state for a project, inchikey; salt-stripped, charge-neutralised is the one that matches a neutral reference set. The one-term form stays in the vocabulary because files really are in that state — a salt-stripped export is a thing you can be handed — and not because it is a sensible thing to ask for.

What the keys can prove, and what they cannot

Charge state is the one part of a representation that is not taken on your word, because it is the one part a key records. The third block of an InChIKey is its protonation layer, so declaring charge-neutralised over keys that are not is a claim the file refutes, and it is refused at load:

Either the declaration names work the table never had, or it had it and the keys were never recomputed; the message says which state these keys are in, so the honest declaration is in front of you.

The other direction cannot be proved, only noticed. A salt-stripped table whose keys are all neutral may be exactly that — a set of compounds with nothing to neutralise keys -N with nothing applied — or it may have been neutralised by a step nobody wrote down. So it is a note rather than a refusal:

library = lg.molecule_set("vendor/catalogue.tsv",
                          representation="inchikey; salt-stripped")
library.notes

('`salt-stripped` was declared without `charge-neutralised`, and all 2314 of
this table\'s InChIKeys are protonation-neutral -- their third block is `N`,
which is what a charge-neutralised set looks like',)

ligase inspect prints the same line under noticed. If the table was neutralised, declaring it is what lets the set join the collections that declare it.

Salt stripping and tautomer normalisation leave no trace of this kind in a key, which is why they are recorded and not checked (D65).

And the one that costs stereochemistry

Asking for tautomer-normalised usually discards stereochemistry, because RDKit removes it at an sp3 centre inside a tautomeric system. The step will not let that happen quietly:

lg.chem.standardize(library, to="inchikey; tautomer-normalised")

ValueError: standardising discarded stereochemistry on 418 of 2314 molecules,
and the state you asked for does not say so. [...] Add `stereo-discarded` to
`to=`, so this set refuses to join one that still knows which enantiomer it
holds, or drop `tautomer-normalised` to keep the difference.

So you choose, with the number in front of you. Declare stereo-discarded and the set is honest about having lost the distinction — and will refuse to join a set that kept it. Drop tautomer-normalised and you keep the distinction and lose the tautomer matching. There is no option that gives you both, because RDKit does not offer one.

The reverse — declaring stereo-discarded when nothing was discarded — is only a note. It can cause a refusal you fix by declaring one thing consistently, which is the safe direction.

Ranking compounds

A virtual screen’s output is a ranking, and its order is the answer:

scored = lg.molecule_ranking("screen/docked.tsv", subject="EGFR")

Everything above applies unchanged. The one difference is the canonical form: a ranking’s row order is part of what it means, so reordering it changes its hash and reruns everything below it, where reordering a set changes nothing. For a screen that is the right way round — the order decides which compounds get bought.

Scoring a ranking against a library is the commonest join of the two kinds, and representation is demanded of it in either direction.

lg.check_step("score", {"library": library, "hits": scored})

What the engine will carry for you

A step over a compound set inherits its representation, so a ranking computed from a tautomer-normalised library carries that state, and anything that tries to score it against a set that kept its stereochemistry is refused the same way — several steps later, on a file that mentions no standardisation anywhere.

Standardising is a conversion, so it is also the one thing that changes the axis: its output declares the state it arrived in, and carries the organism, the subject and the release it was handed.

Molecules that do not standardise

A structure that does not parse, or a mixture with no single parent, is left out of the output and named — never dropped in silence:

the molecule `not a molecule at all` could not be standardised, so it is not
in the output
1 molecules could not be standardised, so the output is 1 molecules smaller
than the input

Twenty are named individually and the rest carried by the count, so the output being smaller than the input is a recorded fact rather than a discrepancy somebody notices six months later. The notes travel into the run report, provenance, a reproducibility bundle and the methods paragraph, the way every note in this engine does.

And molecules that become each other

Standardising a table that holds both a salt and its parent brings the two rows together: one substance, two rows, one key.

two or more structures standardise to `BSYNRYMUTXBXSQ-UHFFFAOYSA-N`, so this
table names one substance more than once

Both rows are kept. Dropping one would make the output quietly smaller, and choosing one would be a guess about whose activity data to carry — which is the argument against coercing identifiers, applied unchanged. So the collision is named and the decision is left with you: drop the row your analysis is not about, or read the table as a ranking, whose rows are not a set.

It is worth knowing why this one matters. A molecule set’s rows are a set, so a table naming one substance twice is not the kind it claims to be: anything counting its compounds double-counts, and a join on the key matches both rows. Nothing else would have told you — both rows are well formed, and the key is correct for each of them.