Fix StopIteration in alignChains with Biopython > 1.79

[Copilot]

Biopython > 1.79's PairwiseAligner can produce alignment strings with more non-gap characters than actual residues in the sequences. This causes getAlignedMapping() to exhaust sequence iterators and raise unhandled StopIteration.

Changes

• Added StopIteration handling in getAlignedMapping() around all next() calls on sequence iterators
• Synchronized list lengths: When chain iterator exhausts after adding target residue, pop from amatch to prevent length mismatch between amatch and bmatch
• Added debug logging: Warning messages when alignment truncation occurs, showing which sequence and position

Key Fix

# Before: Unhandled StopIteration
ares = next(aiter)
amatch.append(ares.getResidue())
if b not in gap_chars:
    bres = next(biter)  # Can raise StopIteration
    bmatch.append(bres.getResidue())

# After: Graceful handling with synchronization
ares = next(aiter)
amatch.append(ares.getResidue())
if b not in gap_chars:
    try:
        bres = next(biter)
    except StopIteration:
        amatch.pop()  # Keep lists synchronized
        LOGGER.warning(...)
        break
    bmatch.append(bres.getResidue())

This maintains backward compatibility while handling alignment format differences across Biopython versions.

Original prompt

---

This section details on the original issue you should resolve

<issue_title>Biopython > 1.79 sometimes breaks alignChains</issue_title>
<issue_description>## Description of the bug.
When I run alignChains on the attached pdbs with Python 3.11 with Biopython > 1.79, I get the following error. With Biopython 1.79, the same code runs without any problems. This also means I can't use Python 3.12 as it doesn't support Biopython 1.79

ref_model.pdb
centered_model.pdb

Error log from code run shown below

@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain A from centered_model:
@>      Mapped: 1132 residues match with 100% sequence identity and 99% overlap.
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain A from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain B from ref_model (len=381) with Chain A from centered_model:
@>      Failed to match chains (seqid=33%, overlap=32%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain B from ref_model (len=381) with Chain A from centered_model:
@>      Failed to match chains (seqid=12%, overlap=6%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain E from ref_model (len=709) with Chain A from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain E from ref_model (len=709) with Chain A from centered_model:
@>      Failed to match chains (seqid=34%, overlap=58%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain E from ref_model (len=709) with Chain A from centered_model:
@>      Failed to match chains (seqid=3%, overlap=6%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain F from ref_model (len=21) with Chain A from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain F from ref_model (len=21) with Chain A from centered_model:
@>      Failed to match chains (seqid=40%, overlap=2%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain F from ref_model (len=21) with Chain A from centered_model:
@>      Failed to match chains (seqid=12%, overlap=1%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain B from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain A from ref_model (len=1132) with Chain B from centered_model:
@>      Failed to match chains (seqid=34%, overlap=29%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain A from ref_model (len=1132) with Chain B from centered_model:
@>      Failed to match chains (seqid=11%, overlap=6%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain B from centered_model:
@>      Mapped: 346 residues match with 100% sequence identity and 99% overlap.
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain E from ref_model (len=709) with Chain B from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain E from ref_model (len=709) with Chain B from centered_model:
@>      Failed to match chains (seqid=29%, overlap=46%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain E from ref_model (len=709) with Chain B from centered_model:
@>      Failed to match chains (seqid=5%, overlap=15%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain F from ref_model (len=21) with Chain B from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain F from ref_model (len=21) with Chain B from centered_model:
@>      Failed to match chains (seqid=28%, overlap=5%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain F from ref_model (len=21) with Chain B from centered_model:
@>      Failed to match chains (seqid=6%, overlap=5%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain E from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain A from ref_model (len=1132) with Chain E from centered_model:
@>      Failed to match chains (seqid=34%, overlap=58%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain A from ref_model (len=1132) with Chain E from centered_model:
@>      Failed to match chains (seqid=5%, overlap=7%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain E from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain B from re...

</details>



<!-- START COPILOT CODING AGENT SUFFIX -->

- Fixes prody/ProDy#2205

<!-- START COPILOT CODING AGENT TIPS -->
---

✨ Let Copilot coding agent [set things up for you](https://github.com/prody/ProDy/issues/new?title=✨+Set+up+Copilot+instructions&body=Configure%20instructions%20for%20this%20repository%20as%20documented%20in%20%5BBest%20practices%20for%20Copilot%20coding%20agent%20in%20your%20repository%5D%28https://gh.io/copilot-coding-agent-tips%29%2E%0A%0A%3COnboard%20this%20repo%3E&assignees=copilot) — coding agent works faster and does higher quality work when set up for your repo.

[AnthonyBogetti]

@jamesmkrieger This looks reasonable to me. Can you test it and see if it works? If so, let's merge it.

[jamesmkrieger]

Will do. I agree it seems reasonable. Thanks

[jamesmkrieger]

No, this doesn't work. I now get another error

In [1]: import prody
/Users/kriegej/scipion3/software/em/prody-github/ProDy/prody/drugui/gui.py:1867: SyntaxWarning: invalid escape sequence '\$'
  """
/Users/kriegej/scipion3/software/em/prody-github/ProDy/prody/drugui/no_gui.py:1160: SyntaxWarning: invalid escape sequence '\$'
  """

In [2]: prody_ref_model = prody.parsePDB("ref_model.pdb")
@> 2243 atoms and 1 coordinate set(s) were parsed in 0.01s.

In [3]: prody_model = prody.parsePDB("centered_model.pdb")
@> 2280 atoms and 1 coordinate set(s) were parsed in 0.01s.

In [4]: prody_ref_amap = prody.alignChains(prody_ref_model, prody_model, overlap=15, matchFunc=prody.sameChid)
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain A from centered_model:
@>      Mapped: 1132 residues match with 100% sequence identity and 99% overlap.
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain A from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain B from ref_model (len=381) with Chain A from centered_model:
@>      Failed to match chains (seqid=33%, overlap=32%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain B from ref_model (len=381) with Chain A from centered_model:
@>      Failed to match chains (seqid=12%, overlap=6%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain E from ref_model (len=709) with Chain A from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain E from ref_model (len=709) with Chain A from centered_model:
@>      Failed to match chains (seqid=34%, overlap=58%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain E from ref_model (len=709) with Chain A from centered_model:
@>      Failed to match chains (seqid=3%, overlap=6%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain F from ref_model (len=21) with Chain A from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain F from ref_model (len=21) with Chain A from centered_model:
@>      Failed to match chains (seqid=40%, overlap=2%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain F from ref_model (len=21) with Chain A from centered_model:
@>      Failed to match chains (seqid=12%, overlap=1%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain B from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain A from ref_model (len=1132) with Chain B from centered_model:
@>      Failed to match chains (seqid=34%, overlap=29%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain A from ref_model (len=1132) with Chain B from centered_model:
@>      Failed to match chains (seqid=11%, overlap=6%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain B from centered_model:
@>      Mapped: 346 residues match with 100% sequence identity and 99% overlap.
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain E from ref_model (len=709) with Chain B from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain E from ref_model (len=709) with Chain B from centered_model:
@>      Failed to match chains (seqid=29%, overlap=46%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain E from ref_model (len=709) with Chain B from centered_model:
@>      Failed to match chains (seqid=5%, overlap=15%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain F from ref_model (len=21) with Chain B from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain F from ref_model (len=21) with Chain B from centered_model:
@>      Failed to match chains (seqid=28%, overlap=5%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain F from ref_model (len=21) with Chain B from centered_model:
@>      Failed to match chains (seqid=6%, overlap=5%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain E from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain A from ref_model (len=1132) with Chain E from centered_model:
@>      Failed to match chains (seqid=34%, overlap=58%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain A from ref_model (len=1132) with Chain E from centered_model:
@>      Failed to match chains (seqid=5%, overlap=7%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain E from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain B from ref_model (len=381) with Chain E from centered_model:
@>      Failed to match chains (seqid=32%, overlap=49%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain B from ref_model (len=381) with Chain E from centered_model:
@>      Failed to match chains (seqid=7%, overlap=14%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain E from ref_model (len=709) with Chain E from centered_model:
@>      Mapped: 709 residues match with 100% sequence identity and 100% overlap.
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain F from ref_model (len=21) with Chain E from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain F from ref_model (len=21) with Chain E from centered_model:
@>      Failed to match chains (seqid=50%, overlap=2%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain F from ref_model (len=21) with Chain E from centered_model:
@>      Failed to match chains (seqid=0%, overlap=2%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain A from ref_model (len=1132) with Chain F from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain A from ref_model (len=1132) with Chain F from centered_model:
@>      Failed to match chains (seqid=38%, overlap=6%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain A from ref_model (len=1132) with Chain F from centered_model:
@>      Failed to match chains (seqid=7%, overlap=5%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain B from ref_model (len=381) with Chain F from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain B from ref_model (len=381) with Chain F from centered_model:
@>      Failed to match chains (seqid=46%, overlap=21%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain B from ref_model (len=381) with Chain F from centered_model:
@>      Failed to match chains (seqid=4%, overlap=13%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain E from ref_model (len=709) with Chain F from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain E from ref_model (len=709) with Chain F from centered_model:
@>      Failed to match chains (seqid=30%, overlap=11%).
@> Trying to map atoms based on CEalign:
@>   Comparing Chain E from ref_model (len=709) with Chain F from centered_model:
@>      Failed to match chains (seqid=4%, overlap=8%).
@> Trying to map atoms based on residue numbers and identities:
@>   Comparing Chain F from ref_model (len=21) with Chain F from centered_model:
@> Trying to map atoms based on local sequence alignment:
@>   Comparing Chain F from ref_model (len=21) with Chain F from centered_model:
@>      Mapped: 20 residues match with 100% sequence identity and 24% overlap.
@> Finding the atommaps based on their coverages...
@> Identified that there exists 1 atommap(s) potentially.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 1
----> 1 prody_ref_amap = prody.alignChains(prody_ref_model, prody_model, overlap=15, matchFunc=prody.sameChid)

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/compare.py:1748, in alignChains(atoms, target, match_func, **kwargs)
   1745     LOGGER.warn('%s has fewer chains than %s'%(atoms.getTitle(), target.getTitle()))
   1746     return []
-> 1748 atommaps = combineAtomMaps(mappings, target, **kwargs)
   1750 return atommaps

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/compare.py:1676, in combineAtomMaps(mappings, target, **kwargs)
   1674 # optimize atommaps based on superposition if target is given
   1675 if target is not None and len(nodes):
-> 1676     atommaps = _optimize(atommaps)
   1677     i = 2
   1679     if len(atommaps) < least_n_atommaps:

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/compare.py:1588, in combineAtomMaps.<locals>._optimize(atommaps)
   1585 def _optimize(atommaps):
   1586     # extract nonoverlaping mappings
   1587     if len(atommaps):
-> 1588         atommaps, rmsds = rankAtomMaps(atommaps, target)
   1590         if rmsds is not None:
   1591             debug['rmsd'] = list(rmsds)

File ~/scipion3/software/em/prody-github/ProDy/prody/proteins/compare.py:1723, in rankAtomMaps(atommaps, target)
   1721 weights = atommap.getFlags('mapped')
   1722 coords = atommap.getCoords()
-> 1723 rcoords, t = superpose(coords, coords0, weights)
   1724 rmsd = calcRMSD(rcoords, coords0, weights)
   1726 rmsds.append(rmsd)

File ~/scipion3/software/em/prody-github/ProDy/prody/measure/transform.py:213, in superpose(mobile, target, weights)
    209 def superpose(mobile, target, weights=None):
    210     """Returns *mobile*, after its RMSD minimizing superposition onto *target*,
    211     and the transformation that minimizes the RMSD."""
--> 213     t = calcTransformation(mobile, target, weights)
    214     result = applyTransformation(t, mobile)
    215     return (result, t)

File ~/scipion3/software/em/prody-github/ProDy/prody/measure/transform.py:117, in calcTransformation(mobile, target, weights)
    114     tar = target
    116 if mob.shape != tar.shape:
--> 117     raise ValueError('reference and target coordinate arrays '
    118                      'must have same number of atoms')
    120 if mob.shape[1] != 3:
    121     raise ValueError('reference and target must be coordinate arrays')

ValueError: reference and target coordinate arrays must have same number of atoms

[jamesmkrieger]

It seems like the problem is that the weights have the wrong shape because the last mapping has the wrong shape. We're missing 62 dummy atoms for the last chain to get it to length 82, so the whole length goes to 2280

In [39]: coords.shape
Out[39]: (2218, 3)

In [40]: coords0.shape
Out[40]: (2280, 3)

In [41]: weights.shape
Out[41]: (2218,)

In [42]: list(prody_model.getHierView())
Out[42]: 
[<Chain: A from centered_model (1140 residues, 1140 atoms)>,
 <Chain: B from centered_model (349 residues, 349 atoms)>,
 <Chain: E from centered_model (709 residues, 709 atoms)>,
 <Chain: F from centered_model (82 residues, 82 atoms)>]

In [43]: len(prody_ref_model.getHierView())
Out[43]: 4

In [44]: list(prody_ref_model.getHierView())
Out[44]: 
[<Chain: A from Segment A from ref_model (1132 residues, 1132 atoms)>,
 <Chain: B from Segment B from ref_model (381 residues, 381 atoms)>,
 <Chain: E from Segment E from ref_model (709 residues, 709 atoms)>,
 <Chain: F from Segment F from ref_model (21 residues, 21 atoms)>]

In [45]: prody_model
Out[45]: <AtomGroup: centered_model (2280 atoms)>

In [46]: prody_ref_model
Out[46]: <AtomGroup: ref_model (2243 atoms)>

In [47]: atommap
Out[47]: <AtomMap: (Chain F from ref_model -> Chain F from centered_model) + (Chain E from ref_model -> Chain E from centered_model) + (Chain B from ref_model -> Chain B from centered_model) + (Chain A from ref_model -> Chain A from centered_model) from ref_model (2218 atoms, 2207 mapped, 11 dummy)>

In [48]: mappings
Out[48]: 
array([[(<AtomMap: Chain A from ref_model -> Chain A from centered_model from ref_model (1140 atoms, 1132 mapped, 8 dummy)>, <AtomMap: Chain A from centered_model -> Chain A from ref_model from centered_model (1140 atoms)>, 100.0, 99.29824561403508),
        None, None, None],
       [None,
        (<AtomMap: Chain B from ref_model -> Chain B from centered_model from ref_model (349 atoms, 346 mapped, 3 dummy)>, <AtomMap: Chain B from centered_model -> Chain B from ref_model from centered_model (349 atoms)>, 100.0, 99.14040114613181),
        None, None],
       [None, None,
        (<AtomMap: Chain E from ref_model -> Chain E from centered_model from ref_model (709 atoms)>, <AtomMap: Chain E from centered_model -> Chain E from ref_model from centered_model (709 atoms)>, 100.0, 100.0),
        None],
       [None, None, None,
        (<AtomMap: Chain F from ref_model -> Chain F from centered_model from ref_model (20 atoms)>, <AtomMap: Chain F from centered_model -> Chain F from ref_model from centered_model (20 atoms)>, 100.0, 24.390243902439025)]],
      dtype=object)

In [49]: 2218+62
Out[49]: 2280

[jamesmkrieger]

@copilot fix this to add enough dummy atoms after breaking the loop

[jamesmkrieger]

@AnthonyBogetti, I don't think this actually fixes the problem at all. The problem is that Biopython pairwise local alignment no longer includes gaps.

Biopython 1.79:

In [1]: import prody

In [2]: import Bio

In [3]: Bio.__version__
Out[3]: '1.79'

In [4]: prody_ref_model = prody.parsePDB("ref_model.pdb")
@> 2243 atoms and 1 coordinate set(s) were parsed in 0.02s.

In [5]: prody_model = prody.parsePDB("centered_model.pdb")
@> 2280 atoms and 1 coordinate set(s) were parsed in 0.01s.

In [6]: target = list(prody_model.getHierView())[3]

In [7]: chain = list(prody_ref_model.getHierView())[3]

In [8]: from prody.utilities import MATCH_SCORE, MISMATCH_SCORE, GAP_PENALTY, GAP_EXT_PENALTY, ALIGNMENT_METHOD, alignBioPairwise

In [9]:             alignments = alignBioPairwise(target.getSequence(),
   ...:                                           chain.getSequence()[3:],
   ...:                                           "local",
   ...:                                           MATCH_SCORE, MISMATCH_SCORE,
   ...:                                           GAP_PENALTY,  GAP_EXT_PENALTY,
   ...:                                           max_alignments=1)

In [10]: alignments
Out[10]: 
[('KKRFEVKKWNAVALWAWDIVNCAICRNHIMDLCIECQANQASATSEECTVAWGVCNHAFHFHCISRWLKTRQVCPLDNREWE',
  '---FEVKKWNAVALWAWDIVV-------------------------------------------------------------',
  17.0,
  3,
  19)]

Biopython 1.80

In [1]: from prody.utilities import MATCH_SCORE, MISMATCH_SCORE, GAP_PENALTY, GAP_EXT_PENALTY, ALIGNMENT_METHOD, alignBioPairwise

In [2]: import prody

In [3]: prody_model = prody.parsePDB("centered_model.pdb")
@> 2280 atoms and 1 coordinate set(s) were parsed in 0.02s.

In [4]: prody_ref_model = prody.parsePDB("ref_model.pdb")
@> 2243 atoms and 1 coordinate set(s) were parsed in 0.02s.

In [5]: target = list(prody_model.getHierView())[3]

In [6]: chain = list(prody_ref_model.getHierView())[3]

In [7]:             alignments = alignBioPairwise(target.getSequence(),
   ...:                                           chain.getSequence()[3:],
   ...:                                           "local",
   ...:                                           MATCH_SCORE, MISMATCH_SCORE,
   ...:                                           GAP_PENALTY,  GAP_EXT_PENALTY,
   ...:                                           max_alignments=1)

In [8]: alignments
Out[8]: [('FEVKKWNAVALWAWDIV', 'FEVKKWNAVALWAWDIV', 17.0, '3', '20')]

I will make a fix to alignBioPairwise, I think