polypeptide

.vscode/settings.json

@@ -0,0 +1 @@
+{}

src/distributions/EvolHMM.jl

@@ -1,5 +1,3 @@
-
-
 struct EvolHMM <: ContinuousMatrixDistribution
     M::Integer                      # Total alignment length M = sum(lengths)
     lengths::AbstractVector{Int}
@@ -10,12 +8,13 @@ struct EvolHMM <: ContinuousMatrixDistribution
     α::Array{Real}
 end
 
-function PairHMM(n, m, A, sub, diff, t)
+function EvolHMM(n, m, A, sub, diff, t)
     α = Array{Real}(undef, n+2, m+2, 4)
     return PairHMM(n, m, A, sub, diff, t, α)
 end
 
-function _logpdf(d::PairHMM, x::AbstractMatrix{<:Real})
+function _logpdf(d::EvolHMM, x::AbstractMatrix{<:Real})
+
     return forward!(d.α, x[1:3, 1:d.n], x[4:6, 1:d.m], d.A, d.sub, d.diff, d.t)
 end
 

src/distributions/SubstitutionModel.jl

@@ -20,13 +20,7 @@ using DelimitedFiles
 # Q = S diag(Π) where S is symmetric exchangeability matrix,
 #        Π is equilibrium frequencies
 
-aminoacids = "ARNDCQEGHILKMFPSTWYV-"
 
-id_to_aa(i) = aminoacids[round(Int, i)]
-
-aminoacid_ids = Dict((aminoacids[i], i) for i ∈ eachindex(aminoacids))
-
-aa_to_id(a) = aminoacid_ids[a]
 
 WAG_S = LowerTriangular(map((x) -> x=="" ? 0.0 : x,
                         readdlm("./data/params/WAG_S.csv"; skipblanks=false, dims=(20,20))))

src/model/Polypeptide.jl

@@ -0,0 +1,50 @@
+using BioSequences
+using BioStructures
+
+aminoacids = "ARNDCQEGHILKMFPSTWYV-"
+id_to_aa(i) = aminoacids[round(Int, i)]
+
+aminoacid_ids = Dict((aminoacids[i], Float(i)) for i ∈ eachindex(aminoacids))
+aa_to_id(a) = aminoacid_ids[a]
+
+# Methods to extract internal coordinates data from BioStructures.Chain object
+sequence(chain::Chain) = aa_to_id.(collect(string(LongAA(chain, standardselector))))
+phi_angles(chain::Chain) = phiangles(chain, standardselector)
+psi_angles(chain::Chain) = phiangles(chain, standardselector)
+omega_angles(chain::Chain) = omegaangles(chain, standardselector)
+calpha_coords(chain::Chain) = coordarray(chain, calphaselector)
+
+# ▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄
+# Polypeptide
+
+# Each column of data represents the internal coordinates of a residue
+# The data in each row is specific by row_names
+struct Polypeptide
+    data::AbstractMatrix{<:Real}
+    row_names::AbstractVector{String}
+    chain::BioStructures.Chain
+end
+
+# Construct Polypeptide from BioStructures.Chain
+function Polypeptide(chain::Chain; primary=true,
+                     phi=true, psi=true, omega=false,
+                     cartesian=false)
+    rows = []
+    row_names = []
+
+    primary && (push!(rows, sequence(chain)); push!(row_names, "primary"))
+    phi && (push!(rows, phi_angles(chain)); push!(row_names, "ϕ angles"))
+    psi && (push!(rows, psi_angles(chain)); push!(row_names, "ψ angles"))
+    omega && (push!(rows, omega_angles(chain)); push!(row_names, "ω angles"))
+    cartesian && (push!(rows, calpha_coords(chain)); push!(row_names, "Cα x coords");
+                  push!(row_names, "Cα y coords"); push!(row_names, "Cα z coords"))
+
+    data = hcat(rows...)'
+
+    return Polypeptide(data, row_names, chain)
+end
+
+num_residues(p::Polypeptide) = size(p.data, 2)
+num_coords(p::Polypeptide) = size(p.data, 1)
+data(p::Polypeptide) = p.data
+chain(p::Polypeptide) = p.chain

src/utils/Backbone.jl

@@ -96,7 +96,7 @@ function build_chain_from_alignment(chain::Chain, alignment, Y)
 
     bond_angles_X, bond_lengths_X = backbone_bond_angles_and_lengths(chain)
 
-    #TODO write this more succintly lol
+    #TODO write this more succintly
     bond_angles = Matrix{Real}(undef, 3, n)
     for i in indicesIY
         bond_angles[1, i] = rand(WrappedNormal(ideal_bond_angle["CA-C-N", aminoacids[i]]...))[1]