Add bonding data for residues (#66)

This provides deterministic bonding between atoms in peptides that
follow standard naming conventions. The bonding data is taken from
OpenMM's ff14SB force field.

To encourage careful work, this is designed to error unless all atoms
(including hydrogens) are present in the structure, and that HIS has
been disambiguated as HIE, HID, or HIP. One can override this behavior
by passing `strict=false` as a keyword argument to `MetaGraph`.

Author: timholy

.gitignore

@@ -8,3 +8,4 @@ benchmark/tune.json
 benchmark/results
 Manifest.toml
 *.swp
+extractdata/protein.ff14SB.*

ext/BioStructuresGraphsExt.jl

@@ -3,11 +3,12 @@ module BioStructuresGraphsExt
 using BioStructures
 using Graphs
 using MetaGraphs
+using BioStructures: findatombyname
 
 """
     MetaGraph(element, contact_distance)
 
-Construct a graph of atoms where edges are contacts.
+Construct a graph of atoms where edges are contacts separated by less than `contact_distance`.
 
 See Graphs.jl and MetaGraphs.jl for more on how to use graphs.
 """
@@ -28,4 +29,92 @@ function MetaGraphs.MetaGraph(el::StructuralElementOrList, contact_dist::Real)
     return mg
 end
 
+"""
+    MetaGraph(chain::Chain; strict::Bool = true)
+
+Construct a graph of atoms where edges are determined by the known bonds of residues in the chain.
+
+By default, the graph is constructed in `strict` mode, which means that:
+
+- residue and atom names must be standard
+- all hydrogens are present
+- HIS must be disambiguated as HIE, HID, or HIP
+
+These constraints can be relaxed by setting `strict = false`, at some risk to accuracy.
+
+See Graphs.jl and MetaGraphs.jl for more on how to use graphs.
+"""
+function MetaGraphs.MetaGraph(chain::Chain; strict::Bool = true)
+    el_list = collectatoms(chain; expand_disordered = true)
+    mg = MetaGraph(length(el_list))
+    for (i, el) in enumerate(el_list)
+        set_prop!(mg, i, :element, el)
+    end
+    set_indexing_prop!(mg, :element)
+
+    prev = nothing
+    for r in chain
+        ishetero(r) && continue
+        if prev !== nothing
+            if resnumber(r) == resnumber(prev) + 1
+                # Add the peptide bond(s) (disordered residues may need multiples)
+                for _r in collectresidues(r; expand_disordered=true), _rp in collectresidues(prev; expand_disordered=true)
+                    for a in _r["N"], ap in _rp["C"]
+                        add_edge!(mg, mg[a, :element], mg[ap, :element])
+                    end
+                end
+            else
+                prev = nothing
+            end
+        end
+        # Add the residue bonds
+        addresiduebonds!(mg, r, strict)
+        prev = r
+        # The "OXT" (C-terminus oxygen) is not connected
+        for _r in collectresidues(r; expand_disordered=true)
+            aj = findatombyname(_r, "OXT"; strict=false)
+            if aj !== nothing
+                ai = findatombyname(_r, "C"; strict)
+                connect_atoms!(mg, ai, aj)
+            end
+        end
+    end
+
+    return mg
+end
+
+function addresiduebonds!(mg, r::Residue, strict::Bool)
+    rname = residuekey(r, strict)
+    strict || haskey(BioStructures.residuedata, rname) || return
+    rd = BioStructures.residuedata[rname]
+    for (ni, nj) in rd.bonds
+        connect_atoms!(mg, r, ni, nj, strict)
+    end
+end
+
+connect_atoms!(mg, r, ni::AbstractString, nj::AbstractString, strict::Bool) =
+    connect_atoms!(mg, findatombyname(r, ni; strict), findatombyname(r, nj; strict))
+
+function connect_atoms!(mg, ai::Union{AbstractAtom,Nothing}, aj::Union{AbstractAtom,Nothing})
+    (ai === nothing || aj === nothing) && return
+    for _ai in ai, _aj in aj  # handle disordered atoms
+        i = mg[_ai, :element]
+        j = mg[_aj, :element]
+        add_edge!(mg, i, j)
+    end
+end
+
+addresiduebonds!(mg, dr::DisorderedResidue, strict::Bool) = for r in values(dr.names)
+    addresiduebonds!(mg, r, strict)
+end
+
+function residuekey(r::Residue, strict::Bool)
+    rname = resname(r)
+    if rname == "HIS"
+        strict && error("HIS is not a valid residue name in strict mode; use HIE, HID, or HIP")
+        rname = "HIE"
+    end
+    return rname
+end
+
 end # BioStructuresGraphsExt

extractdata/bonding.jl

@@ -0,0 +1,121 @@
+# Running this file generates src/bonding.jl, which contains the atomtypes and residuedata dictionaries.
+# Thanks to OpenMM for the ff14SB force field XML file.
+
+using Downloads
+
+if !isfile(joinpath(@__DIR__, "protein.ff14SB.xml"))
+    Downloads.download("https://raw.githubusercontent.com/openmm/openmm/refs/heads/master/wrappers/python/openmm/app/data/amber14/protein.ff14SB.xml", "protein.ff14SB.xml")
+end
+
+function parsexmblock(f, io::IO, key)
+    while !eof(io)
+        line = strip(readline(io))
+        line == key && return nothing
+        f(line)
+    end
+end
+
+function parsestring(str)
+    @assert startswith(str, '"')
+    @assert endswith(str, '"')
+    return String(str[2:end-1])
+end
+
+function parsexmlline(f, line, tag, keyname)
+    @assert startswith(line, "<$tag ")
+    @assert endswith(line, "/>")
+    kv = split(strip(line[length(tag)+2:end-2]), ' ')
+    key = ""
+    vals = Pair{Symbol,Any}[]
+    for kvp in kv
+        k, v = split(kvp, '=')
+        if k == keyname
+            key = parsestring(v)
+        else
+            push!(vals, Symbol(k) => f(k, v))
+        end
+    end
+    return key => (; vals...)
+end
+
+
+atomtypes, residues = open("protein.ff14SB.xml", "r") do io
+    line = readline(io)
+    @assert line == "<ForceField>"
+    atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String}}()
+    residues = Dict{String, @NamedTuple{atoms::Dict{String, @NamedTuple{charge::Float32, type::String}}, bonds::Vector{Tuple{String,String}}, externalbonds::Vector{String}}}()
+    parsexmblock(io, "</ForceField>") do line
+        if line == "<AtomTypes>"
+            parsexmblock(io, "</AtomTypes>") do line
+                push!(atomtypes, parsexmlline(line, "Type", "class") do k, v
+                    if k == "element"
+                        return parsestring(v)
+                    elseif k == "mass"
+                        return parse(Float32, v[2:end-1])  # strip the quotes
+                    elseif k == "name"
+                        return parsestring(v)
+                    else
+                        error("Unknown AtomType key $k")
+                    end
+                end)
+            end
+        elseif line == "<Residues>"
+            parsexmblock(io, "</Residues>") do line
+                if startswith(line, "<Residue name=")
+                    resname = parsestring(line[15:end-1])
+                    atoms = Dict{String, @NamedTuple{charge::Float32, type::String}}()
+                    bonds = Vector{Tuple{String,String}}()
+                    externalbonds = Vector{String}()
+                    parsexmblock(io, "</Residue>") do line
+                        if startswith(line, "<Atom")
+                            push!(atoms, parsexmlline(line, "Atom", "name") do k, v
+                                if k == "charge"
+                                    return parse(Float32, v[2:end-1])  # strip the quotes
+                                elseif k == "type"
+                                    return parsestring(v)
+                                else
+                                    error("Unknown Atom key $k")
+                                end
+                            end)
+                        elseif startswith(line, "<Bond")
+                            line = line[6:end-2]
+                            a1, a2 = split(strip(line), ' ')
+                            push!(bonds, (only(match(r"atomName1=\"(.*)\"", a1).captures), only(match(r"atomName2=\"(.*)\"", a2).captures)))
+                        elseif startswith(line, "<ExternalBond")
+                            line = line[14:end-2]
+                            push!(externalbonds, only(match(r"atomName=\"(.*)\"", line).captures))
+                        else
+                            error("Unknown Residue line $line")
+                        end
+                    end
+                    residues[resname] = (; atoms, bonds, externalbonds)
+                else
+                    error("Unknown Residues line $line")
+                end
+            end
+        end
+    end
+    atomtypes, residues
+end
+
+open(joinpath(dirname(@__DIR__), "src", "bonding.jl"), "w") do io
+    println(io, "const atomtypes = Dict{String, @NamedTuple{element::String, mass::Float32, name::String}}(")
+    at = sort!(collect(atomtypes); by=first)
+    for pr in at
+        println(io, "    ", pr, ',')
+    end
+    println(io, ")\n")
+
+    println(io, "const RDADict = Dict{String, @NamedTuple{charge::Float32, type::String}}")
+
+    println(io, "const residuedata = Dict{String, @NamedTuple{atoms::RDADict, bonds::Vector{Tuple{String,String}}, externalbonds::Vector{String}}}(")
+    rd = sort!(collect(residues); by=first)
+    for (k, v) in rd
+        print(io, "    "*" "^(4-length(k)))
+        show(io, k)
+        println(io, " => (atoms = ", replace(sprint(show, v.atoms), "Dict{String, @NamedTuple{charge::Float32, type::String}}" => "RDADict"), ',')
+        println(io, "               bonds = ", v.bonds, ",")
+        println(io, "               externalbonds = ", v.externalbonds, "),")
+    end
+    println(io, ")")
+end

src/BioStructures.jl

@@ -26,6 +26,7 @@ include("pdb.jl")
 include("mmcif.jl")
 include("download.jl")
 include("spatial.jl")
+include("bonding.jl")
 
 @compile_workload begin
     let