diff --git a/src/backbone/backbone.jl b/src/backbone/backbone.jl index c31eb24..36f1120 100644 --- a/src/backbone/backbone.jl +++ b/src/backbone/backbone.jl @@ -34,5 +34,4 @@ Returns the coordinates of specific columns of atoms in a backbone. @inline atom_coords(backbone::Backbone, i) = view(backbone.coords, :, i, :) include("rotations.jl") -include("bonds.jl") -include("dihedrals.jl") \ No newline at end of file +include("bonds.jl") \ No newline at end of file diff --git a/src/backbone/bonds.jl b/src/backbone/bonds.jl index 6c5014a..d20332f 100644 --- a/src/backbone/bonds.jl +++ b/src/backbone/bonds.jl @@ -1,28 +1,25 @@ -export get_atom_displacements, get_atom_distances, get_bond_vectors, get_bond_angles +export get_atom_displacements, get_atom_distances +export get_bond_vectors, get_bond_lengths +export get_bond_angles, get_dihedrals +export ChainedBonds -function get_atom_displacements(backbone::Backbone{A}, atom1::Integer, atom2::Integer, residue_offset::Integer=0) where A +column_norms(vectors::AbstractMatrix) = reshape(mapslices(norm, vectors, dims=1), :) + +function get_atom_displacements( + backbone::Backbone{A}, atom1::Integer, atom2::Integer, residue_offset::Integer = 0, +) where A @assert 1 <= atom1 <= A && 1 <= atom2 <= A "Backbone{$N} does not have atoms $atom1 and $atom2" coords = backbone.coords displacements = @view(coords[:, atom2, 1+residue_offset:end]) .- @view(coords[:, atom1, 1:end-residue_offset]) return displacements end -norms(vectors::AbstractMatrix) = reshape(mapslices(norm, vectors, dims=1), :) - -""" - get_atom_distances(backbone::Backbone, atom1::Integer, atom2::Integer, residue_offset::Integer) - -Calculate the distances between all pairs of two types atoms in a backbone, e.g. -the distances between all pairs of contiguous carbonyl and nitrogen atoms. -atom1 and atom2 are the indices of the atoms in the backbone, and residue_offset -is the number of residues between the atoms (0 by default). - -Returns a vector of distances of length (length(backbone) - residue_offset). -""" -function get_atom_distances(backbone::Backbone{A}, atom1::Integer, atom2::Integer, residue_offset::Integer=0) where A +function get_atom_distances( + backbone::Backbone{A}, atom1::Integer, atom2::Integer, residue_offset::Integer = 0, +) where A @assert 1 <= atom1 <= A && 1 <= atom2 <= A "Backbone{$N} does not have atoms $atom1 and $atom2" displacements = get_atom_displacements(backbone, atom1, atom2, residue_offset) - distances = norms(displacements) + distances = column_norms(displacements) return distances end @@ -34,15 +31,15 @@ function get_bond_vectors(backbone::Backbone{A}) where A end -function get_bond_lengths(bond_vectors::AbstractVector{<:AbstractVector{T}}) where T - bond_lengths = norms(bond_vectors) +function get_bond_lengths(bond_vectors::AbstractMatrix{T}) where T + bond_lengths = column_norms(bond_vectors) return bond_lengths end get_bond_lengths(backbone::Backbone) = get_bond_lengths(get_bond_vectors(backbone)) -get_bond_angle(v1::T, v2::T) where T <: AbstractVector = acos(dot(v1, v2) / (norm(v1) * norm(v2))) +get_bond_angle(v1::V, v2::V) where V <: AbstractVector{<:Real} = acos((-dot(v1, v2)) / (norm(v1) * norm(v2))) function get_bond_angles(bond_vectors::AbstractMatrix{T}) where T bond_vector_pairs = zip(eachcol(bond_vectors), Iterators.drop(eachcol(bond_vectors), 1)) @@ -52,28 +49,124 @@ end get_bond_angles(backbone::Backbone) = get_bond_angles(get_bond_vectors(backbone)) +# source: en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics +function dihedral_angle(u1::V, u2::V, u3::V) where V <: AbstractVector{<:Real} + c12, c23 = cross(u1, u2), cross(u2, u3) + atan(dot(u2, cross(c12, c23)), norm(u2) * dot(c12, c23)) +end + +function get_dihedrals(vectors::AbstractMatrix{T}) where T + dihedrals = Vector{T}(undef, size(vectors, 2) - 2) + u1s = eachcol(vectors) + u2s = Iterators.drop(u1s, 1) + u3s = Iterators.drop(u1s, 2) + for (i, u1, u2, u3) in zip(eachindex(dihedrals), u1s, u2s, u3s) + dihedrals[i] = dihedral_angle(u1, u2, u3) + end + return dihedrals +end + +get_dihedrals(backbone::Backbone) = get_dihedrals(get_bond_vectors(backbone)) + + +""" + ChainedBonds{T <: Real} -struct ContinuousBonds{T} <: AbstractVector{AbstractVector{T}} - vectors::AbstractMatrix{T} +A lazy way to store a backbone as a series of bond lengths, angles, and dihedrals. + +Can be instantiated from a Backbone or a matrix of bond vectors. + +Can be used to instantiate a Backbone using the `Backbone{A}(bonds::ChainedBonds)` constructor. +""" +struct ChainedBonds{T <: Real} lengths::AbstractVector{T} angles::AbstractVector{T} + dihedrals::AbstractVector{T} - function ContinuousBonds(vectors::AbstractMatrix{T}) where T + function ChainedBonds(vectors::AbstractMatrix{T}) where T lengths = get_bond_lengths(vectors) angles = get_bond_angles(vectors) + dihedrals = get_dihedrals(vectors) - L = size(vectors, 2) + 1 - @assert size(lengths) == L-1 - @assert size(angles) == L-2 + return new{T}(lengths, angles, dihedrals) + end - return new{T}(vectors, lengths, angles) + function ChainedBonds(backbone::Backbone{A, T}) where {A, T} + return ChainedBonds(get_bond_vectors(backbone)) end +end + +Base.:(==)(b1::ChainedBonds, b2::ChainedBonds) = b1.lengths == b2.lengths && b1.angles == b2.angles && b1.dihedrals == b2.dihedrals +Base.:(≈)(b1::ChainedBonds, b2::ChainedBonds) = b1.lengths ≈ b2.lengths && b1.angles ≈ b2.angles && b1.dihedrals ≈ b2.dihedrals +Base.length(bonds::ChainedBonds) = length(bonds.lengths) +Base.size(bonds::ChainedBonds) = Tuple(length(bonds)) + +function get_first_points(bonds::ChainedBonds{T}) where T + L = length(bonds) + 1 + l = min(3, L) + coords = Matrix{T}(undef, 3, l) + + coords[:, 1] = [0, 0, 0] + if l >= 2 + coords[:, 2] = [bonds.lengths[1], 0, 0] + if l == 3 + N = normalize([0, 0, 1]) + bond_angle_rotation = Rotations.AngleAxis(π - bonds.angles[1], N...) + coords[:, 3] = coords[:, 2] + bond_angle_rotation * [bonds.lengths[2], 0, 0] + end + end + + return coords +end + +# first points currently don't get adjusted to fit the bonds +function Backbone{ATOMS_PER_RESIDUE}( + bonds::ChainedBonds{T}, + first_points::AbstractMatrix{T} = get_first_points(bonds), +) where {ATOMS_PER_RESIDUE, T} + @assert (length(bonds) + 1) % ATOMS_PER_RESIDUE == 0 "Invalid number of atoms per residue in backbone" + + L = length(bonds) + 1 + coords = fill(T(NaN), 3, L) + + l = size(first_points, 2) + @assert l <= L + coords[:, 1:l] = first_points + + for i in l+1:L + A, B, C = eachcol(coords[:, i-3:i-1]) + AB, BC = B - A, C - B + n_AB, n_BC = normalize(AB), normalize(BC) + CD_init = n_BC * bonds.lengths[i-1] + + N = normalize(cross(n_AB, n_BC)) # wont work if AB == BC + bond_angle_rotation = Rotations.AngleAxis(π - bonds.angles[i-2], N...) + CD_rot1 = bond_angle_rotation * CD_init + + dihedral_rotation = Rotations.AngleAxis(bonds.dihedrals[i-3], n_BC...) + CD_rot2 = dihedral_rotation * CD_rot1 + + D = C + CD_rot2 + coords[:, i] = D + end + + backbone = Backbone(reshape(coords, 3, ATOMS_PER_RESIDUE, :)) + return backbone +end + +Backbone(bonds::ChainedBonds; kwargs...) = Backbone{1}(bonds; kwargs...) + +get_skip_length(L1, L2, θ) = sqrt(L1^2 + L2^2 - 2*L1*L2*cos(θ)) - function ContinuousBonds(backbone::Backbone{A, T}) where {A, T} - return ContinuousBonds(get_bond_vectors(backbone)) +function get_skip_lengths(lengths::V, angles::V) where V <: AbstractVector{<:Real} + @assert length(lengths) == length(angles) + 1 + L1s = lengths + L2s = Iterators.drop(lengths, 1) + skip_lengths = similar(angles) + for (L1, L2, (i, θ)) in zip(L1s, L2s, enumerate(angles)) + skip_lengths[i] = get_skip_length(L1, L2, θ) end + return skip_lengths end -Base.length(bonds::ContinuousBonds) = size(bonds.vectors, 2) -Base.size(bonds::ContinuousBonds) = Tuple(length(bonds)) -Base.getindex(bonds::ContinuousBonds, i) = @view(bonds.vectors[:, i]) \ No newline at end of file +get_skip_lengths(bonds::ChainedBonds{T}) where T = get_skip_lengths(bonds.lengths, bonds.angles) diff --git a/src/backbone/d2.jl b/src/backbone/d2.jl deleted file mode 100644 index 654ddc6..0000000 --- a/src/backbone/d2.jl +++ /dev/null @@ -1,295 +0,0 @@ -using LinearAlgebra -using NNlib: batched_mul -using Rotations - -export dihedrals2xyz, dihedrals2xyz_exact, get_dihedrals, get_ks_ls, get_ks_ls_dihs, idealize_lengths_angles, new_frame_dihedrals - -# Mean value of bond lengths and bond angles, from approx. 160k PDB values -const MEAN_BOND_LENGTHS = (1.47, 1.53, 1.32) -const MEAN_BOND_ANGLES = (deg2rad(110), deg2rad(117), deg2rad(120)) - -#= -TODO: think about how to better integrate this with the Backbone type. -could also probably distill the set of functions to be more concise and general, -e.g. -=# - -""" -Returns the vectors and lengths connecting each pair of adjacent atoms in the backbone -""" -function bonds_vecs_and_lens(backbone::Backbone{A}) where A - @assert A >= 3 "backbone needs at least the N, Cα and C atoms to calculate bond vectors" - bond_vectors = bond_vectors(backbone) - lengths = reshape(mapslices(norm, bond_vectors, dims=1), :) - return bond_vectors, lengths -end - -""" -Turns a list of vectors into a set of points starting at the origin, where p1 = 0, p2 = v1 + p1, p3 = v2 + p2, etc. -""" -function coords_from_vecs(vecs) - num_atoms = size(vecs, 2) + 1 - coords = zeros(3, num_atoms) - for i in 2:num_atoms - coords[:, i] = coords[:, i-1] .+ vecs[:, i-1] - end - return coords -end - -""" - get_dihedrals(vecs::AbstractMatrix, lengths::AbstractVector) - -Getting dihedral angles from vectors stored as a 3xN matrix and their lengths. -""" -function get_dihedrals(vecs::AbstractMatrix, lengths::AbstractVector) - len_prot = size(vecs, 2) + 1 - - crosses = stack(cross.(eachcol(vecs[:, 1:len_prot-2]), eachcol(vecs[:, 2:len_prot-1])), dims=2) - lengths_cross = reshape(sqrt.(crosses[1, :].^2 .+ crosses[2, :].^2 .+ crosses[3, :].^2), 1, :) - normalized_crosses = crosses ./ lengths_cross - cos_theta = dot.(eachcol(normalized_crosses[:, 1:len_prot-3]), eachcol(normalized_crosses[:, 2:len_prot-2])) - - cross_crosses = stack(cross.(eachcol(normalized_crosses[:, 1:len_prot-3]), eachcol(normalized_crosses[:, 2:len_prot-2])), dims=2) - normalized_vecs = (vecs ./ reshape(lengths, 1, :))[:, 2:len_prot-2] - - sin_theta = dot.(eachcol(cross_crosses), eachcol(normalized_vecs)) - - thetas = atan.(sin_theta, cos_theta) - return thetas -end - -@inline get_dihedrals(backbone::Backbone) = get_dihedrals(bonds_vecs_and_lens(backbone)...) -@inline get_dihedrals(coords::AbstractArray{T, 3}) where T = get_dihedrals(Backbone(coords)) - -""" -Fix the bond lengths of respective bonds to NCa, CaC, and CN, and return the new coordinates. -""" -function fixed_bond_lengths(coords::AbstractArray{T, 3}, NCa, CaC, CN) where T - vecs, lengths = bonds_vecs_and_lens(Backbone(coords)) - Nvecs = vecs ./ reshape(lengths, 1, :) - new_coords = reshape(zeros(size(coords)), 3, :) - new_coords[:, 1] = coords[:, 1, 1] - for i in axes(Nvecs, 2) - if i % 3 == 1 - new_coords[:, i+1] = new_coords[:, i] .+ NCa*Nvecs[:, i] - elseif i % 3 == 2 - new_coords[:, i+1] = new_coords[:, i] .+ CaC*Nvecs[:, i] - else - new_coords[:, i+1] = new_coords[:, i] .+ CN*Nvecs[:, i] - end - end - return reshape(new_coords, 3, 3, :) -end - -# Finds the length of the third side of a triangle given two lengths and their intermediate angle. -@inline law_of_cosines(a, b, C) = sqrt(a^2 + b^2 - 2*a*b*cos(C)) - -""" - idealize_lengths_angles(coords_vector::AbstractVector{<:AbstractArray{T, 3}}; bond_lengths=MEAN_BOND_LENGTHS, bond_angles=MEAN_BOND_ANGLES) where T - -Idealizes the bond lengths and angles of coords_vector while maintaining the same overall structure. coords_vector can be a single 3x3xN matrix or a vector of 3x3xN matrices. -""" -function idealize_lengths_angles( - coords_vector::AbstractVector{<:AbstractArray{T, 3}}; - bond_lengths=MEAN_BOND_LENGTHS, - bond_angles=MEAN_BOND_ANGLES, -) where T - # Code is general as to allow for different bond lengths and angles for different bonds, however it is currently not used. - NCa, CaC, CN = bond_lengths - NCaC, CaCN, CNCa = bond_angles - - NC = law_of_cosines(NCa, CaC, NCaC) - CaN = law_of_cosines(CaC, CN, CaCN) - CCa = law_of_cosines(CN, NCa, CNCa) - - prots_prepped = Array{T, 3}[] - - for i in axes(coords_vector, 1) - org = fixed_bond_lengths(coords_vector[i], NCa, CaC, CN) - orgp = coords_vector[i] - points = reshape(org[:, :, :], :) - - # ks and ls are lengths used for fixing bond angles. - ks = repeat([CaC, CN, NCa], size(orgp, 3)) - ls = repeat([CaN, CCa, NC], size(orgp, 3)) - p = reshape(fix_sequence_of_points(points, ks, ls), 3, 3, :) - push!(prots_prepped, p) - end - return prots_prepped -end - -function idealize_lengths_angles( - coords::AbstractArray{T, 3}; - bond_lengths=MEAN_BOND_LENGTHS, - bond_angles=MEAN_BOND_ANGLES -) where T - return idealize_lengths_angles([coords]; bond_lengths=bond_lengths, bond_angles=bond_angles)[1] -end - -""" -Fixes P3 to the correct bond angle given its bond length and skip length, while keeping it in the same plane defined by the initial P1-P2-P3. -""" -function fix_bond_angle(P1, P2, P3, k, l) - v1 = P1 - P2 - v2 = P3 - P2 - n = cross(v1, v2) - p = P2 - d = norm(P1 - P2) - - # circle of intersection center - h = (k^2-l^2+d^2) / (2*d) - c = P2 + (h/d) * (P1-P2) - - # circle of intersection radius - r = sqrt(k^2 - h^2) - - # plane's point closest to circle - t = dot(n, (p-c)) / dot(n, n) - proj = c + t * n - - # distance from the circle's center to the projection - dist_to_proj = norm(proj - c) - - d = sqrt(r^2 - dist_to_proj^2) # distance from projection to intersection points - dir = cross(n, (P1 - P2)) # direction from proj to the intersection points - dir /= norm(dir) # normalize the direction - - # 2 intersection points - P_1st = proj + d*dir - P_2nd = proj - d*dir - - # choose one with right chirality - i = argmin([norm(P_1st - P3), norm(P_2nd - P3 )]) - return [P_1st, P_2nd][i] -end - -""" -Fix a sequence of points to their bond angles determined by ks and ls (bond lenths and skip lengths). -""" -function fix_sequence_of_points(points, ks, ls) - # Initial checks - new_points = deepcopy(points) - for i in 3:size(points, 2) - new_points[:, i] = fix_bond_angle(new_points[:, i-2], new_points[:, i-1], new_points[:, i], ks[i-2], ls[i-2]) - end - - return new_points -end - -""" -Edits the dihedrals of the given coords to a list or 3xN matrix of new dihedrals. -""" -function dihedrals_to_vecs_respect_bond_angles( - coords::AbstractArray{T, 3}, - new_dihedrals::AbstractVecOrMat -) where T - # Start at v_1. - # Compute dihedral angle v_1 - v2 - v_3 - # Update v_3 so that v_1 - v_2 - v_3 has the new dihedral angle. - # Continue for v_2 - v_3 - v_4 and so on... - new_dihedrals = reshape(new_dihedrals, :) - vectahs, lens = bonds_vecs_and_lens(Backbone(coords)) - for i in 1:size(vectahs, 2)-2 # Seems like there should exist a smart way to vectorize this - - # In every iteration we consider v_{i-1} - v_i - v_{i+1} - curr_vecs = vectahs[:, i:i+2] - lengths = reshape(mapslices(norm, curr_vecs, dims=1), :) - - # Get the rotational update needed, in radians and about the v_i axis - theta = mod.(new_dihedrals[i] .- get_dihedrals(curr_vecs, lengths) .+ π, 2π) .- π - - # Compute the rotation matrix corresponding to that rotation - current_rotation = AngleAxis(theta..., curr_vecs[:, 2]...) - - # Apply the rotation. - vectahs[:, i+2:end] = batched_mul(current_rotation, vectahs[:, i+2:end]) - end - return vectahs, lens -end - -""" - dihedrals2xyz(dihedrals::AbstractVecOrMat, start_res::AbstractMatrix; bond_lengths=MEAN_BOND_LENGTHS, bond_angles=MEAN_BOND_ANGLES) - -Takes an array or a 3xN matrix of dihedrals and a starting residue and returns the xyz coordinates determined by the dihedrals, bond lengths and bond angles. -""" -function dihedrals2xyz( - dihedrals::AbstractVecOrMat, - start_res::AbstractMatrix; - bond_lengths=MEAN_BOND_LENGTHS, - bond_angles=MEAN_BOND_ANGLES, -) - dihedrals3xL = reshape(dihedrals, 3, :) # note there are N-1 sets of three dihedrals for N residues. - init_points = cat(start_res, randn(3, 3, size(dihedrals3xL, 2)), dims = 3) - st = idealize_lengths_angles(init_points, bond_lengths=bond_lengths, bond_angles=bond_angles) - return reshape(coords_from_vecs(dihedrals_to_vecs_respect_bond_angles(st, dihedrals3xL)[1]) .+ start_res[:, 1], 3, 3, :) -end - -""" - get_ks_ls(coords::AbstractArray{T, 3}) where T - -Returns the bond lengths between adjacent atoms (ks) and the skip lengths (ls). -""" -function get_ks_ls(coords::AbstractArray{T, 3}) where T - _, ks = bonds_vecs_and_lens(Backbone(coords)) - flat_coords = reshape(coords, 3, :) - ls = [norm(@view(flat_coords[:, i]) .- @view(flat_coords[:, i+2])) for i in 1:size(flat_coords, 2)-2] - return ks, ls -end - -""" - get_ks_ls_dihs(coords::AbstractArray{T, 3}) where T - -Gets bond lengths, skip lengths, and dihedrals from coords. -""" -function get_ks_ls_dihs(coords::AbstractArray{T, 3}) where T - ks, ls = get_ks_ls(coords) - dihedrals = get_dihedrals(coords) - return ks, ls, dihedrals -end - -""" -Fix all bond lengths in coords to the bond lengths given in l. -""" -function fix_bond_lengths_sequence(coords::AbstractArray{T, 3}, ks::AbstractVector) where T - vecs, lengths = bonds_vecs_and_lens(Backbone(coords)) - Nvecs = vecs ./ reshape(lengths, 1, :) - new_coords = reshape(zeros(size(coords)), 3, :) - new_coords[:, 1] = coords[:, 1, 1] - for i in axes(Nvecs, 2) - new_coords[:, i+1] = new_coords[:, i] .+ ks[i]*Nvecs[:, i] - end - return new_coords -end - -""" -Fixes the coords_vector to the exact sequence of lengths and angles (angles parametrized by skip lengths "ls"), where ks = [NCa1, CaC1, CN1, ...] and ls = [NC, CaN, CCa, ...]. -""" -function fix_sequence_lengths_angles(coords::AbstractArray{T, 3}, ks::AbstractVector, ls::AbstractVector) where T - org = fix_bond_lengths_sequence(coords, ks) - points = reshape(org[:, :, :], :) - p = reshape(fix_sequence_of_points(points, ks[1:end], ls), 3, 3, :) - return p -end - -""" - dihedrals2xyz_exact(dihedrals::AbstractVecOrMat, start_res::AbstractMatrix, ks::AbstractVector, ls::AbstractVector) - -Maps dihedrals, adjacent bond lengths, skips lengths, and a starting residue to xyz coordinates. -""" -function dihedrals2xyz_exact(dihedrals::AbstractVecOrMat, start_res::AbstractMatrix, ks::AbstractVector, ls::AbstractVector) - dihedrals3xL = reshape(dihedrals, 3, :) - init_points = cat(start_res, randn(3, 3, size(dihedrals3xL, 2)), dims = 3) - st = fix_sequence_lengths_angles(init_points, ks, ls) - return reshape(coords_from_vecs(dihedrals_to_vecs_respect_bond_angles(st, dihedrals3xL)[1]) .+ start_res[:, 1], 3, 3, :) -end - -""" - new_frame_dihedrals(frames_prev::AbstractArray{T, 4}, dihedrals::AbstractMatrix) where T - -Takes protxyz and a singular set of dihedrals to place the next frame with idealized bond lengths and angles -""" -function new_frame_dihedrals(frames_prev::AbstractArray{T, 3}, dihedrals::AbstractMatrix) where T - new_frame = dihedrals2xyz(dihedrals[:, end], frames_prev[:, :, end]) - pxyz = cat(frames_prev[:, :, :], new_frame[:, :, end], dims=3) - return pxyz -end diff --git a/src/backbone/dihedrals.jl b/src/backbone/dihedrals.jl deleted file mode 100644 index c0866b5..0000000 --- a/src/backbone/dihedrals.jl +++ /dev/null @@ -1,32 +0,0 @@ -export Dihedrals - -# source: en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics -function dihedral_angle( - u1::T, u2::T, u3::T, -) where T <: AbstractVector - c12 = cross(u1, u2) - c23 = cross(u2, u3) - atan2( - dot(u2, cross(c12, c23)), - abs(u2) * dot(c12, c23) - ) -end - -struct Dihedrals{T} <: AbstractVector{T} - angles::AbstractVector{T} - - function Dihedrals(bonds::ContinuousBonds{T}) where T - angles = Vector{T}(undef, length(bonds)-2) - vectors = eachcol(bonds.vectors) - for (i, j) in zip(eachindex(angles), eachindex(vectors)) - angles[i] = dihedral_angle(vectors[j], vectors[j+1], vectors[j+2]) - end - return new{T}(angles) - end - - Dihedrals(backbone::Backbone) = Dihedrals(ContinuousBonds(backbone)) -end - -@inline Base.size(dihedrals::Dihedrals) = size(dihedrals.angles) -@inline Base.length(dihedrals::Dihedrals) = size(dihedrals, 2) -@inline Base.getindex(dihedrals::Dihedrals, i) = dihedrals.angles[i] \ No newline at end of file diff --git a/src/protein/chain.jl b/src/protein/chain.jl index 020a966..699c10e 100644 --- a/src/protein/chain.jl +++ b/src/protein/chain.jl @@ -34,7 +34,8 @@ end Base.summary(chain::ProteinChain) = "ProteinChain $(chain.id) with $(length(chain)) residue$(length(chain) == 1 ? "" : "s")" Base.show(io::IO, chain::ProteinChain) = print(io, summary(chain)) -@inline Base.getindex(protein::AbstractVector{ProteinChain}, id::AbstractString) = protein[findfirst(c -> c.id == id, protein)] +@inline Base.getindex(protein::AbstractVector{ProteinChain}, i::AbstractString) = protein[findfirst(c -> c.id == i, protein)] +@inline Base.getindex(protein::AbstractVector{ProteinChain}, i::Symbol) = protein[String(i)] export nitrogen_alphacarbon_distances, alphacarbon_carbonyl_distances, carbonyl_nitrogen_distances diff --git a/test/backbone/backbone.jl b/test/backbone/backbone.jl index db7f87a..b08ae8b 100644 --- a/test/backbone/backbone.jl +++ b/test/backbone/backbone.jl @@ -4,6 +4,7 @@ coords = randn(3, 4, 5) backbone = Backbone(coords) + @test backbone isa Backbone{4} @test size(backbone) == (3, 4, 5) @test length(backbone) == 5 @test backbone[1] == coords[:, :, 1] @@ -12,6 +13,5 @@ include("rotations.jl") include("bonds.jl") - include("dihedrals.jl") end \ No newline at end of file diff --git a/test/backbone/bonds.jl b/test/backbone/bonds.jl index 5aba1fe..850430f 100644 --- a/test/backbone/bonds.jl +++ b/test/backbone/bonds.jl @@ -1,9 +1,69 @@ -@testset "distance.jl" begin +@testset "bonds.jl" begin - protein = pdb_to_protein("data/1ASS.pdb") - backbone = protein["A"].backbone - cn_distances = carbonyl_nitrogen_distances(backbone) - @test length(cn_distances) == length(backbone) - 1 - @test 1.32 <= sum(cn_distances) / length(cn_distances) <= 1.34 + @testset "vectors and lengths" begin + + backbone = Backbone{2}([ + 0.0 1.0; + 0.0 2.0; + 0.0 2.0;;; + + 2.0 5.0; + 4.0 8.0; + 4.0 4.0 + ]) + + @testset "get_atom_displacements" begin + @test get_atom_displacements(backbone, 1, 2, 0) == [1.0 3.0; 2.0 4.0; 2.0 0.0] + @test get_atom_displacements(backbone, 2, 1, 1) == [1.0; 2.0; 2.0;;] + end + + @testset "get_atom_distance" begin + @test get_atom_distances(backbone, 1, 2, 0) == [3.0, 5.0] + @test get_atom_distances(backbone, 2, 1, 1) == [3.0] + end + + @testset "get_bond_vectors" begin + @test get_bond_vectors(backbone) == [ + 1.0 1.0 3.0; + 2.0 2.0 4.0; + 2.0 2.0 0.0 + ] + end + + @testset "get_bond_lengths" begin + @test get_bond_lengths(backbone) == [3.0, 3.0, 5.0] + end + + end + + @testset "angles" begin + + coords = [ + 0.0 1.0 1.0 1.0 1.0 2.0 2.0 2.0; + 0.0 0.0 1.0 1.0 2.0 2.0 1.0 0.0; + 0.0 0.0 0.0 1.0 1.0 1.0 0.0 -1.0 + ] + + backbone = Backbone(reshape(coords, 3, 1, :)) + + @testset "bond angles" begin + @test get_bond_angles(backbone) ≈ [π/2, π/2, π/2, π/2, π/2, π] + end + + @testset "dihedrals" begin + @test get_dihedrals(backbone) ≈ [π/2, π, -π/2, π/4, 0] + end + + end + + @testset "ChainedBonds" begin + + @testset "invertibility" begin + backbone = pdb_to_protein("data/1ASS.pdb")["A"].backbone + @test ChainedBonds(Backbone{3}(ChainedBonds(backbone))) ≈ ChainedBonds(backbone) + @test ChainedBonds(Backbone(ChainedBonds(backbone))) ≈ ChainedBonds(backbone) + end + + end end \ No newline at end of file diff --git a/test/backbone/dihedrals.jl b/test/backbone/dihedrals.jl deleted file mode 100644 index 0e014e8..0000000 --- a/test/backbone/dihedrals.jl +++ /dev/null @@ -1,17 +0,0 @@ -@testset "dihedrals.jl" begin - - protein = pdb_to_protein("data/1ASS.pdb") - chain_A = protein["A"] - backbone = chain_A.backbone - backbone3 = Backbone(atom_coords(backbone, 1:3)) - - @testset "idealized lengths and angles" begin - - ideal_backbone_coords = idealize_lengths_angles(backbone3.coords) - ks, ls, dihs = get_ks_ls_dihs(backbone3.coords) - fixed_ideal = Backboner.fix_sequence_lengths_angles(ideal_backbone_coords, ks, ls) - @test all(mapslices(norm, fixed_ideal .- backbone3.coords, dims=1) .< 0.1) - - end - -end \ No newline at end of file diff --git a/test/protein/chain.jl b/test/protein/chain.jl index c1f6ed7..3c2fe3e 100644 --- a/test/protein/chain.jl +++ b/test/protein/chain.jl @@ -13,6 +13,9 @@ @test length(chain) == 5 @test size(chain) == (5,) @test ProteinChain(backbone).id == "_" + + protein = [chain] + @test protein["A"] == protein[:A] == protein[1] @test chain[1] == Residue(1, backbone, 'G', ' ') @@ -24,4 +27,14 @@ end + @testset "atom_distances" begin + + protein = pdb_to_protein("data/1ASS.pdb") + backbone = protein["A"].backbone + cn_distances = carbonyl_nitrogen_distances(backbone) + @test length(cn_distances) == length(backbone) - 1 + @test 1.32 <= sum(cn_distances) / length(cn_distances) <= 1.34 + + end + end \ No newline at end of file