Skip to content

ResidueGroup silently drops residues when one ResidueId maps to several disconnected fragments #366

Description

@alex-sbaq

What I expected, and what happens instead

split_mol_by_residues emits one Residue per disconnected fragment. When a single PDB
residue arrives as two or more fragments, they all carry the same ResidueId, and
ResidueGroup.__init__ builds its dict from a list of pairs:

# prolif/residue.py:223
super().__init__([(r.resid, r) for r in self._residues])

Duplicate keys mean every fragment but the last is discarded. I would expect either a merge,
distinct keys, or a warning naming the affected residues. What happens is that atoms vanish
from the fingerprint with nothing logged.

The class also ends up internally inconsistent, which may be the easier half to fix.
self._residues and the parallel name / number / chain arrays keep every fragment while
self.data keeps one, so len(rg.name) != rg.n_residues, and rg[i] can reach a fragment
that rg[resid] cannot.

Why it is worth reporting rather than a curiosity: on a trimmed multi-chain receptor loaded
through MDAnalysis, connectivity gaps produced fragments including a residue consisting of a
single hydrogen atom. ProLIF returned almost no interactions for a pose that plainly had many,
and raised nothing. Nothing in the output separates "this pose makes few contacts" from "most
of the binding site was dropped at load time," so the wrong answer is entirely plausible-looking.

Full error message

There is none, and that is the substance of the report. The collapse is silent. I grepped
prolif/**/*.py on master for duplicat|already exists|overwrit|collision|clobber and the
only hits are unrelated water dedup in water_bridge.py.

The one place an error does surface is select(), indirectly. It does self._residues[mask],
which is fragment-length, while its docstring says the mask should have "the same length as the
number of residues in the ResidueGroup", and both len(rg) and n_residues report the
collapsed dict length. So on an affected molecule a mask built the documented way raises:

IndexError: boolean index did not match indexed array along axis 0; size of axis is 3 but
size of corresponding boolean axis is 2

Code snippet

import numpy as np
import prolif as plf
from rdkit import Chem

# A/ALA1 written as two atoms far enough apart that no bond is inferred, so it
# becomes two disconnected fragments carrying the same ResidueId.
pdb = """\
ATOM      1  N   ALA A   1       0.000   0.000   0.000  1.00  0.00           N
ATOM      2  CA  ALA A   1      10.000  10.000  10.000  1.00  0.00           C
ATOM      3  N   GLY A   2      20.000  20.000  20.000  1.00  0.00           N
END
"""
mol = Chem.MolFromPDBBlock(pdb, removeHs=False, sanitize=False)
rg = plf.Molecule(mol).residues

print(len(rg._residues), [str(r.resid) for r in rg._residues])
print(len(rg), [str(k) for k in rg.data])
print(rg["ALA1.A"].GetNumAtoms())
print(len(rg.name), rg.n_residues)
rg.select(np.ones(len(rg), dtype=bool))

Output on 2.2.1:

3 ['ALA1.A', 'ALA1.A', 'GLY2.A']
2 ['ALA1.A', 'GLY2.A']
1
3 2
IndexError: boolean index did not match indexed array along axis 0; size of axis is 3 but
size of corresponding boolean axis is 2

The N atom of A/ALA1 is gone. Last fragment wins, so the CA is what survives.

Workarounds tried

Routing the receptor through RDKit's PDB reader rather than
prolif.Molecule.from_mda(mda.Universe(...)) avoided it on my system, because the RDKit
reader inferred the connectivity that MDAnalysis did not. That is a workaround for one source
of fragmentation, not for the collapse itself. guess_bonds() is the equivalent fix on the
MDAnalysis side and is already in the PDB tutorial's troubleshooting note.

The only reliable detection I found is asserting len(rg.name) == rg.n_residues after load,
which is not something a user would think to do.

Notes on related issues

I searched open and closed issues, discussions, and the changelog before filing, and did not
find this mechanism reported. Three adjacent items, none of which I think is a duplicate:

  • Incorrect detection of interactions #358 looks like the same root cause with a different symptom. The note there about the
    proper fix being substructure matching on the unfragmented molecule, and needing a
    substantial refactor, would presumably cover this too. Happy for this to be folded in there
    if you see it the same way.
  • ProLif not showing hydrophobic interaction (and some others) on Ligand-Protein #242 may be this bug misdiagnosed. The reported symptom is missing hydrophobic
    interactions after a spatial selection without byres, and the explanation given is
    incomplete valence causing atoms to be perceived as charged. Fragmented residues would also
    have been silently collapsed in that case, so both effects could have been in play.
  • The v2.1.0 use_segid changelog entry is close but does not help here. It disambiguates
    genuinely distinct residues that share resname, resnumber and chain. Two fragments of one
    residue share the segid as well, so they still collide.

I also could not find a test covering the duplicate case. tests/test_residues.py:220 asserts
rg.n_residues == len(rg), which holds only when no collapse has occurred.

Environment

  • ProLIF 2.2.1 (also reproduces on 2.1.0, identical output)
  • RDKit 2026.03.5
  • Python 3.12, Linux
  • prolif/residue.py:223 is byte-identical on master and develop as of today

Happy to help with a patch or a test if that is useful, though given #358 you may prefer to
fold this into the refactor. Either way, tell me what would help.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions