Alignment by a subset of atoms produces incorrect coordinates

[franciscoadasme]

It's very common to align two structures by a subset of atoms, e.g., aligning two proteins by the backbone. While the alignment between two atom sets works, computing the transformation of a subset of atoms and applying it to a superset produces incorrect coordinates.

For instance, the following code aligns the alpha carbons of chain A (sorted by their position along Z):

structure = Chem::Structure.from_pdb("actual.pdb")
ref_structure = Chem::Structure.from_pdb("ref.pdb")

atoms, ref_atoms = {structure, ref_structure}.map &.atoms[[52, 2, 12, 22, 32, 42]]
transformation = AffineTransform.aligning(atoms, to: ref_atoms)
atoms.coords.transform! transformation

Chem::Spatial.rmsd(atoms.coords, ref_atoms.coords) # => 0.09093355386835195

The resulting RMSD equals the minimized RMSD (Chem::Spatial.rmsd(atoms.coords, ref_atoms.coords, minimize: true)), which is correct. However, if the same transformation is applied to the entire structure instead produces different coordinates:

transformation = AffineTransform.aligning(atoms, to: ref_atoms)
structure.coords.transform! transformation
Chem::Spatial.rmsd(atoms.coords, ref_atoms.coords) # => 1.5788702192923796

There seems to be an issue when applying a transformation to an atom set due to intermediate translations.

actual.pdb.txt ref.pdb.txt

[franciscoadasme]

Another example of applying a transformation to a set of atoms. The idea is to place the protein chain (using the four carbonyl oxygens as reference) along the X axis.

structure = Chem::Structure.from_pdb "7M2I_sf-helice--aligned.pdb"

ks = structure.atoms.select(&.potassium?)
os = (75..78).map { |i| structure.chains.first.dig(i, "O") }

s.coords.translate! -ks.mean(&.coords)

# compute center of selected oxygens and set Y component to 0
oc = os.mean(&.coords).reject(Y).normalize
transform = Chem::Spatial::AffineTransform.aligning oc, to: Chem::Spatial::Vec3[1, 0, 0]
s.coords.transform! transform # produces incorrect coordinates

Note that applying the transformation to the original vector oc does work.

oc                      # => Vec3[-0.6698114, 0, 0.7425313]
oc.transform(transform) # => Vec3[ 1  0 -2.220446e-16 ]

7M2I_sf-helice--aligned.pdb.txt 7M2I_sf-helice--aligned-bad.pdb.txt

[franciscoadasme]

It seems the issue is related to the CoordinatesProxy#transform method because atoms.each { |atom| atom.coords = atom.coords.transform(transform) } works as expected.