Add option to query CCD for non-standard residues

Author: benmwebb

ihm/util/make_mmcif.py

@@ -27,10 +27,13 @@
 import ihm.model
 import ihm.protocol
 import ihm.util
+import ihm.format
+import urllib.request
 import os
 import argparse
 import collections
 import operator
+import warnings
 
 
 # All canonical atom names for each standard residue type, as per CCD.
@@ -148,7 +151,7 @@ def add_ihm_info(s, fix_histidines, check_atom_names):
                     if fix_histidines:
                         _fix_histidine_het_atoms(model, histidines)
                     if check_atom_names != 'no':
-                        _check_atom_names(model)
+                        _check_atom_names(model, check_atom_names == 'all')
     if fix_histidines:
         _fix_histidine_chem_comps(s, histidines)
     return s
@@ -165,26 +168,53 @@ def _fix_histidine_het_atoms(model, histidines):
             atom.het = False
 
 
-def _get_non_canon(seen_atom_names):
+class _ChemCompAtomHandler:
+    not_in_file = omitted = unknown = None
+
+    def __init__(self):
+        super().__init__()
+        self.atoms = collections.defaultdict(set)
+
+    def __call__(self, comp_id, atom_id):
+        self.atoms[comp_id].add(atom_id)
+
+
+def _get_non_std_restyp(restyp):
+    """Return CCD info for the given residue type"""
+    url_pattern = 'http://ligand-expo.rcsb.org/reports/%s/%s/%s.cif'
+    url = url_pattern % (restyp[:1], restyp, restyp)
+    cca = _ChemCompAtomHandler()
+    try:
+        with urllib.request.urlopen(url) as fh:
+            c = ihm.format.CifReader(fh,
+                                     category_handler={'_chem_comp_atom': cca})
+            c.read_file()
+    except urllib.error.URLError as exc:
+        warnings.warn(
+            "Component %s could not be found in CCD: %s" % (restyp, exc))
+    return cca.atoms
+
+
+def _get_non_canon(seen_atom_names, check_all):
     """Get all non-canonical atom names for each residue type"""
     for restyp, atoms in seen_atom_names.items():
-        # todo: if restyp not known, query Ligand Expo and parse the
-        # resulting mmCIF
+        if check_all and restyp not in KNOWN_ATOM_NAMES:
+            KNOWN_ATOM_NAMES.update(_get_non_std_restyp(restyp))
         if restyp in KNOWN_ATOM_NAMES:
             non_canon_atoms = atoms - KNOWN_ATOM_NAMES[restyp]
             if non_canon_atoms:
                 yield restyp, non_canon_atoms
 
 
-def _check_atom_names(model):
+def _check_atom_names(model, check_all):
     """Check that only standard atom names are used for known
        residue types"""
     seen_atom_names = collections.defaultdict(set)
     for atom in model._atoms:
         seq_id = 1 if atom.seq_id is None else atom.seq_id
         comp = atom.asym_unit.sequence[seq_id - 1]
         seen_atom_names[comp.id].add(atom.atom_id)
-    non_canon = sorted(_get_non_canon(seen_atom_names),
+    non_canon = sorted(_get_non_canon(seen_atom_names, check_all),
                        key=operator.itemgetter(0))
     if non_canon:
         raise ValueError(
@@ -370,11 +400,12 @@ def get_args():
     p.add_argument("--histidines", action='store_true', dest="fix_histidines",
                    help="Convert any non-standard histidine names (HIP, HID, "
                         "HIE, for different protonation states) to HIS")
-    p.add_argument('--check_atom_names', choices=['no', 'standard'],
+    p.add_argument('--check_atom_names', choices=['no', 'standard', 'all'],
                    dest="check_atom_names", default='no',
                    help="If 'standard', check for non-canonical atom names "
                         "in standard amino acid and nucleic acid chemical "
-                        "components")
+                        "components; if 'all', also check non-standard "
+                        "residue types by querying CCD (needs network access)")
     return p.parse_args()
 
 

test/input/non_canon_atom.cif

@@ -5,10 +5,20 @@ _exptl.method 'model, MODELLER Version 9.24 2020/08/21 11:54:31'
 _modeller.version 9.24
 #
 loop_
+_entity.id
+_entity.type
+1 polymer
+2 non-polymer
+3 non-polymer
+#
+loop_
 _struct_asym.id
 _struct_asym.entity_id
 _struct_asym.details
 A 1 ?
+B 2 ?
+C 2 ?
+D 3 ?
 #
 loop_
 _entity_poly_seq.entity_id
@@ -17,7 +27,13 @@ _entity_poly_seq.mon_id
 1 1 VAL
 1 2 GLY
 1 3 GLN
-1 4 MSE
+#
+loop_
+_pdbx_entity_nonpoly.entity_id
+_pdbx_entity_nonpoly.name
+_pdbx_entity_nonpoly.comp_id
+2 ? ZN
+3 ? invalid-comp-name
 #
 loop_
 _atom_site.group_PDB
@@ -38,8 +54,10 @@ _atom_site.B_iso_or_equiv
 _atom_site.label_entity_id
 _atom_site.id
 _atom_site.pdbx_PDB_model_num
-ATOM C bad1 . VAL A A 1 2 ? 114.370 27.980 -26.088 1.000 143.490 1 1 1
-ATOM C bad2 . VAL A A 1 2 ? 114.370 27.980 -26.088 1.000 143.490 1 2 1
-ATOM C CA   . GLY A A 2 3 ? 111.506 26.368 -28.075 1.000 137.530 1 3 1
-ATOM C bad3 . GLN A A 3 4 ? 113.468 23.113 -28.639 1.000 128.420 1 4 1
-ATOM C ig1  . MSE A A 4 5 ? 113.808 21.534 -32.168 1.000 117.620 1 5 1
+ATOM   C  bad1 . VAL A A 1 2 ? 114.370 27.980 -26.088 1.000 143.490 1 1 1
+ATOM   C  bad2 . VAL A A 1 2 ? 114.370 27.980 -26.088 1.000 143.490 1 2 1
+ATOM   C  CA   . GLY A A 2 3 ? 111.506 26.368 -28.075 1.000 137.530 1 3 1
+ATOM   C  bad3 . GLN A A 3 4 ? 113.468 23.113 -28.639 1.000 128.420 1 4 1
+HETATM ZN ZN   . ZN  B A . 5 ? 113.808 21.534 -32.168 1.000 117.620 2 5 1
+HETATM ZN bad4 . ZN  C A . 6 ? 113.808 21.534 -32.168 1.000 117.620 2 6 1
+HETATM ZN ZN   . invalid-comp-name D A . 7 ? 113.808 21.534 -32.168 3.000 117.620 3 7 1

test/mock/non_canon_atom/urllib/__init__.py

@@ -0,0 +1 @@
+# noop

test/mock/non_canon_atom/urllib/error.py

@@ -0,0 +1,6 @@
+class URLError(Exception):
+    pass
+
+
+class HTTPError(URLError):
+    pass

test/mock/non_canon_atom/urllib/request.py

@@ -0,0 +1,10 @@
+from io import BytesIO
+import urllib.error
+
+zinc_atoms = b'_chem_comp_atom.comp_id ZN\n_chem_comp_atom.atom_id ZN\n'
+
+
+def urlopen(url):
+    if 'invalid' in url:
+        raise urllib.error.HTTPError("404")
+    return BytesIO(zinc_atoms)

test/test_make_mmcif.py

@@ -257,8 +257,8 @@ def test_histidine(self):
                                          [False, False, False, False, True])
         os.unlink('output.cif')
 
-    def test_check_non_canon_atom(self):
-        """Test check for non-canonical atom names"""
+    def test_check_non_canon_atom_standard(self):
+        """Test check for non-canonical atom names, standard restypes"""
         incif = utils.get_input_file_name(TOPDIR, 'non_canon_atom.cif')
 
         # Should work fine without check
@@ -273,10 +273,37 @@ def test_check_non_canon_atom(self):
         out, err = r.communicate()
         self.assertEqual(r.returncode, 1)
         # Non-canonical atoms in standard residues should be reported
-        # Non-standard residues (MSE) are not checked
+        # Non-standard residues (ZN, ...) are not checked
         self.assertIn("Non-canonical atom names found in the following "
                       "residues: GLN: ['bad3']; VAL: ['bad1', 'bad2']",
                       err)
+        os.unlink('output.cif')
+
+    def test_check_non_canon_atom_all(self):
+        """Test check for non-canonical atom names, all restypes"""
+        incif = utils.get_input_file_name(TOPDIR, 'non_canon_atom.cif')
+
+        # Use mock urllib so we don't hit the network during this test
+        env = os.environ.copy()
+        mockdir = os.path.join(TOPDIR, 'test', 'mock', 'non_canon_atom')
+        env['PYTHONPATH'] = mockdir + os.pathsep + env['PYTHONPATH']
+
+        r = subprocess.Popen([sys.executable, MAKE_MMCIF,
+                             "--check_atom_names=all", incif],
+                             stdout=subprocess.PIPE,
+                             stderr=subprocess.PIPE,
+                             universal_newlines=True, env=env)
+        out, err = r.communicate()
+        self.assertEqual(r.returncode, 1)
+        # Non-canonical atoms in standard residues should be reported
+        # Non-standard residue (ZN) should also be checked
+        self.assertIn("Non-canonical atom names found in the following "
+                      "residues: GLN: ['bad3']; VAL: ['bad1', 'bad2']; "
+                      "ZN: ['bad4']", err)
+        # Residues not in CCD should give a warning
+        self.assertIn("Component invalid-comp-name could not be found in CCD",
+                      err)
+        os.unlink('output.cif')
 
 
 if __name__ == '__main__':