`_compute_common_indices` fails when homologous structures have different residue offsets

[minhuanli]

Currently, _compute_common_indices relies on exact string matching between moving_cra and reference_cra to find overlapping atoms, which is to say it only works when the two pdb objects has the exact same sequence.

This can be an issue in practice when there are missing residues or just truncation residues in the N or C terminal.

For example, in this simple 4lZT example (I attached the two pdb files here), where the second one is just a truncation of 4 residues on the N terminus, we would want the methods to help us find all the remaining residues as overlapping

4LZT_truncated.pdb
4LZT.pdb

# Sequence 1
KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNTDGSTDYGILQINSRWWCNDGRTPGSRNLCNIPCSALLSSDITASVNCAKKIVSDGNGMNAWVAWRNRCKGTDVQAWIRGCRL
# Sequence 2
RCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNTDGSTDYGILQINSRWWCNDGRTPGSRNLCNIPCSALLSSDITASVNCAKKIVSDGNGMNAWVAWRNRCKGTDVQAWIRGCRL

But as the function requires exact string matching, it will only find something occasionally because of the repetition.

from SFC_Torch import PDBParser
pdb1 = PDBParser("4LZT.pdb")
pdb2 = PDBParser("4LZT_truncated.pdb")
ind1, ind2 = compute_common_indices(pdb1.cra_name, pdb2.cra_name, "CA")

>>> print(ind1), print(ind2)
[172 262 336 373 380 539 605 640]
[171 263 340 367 374 535 601 634]

Solution

To make it more robust, we should use what we had in the ROCKET codebase, using the Smith Waterman sequence alignment:

https://github.com/rs-station/ROCKET/blob/de3041a08265a7e4933960aa969974e3453f146c/rocket/refinement_utils.py#L72C1-L73C1