Merge pull request #42 from xchem/assembly_based_alignment

Assembly based alignment

data/assemblies.yaml

@@ -0,0 +1,36 @@
+assemblies:
+    dimer:
+        reference: Mpro-IBM0045
+        biomol: A,B
+        chains: A,A(-x,y,-z)
+    monomer:
+        reference: 7ql8
+        biomol: B
+        chains: A
+    fake_dimer:
+        reference: 8e1y
+        biomol: A,C
+        chains: A,B
+    fake_tetramer:
+        reference: 8dz9
+        biomol: D,C,E,F
+        chains: A,B,C,D
+crystalforms:
+    xtalform1:
+        reference: Mpro-IBM0045
+        assemblies:
+            1:
+                assembly: dimer
+                chains: A,A(-x,y,-z)
+    xtalform2:
+        reference: 7ql8
+        assemblies:
+            1:
+                assembly: monomer
+                chains: A
+    xtalform3:
+        reference: 8e1y
+        assemblies:
+            1:
+                assembly: dimer
+                chains: A,B

pyproject.toml

@@ -79,16 +79,19 @@ exclude = [".tox", "venv"]
 [tool.black]
 line-length = 119
 
+#[tool.pytest.ini_options]
+## Run pytest with all our checkers, and don't spam us with massive tracebacks on error
+#addopts = """
+#    --tb=native -vv --doctest-modules --doctest-glob="*.rst"
+#    --cov=ligand_neighbourhood_alignment --cov-report term --cov-report xml:cov.xml
+#    """
+## https://iscinumpy.gitlab.io/post/bound-version-constraints/#watch-for-warnings
+#filterwarnings = "error"
+## Doctest python code in docs, python code in src docstrings, test functions in tests
+#testpaths = "docs src tests"
+
 [tool.pytest.ini_options]
-# Run pytest with all our checkers, and don't spam us with massive tracebacks on error
-addopts = """
-    --tb=native -vv --doctest-modules --doctest-glob="*.rst"
-    --cov=ligand_neighbourhood_alignment --cov-report term --cov-report xml:cov.xml
-    """
-# https://iscinumpy.gitlab.io/post/bound-version-constraints/#watch-for-warnings
-filterwarnings = "error"
-# Doctest python code in docs, python code in src docstrings, test functions in tests
-testpaths = "docs src tests"
+testpaths = ['tests']
 
 [tool.coverage.run]
 data_file = "/tmp/ligand_neighbourhood_alignment.coverage"

scripts/generate_p1_version.py

@@ -0,0 +1,46 @@
+import gemmi
+from ligand_neighbourhood_alignment import alignment_heirarchy, dt
+
+
+def main():
+    dtag = 'Mpro-x0107'
+    file = 'data/refine_8.split.bound-state.pdb'
+    st = gemmi.read_structure(file)
+    generators = [
+        dt.Generator('A', 'A', 'x,y,z'),
+        dt.Generator('B', 'A', '-x,y,-z'),
+    ]
+
+    # Setup new structure to add biochains to
+    new_st = gemmi.Structure()
+    new_st.spacegroup_hm = gemmi.SpaceGroup('P 1').xhm()
+    new_st.cell = st.cell
+    new_model = gemmi.Model("0")
+    new_st.add_model(new_model)
+
+    # Iterate over chain, biochain, transform tuples in the assembly
+    for generator in generators:
+        # Generate the symmetry operation
+        op = gemmi.Op(generator.triplet)
+
+        # Create a clone of base chain to transform
+        new_chain = st[0][generator.chain].clone()
+
+        # Transform the residues
+        for residue in new_chain:
+            for atom in residue:
+                atom_frac = st.cell.fractionalize(atom.pos)
+                new_pos_frac = op.apply_to_xyz([atom_frac.x, atom_frac.y, atom_frac.z])
+                new_pos_orth = st.cell.orthogonalize(gemmi.Fractional(*new_pos_frac))
+                atom.pos = gemmi.Position(*new_pos_orth)
+        new_chain.name = generator.biomol
+        new_st[0].add_chain(new_chain)
+
+    new_st.setup_entities()
+    print(new_st.cell)
+    print(new_st.spacegroup_hm)
+    new_st.write_pdb('Mpro-x0107_fake_P1.pdb')
+
+
+if __name__ == "__main__":
+    main()

src/ligand_neighbourhood_alignment/align_xmaps.py

@@ -7,6 +7,7 @@
 import numpy as np
 from loguru import logger
 
+from ligand_neighbourhood_alignment import alignment_heirarchy
 from ligand_neighbourhood_alignment.data import (
     Block,
     CanonicalSites,
@@ -529,6 +530,8 @@ def __align_xmap(
     canonical_site_id,
     output_path: Path,
     aligned_res,
+chain_to_assembly_transform,
+assembly_transform
 ):
     # Get the ligand neighbourhood
     # neighbourhood: LigandNeighbourhood = neighbourhoods.get_neighbourhood(lid)
@@ -557,13 +560,20 @@ def __align_xmap(
 
     # Get the subsite transform
     # print(conformer_site_transforms)
-    conformer_site_transform = transform_to_gemmi(conformer_site_transforms[(canonical_site_id, conformer_site_id)])
-    running_transform = conformer_site_transform.combine(running_transform)
+    # conformer_site_transform = transform_to_gemmi(conformer_site_transforms[(canonical_site_id, conformer_site_id)])
+    # running_transform = conformer_site_transform.combine(running_transform)
 
     # Get the site transform
     # canonical_site_transform = transform_to_gemmi(canonical_site_transforms[canonical_site_id])
     # running_transform = canonical_site_transform.combine()
 
+    # Update the transform with the assembly alignment
+    # # Get the xtalform to assembly transform
+    running_transform = alignment_heirarchy._transform_to_gemmi(chain_to_assembly_transform).combine(running_transform)
+
+    # # Get the assembly alignment transform
+    running_transform = alignment_heirarchy._transform_to_gemmi(assembly_transform).combine(running_transform)
+
     logger.debug(
         f"Transform from native frame to subsite frame to site frame is: {gemmi_to_transform(running_transform)}"
     )