Non-consecutive atom index numbers break mapping

[francosimonetti]

Hi Hami and Alex,

I tried the package (thanks a lot!) and I had an issue with a PDB file with non-consecutive atom index (some had been deleted due to alternative atomic positions). This produced an invalid mapping of residue number before and after loading the pdb into mdtraj. The issue is that mdtraj re-indexes atoms by topological order and so the atom_res_map no longer holds the right atom_idx->res_idx mapping in the lines variable after re-reading the file. It is a weird behaviour. You can reproduce by just deleting one or two atoms in a PDB and breaking the consecutive atom indices.
I fixed it with a small change in the _save_clean_topology() function:

def _save_clean_topology(self, frame, output_path):
       
        pdb_res_nums = [res.resSeq for res in frame.topology.residues]
        min_residue_number = min(pdb_res_nums)
        shift = 1 - min_residue_number
        
        import tempfile
        with tempfile.NamedTemporaryFile(suffix=".pdb", delete=False) as temp_pdb:
            frame.save(temp_pdb.name)
            with open(temp_pdb.name, 'r') as f:
                lines = f.readlines()

        topology_atoms = list(frame.topology.atoms)
        atom_line_index = 0

        with open(output_path, 'w') as f:
            for line in lines:
                if line.startswith(("ATOM", "HETATM")):
                    atom = topology_atoms[atom_line_index]
                    res_idx = atom.residue.index
                    atom_line_index += 1
                    original_resnum = pdb_res_nums[res_idx]
                    corrected_resnum = original_resnum + shift
                    line = line[:22] + f"{corrected_resnum:4d}" + line[26:]
                f.write(line)

I hope you can propagate this fix into the main code and thanks for this great package!
Best regards,
Franco

[hamidrgh]

Hello Franco,

Thank you very much for your valuable feedback. I'll check it again and resolve the issue.