Specialize residue names (#80)

If you've added hydrogens to a structure, this supports the full
information present in the `residuedata` tables imported from ff19SB:
- renames HIS to one of HID/HIE/HIP (using the protonation state from
  the added hydrogens)
- renames the amino- and carboxyl-terminal residues by prepending
  `N`/`C`

Author: timholy

Project.toml

@@ -1,7 +1,7 @@
 name = "BioStructures"
 uuid = "de9282ab-8554-53be-b2d6-f6c222edabfc"
 authors = ["Joe G Greener <jgreener@hotmail.co.uk>"]
-version = "4.6.1"
+version = "4.7.0"
 
 [deps]
 BioGenerics = "47718e42-2ac5-11e9-14af-e5595289c2ea"

docs/src/documentation.md

@@ -568,7 +568,34 @@ Similar to [`ContactMap`](@ref), contacts are found between any element type pas
 So if you wanted the graph of chain contacts in a protein complex you could give a [`Model`](@ref) as the first argument.
 
 Calling `MetaGraph` on a [`Chain`](@ref) without a contact distance does something different: it constructs a graph of atoms where edges are determined by the known bonds of standard amino acids in the chain.
-Hydrogens should be present and atom names should match those in OpenMM.
+Hydrogens should be present and atom names should match those in OpenMM. See [`specialize_resnames!`](@ref) to
+ensure that residue names have been specialized to their actual configuration in the protein:
+
+```julia-repl
+julia> using BioStructures, MetaGraphs
+
+julia> struc_M3YHX5 = read(joinpath(pkgdir(BioStructures), "test", "data", "AF-M3YHX5-F1-model_v4_hydrogens.cif"), MMCIFFormat)
+MolecularStructure AF-M3YHX5-F1-model_v4_hydrogens.cif with 1 models, 1 chains (A), 112 residues, 1833 atoms
+
+julia> MetaGraph(struc_M3YHX5["A"])
+r = Residue 1:A with name MET, 19 atoms
+ERROR: KeyError: key "H" not found
+[...]
+
+julia> mg = MetaGraph(struc_M3YHX5["A"]; strict=false)
+{1833, 1850} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)
+```
+
+but if you really want to get all the bonds, you'd better use the default `strict=true`.
+Fortunately, this example only requires residue renaming:
+
+```julia-repl
+julia> specialize_resnames!(struc_M3YHX5)
+MolecularStructure AF-M3YHX5-F1-model_v4_hydrogens.cif with 1 models, 1 chains (A), 112 residues, 1833 atoms
+
+julia> mg = MetaGraph(struc_M3YHX5["A"])
+{1833, 1853} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)
+```
 
 ## Assigning secondary structure
 

src/model.jl

@@ -39,6 +39,7 @@ export
     defaultresname,
     defaultresidue,
     resnames,
+    specialize_resnames!,
     sscode,
     sscode!,
     chain,
@@ -1097,6 +1098,54 @@ This is the `Model` with the lowest model number.
 """
 defaultmodel(struc::MolecularStructure) = first(sort(collect(values(models(struc)))))
 
+"""
+    specialize_resnames!(r)
+
+If `r` is a histidine residue, rename it to reflect a specific protonation configuration (HID,
+HIE, or HIP) based on the presence of HD1 and HE2 atoms. For any other residue,
+`r` will be unchanged.
+
+Hydrogens must have been assigned prior to calling this function.
+"""
+function specialize_resnames!(r::Residue)
+    rname = resname(r)
+    if rname == "HIS"
+        hd1 = findatombyname(r, "HD1"; strict=false)
+        he2 = findatombyname(r, "HE2"; strict=false)
+        if !isnothing(hd1) && !isnothing(he2)
+            r.name = "HIP"
+        elseif !isnothing(hd1)
+            r.name = "HID"
+        elseif !isnothing(he2)
+            r.name = "HIE"
+        else
+            error("at least one of HD1 or HE2 must be present")
+        end
+    end
+    # N and C termini
+    if haskey(r.atoms, "H3") && !haskey(residuedata[rname].atoms, "H3") && haskey(residuedata, "N" * rname)
+        r.name = "N" * rname
+    elseif haskey(r.atoms, "OXT") && !haskey(residuedata[rname].atoms, "OXT") && haskey(residuedata, "C" * rname)
+        r.name = "C" * rname
+    end
+    return r
+end
+specialize_resnames!(r::DisorderedResidue) = (foreach(specialize_resnames!, collectresidues(r; expand_disordered=true)); r)
+function specialize_resnames!(c::Chain)
+    for (key, res) in c.residues
+        if isa(res, DisorderedResidue)
+            # Create a new `names` dictionary with specialized residue names
+            names = Dict((specialize_resnames!(r); resname(r, strip=false) => r) for r in values(res.names))
+            c.residues[key] = DisorderedResidue(names, resname(res.names[defaultresname(res)]))
+        else
+            specialize_resnames!(res)
+        end
+    end
+    return c
+end
+specialize_resnames!(m::Model) = (foreach(specialize_resnames!, m); m)
+specialize_resnames!(s::MolecularStructure) = (foreach(specialize_resnames!, s); s)
+
 # Sort lists of elements
 
 # Sort atoms by serial