Add chiangles (#67)

This computes the standard `chi` angles for standard residues.
The table was taken from http://www.mlb.co.jp/linux/science/garlic/doc/commands/dihedrals.html
Cases where there are 180° symmetries in atom-naming are handled by ensuring that the
corresponding angle is positive.

Author: timholy

docs/src/documentation.md

@@ -411,9 +411,11 @@ Various functions are provided to calculate spatial quantities for proteins:
 | [`omegaangle`](@ref)         | Omega dihedral angle between a residue and the previous residue                                  |
 | [`phiangle`](@ref)           | Phi dihedral angle between a residue and the previous residue                                    |
 | [`psiangle`](@ref)           | Psi dihedral angle between a residue and the next residue                                        |
+| [`chiangle`](@ref)           | Chi dihedral angle within a residue                                                              |
 | [`omegaangles`](@ref)        | `Vector` of omega dihedral angles of an element                                                  |
 | [`phiangles`](@ref)          | `Vector` of phi dihedral angles of an element                                                    |
 | [`psiangles`](@ref)          | `Vector` of psi dihedral angles of an element                                                    |
+| [`chiangles`](@ref)          | All chi dihedral angles for a residue or element                                                 |
 | [`ramachandranangles`](@ref) | `Vector`s of phi and psi angles of an element                                                    |
 | [`ContactMap`](@ref)         | `ContactMap` of two elements, or one element with itself                                         |
 | [`DistanceMap`](@ref)        | `DistanceMap` of two elements, or one element with itself                                        |
@@ -446,6 +448,17 @@ julia> rad2deg(psiangle(struc['A'][50], struc['A'][51]))
 
 julia> rad2deg(psiangle(struc['A'], 50))
 -177.38288114072924
+
+julia> rad2deg.(chiangles(struc['A'][2]))    # ASP, χ₁...χ₅
+5-element Vector{Float64}:
+  -55.5708140556466
+ -118.03264320054458
+   56.50403552503829
+ -176.37208357380604
+    0.04991054795811607
+
+julia> rad2deg.(chiangles(struc['A'][4]))    # GLY, no χ angles
+Float64[]
 ```
 
 [`ContactMap`](@ref) takes in a structural element or a list, such as a `Chain` or `Vector{Atom}`, and returns a [`ContactMap`](@ref) object showing the contacts between the elements for a specified distance.

src/spatial.jl

@@ -16,6 +16,8 @@ export
     omegaangles,
     phiangles,
     psiangles,
+    chiangle,
+    chiangles,
     ramachandranangles,
     SpatialMap,
     ContactMap,
@@ -341,7 +343,7 @@ function dihedralangle(vec_a::AbstractVector{<:Real},
                     vec_b::AbstractVector{<:Real},
                     vec_c::AbstractVector{<:Real})
     return atan(
-        dot(cross(cross(vec_a, vec_b), cross(vec_b, vec_c)), vec_b / norm(vec_b)),
+        norm(vec_b) * dot(vec_a, cross(vec_b, vec_c)),
         dot(cross(vec_a, vec_b), cross(vec_b, vec_c)))
 end
 
@@ -462,6 +464,126 @@ function psiangle(chain::Chain, res_id::Union{Integer, AbstractString})
     return psiangle(chain[resids(chain)[i]], chain[resids(chain)[i + 1]])
 end
 
+# source: http://www.mlb.co.jp/linux/science/garlic/doc/commands/dihedrals.html
+const chitables = [
+    Dict{String,NTuple{4,String}}(
+        "ARG" => ("N", "CA", "CB", "CG"),
+        "ASN" => ("N", "CA", "CB", "CG"),
+        "ASP" => ("N", "CA", "CB", "CG"),
+        "CYS" => ("N", "CA", "CB", "SG"),
+        "GLN" => ("N", "CA", "CB", "CG"),
+        "GLU" => ("N", "CA", "CB", "CG"),
+        "HIS" => ("N", "CA", "CB", "CG"),
+        "ILE" => ("N", "CA", "CB", "CG1"),
+        "LEU" => ("N", "CA", "CB", "CG"),
+        "LYS" => ("N", "CA", "CB", "CG"),
+        "MET" => ("N", "CA", "CB", "CG"),
+        "PHE" => ("N", "CA", "CB", "CG"),
+        "PRO" => ("N", "CA", "CB", "CG"),
+        "SER" => ("N", "CA", "CB", "OG"),
+        "THR" => ("N", "CA", "CB", "OG1"),
+        "TRP" => ("N", "CA", "CB", "CG"),
+        "TYR" => ("N", "CA", "CB", "CG"),
+        "VAL" => ("N", "CA", "CB", "CG1"),
+    ),
+    Dict{String,NTuple{4,String}}(
+        "ARG" => ("CA", "CB", "CG", "CD"),
+        "ASN" => ("CA", "CB", "CG", "OD1"),
+        "ASP" => ("CA", "CB", "CG", "OD1"),
+        "GLN" => ("CA", "CB", "CG", "CD"),
+        "GLU" => ("CA", "CB", "CG", "CD"),
+        "HIS" => ("CA", "CB", "CG", "ND1"),
+        "ILE" => ("CA", "CB", "CG1", "CD1"),
+        "LEU" => ("CA", "CB", "CG", "CD1"),
+        "LYS" => ("CA", "CB", "CG", "CD"),
+        "MET" => ("CA", "CB", "CG", "SD"),
+        "PHE" => ("CA", "CB", "CG", "CD1"),
+        "PRO" => ("CA", "CB", "CG", "CD"),
+        "TRP" => ("CA", "CB", "CG", "CD1"),
+        "TYR" => ("CA", "CB", "CG", "CD1"),
+    ),
+    Dict{String,NTuple{4,String}}(
+        "ARG" => ("CB", "CG", "CD", "NE"),
+        "GLN" => ("CB", "CG", "CD", "OE1"),
+        "GLU" => ("CB", "CG", "CD", "OE1"),
+        "LYS" => ("CB", "CG", "CD", "CE"),
+        "MET" => ("CB", "CG", "SD", "CE"),
+    ),
+    Dict{String,NTuple{4,String}}(
+        "ARG" => ("CG", "CD", "NE", "CZ"),
+        "LYS" => ("CG", "CD", "CE", "NZ"),
+    ),
+    Dict{String,NTuple{4,String}}(
+        "ARG" => ("CD", "NE", "CZ", "NH1"),
+    ),
+]
+
+# source: Jumper et al 2021, Supp. Table 2
+const chisymmetries = [
+    ("ASP", 2),
+    ("GLU", 3),
+    ("PHE", 2),
+    ("TYR", 2),
+]
+
+function chiangle(a1, a2, a3, a4, sym::Bool)
+    χ = dihedralangle(a1, a2, a3, a4)
+    if sym && χ < 0
+        χ += π
+    end
+    return χ
+end
+
+"""
+    chiangle(res, i)
+
+Calculate the χᵢ angle in radians for a standard `AbstractResidue` with standard atom names.
+The angle is in the range -π to π.
+"""
+function chiangle(res::AbstractResidue, i::Integer)
+    1 <= i <= 5 || throw(ArgumentError("χᵢ index `i` must be between 1 and 5"))
+    ct = chitables[i]
+    rname = resname(res)
+    t = get(ct, rname, nothing)
+    if t === nothing
+        # is it because this isn't a standard residue?
+        (rname == "GLY" || rname == "ALA") && throw(ArgumentError("$rname does not have any χ angles"))
+        haskey(chitables[1], rname) || throw(ArgumentError("no χ angles are defined for residues with name $rname"))
+        # it must be missing specifically for `i`
+        j = 1
+        while haskey(chitables[j+1], rname)
+            j += 1
+        end
+        throw(ArgumentError("χ angle with index $i does not exist for residue $rname (max is $j)"))
+    end
+    return chiangle(res[t[1]], res[t[2]], res[t[3]], res[t[4]], (rname, i) in chisymmetries)
+end
+
+"""
+    chiangles(res)
+
+Calculate the `Vector` of standard χ angles for a standard `AbstractResidue` with standard atom names.
+The length of the vector ranges from 0 (GLY, ALA) to 5 (ARG).
+The angles are each in the range -π to π.
+"""
+function chiangles(res::AbstractResidue)
+    chi = Float64[]
+    rname = resname(res)
+    for i = 1:5
+        ct = chitables[i]
+        t = get(ct, rname, nothing)
+        if t === nothing
+            if i == 1
+                (rname == "GLY" || rname == "ALA") && break
+                throw(ArgumentError("no χ angles are defined for residues with name $rname"))
+            end
+            break
+        end
+        push!(chi, chiangle(res[t[1]], res[t[2]], res[t[3]], res[t[4]], (rname, i) in chisymmetries))
+    end
+    return chi
+end
+
 """
     omegaangles(element, residue_selectors...)
 
@@ -598,6 +720,17 @@ function ramachandranangles(el::StructuralElementOrList,
     return phiangles(el, residue_selectors...), psiangles(el, residue_selectors...)
 end
 
+"""
+    chiangles(element, residue_selectors...)
+
+Calculate the `Vector` of standard χ angles for each residue in a `StructuralElementOrList`.
+This returns a `Vector` of `Vector`s, where all angles are in the range -π to π.
+Additional arguments are residue selector functions - only residues that return
+`true` from the functions are retained.
+"""
+chiangles(el::StructuralElementOrList, residue_selectors::Function...) =
+    [chiangles(r) for r in collectresidues(el, residue_selectors...)]
+
 "A map of a structural property, e.g. a `ContactMap` or a `DistanceMap`."
 abstract type SpatialMap end
 

test/runtests.jl

@@ -3191,6 +3191,59 @@ end
     @test isapprox(phis[10], phiangle(struc_1AKE['A'], 10), atol=1e-5)
     @test isapprox(omegas[10], omegaangle(struc_1AKE['A'], 10), atol=1e-5)
 
+    # Test that the entries in `chitables` are bonded
+    sortt((a, b)) = a < b ? (a, b) : (b, a)
+    rd = BioStructures.residuedata
+    for ct in BioStructures.chitables
+        for (rname, alist) in ct
+            if rname == "HIS"
+                rname = "HID"
+            end
+            sb = sortt.(rd[rname].bonds)
+            for i = 1:3
+                @test sortt((alist[i], alist[i+1])) ∈ sb
+            end
+        end
+    end
+    chis = chiangles(struc_1AKE['A'], standardselector)
+    @test length(chis) == countresidues(struc_1AKE['A'], standardselector)
+    @test chis[1] ≈ deg2rad.([-176.231, 172.056, 56.069]) rtol=1e-5   # MET
+    @test chis[2] ≈ deg2rad.([-76.551, 171.696, 171.162, -175.969, 0.424]) rtol=1e-5   # ARG
+    @test isempty(chis[7])       # GLY
+    @test chiangle(struc_1AKE['A'][2], 3) == chis[2][3]
+    @test_throws "GLY does not have any χ angles" chiangle(struc_1AKE['A'][7], 1)
+    @test_throws "χ angle with index 4 does not exist for residue MET (max is 3)" chiangle(struc_1AKE['A'][1], 4)
+    fakeres = Residue("UKN", 10, ' ', false, struc_1AKE['A'])
+    @test_throws "no χ angles are defined for residues with name UKN" chiangle(fakeres, 1)
+    @test_throws "no χ angles are defined for residues with name UKN" chiangles(fakeres)
+    # Test symmetries: two ASPs with OD1 and OD2 swapped
+    io = IOBuffer("""
+        ATOM    534  N   ASP A   1      -8.068   7.150  -2.008  1.00 95.46           N
+        ATOM    535  CA  ASP A   1      -8.464   6.572  -3.297  1.00 95.46           C
+        ATOM    536  C   ASP A   1      -8.522   7.636  -4.409  1.00 95.46           C
+        ATOM    537  CB  ASP A   1      -9.844   5.896  -3.142  1.00 95.46           C
+        ATOM    538  O   ASP A   1      -8.015   7.431  -5.518  1.00 95.46           O
+        ATOM    539  CG  ASP A   1      -9.818   4.542  -2.418  1.00 95.46           C
+        ATOM    540  OD1 ASP A   1      -8.738   3.930  -2.377  1.00 95.46           O
+        ATOM    541  OD2 ASP A   1     -10.899   4.073  -1.981  1.00 95.46           O
+        """)
+    r = only(only(only(read(io, PDBFormat))))
+    χs1 = chiangles(r)
+    io = IOBuffer("""
+        ATOM    534  N   ASP A   1      -8.068   7.150  -2.008  1.00 95.46           N
+        ATOM    535  CA  ASP A   1      -8.464   6.572  -3.297  1.00 95.46           C
+        ATOM    536  C   ASP A   1      -8.522   7.636  -4.409  1.00 95.46           C
+        ATOM    537  CB  ASP A   1      -9.844   5.896  -3.142  1.00 95.46           C
+        ATOM    538  O   ASP A   1      -8.015   7.431  -5.518  1.00 95.46           O
+        ATOM    539  CG  ASP A   1      -9.818   4.542  -2.418  1.00 95.46           C
+        ATOM    540  OD1 ASP A   1     -10.899   4.073  -1.981  1.00 95.46           O
+        ATOM    541  OD2 ASP A   1      -8.738   3.930  -2.377  1.00 95.46           O
+        """)
+    r = only(only(only(read(io, PDBFormat))))
+    χs2 = chiangles(r)
+    @test χs1[1] ≈ χs2[1]
+    @test χs1[2] ≈ χs2[2] rtol=0.05  # the bonds are not exactly symmetric
+
     # Test ContactMap
     cas = collectatoms(struc_1AKE, calphaselector)[1:10]
     @test isa(ContactMap(cas, 10).data, BitArray{2})