Fill in missing chem_comp data in make_mmcif

Add a new --fix_chem_comp flag to make_mmcif
which, if specified, will fill in any missing
information in the chem_comp table (type, name,
formula) by querying CCD. Closes #179.

Author: benmwebb

ihm/util/make_mmcif.py

@@ -34,6 +34,7 @@
 import collections
 import operator
 import warnings
+import inspect
 
 
 # All canonical atom names for each standard residue type, as per CCD.
@@ -118,7 +119,7 @@
 }
 
 
-def add_ihm_info(s, fix_histidines, check_atom_names):
+def add_ihm_info(s, fix_histidines, check_atom_names, fix_chem_comp):
     # Non-standard histidine names (protonation states)
     histidines = frozenset(('HIP', 'HID', 'HIE'))
 
@@ -155,6 +156,8 @@ def add_ihm_info(s, fix_histidines, check_atom_names):
     if fix_histidines:
         _fix_histidine_chem_comps(s, histidines)
     _fix_empty_assemblies(s)
+    if fix_chem_comp:
+        _fix_incomplete_chem_comps(s)
     return s
 
 
@@ -180,20 +183,25 @@ 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"""
+def _read_ccd(restyp, category_handler):
+    """Read the given residue type from CCD using the given CIF handler(s)"""
     url_top = 'https://files.rcsb.org'
     url_pattern = url_top + '/pub/pdb/refdata/chem_comp/%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 = ihm.format.CifReader(
+                fh, category_handler=category_handler)
             c.read_file()
     except urllib.error.URLError as exc:
         warnings.warn(
             "Component %s could not be found in CCD: %s" % (restyp, exc))
+
+
+def _get_non_std_restyp(restyp):
+    """Return CCD atom info for the given residue type"""
+    cca = _ChemCompAtomHandler()
+    _read_ccd(restyp, {'_chem_comp_atom': cca})
     return cca.atoms
 
 
@@ -290,15 +298,48 @@ def _fix_empty_assemblies(s):
                 s.complete_assembly.description = asmb.description
 
 
-def add_ihm_info_one_system(fname, fix_histidines, check_atom_names):
+class _ChemCompHandler:
+    """Read the _chem_comp table from a CCD entry"""
+    not_in_file = omitted = unknown = None
+
+    def __call__(self, name, type, formula):
+        self.name, self.type, self.formula = name, type, formula
+
+
+def _fix_incomplete_chem_comps(s):
+    """Add any missing information to ChemComps using CCD"""
+    # Map type to ChemComp subclass. Map nonpolymer to NonPolyChemComp,
+    # not WaterChemComp
+    typmap = {x[1].type.lower(): x[1]
+              for x in inspect.getmembers(ihm, inspect.isclass)
+              if issubclass(x[1], ihm.ChemComp)
+              and x[1] is not ihm.WaterChemComp}
+    for cc in s._orphan_chem_comps:
+        if cc.type == 'other' or cc.name is None or cc.formula is None:
+            _fix_chem_comp(cc, typmap)
+
+
+def _fix_chem_comp(cc, typmap):
+    """Add missing information to a single ChemComp from CCD"""
+    h = _ChemCompHandler()
+    _read_ccd(cc.id.upper(), {'_chem_comp': h})
+    if hasattr(h, 'name') and h.name is not None:
+        cc.name = h.name
+        cc.formula = h.formula
+        cc.__class__ = typmap.get(h.type.lower(), ihm.ChemComp)
+
+
+def add_ihm_info_one_system(fname, fix_histidines, check_atom_names,
+                            fix_chem_comp):
     """Read mmCIF file `fname`, which must contain a single System, and
        return it with any missing IHM data added."""
     with open(fname) as fh:
         systems = ihm.reader.read(fh)
     if len(systems) != 1:
         raise ValueError("mmCIF file %s must contain exactly 1 data block "
                          "(%d found)" % (fname, len(systems)))
-    return add_ihm_info(systems[0], fix_histidines, check_atom_names)
+    return add_ihm_info(systems[0], fix_histidines, check_atom_names,
+                        fix_chem_comp)
 
 
 def combine(s, other_s):
@@ -422,6 +463,10 @@ def get_args():
                         "in standard amino acid and nucleic acid chemical "
                         "components; if 'all', also check non-standard "
                         "residue types by querying CCD (needs network access)")
+    p.add_argument("--fix_chem_comp", action='store_true',
+                   dest="fix_chem_comp",
+                   help="Add any missing data to the chem_comp table by"
+                        "querying CCD (needs network access)")
     return p.parse_args()
 
 
@@ -434,10 +479,12 @@ def main():
 
     if args.add:
         s = add_ihm_info_one_system(args.input, args.fix_histidines,
-                                    args.check_atom_names)
+                                    args.check_atom_names,
+                                    args.fix_chem_comp)
         for other in args.add:
             other_s = add_ihm_info_one_system(other, args.fix_histidines,
-                                              args.check_atom_names)
+                                              args.check_atom_names,
+                                              args.fix_chem_comp)
             combine(s, other_s)
         with open(args.output, 'w') as fhout:
             ihm.dumper.write(
@@ -448,7 +495,8 @@ def main():
             with open(args.output, 'w') as fhout:
                 ihm.dumper.write(
                     fhout, [add_ihm_info(s, args.fix_histidines,
-                                         args.check_atom_names)
+                                         args.check_atom_names,
+                                         args.fix_chem_comp)
                             for s in ihm.reader.read(fh)],
                     variant=ihm.dumper.IgnoreVariant(['_audit_conform']))
 

test/input/missing_chem_comp.cif

@@ -0,0 +1,10 @@
+loop_
+_chem_comp.id
+_chem_comp.type
+_chem_comp.name
+_chem_comp.formula
+_chem_comp.formula_weight
+ALA 'L-peptide linking' ALANINE 'C3 H7 N O2' 89.094
+MG . . . .
+ZN . . . .
+invalid-chem-comp . . . .

test/mock/non_canon_atom/urllib/request.py

@@ -3,8 +3,16 @@
 
 zinc_atoms = b'_chem_comp_atom.comp_id ZN\n_chem_comp_atom.atom_id ZN\n'
 
+mg = b"""
+_chem_comp.id                                   MG
+_chem_comp.name                                 "MAGNESIUM ION"
+_chem_comp.type                                 NON-POLYMER
+_chem_comp.pdbx_type                            HETAI
+_chem_comp.formula                              Mg
+"""
+
 
 def urlopen(url):
-    if 'invalid' in url:
+    if 'invalid' in url or 'INVALID' in url:
         raise urllib.error.HTTPError("404")
-    return BytesIO(zinc_atoms)
+    return BytesIO(zinc_atoms if url.endswith('ZN.cif') else mg)

test/test_make_mmcif.py

@@ -352,6 +352,39 @@ def test_empty_assembly_no_description(self):
         os.unlink(incif_new)
         os.unlink('output.cif')
 
+    def test_missing_chem_comp(self):
+        """Test fix of incomplete chem_comp table"""
+        incif = utils.get_input_file_name(TOPDIR, 'missing_chem_comp.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,
+                             "--fix_chem_comp", incif],
+                             stdout=subprocess.PIPE,
+                             stderr=subprocess.PIPE,
+                             universal_newlines=True, env=env)
+        out, err = r.communicate()
+        self.assertEqual(r.returncode, 0)
+        # Residues not in CCD should give a warning
+        self.assertIn("Component INVALID-CHEM-COMP could not be found in CCD",
+                      err)
+        # ALA should be left unchanged (already present in the input file);
+        # MG should be filled in with CCD info;
+        # ZN should be left unchanged (no chem_comp table in our mock CCD)
+        with open('output.cif') as fh:
+            contents = fh.readlines()
+        ind = contents.index('_chem_comp.formula_weight\n')
+        self.assertEqual(
+            contents[ind + 1:ind + 5],
+            ["ALA 'L-peptide linking' ALANINE 'C3 H7 N O2' 89.094\n",
+             "MG non-polymer 'MAGNESIUM ION' Mg 24.305\n",
+             'ZN other . . .\n',
+             'invalid-chem-comp other . . .\n'])
+        os.unlink('output.cif')
+
 
 if __name__ == '__main__':
     unittest.main()