Add atom-in-cyclinder features

Project.toml

@@ -1,6 +1,6 @@
 name = "GPCRAnalysis"
 uuid = "c1d73f9e-d42a-418a-8d5b-c7b00ec0358f"
-version = "0.6.1"
+version = "0.6.2"
 authors = ["Tim Holy <tim.holy@gmail.com> and contributors"]
 
 [deps]
@@ -15,6 +15,7 @@ GaussianMixtureAlignment = "f2431ed1-b9c2-4fdb-af1b-a74d6c93b3b3"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 Hungarian = "e91730f6-4275-51fb-a7a0-7064cfbd3b39"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411"
@@ -51,6 +52,7 @@ HTTP = "1"
 HiGHS = "1"
 Hungarian = "0.6, 0.7"
 Interpolations = "0.15, 0.16"
+IntervalSets = "0.7.11"
 JSON3 = "1"
 JuMP = "1"
 LinearAlgebra = "1"

src/GPCRAnalysis.jl

@@ -15,6 +15,7 @@ using Hungarian
 using ProgressMeter
 using StaticArrays
 using Rotations
+using IntervalSets
 using TravelingSalesmanHeuristics
 
 using Interpolations

src/align.jl

@@ -88,13 +88,14 @@ and false otherwise.
 
 `applytransform!(chain, tform)` (or the `model` that includes `chain`) will
 re-orient `chain` so that the center of the membrane is `z=0` and extracellular
-is positive.
+is positive. Moreover, the mean `x` and `y` position of atoms in the
+transmembrane residues will be zero.
 
 The algorithm finds the membrane normal `u` by maximizing the ratio
 
-    sumᵢ (u ⋅ vᵢ)²
-  -----------------
-    sumᵢ (u ⋅ δᵢ)²
+    Σᵢ (u ⋅ vᵢ)²
+    ------------
+    Σᵢ (u ⋅ δᵢ)²
 
 where `vᵢ` is a vector parallel to the `i`th TM helix, and `δᵢ` is a
 within-leaflet atomic displacement.
@@ -111,8 +112,8 @@ function align_to_membrane(chain::ChainLike, tms; extracellular::Bool=true)
     # Collect the atom displacements in each leaflet
     delta_e = leaflet_displacements(leaflet_e)
     delta_i = leaflet_displacements(leaflet_i)
-    # Also collect the vectors for the TM helices
-    tm_vectors = [residue_centroid(chain[first(r)]) - residue_centroid(chain[last(r)]) for r in tms]
+    # Also collect the vectors for the TM helices. Use the CA atoms since these are part of the backbone.
+    tm_vectors = [coords(chain[first(r)]["CA"]) - coords(chain[last(r)]["CA"]) for r in tms]
     # Find the membrane normal u, maximizing the ratio
     #         sumᵢ (u ⋅ vᵢ)²
     #       -----------------
@@ -133,7 +134,7 @@ function align_to_membrane(chain::ChainLike, tms; extracellular::Bool=true)
     t = [0, 0, mean(proj_e) > mean(proj_i) ? 1 : -1]   # extracellular is positive
     q = QuatRotation(1 + u'*t, cross(u, t)...)
     # Apply the rotation to the coordinates, and offset the z-coordinate to the mean of the leaflets
-    ccoords = coordarray(chain)
+    ccoords = coordarray(collectresidues(chain)[reduce(vcat, tms)])
     μ = vec(mean(ccoords, dims=2))
     return Transformation(μ, [0, 0, (median(proj_e) + median(proj_i))/2 - μ'*u], RotMatrix(q))
 end

src/features.jl

@@ -163,3 +163,36 @@ function features_from_structure(seq, idxs=1:length(seq); kwargs...)
     end
     return mgmm
 end
+
+"""
+    mgmm = features_from_structure(seq::ChainLike, ρmax::Real, zi::AbstractInterval)
+
+Construct an `IsotropicMultiGMM` from `seq` including all atoms that lie within
+the cylinder
+
+    x^2 + y^2 <= ρmax^2
+    z ∈ zi
+
+`zi` is an `AbstractInterval`, e.g. `0..30` or `-15..15` (see IntervalSets.jl).
+
+This implicitly assumes that you've aligned `seq` to the membrane, or aligned
+`seq` to a homolog that is membrane-aligned. See [`align_to_membrane`](@ref),
+[`align`](@ref).
+"""
+function features_from_structure(
+        seq, ρmax::Real, zi::AbstractInterval;
+        σfun = atomic_σfun,
+        ϕfun = atomic_ϕfun,
+    )
+    mgmm = IsotropicMultiGMM(Dict{Symbol,IsotropicGMM{3,Float64}}())
+    for r in seq
+        for a in r
+            c = coords(a)
+            ρ = hypot(c[1], c[2])
+            if ρ < ρmax && c[3] ∈ zi
+                add_features_from_atom!(mgmm, a, r, σfun, ϕfun)
+            end
+        end
+    end
+    return mgmm
+end

test/runtests.jl

@@ -9,6 +9,7 @@ using BioStructures
 using FASTX
 using GaussianMixtureAlignment
 using InvertedIndices
+using IntervalSets
 using Statistics
 using LinearAlgebra
 using JuMP, HiGHS
@@ -240,6 +241,8 @@ using Test
         zi = median(ci[3,:])
         @test ze > zi    # extracellular is positive
         @test abs(ze + zi) < 1e-6 * (ze - zi) # center of membrane is at z=0
+        ccoords = coordarray(collectresidues(opsd)[reduce(vcat, opsd_tms)])
+        @test norm(mean(ccoords, dims=2)[1:2]) < 1000 * eps(eltype(ccoords))  # mean x and y is zero
 
         tm_res = inward_tm_residues(opsd, opsd_tms[[2,3,5,6,7]])
         @test length(tm_res) == 5
@@ -285,6 +288,14 @@ using Test
                 @test g1.ϕ >= g2.ϕ
             end
         end
+
+        cyl_mgmm = features_from_structure(opsd, 12, 0 .. 8)
+        for (k, gmm) in cyl_mgmm
+            for g in gmm
+                @test g.μ[1]^2 + g.μ[2]^2 <= 12^2
+                @test g.μ[3] ∈ 0 .. 8
+            end
+        end
     end
 
     @testset "Forces" begin