Support pdbx_poly_seq_scheme heterogeneous composition

The _pdbx_poly_seq_scheme table now requires
a new ihm_model_id_list data item that lists all
models that have coordinates for a particular residue.
This is to enable PDB tools to fully support
heterogeneous composition. Use our existing
not-modeled-residue ranges to populate this data
item on output.

Author: benmwebb

ihm/dumper.py

@@ -914,21 +914,22 @@ def dump(self, system, writer):
                          ["asym_id", "entity_id", "seq_id", "mon_id",
                           "pdb_seq_num", "auth_seq_num", "pdb_mon_id",
                           "auth_mon_id", "pdb_strand_id",
-                          "pdb_ins_code"]) as lp:
+                          "pdb_ins_code", "ihm_model_id_list"]) as lp:
             for asym in system.asym_units:
                 entity = asym.entity
                 if not entity.is_polymeric():
                     continue
-                for start, end, modeled in self._get_ranges(system, asym):
+                for (start, end, modeled,
+                     model_ids) in self._get_ranges(system, asym):
                     for num in range(start, end + 1):
                         comp = entity.sequence[num - 1]
                         auth_comp_id = comp.id
                         pdb_seq_num, auth_seq_num, ins = \
                             asym._get_pdb_auth_seq_id_ins_code(num)
                         if not modeled:
-                            # If a residue wasn't modeled, PDB convention is
-                            # to state ? for auth_seq_num, pdb_mon_id,
-                            # auth_mon_id.
+                            # If a residue wasn't modeled in any models,
+                            # PDB convention is to state ? for auth_seq_num,
+                            # pdb_mon_id, auth_mon_id.
                             # See, e.g., https://files.rcsb.org/view/8QB4.cif
                             auth_comp_id = ihm.unknown
                             auth_seq_num = ihm.unknown
@@ -942,13 +943,16 @@ def dump(self, system, writer):
                                  pdb_seq_num=pdb_seq_num,
                                  auth_seq_num=auth_seq_num, mon_id=comp.id,
                                  pdb_mon_id=auth_comp_id,
-                                 auth_mon_id=auth_comp_id, pdb_ins_code=ins)
+                                 auth_mon_id=auth_comp_id, pdb_ins_code=ins,
+                                 ihm_model_id_list=model_ids)
 
     def _get_ranges(self, system, asym):
-        """Get a list of (seq_id_begin, seq_id_end, modeled) residue ranges
-           for the given asym. The list is guaranteed to be sorted and to cover
-           all residues in the asym. `modeled` is True if no Model has any
-           residue in that range in a NotModeledResidueRange."""
+        """Get a list of (seq_id_begin, seq_id_end, modeled, model_ids)
+           residue ranges for the given asym. The list is guaranteed to be
+           sorted and to cover all residues in the asym. `model_ids` is a list
+           of all Models (as a comma-separated string of model IDs) that have
+           no residue in that range in a NotModeledResidueRange, or None if
+           the residues are not modeled in any Model."""
         _all_modeled = []
         num_models = 0
         for group, model in system._all_models():
@@ -964,14 +968,15 @@ def _get_ranges(self, system, asym):
                 (rr.seq_id_begin, rr.seq_id_end)
                 for rr in ranges if rr.asym_unit is asym)
             # Invert to get a list of modeled ranges for this model
-            _all_modeled.extend(util._invert_ranges(_all_not_modeled,
-                                                    len(asym.entity.sequence)))
-        # If no models, there are no "not modeled residues", so say everything
-        # was modeled
+            for rng in util._invert_ranges(_all_not_modeled,
+                                           len(asym.entity.sequence)):
+                _all_modeled.append((rng[0], rng[1], model._id))
+        # If no models, there are no explicitly "not modeled residues",
+        # so say everything was modeled (even though there are no model IDs)
         if num_models == 0:
-            _all_modeled = [(1, len(asym.entity.sequence))]
-        return util._pred_ranges(util._combine_ranges(_all_modeled),
-                                 len(asym.entity.sequence))
+            return [(1, len(asym.entity.sequence), True, None)]
+        return util._pred_id_ranges(util._combine_id_ranges(_all_modeled),
+                                    len(asym.entity.sequence))
 
 
 class _NonPolySchemeDumper(Dumper):

ihm/util/__init__.py

@@ -223,20 +223,22 @@ def _invert_ranges(ranges, end, start=1):
         yield (start, end)
 
 
-def _pred_ranges(ranges, end):
-    """Given a sorted list of non-overlapping ranges, yield a new list which
-       covers the range 1-end. Each element in the new list contains a new
-       third bool member which is True iff the element was in the original
-       list. For example, if end=4,
-       [(2, 3)] -> [(1, 1, False), (2, 3, True), (4, 4, False)]"""
+def _pred_id_ranges(id_ranges, end):
+    """Given a sorted list of ranges with IDs, yield a new list which covers
+       the range 1-end. The IDs in the new list are the IDs from the original
+       list, or None. In addition, a new bool item is returned, True only if
+       there is at least one model ID.
+       For example, if end=4,
+       [(2, 3, 'x')] -> [(1, 1, False, None), (2, 3, True, 'x'),
+                         (4, 4, False, None)]"""
     start = 1
-    for r in ranges:
+    for r in id_ranges:
         if r[0] > start:
-            yield (start, r[0] - 1, False)
-        yield (r[0], r[1], True)
+            yield (start, r[0] - 1, False, None)
+        yield (r[0], r[1], True, r[2])
         start = r[1] + 1
     if end >= start:
-        yield (start, end, False)
+        yield (start, end, False, None)
 
 
 def _combine_ranges(ranges):
@@ -255,6 +257,58 @@ def _combine_ranges(ranges):
     yield current
 
 
+def _sorted_nonoverlap_id_ranges(id_ranges):
+    """Sort the input ranges with IDs and remove any overlaps with the
+       same ID; yield the result.
+       For example, [(8, 10, 'x'), (1, 2, 'y'), (3, 4, 'y'), (5, 6, 'z')]
+                    -> [(1, 4, 'y'), (5, 6, 'z'), (8, 10, 'x')]"""
+    id_ranges = sorted(id_ranges)
+    if not id_ranges:
+        return
+    current = id_ranges[0]
+    for r in id_ranges[1:]:
+        if current[1] + 1 >= r[0] and r[2] == current[2]:
+            current = (current[0], max(r[1], current[1]), current[2])
+        else:
+            yield current
+            current = r
+    yield current
+
+
+def _combine_id_ranges(id_ranges):
+    """Sort the input ranges with IDs into non-overlapping ranges per ID;
+       yield the result.
+       For example, [(8, 10, 'x'), (1, 2, 'y'), (3, 4, 'z'), (2, 3, 'z')]
+                    -> [(1, 1, 'y'), (2, 2, 'y,z'), (3, 4, 'z'),
+                        (8, 10, 'x')]"""
+    id_ranges = list(_sorted_nonoverlap_id_ranges(id_ranges))
+    if not id_ranges:
+        return
+    yield_end = 0
+    for i_r, r in enumerate(id_ranges):
+        start, end, r_id = r
+        start = max(start, yield_end + 1)
+        # Get subranges of r and their IDs
+        while start <= end:
+            ids = set([r_id])
+            new_end = end
+            for r2 in id_ranges[i_r + 1:]:
+                start2, end2, id2 = r2
+                # If r2 has same start as r1, adjust end to cover overlap
+                if start2 <= start:
+                    if end2 >= start:
+                        new_end = min(new_end, end2)
+                        ids.add(id2)
+                # If only partial overlap, exclude this range; we will
+                # look at it again later
+                elif start2 <= new_end:
+                    new_end = start2 - 1
+                # Otherwise, no overlap
+            yield start, new_end, ",".join("%s" % d for d in sorted(ids))
+            yield_end = new_end
+            start = new_end + 1
+
+
 def _make_range_from_list(rr):
     """Yield a list of ranges given a sorted list of values.
        For example, [1, 2, 5, 6] -> [[1, 2], [5, 6]]"""

test/test_dumper.py

@@ -1303,17 +1303,18 @@ def test_poly_seq_scheme_dumper(self):
 _pdbx_poly_seq_scheme.auth_mon_id
 _pdbx_poly_seq_scheme.pdb_strand_id
 _pdbx_poly_seq_scheme.pdb_ins_code
-A 1 1 ALA 1 1 ALA ALA A .
-A 1 2 CYS 2 2 CYS CYS A .
-A 1 3 GLY 3 3 GLY GLY A .
-A 1 4 THR 4 4 THR THR A .
-B 2 1 ALA 6 6 ALA ALA B .
-B 2 2 CYS 7 7 CYS CYS B .
-B 2 3 CYS 8 8 CYS CYS B .
-C 3 1 A 1 1 A A C .
-C 3 2 C 2 2 C C C .
-D 4 1 DA 1 1 DA DA X A
-D 4 2 DC 1 1 DC DC X B
+_pdbx_poly_seq_scheme.ihm_model_id_list
+A 1 1 ALA 1 1 ALA ALA A . .
+A 1 2 CYS 2 2 CYS CYS A . .
+A 1 3 GLY 3 3 GLY GLY A . .
+A 1 4 THR 4 4 THR THR A . .
+B 2 1 ALA 6 6 ALA ALA B . .
+B 2 2 CYS 7 7 CYS CYS B . .
+B 2 3 CYS 8 8 CYS CYS B . .
+C 3 1 A 1 1 A A C . .
+C 3 2 C 2 2 C C C . .
+D 4 1 DA 1 1 DA DA X A .
+D 4 2 DC 1 1 DC DC X B .
 #
 """)
 
@@ -1344,29 +1345,33 @@ def test_poly_seq_scheme_unknown_auth_seq(self):
 _pdbx_poly_seq_scheme.auth_mon_id
 _pdbx_poly_seq_scheme.pdb_strand_id
 _pdbx_poly_seq_scheme.pdb_ins_code
-A 1 1 ALA 1 3 ALA ALA A .
-A 1 2 CYS 2 4 CYS CYS A .
-A 1 3 GLY 3 ? ? ? A .
-A 1 4 THR 4 6 THR THR A .
+_pdbx_poly_seq_scheme.ihm_model_id_list
+A 1 1 ALA 1 3 ALA ALA A . .
+A 1 2 CYS 2 4 CYS CYS A . .
+A 1 3 GLY 3 ? ? ? A . .
+A 1 4 THR 4 6 THR THR A . .
 #
 """)
 
     def test_poly_seq_scheme_dumper_not_modeled(self):
         """Test PolySeqSchemeDumper with not-modeled residues"""
         system, m1, asym = self._make_test_model()
+        m1._id = 1
         del asym.entity._id
         rr = ihm.model.NotModeledResidueRange(asym, 1, 2)
         m1.not_modeled_residue_ranges.append(rr)
 
         m2 = ihm.model.Model(assembly=m1.assembly, protocol=m1.protocol,
                              representation=m1.representation,
                              name='2nd test model')
+        m2._id = 2
         rr = ihm.model.NotModeledResidueRange(asym, 2, 4)
         m2.not_modeled_residue_ranges.append(rr)
 
         m3 = ihm.model.Model(assembly=m1.assembly, protocol=m1.protocol,
                              representation=m1.representation,
                              name='3rd test model')
+        m3._id = 3
         rr = ihm.model.NotModeledResidueRange(asym, 2, 3)
         m3.not_modeled_residue_ranges.append(rr)
 
@@ -1390,10 +1395,11 @@ def test_poly_seq_scheme_dumper_not_modeled(self):
 _pdbx_poly_seq_scheme.auth_mon_id
 _pdbx_poly_seq_scheme.pdb_strand_id
 _pdbx_poly_seq_scheme.pdb_ins_code
-A 1 1 ALA 1 1 ALA ALA A .
-A 1 2 CYS 2 ? ? ? A .
-A 1 3 GLY 3 3 GLY GLY A .
-A 1 4 THR 4 4 THR THR A .
+_pdbx_poly_seq_scheme.ihm_model_id_list
+A 1 1 ALA 1 1 ALA ALA A . 2,3
+A 1 2 CYS 2 ? ? ? A . .
+A 1 3 GLY 3 3 GLY GLY A . 1
+A 1 4 THR 4 4 THR THR A . 1,3
 #
 """)
 
@@ -1404,6 +1410,7 @@ def test_poly_seq_scheme_dumper_no_not_modeled(self):
         system, m1, asym = self._make_test_model()
         del asym.entity._id
         del m1.not_modeled_residue_ranges
+        m1._id = 1
 
         mg = system.state_groups[0][0][0]
         mg.append(m1)
@@ -1425,10 +1432,11 @@ def test_poly_seq_scheme_dumper_no_not_modeled(self):
 _pdbx_poly_seq_scheme.auth_mon_id
 _pdbx_poly_seq_scheme.pdb_strand_id
 _pdbx_poly_seq_scheme.pdb_ins_code
-A 1 1 ALA 1 1 ALA ALA A .
-A 1 2 CYS 2 2 CYS CYS A .
-A 1 3 GLY 3 3 GLY GLY A .
-A 1 4 THR 4 4 THR THR A .
+_pdbx_poly_seq_scheme.ihm_model_id_list
+A 1 1 ALA 1 1 ALA ALA A . 1
+A 1 2 CYS 2 2 CYS CYS A . 1
+A 1 3 GLY 3 3 GLY GLY A . 1
+A 1 4 THR 4 4 THR THR A . 1
 #
 """)
 

test/test_examples.py

@@ -41,7 +41,7 @@ def test_simple_docking_example(self):
             # can read it
             with open(os.path.join(tmpdir, 'output.cif')) as fh:
                 contents = fh.readlines()
-            self.assertEqual(len(contents), 321)
+            self.assertEqual(len(contents), 322)
             with open(os.path.join(tmpdir, 'output.cif')) as fh:
                 s, = ihm.reader.read(fh)
 
@@ -70,7 +70,7 @@ def test_ligands_water_example(self):
         # can read it
         with open(out) as fh:
             contents = fh.readlines()
-        self.assertEqual(len(contents), 255)
+        self.assertEqual(len(contents), 256)
         with open(out) as fh:
             s, = ihm.reader.read(fh)
         # Make sure that resulting Python objects are picklable

test/test_util.py

@@ -199,19 +199,20 @@ def test_invert_ranges(self):
         self.assertEqual(list(ihm.util._invert_ranges(inrng, 4, start=1)),
                          [(1, 4)])
 
-    def test_pred_ranges(self):
-        """Test _pred_ranges function"""
-        inrng = [(2, 3)]
-        self.assertEqual(list(ihm.util._pred_ranges(inrng, 4)),
-                         [(1, 1, False), (2, 3, True), (4, 4, False)])
-        inrng = [(1, 1), (4, 7)]
-        self.assertEqual(list(ihm.util._pred_ranges(inrng, 8)),
-                         [(1, 1, True), (2, 3, False), (4, 7, True),
-                          (8, 8, False)])
-        inrng = [(2, 2), (4, 7)]
-        self.assertEqual(list(ihm.util._pred_ranges(inrng, 7)),
-                         [(1, 1, False), (2, 2, True), (3, 3, False),
-                          (4, 7, True)])
+    def test_pred_id_ranges(self):
+        """Test _pred_id_ranges function"""
+        inrng = [(2, 3, 'a')]
+        self.assertEqual(list(ihm.util._pred_id_ranges(inrng, 4)),
+                         [(1, 1, False, None), (2, 3, True, 'a'),
+                          (4, 4, False, None)])
+        inrng = [(1, 1, 'a'), (4, 7, 'b')]
+        self.assertEqual(list(ihm.util._pred_id_ranges(inrng, 8)),
+                         [(1, 1, True, 'a'), (2, 3, False, None),
+                          (4, 7, True, 'b'), (8, 8, False, None)])
+        inrng = [(2, 2, 'a'), (4, 7, 'b')]
+        self.assertEqual(list(ihm.util._pred_id_ranges(inrng, 7)),
+                         [(1, 1, False, None), (2, 2, True, 'a'),
+                          (3, 3, False, None), (4, 7, True, 'b')])
 
     def test_combine_ranges(self):
         """Test _combine_ranges function"""
@@ -227,6 +228,36 @@ def test_combine_ranges(self):
                          [(1, 2), (4, 4)])
         self.assertEqual(list(ihm.util._combine_ranges([])), [])
 
+    def test_sorted_nonoverlap_id_ranges(self):
+        """Test _sorted_nonoverlap_id_ranges function"""
+        inrng = [(8, 10, 'a'), (1, 2, 'a'), (3, 4, 'a')]
+        self.assertEqual(list(ihm.util._sorted_nonoverlap_id_ranges(inrng)),
+                         [(1, 4, 'a'), (8, 10, 'a')])
+        inrng = [(8, 10, 'a'), (1, 2, 'a'), (3, 4, 'b')]
+        self.assertEqual(list(ihm.util._sorted_nonoverlap_id_ranges(inrng)),
+                         [(1, 2, 'a'), (3, 4, 'b'), (8, 10, 'a')])
+        inrng = [(1, 10, 'a'), (3, 4, 'a')]
+        self.assertEqual(list(ihm.util._sorted_nonoverlap_id_ranges(inrng)),
+                         [(1, 10, 'a')])
+        inrng = [(1, 10, 'a'), (3, 4, 'b')]
+        self.assertEqual(list(ihm.util._sorted_nonoverlap_id_ranges(inrng)),
+                         [(1, 10, 'a'), (3, 4, 'b')])
+        self.assertEqual(list(ihm.util._sorted_nonoverlap_id_ranges([])), [])
+
+    def test_combine_id_ranges(self):
+        """Test _combine_id_ranges function"""
+        inrng = [(1, 30, 'd'), (1, 50, 'a'), (1, 60, 'c'), (2, 40, 'b'),
+                 (80, 100, 'e')]
+        self.assertEqual(list(ihm.util._combine_id_ranges(inrng)),
+                         [(1, 1, 'a,c,d'), (2, 30, 'a,b,c,d'),
+                          (31, 40, 'a,b,c'), (41, 50, 'a,c'),
+                          (51, 60, 'c'), (80, 100, 'e')])
+        inrng = [(8, 10, 'x'), (1, 2, 'y'), (3, 4, 'z'), (2, 3, 'z')]
+        self.assertEqual(list(ihm.util._combine_id_ranges(inrng)),
+                         [(1, 1, 'y'), (2, 2, 'y,z'),
+                          (3, 4, 'z'), (8, 10, 'x')])
+        self.assertEqual(list(ihm.util._combine_id_ranges([])), [])
+
     def test_make_range_from_list(self):
         """Test _make_range_from_list function"""
         rr = []