Clustering pockets by their interactions, and a gotcha in PLINDER 2024-06/v2 pli_unique_qcov

Lately I’ve been training generative models for proteins. While diving into this area, I went down a rabbit hole on how to cluster pockets by similarity. This matters both for building train/test splits without leakage and sampling training data in a balanced way.

While clustering based on protein sequence similarity or ligand similarity alone seems like a decent way to start with, my main concern was that this might not capture the actual pocket interactions that we are interested in. This means that there could be potentially two ligands that are somewhat different in similarity, yet they interact with the protein in the same way. The other way, of course could be true as well. Same ligand interacting with a several different proteins, yet forming similar interactions. ¹¹ Just a heads up: if we compare interactions only at residue pairs that a whole-chain MMseqs2/Foldseek alignment lines up, two unrelated proteins that hold a ligand the same way may never get compared, since the protein is so different. That’s not the focus of this post. Ideally, a model trained on these complexes should see the full range of interaction types, rather than leaving some interaction modes behind.

When I started pondering about this question, I thought to look at the PLINDER dataset [1]. They have done an amazing job in processing the PDB and annotating protein-ligand interactions. So my initial motivation was to use their own clustering for my training sampling.²² Because they release pli_unique_qcov__50__community. Let me lay out the key idea in clustering pockets, by leaving some detailed information later in the post.

How PLINDER compares two pockets (pli_unique_qcov)

First, for a protein complex A and B, we need to figure out which residue in A corresponds to residue in complex B. Plinder does this through MMseqs [2], and Foldseek [3]. Foldseek matters beause even if the sequences have diverged, they can share a fold ³³ I wonder if the counter example i stated in the above footnote would be considered. PLINDER computes each similarity score from both alignments and keeps the higher one. Put simply, the key idea in PLINDER is to figure out :

1. Align residues of protein A and B
2. keep only the aligned pairs where both residues are in their binding pocket (within 6 Å of the ligand, or interacting with it),
3. compare the interaction fingerprints of each such pair, add up the matches, and divide by the size of A’s fingerprint, and
4. connect A and B in a graph when that score passes a threshold, then take clusters of the graph.

After step 1, say we have an ARG176 and a ARG98, and assume that they are aligned. Now, say both of them have a fingerprint like ⁴⁴ This is a simplified example. For a real one from the PLINDER release, see the Arg224 example below.:

ARG176 = Counter({"hydrogen_bond": 2, "hydrophobic": 1})
ARG98  = Counter({"hydrogen_bond": 1, "hydrophobic": 1})

Now in a way, if we disregard the fact that counts are different, both these residues have the same two kinds of interactions, a hydrogen bond and a hydrophobic contact. So these two residues should count as a full match. Note that we have to do this for each pair. That is the idea behind PLINDER’s pli_unique_qcov. ⁵⁵ Note that its sibling score, pli_qcov is different, it counts every interaction, repeats included, and asks what fraction of A’s interactions B reproduces.

Let’s get back to pli_unique_qcov , which is used for several the downstream graph making (pli_unique_qcov__50__community) in PLINDER 2024-06/v2. To define pli_unique_qcov more precisely, let be the aligned pairs where both residues are in their pockets, and let be the set of interaction types at residue (the keys of its Counter). Then ⁶⁶ For simplicity, this formula assumes a single ligand and no artifact ligands. PLINDER counts unique fingerprints separately per ligand in the denominator, but merges them across non-artifact ligands at each residue in the numerator. Even after fixing the bug, a multi-ligand system can therefore score below 100% against itself.

where the bottom sum runs over A’s interacting residues. Note that the score is divided by A’s total, so and can differ.

So what’s the gotcha ?

As I mentioned before, the reason I turned to PLINDER was to reuse their clustering for my training. However, in the process of visualizing the clusters, and understanding the, pli_unique_qcov metric, codex noticed a bug in lines 796–801 of get_similarity_scores.py. . Even though pli_unique_qcov is supposed to look at the interaction type of a paired residue, it was comparing the interaction counts.

len(set(q_interactions[q_n].values()) & set(t_interactions[t_n].values()))

Because of this, going back to our previous example,

ARG176 = Counter({"hydrogen_bond": 2, "hydrophobic": 1})
ARG98  = Counter({"hydrogen_bond": 1, "hydrophobic": 1})

Suppose these are the only interacting residues in their pockets, so the whole score comes from this one pair. For ARG176 the code sees the counts {2, 1}, and for ARG98 it sees {1}. They share one count, so the pair scores 0.5 instead of 1.0. Even ARG98 compared with itself scores only 0.5. A more prominent mistake would be: a residue with one hydrophobic contact and a residue with one salt bridge both have the count 1, so they score as a full match despite sharing no interaction at all. Even though this seems like a big issue, I think partly its watered down since we align the two proteins first. However, for my train set I did see clusters drastically change. So the aligned residues anyway are somewhat similar interaction forming residues. But it does compare the incorrect thing. For instance:

A PDB example: 7FK1 (W3H), 7FLE (VB9), 7FL1 (VO9)

The three PLINDER systems below are crystal structures of Aar2/RNaseH from a fragment screen [4], and in each one PLINDER annotates one interacting residue, Arg224:

• 7fk1__1__1.A__1.C
• 7fle__1__1.A__1.C
• 7fl1__1__1.A__1.D

These are PLINDER’s annotations for Arg224 in each:

arg224_7FK1 = Counter({"type:salt_bridges__protispos:True": 1})
arg224_7FLE = Counter({"type:hydrogen_bonds__protisdon:True__sidechain:True": 1})
arg224_7FL1 = Counter({"type:hydrogen_bonds__protisdon:True__sidechain:True": 2})

Because each pocket has only this one interacting residue, the score between two of these systems is just this residue pair, this makes it easier to show the result of the equation:

So I fixed this locally, which took some groundwork. The released data don’t include residue-level alignments, so this meant re-running MMseqs2 and Foldseek on the system pairs PLINDER had already matched. Afterwards the clusters looked bit better. In a later post I’d like to share that!

Summary

I think PLINDER authors might be actively working on new releases, and better clustering approaches. This article hopes to just document a potential bug in PLINDER I got stuck in and document it for anyone who else might be going to train with clustered pli systems. On the other hand, this shows the human-AI collaboration aspect 😅, if I hadn’t tried to understand how clustering was done with codex, or if I tried to do this alone, or had used pli_unique_qcov__50__community directly, I wouldn’t have noticed this myself.

I’m really grateful for the PLINDER paper for their work, diving into these metrics made me realize how complicated it is to define metrics, and cluster things. Even if you have the right metric, you have think through lots of different views on what it means to be included in a cluster. I’ve been postponing of writing this for around 2 months, I’m glad I documented this finally and I’d like to learn from you as well!