Skip to content

Add chiangles #67

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
May 7, 2025
Merged

Add chiangles #67

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
72b1a2d
Add `chiangles`
timholy May 4, 2025
f134295
Handle symmetries
timholy May 4, 2025
a85ffd4
Improve error messages
timholy May 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 13 additions & 0 deletions docs/src/documentation.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 |
Expand Down Expand Up @@ -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.
Expand Down
135 changes: 134 additions & 1 deletion src/spatial.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,8 @@ export
omegaangles,
phiangles,
psiangles,
chiangle,
chiangles,
ramachandranangles,
SpatialMap,
ContactMap,
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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...)

Expand Down Expand Up @@ -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

Expand Down
53 changes: 53 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3184,6 +3184,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})
Expand Down
Loading