diff --git a/ihm/dumper.py b/ihm/dumper.py index cd892ca..d34dba8 100644 --- a/ihm/dumper.py +++ b/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): diff --git a/ihm/util/__init__.py b/ihm/util/__init__.py index 77c5693..ea7fa35 100644 --- a/ihm/util/__init__.py +++ b/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,60 @@ 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 (and no further ranges will overlap) + else: + break + 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]]""" diff --git a/test/test_dumper.py b/test/test_dumper.py index a2e8e49..671318c 100644 --- a/test/test_dumper.py +++ b/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,16 +1345,18 @@ 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) @@ -1361,12 +1364,14 @@ def test_poly_seq_scheme_dumper_not_modeled(self): 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 # """) diff --git a/test/test_examples.py b/test/test_examples.py index 03cecfc..3aa5ba1 100644 --- a/test/test_examples.py +++ b/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 diff --git a/test/test_util.py b/test/test_util.py index f353936..19c1026 100644 --- a/test/test_util.py +++ b/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 = []