Account for gaps in entity-reference alignment

When checking that the entity sequence matches
that in the reference database, consider gaps
in the alignment resulting from SeqDif records
that correspond to insertions or deletions.
Closes #181.

Author: benmwebb

ihm/dumper.py

@@ -552,9 +552,10 @@ def _get_sequence(self, reference):
             (code if len(code) == 1 else '(%s)' % code
              for code in fullrefseq[db_begin - 1:db_end]), 70))
 
-    def _check_seq_dif(self, entity, refseq, align):
+    def _check_seq_dif(self, entity, ref, refseq, align):
         """Check all SeqDif objects for the Entity sequence. Return the mutated
-           sequence (to match the reference)."""
+           sequence (to match the reference) plus sets of all insertion
+           and deletion positions."""
         entseq = [comp.code_canonical for comp in entity.sequence]
         mutations = []
         insertions = []
@@ -569,14 +570,15 @@ def _check_seq_dif(self, entity, refseq, align):
         for sd in insertions:
             self._check_seq_dif_entity(entity, entseq, sd)
         for sd in deletions:
-            self._check_seq_dif_reference(entity, refseq, sd)
+            self._check_seq_dif_reference(ref, refseq, sd)
         for sd in mutations:
             self._check_seq_dif_entity(entity, entseq, sd)
-            # todo: db_seq_id is often missing so we can't check reference
-            # info - unless we determine it automatically from seq_id
+            # we don't check db_monomer against the reference sequence here,
+            # but this will be done later when we check the alignment
             if sd.db_monomer:
                 entseq[sd.seq_id - 1] = sd.db_monomer.code_canonical
-        return entseq
+        return (entseq, frozenset(sd.seq_id for sd in insertions),
+                frozenset(sd.db_seq_id for sd in deletions))
 
     def _check_seq_dif_entity(self, entity, entseq, sd):
         """Make sure the entity information in a SeqDif record matches"""
@@ -591,18 +593,18 @@ def _check_seq_dif_entity(self, entity, entseq, sd):
                              % (sd.monomer.code_canonical, entity,
                                 entseq[sd.seq_id - 1], sd.seq_id))
 
-    def _check_seq_dif_reference(self, entity, refseq, sd):
+    def _check_seq_dif_reference(self, ref, refseq, sd):
         """Make sure the reference information in a SeqDif record matches"""
         if sd.db_seq_id < 1 or sd.db_seq_id > len(refseq):
             raise IndexError("SeqDif.db_seq_id for %s is %d, out of "
                              "range 1-%d"
-                             % (entity, sd.db_seq_id, len(refseq)))
+                             % (ref, sd.db_seq_id, len(refseq)))
         if (sd.db_monomer and len(refseq[sd.db_seq_id - 1]) == 1
                 and sd.db_monomer.code_canonical != refseq[sd.db_seq_id - 1]):
             raise ValueError("SeqDif.db_monomer one-letter code (%s) does "
                              "not match that in %s (%s at position %d)"
-                             % (sd.db_monomer.code_canonical, entity,
-                                refseq[sd.seq_id - 1], sd.db_seq_id))
+                             % (sd.db_monomer.code_canonical, ref,
+                                refseq[sd.db_seq_id - 1], sd.db_seq_id))
 
     def _get_ranges(self, entity, fullrefseq, align):
         """Get the sequence ranges for an Entity and Reference"""
@@ -612,6 +614,31 @@ def _get_ranges(self, entity, fullrefseq, align):
                 (align.db_begin,
                  len(fullrefseq) if align.db_end is None else align.db_end))
 
+    def _get_gapped_alignment(self, entity_rng, db_rng, entseq, refseq,
+                              insertions, deletions):
+        """Get the given ranges from the entity and reference sequences,
+           with gaps added to account for any insertion or deletion SeqDif
+           records"""
+        gapentseq = []
+        gaprefseq = []
+        ent_i = entity_rng[0]
+        ref_i = db_rng[0]
+        while ent_i <= entity_rng[1] and ref_i <= db_rng[1]:
+            if ent_i in insertions:
+                gapentseq.append(entseq[ent_i - 1])
+                ent_i += 1
+                gaprefseq.append('-')
+            elif ref_i in deletions:
+                gapentseq.append('-')
+                gaprefseq.append(refseq[ref_i - 1])
+                ref_i += 1
+            else:
+                gapentseq.append(entseq[ent_i - 1])
+                gaprefseq.append(refseq[ref_i - 1])
+                ref_i += 1
+                ent_i += 1
+        return gapentseq, gaprefseq
+
     def _check_reference_sequence(self, entity, ref):
         """Make sure that the Entity and Reference sequences match"""
         for align in ref._get_alignments():
@@ -625,7 +652,11 @@ def _check_alignment(self, entity, ref, align):
         # Reference sequence may contain non-standard residues, so parse them
         # out; e.g. "FLGHGGN(WP9)LHFVQLAS"
         fullrefseq = list(util._get_codes(ref.sequence))
-        entseq = self._check_seq_dif(entity, fullrefseq, align)
+
+        # Get mutated entity sequence, plus sets of insertions and deletions,
+        # from SeqDif records
+        entseq, insertions, deletions = self._check_seq_dif(
+            entity, ref, fullrefseq, align)
 
         def check_rng(rng, seq, rngstr, obj):
             if any(r < 1 or r > len(seq) for r in rng):
@@ -636,14 +667,13 @@ def check_rng(rng, seq, rngstr, obj):
         check_rng(entity_rng, entseq, "entity_begin,entity_end", entity)
         check_rng(db_rng, fullrefseq, "db_begin,db_end", ref)
 
-        matchlen = min(entity_rng[1] - entity_rng[0], db_rng[1] - db_rng[0])
-        entseq = entseq[entity_rng[0] - 1:entity_rng[0] + matchlen]
-        refseq = fullrefseq[db_rng[0] - 1:db_rng[0] + matchlen]
+        entseq, refseq = self._get_gapped_alignment(
+            entity_rng, db_rng, entseq, fullrefseq, insertions, deletions)
 
         # Entity sequence is canonical so likely won't match any non-standard
         # residue (anything of length > 1), so just skip checks of these
         def matchseq(a, b):
-            return a == b or len(a) > 1 or len(b) > 1
+            return a == b or len(a) > 1 or len(b) > 1 or a == '-' or b == '-'
         if (len(refseq) != len(entseq)
                 or not all(matchseq(a, b) for (a, b) in zip(refseq, entseq))):
             raise ValueError(

test/test_dumper.py

@@ -950,27 +950,105 @@ def test_struct_ref_seq_dif_mismatch(self):
         # Should work with checks disabled
         _ = _get_dumper_output(dumper, system, check=False)
 
+    def test_struct_ref_deletion_mismatch(self):
+        """Test StructRefDumper with deletion residue type mismatch"""
+        system = ihm.System()
+        lpep = ihm.LPeptideAlphabet()
+        sd = ihm.reference.DeletionSeqDif(db_seq_id=6, db_monomer=lpep['W'])
+        r = ihm.reference.UniProtSequence(
+            db_code='NUP84_YEAST', accession='P52891', sequence='MELWPTYQT',
+            details='test sequence')
+        r.alignments.append(ihm.reference.Alignment(seq_dif=[sd]))
+        system.entities.append(ihm.Entity('LSPT', references=[r]))
+        dumper = ihm.dumper._EntityDumper()
+        dumper.finalize(system)  # Assign entity IDs
+
+        dumper = ihm.dumper._StructRefDumper()
+        dumper.finalize(system)  # Assign IDs
+        with self.assertRaises(ValueError) as cm:
+            _get_dumper_output(dumper, system)
+        self.assertIn('one-letter code (W) does not match', str(cm.exception))
+        self.assertIn('(T at position 6)', str(cm.exception))
+        # Should work with checks disabled
+        _ = _get_dumper_output(dumper, system, check=False)
+
+    def test_struct_ref_deletion_out_of_range(self):
+        """Test StructRefDumper with deletion out of range"""
+        system = ihm.System()
+        lpep = ihm.LPeptideAlphabet()
+        sd = ihm.reference.DeletionSeqDif(db_seq_id=12, db_monomer=lpep['W'])
+        r = ihm.reference.UniProtSequence(
+            db_code='NUP84_YEAST', accession='P52891', sequence='MELWPTYQT',
+            details='test sequence')
+        r.alignments.append(ihm.reference.Alignment(seq_dif=[sd]))
+        system.entities.append(ihm.Entity('LSPT', references=[r]))
+        dumper = ihm.dumper._EntityDumper()
+        dumper.finalize(system)  # Assign entity IDs
+
+        dumper = ihm.dumper._StructRefDumper()
+        dumper.finalize(system)  # Assign IDs
+        with self.assertRaises(IndexError) as cm:
+            _get_dumper_output(dumper, system)
+        self.assertIn('SeqDif.db_seq_id for ', str(cm.exception))
+        self.assertIn('is 12, out of range 1-9', str(cm.exception))
+        # Should work with checks disabled
+        _ = _get_dumper_output(dumper, system, check=False)
+
     def test_struct_ref_seq_dif_ins_del(self):
         """Test StructRefDumper with SeqDif insertions and deletions"""
         system = ihm.System()
         lpep = ihm.LPeptideAlphabet()
-        sd1 = ihm.reference.SeqDif(seq_id=2, db_monomer=lpep['G'],
-                                   monomer=None, details='deletion')
-        sd2 = ihm.reference.SeqDif(seq_id=3, db_monomer=lpep['C'],
-                                   monomer=None, details='insertion')
+        sd1 = ihm.reference.DeletionSeqDif(db_seq_id=7, db_monomer=lpep['Q'])
+        sd2 = ihm.reference.InsertionSeqDif(
+            seq_id=1, monomer=lpep['H'], details='expression tag')
         r = ihm.reference.UniProtSequence(
             db_code='NUP84_YEAST', accession='P52891', sequence='MEWPTYQT',
             details='test sequence')
         r.alignments.append(ihm.reference.Alignment(seq_dif=[sd1, sd2]))
-        system.entities.append(ihm.Entity('MEWPTYQT', references=[r]))
+        system.entities.append(ihm.Entity('HMEWPTYT', references=[r]))
         dumper = ihm.dumper._EntityDumper()
         dumper.finalize(system)  # Assign entity IDs
 
         dumper = ihm.dumper._StructRefDumper()
         dumper.finalize(system)  # Assign IDs
         # Insertions and deletions are not currently checked, so
         # this should pass
-        _ = _get_dumper_output(dumper, system)
+        out = _get_dumper_output(dumper, system)
+        self.assertEqual(out, """#
+loop_
+_struct_ref.id
+_struct_ref.entity_id
+_struct_ref.db_name
+_struct_ref.db_code
+_struct_ref.pdbx_db_accession
+_struct_ref.pdbx_align_begin
+_struct_ref.pdbx_seq_one_letter_code
+_struct_ref.details
+1 1 UNP NUP84_YEAST P52891 1 MEWPTYQT 'test sequence'
+#
+#
+loop_
+_struct_ref_seq.align_id
+_struct_ref_seq.ref_id
+_struct_ref_seq.seq_align_beg
+_struct_ref_seq.seq_align_end
+_struct_ref_seq.db_align_beg
+_struct_ref_seq.db_align_end
+1 1 1 8 1 8
+#
+#
+loop_
+_struct_ref_seq_dif.pdbx_ordinal
+_struct_ref_seq_dif.align_id
+_struct_ref_seq_dif.db_mon_id
+_struct_ref_seq_dif.pdbx_seq_db_seq_num
+_struct_ref_seq_dif.mon_id
+_struct_ref_seq_dif.seq_num
+_struct_ref_seq_dif.details
+1 1 GLN 7 ? ? deletion
+2 1 ? ? HIS 1 'expression tag'
+#
+""")
 
     def test_chem_comp_dumper(self):
         """Test ChemCompDumper"""