Error loading `*.pdb` files generated by OpenBabel or `AtomsBase.jl`-compatible packages

[hellmrf]

I'm getting an error while reading basically any file generated by OpenBabel.

ERROR: Two copies of the same atom have the same alternative location ID.

The same error is also thrown for some *.pdb files generated by Chemfiles.jl, PDBTools.jl or AtomsIO.jl (AtomsBase.jl-compatible).

Open Babel

Let's use obabel to generate some files:

obabel -:"CC(=O)OC1=CC=CC=C1C(O)=O" -O aspirin.pdb --gen3d
obabel -:"CC(=O)OC1=CC=CC=C1C(O)=O" -O aspirin.mmcif --gen3d
obabel -:"CC(=O)NC1=CC=C(C=C1)O" -O paracetamol.pdb --gen3d
obabel *.pdb -O multi.pdb

And now we have aspirin.pdb, aspirin.mmcif, paracetamol.pdb, and multi.pdb (or download.zip).

Chemfiles.jl

For example, with Chemfiles.jl, we can load them:

julia> using Chemfiles

julia> trajectory = Trajectory("multi.pdb")
Trajectory(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_TRAJECTORY}(Ptr{Chemfiles.lib.CHFL_TRAJECTORY} @0x0000000009294110, false), nothing)

julia> trajectory |> length |> Int
2

julia> frame = read(trajectory)
Frame(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_FRAME}(Ptr{Chemfiles.lib.CHFL_FRAME} @0x0000000002417e60, false))

julia> positions(frame)
3×21 Chemfiles.ChemfilesArray:
 1.071  2.533  2.977   3.273   4.679  5.417  6.803   7.457   6.728   5.335   4.64    3.521   5.085  0.494   0.787  0.852  4.928  7.37    8.536   7.252   3.312
 0.277  0.556  1.556  -0.507  -0.5    0.658  0.652  -0.51   -1.679  -1.69   -2.983  -3.163  -3.843  1.136  -0.601  0.128  1.576  1.563  -0.507  -2.585  -2.419
 1.264  1.454  2.0     0.924   1.053  0.783  0.939   1.343   1.57    1.409   1.648   0.915   2.389  1.62    1.849  0.206  0.47   0.761   1.48    1.873   0.319

julia> f0 = read_step(trajectory, 0)
Frame(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_FRAME}(Ptr{Chemfiles.lib.CHFL_FRAME} @0x0000000009cde290, false))

julia> f1 = read_step(trajectory, 1)
Frame(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_FRAME}(Ptr{Chemfiles.lib.CHFL_FRAME} @0x0000000009c23c10, false))

julia> positions(f0)
3×21 Chemfiles.ChemfilesArray:
 1.071  2.533  2.977   3.273   4.679  5.417  6.803   7.457   6.728   5.335   4.64    3.521   5.085  0.494   0.787  0.852  4.928  7.37    8.536   7.252   3.312
 0.277  0.556  1.556  -0.507  -0.5    0.658  0.652  -0.51   -1.679  -1.69   -2.983  -3.163  -3.843  1.136  -0.601  0.128  1.576  1.563  -0.507  -2.585  -2.419
 1.264  1.454  2.0     0.924   1.053  0.783  0.939   1.343   1.57    1.409   1.648   0.915   2.389  1.62    1.849  0.206  0.47   0.761   1.48    1.873   0.319

julia> positions(f1)
3×20 Chemfiles.ChemfilesArray:
  1.078  2.546  2.933   3.328   4.731   5.51    6.902   7.515   6.757   5.364   8.87   0.545   0.684   0.914   2.864  5.061  7.48    7.246   4.79    9.26
 -0.099  0.135  1.177  -0.932  -1.036  -0.175  -0.304  -1.3    -2.184  -2.059  -1.44   0.855  -0.775  -0.523  -1.669  0.606  0.383  -2.962  -2.754  -0.721
  0.096  0.344  0.857  -0.067   0.04    0.816   0.848   0.099  -0.656  -0.678   0.077  0.152   0.859  -0.899  -0.582  1.422  1.457  -1.236  -1.285   0.601

julia> positions(frame) == positions(f0)
true

PDBTools.jl

With PDBTools.jl (which is AtomsBase.jl/AtomsIO.jl-compatible), everything runs smooth:

julia> using PDBTools

julia> asp = read_pdb("aspirin.pdb")
   Vector{Atom{Nothing}} with 21 atoms with fields:
   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
       1    C     UNL     X        1        1    1.071    0.277    1.264  1.00  0.00     1       X         1
       2    C     UNL     X        1        1    2.533    0.556    1.454  1.00  0.00     1       X         2
       3    O     UNL     X        1        1    2.977    1.556    2.000  1.00  0.00     1       X         3
       4    O     UNL     X        1        1    3.273   -0.507    0.924  1.00  0.00     1       X         4
       5    C     UNL     X        1        1    4.679   -0.500    1.053  1.00  0.00     1       X         5
       6    C     UNL     X        1        1    5.417    0.658    0.783  1.00  0.00     1       X         6
       7    C     UNL     X        1        1    6.803    0.652    0.939  1.00  0.00     1       X         7
       8    C     UNL     X        1        1    7.457   -0.510    1.343  1.00  0.00     1       X         8
⋮
      14    H     UNL     X        1        1    0.494    1.136    1.620  1.00  0.00     1       X        14
      15    H     UNL     X        1        1    0.787   -0.601    1.849  1.00  0.00     1       X        15
      16    H     UNL     X        1        1    0.852    0.128    0.206  1.00  0.00     1       X        16
      17    H     UNL     X        1        1    4.928    1.576    0.470  1.00  0.00     1       X        17
      18    H     UNL     X        1        1    7.370    1.563    0.761  1.00  0.00     1       X        18
      19    H     UNL     X        1        1    8.536   -0.507    1.480  1.00  0.00     1       X        19
      20    H     UNL     X        1        1    7.252   -2.585    1.873  1.00  0.00     1       X        20
      21    H     UNL     X        1        1    3.312   -2.419    0.319  1.00  0.00     1       X        21

julia> traj = read_pdb("multi.pdb")
   Vector{Atom{Nothing}} with 41 atoms with fields:
   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
       1    C     UNL     X        1        1    1.071    0.277    1.264  1.00  0.00     1       X         1
       2    C     UNL     X        1        1    2.533    0.556    1.454  1.00  0.00     1       X         2
       3    O     UNL     X        1        1    2.977    1.556    2.000  1.00  0.00     1       X         3
       4    O     UNL     X        1        1    3.273   -0.507    0.924  1.00  0.00     1       X         4
       5    C     UNL     X        1        1    4.679   -0.500    1.053  1.00  0.00     1       X         5
       6    C     UNL     X        1        1    5.417    0.658    0.783  1.00  0.00     1       X         6
       7    C     UNL     X        1        1    6.803    0.652    0.939  1.00  0.00     1       X         7
       8    C     UNL     X        1        1    7.457   -0.510    1.343  1.00  0.00     1       X         8
⋮
      35    H     UNL     X        1        2    0.914   -0.523   -0.899  1.00  0.00     2       X        14
      36    H     UNL     X        1        2    2.864   -1.669   -0.582  1.00  0.00     2       X        15
      37    H     UNL     X        1        2    5.061    0.606    1.422  1.00  0.00     2       X        16
      38    H     UNL     X        1        2    7.480    0.383    1.457  1.00  0.00     2       X        17
      39    H     UNL     X        1        2    7.246   -2.962   -1.236  1.00  0.00     2       X        18
      40    H     UNL     X        1        2    4.790   -2.754   -1.285  1.00  0.00     2       X        19
      41    H     UNL     X        1        2    9.260   -0.721    0.601  1.00  0.00     2       X        20

julia> read_pdb("multi.pdb", "residue 1")
   Vector{Atom{Nothing}} with 21 atoms with fields:
   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
       1    C     UNL     X        1        1    1.071    0.277    1.264  1.00  0.00     1       X         1
       2    C     UNL     X        1        1    2.533    0.556    1.454  1.00  0.00     1       X         2
       3    O     UNL     X        1        1    2.977    1.556    2.000  1.00  0.00     1       X         3
       4    O     UNL     X        1        1    3.273   -0.507    0.924  1.00  0.00     1       X         4
       5    C     UNL     X        1        1    4.679   -0.500    1.053  1.00  0.00     1       X         5
       6    C     UNL     X        1        1    5.417    0.658    0.783  1.00  0.00     1       X         6
       7    C     UNL     X        1        1    6.803    0.652    0.939  1.00  0.00     1       X         7
       8    C     UNL     X        1        1    7.457   -0.510    1.343  1.00  0.00     1       X         8
       9    C     UNL     X        1        1    6.728   -1.679    1.570  1.00  0.00     1       X         9
      10    C     UNL     X        1        1    5.335   -1.690    1.409  1.00  0.00     1       X        10
      11    C     UNL     X        1        1    4.640   -2.983    1.648  1.00  0.00     1       X        11
      12    O     UNL     X        1        1    3.521   -3.163    0.915  1.00  0.00     1       X        12
      13    O     UNL     X        1        1    5.085   -3.843    2.389  1.00  0.00     1       X        13
      14    H     UNL     X        1        1    0.494    1.136    1.620  1.00  0.00     1       X        14
      15    H     UNL     X        1        1    0.787   -0.601    1.849  1.00  0.00     1       X        15
      16    H     UNL     X        1        1    0.852    0.128    0.206  1.00  0.00     1       X        16
      17    H     UNL     X        1        1    4.928    1.576    0.470  1.00  0.00     1       X        17
      18    H     UNL     X        1        1    7.370    1.563    0.761  1.00  0.00     1       X        18
      19    H     UNL     X        1        1    8.536   -0.507    1.480  1.00  0.00     1       X        19
      20    H     UNL     X        1        1    7.252   -2.585    1.873  1.00  0.00     1       X        20
      21    H     UNL     X        1        1    3.312   -2.419    0.319  1.00  0.00     1       X        21

julia> read_pdb("multi.pdb", "residue 2")
   Vector{Atom{Nothing}} with 20 atoms with fields:
   index name resname chain   resnum  residue        x        y        z occup  beta model segname index_pdb
      22    C     UNL     X        1        2    1.078   -0.099    0.096  1.00  0.00     2       X         1
      23    C     UNL     X        1        2    2.546    0.135    0.344  1.00  0.00     2       X         2
      24    O     UNL     X        1        2    2.933    1.177    0.857  1.00  0.00     2       X         3
      25    N     UNL     X        1        2    3.328   -0.932   -0.067  1.00  0.00     2       X         4
      26    C     UNL     X        1        2    4.731   -1.036    0.040  1.00  0.00     2       X         5
      27    C     UNL     X        1        2    5.510   -0.175    0.816  1.00  0.00     2       X         6
      28    C     UNL     X        1        2    6.902   -0.304    0.848  1.00  0.00     2       X         7
      29    C     UNL     X        1        2    7.515   -1.300    0.099  1.00  0.00     2       X         8
      30    C     UNL     X        1        2    6.757   -2.184   -0.656  1.00  0.00     2       X         9
      31    C     UNL     X        1        2    5.364   -2.059   -0.678  1.00  0.00     2       X        10
      32    O     UNL     X        1        2    8.870   -1.440    0.077  1.00  0.00     2       X        11
      33    H     UNL     X        1        2    0.545    0.855    0.152  1.00  0.00     2       X        12
      34    H     UNL     X        1        2    0.684   -0.775    0.859  1.00  0.00     2       X        13
      35    H     UNL     X        1        2    0.914   -0.523   -0.899  1.00  0.00     2       X        14
      36    H     UNL     X        1        2    2.864   -1.669   -0.582  1.00  0.00     2       X        15
      37    H     UNL     X        1        2    5.061    0.606    1.422  1.00  0.00     2       X        16
      38    H     UNL     X        1        2    7.480    0.383    1.457  1.00  0.00     2       X        17
      39    H     UNL     X        1        2    7.246   -2.962   -1.236  1.00  0.00     2       X        18
      40    H     UNL     X        1        2    4.790   -2.754   -1.285  1.00  0.00     2       X        19
      41    H     UNL     X        1        2    9.260   -0.721    0.601  1.00  0.00     2       X        20

BioStructures.jl

But with BioStructures.jl and PDBFormat, we get some errors:

julia> using BioStructures

julia> struc = read("aspirin.pdb", PDBFormat)
ERROR: Two copies of the same atom have the same alternative location ID. Existing atom:
Atom C with serial 1, coordinates [1.071, 0.277, 1.264]
New atom record to add:
AtomRecord C with serial 2, coordinates [2.533, 0.556, 1.454]
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] unsafe_addatomtomodel!(mo::Model, atom_rec::AtomRecord; remove_disorder::Bool)
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/model.jl:1652
  [3] unsafe_addatomtomodel!
    @ ~/.julia/packages/BioStructures/ZBun0/src/model.jl:1559 [inlined]
  [4] read(input::IOStream, ::Type{…}; structure_name::String, remove_disorder::Bool, read_std_atoms::Bool, read_het_atoms::Bool, run_dssp::Bool, run_stride::Bool)
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:66
  [5] read
    @ ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:47 [inlined]
  [6] (::BioStructures.var"#70#71"{String, @Kwargs{}, DataType})(input::IOStream)
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:113
  [7] open(f::BioStructures.var"#70#71"{String, @Kwargs{}, DataType}, args::String; kwargs::@Kwargs{})
    @ Base ./io.jl:410
  [8] open
    @ ./io.jl:407 [inlined]
  [9] #read#69
    @ ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:112 [inlined]
 [10] read(filepath::String, t::Type{PDBFormat})
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:108
 [11] top-level scope
    @ REPL[7]:1
Some type information was truncated. Use `show(err)` to see complete types.

(Same behavior for paracetamol.pdb or multi.pdb).

Same for MMCIFDict (now with aspirin.mmcif):

julia> struc = "aspirin.mmcif" |> MMCIFDict |> MolecularStructure
ERROR: KeyError: key "_atom_site.group_PDB" not found
Stacktrace:
 [1] getindex
   @ ./dict.jl:477 [inlined]
 [2] getindex
   @ ~/.julia/packages/BioStructures/ZBun0/src/mmcif.jl:84 [inlined]
 [3]
   @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/mmcif.jl:332
 [4] MolecularStructure
   @ ~/.julia/packages/BioStructures/ZBun0/src/mmcif.jl:321 [inlined]
 [5] |>(x::MMCIFDict, f::Type{MolecularStructure})
   @ Base ./operators.jl:926
 [6] top-level scope
   @ REPL[9]:1
Some type information was truncated. Use `show(err)` to see complete types.

Non-Open Babel

Maybe the problem is with Open Babel, right? So let's download Aspirin's SDF from PubChem and use Chemfiles.jl to convert it to PDB.

julia> using Downloads; Downloads.download("https://pubchem.ncbi.nlm.nih.gov/rest/pug/conformers/000008C400000001/SDF?response_type=save&response_basename=2244", "./2244.sdf")
"./2244.sdf"

julia> using Chemfiles

julia> traj = Trajectory("2244.sdf")
Trajectory(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_TRAJECTORY}(Ptr{Chemfiles.lib.CHFL_TRAJECTORY} @0x0000000002315a30, false), nothing)

julia> frame = read(traj)
Frame(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_FRAME}(Ptr{Chemfiles.lib.CHFL_FRAME} @0x00000000030c9fd0, false))

julia> Trajectory("chfl.pdb", 'w') do traj
           write(traj, frame)
       end

julia> using BioStructures

julia> struc = read("chfl.pdb", PDBFormat)
ERROR: PDBParseError: could not read model serial number at line 1 of file:
MODEL    1
Stacktrace:
 [1] read(input::IOStream, ::Type{…}; structure_name::String, remove_disorder::Bool, read_std_atoms::Bool, read_het_atoms::Bool, run_dssp::Bool, run_stride::Bool)
   @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:74
 [2] read
   @ ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:47 [inlined]
 [3] (::BioStructures.var"#70#71"{String, @Kwargs{}, DataType})(input::IOStream)
   @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:113
 [4] open(f::BioStructures.var"#70#71"{String, @Kwargs{}, DataType}, args::String; kwargs::@Kwargs{})
   @ Base ./io.jl:410
 [5] open
   @ ./io.jl:407 [inlined]
 [6] #read#69
   @ ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:112 [inlined]
 [7] read(filepath::String, t::Type{PDBFormat})
   @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:108
 [8] top-level scope
   @ REPL[9]:1

caused by: ArgumentError: input string is empty or only contains whitespace
Stacktrace:
  [1] tryparse_internal(::Type{Int64}, s::String, startpos::Int64, endpos::Int64, base_::Int64, raise::Bool)
    @ Base ./parse.jl:115
  [2] #parse#619
    @ ./parse.jl:254 [inlined]
  [3] parse
    @ ./parse.jl:253 [inlined]
  [4] read(input::IOStream, ::Type{…}; structure_name::String, remove_disorder::Bool, read_std_atoms::Bool, read_het_atoms::Bool, run_dssp::Bool, run_stride::Bool)
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:72
  [5] read
    @ ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:47 [inlined]
  [6] (::BioStructures.var"#70#71"{String, @Kwargs{}, DataType})(input::IOStream)
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:113
  [7] open(f::BioStructures.var"#70#71"{String, @Kwargs{}, DataType}, args::String; kwargs::@Kwargs{})
    @ Base ./io.jl:410
  [8] open
    @ ./io.jl:407 [inlined]
  [9] #read#69
    @ ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:112 [inlined]
 [10] read(filepath::String, t::Type{PDBFormat})
    @ BioStructures ~/.julia/packages/BioStructures/ZBun0/src/pdb.jl:108
 [11] top-level scope
    @ REPL[9]:1
Some type information was truncated. Use `show(err)` to see complete types.

Finals

Also, this seem to be related to #48.

download.zip

[jgreener64]

This is indeed related to #48. The issue is lines like this:

HETATM    1  C   UNL     1       1.071   0.277   1.264  1.00  0.00           C  
HETATM    2  C   UNL     1       2.533   0.556   1.454  1.00  0.00           C  

If the file is being read as a series of atoms then this isn't a problem. If a molecule is being constructed from the contents of the file, then this creates ambiguity. You either end up with two different atoms and no way to distinguish them, or you have to make a choice about which atom to retain. BioStructures.jl tries to construct the underlying molecule, hence the error.

Whilst the PDB format is somewhat underspecified, the PDB itself doesn't have duplicate lines like the above (except when distinguished by alternative location IDs).

The chfl.pdb file written by Chemfiles has a different problem, with lines like:

MODEL    1

This puts the model number in column 10 when it should be in columns 11-14 (https://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#MODEL), causing the error seen.

HETATM    1 O            1       1.233   0.554   0.779  1.00  0.00           O

This lacks a residue name and chain ID. Again it's hard to know what is required by the format, but the PDB always has residue names with placeholders like UNK possible. https://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM mentions that the chain ID should be non-blank, though BioStructures does accept blank chain IDs. You could raise these issues at Chemfiles.

I'll add some of this discussion to the docs.

[jgreener64]

I expanded the docs section at https://biojulia.dev/BioStructures.jl/dev/documentation/#Design-decisions.