From 29cd5cff40dd714e1a21bc1a6709e9efa8d7e771 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 6 Dec 2022 07:43:29 -0600 Subject: [PATCH 1/3] Compute PDB "fields" Currently supports electrostatic potential and steric hindrance. --- Project.toml | 1 + src/GPCRAnalysis.jl | 4 ++- src/analyze.jl | 78 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 82 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3a0790f..2b8b598 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411" Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" TravelingSalesmanHeuristics = "8c8f4381-2cdd-507c-846c-be2bcff6f45f" [compat] diff --git a/src/GPCRAnalysis.jl b/src/GPCRAnalysis.jl index 56ae8e1..1c635a9 100644 --- a/src/GPCRAnalysis.jl +++ b/src/GPCRAnalysis.jl @@ -13,6 +13,7 @@ using MultivariateStats using Distances using Hungarian using StaticArrays +using StructArrays using TravelingSalesmanHeuristics export @res_str @@ -21,7 +22,8 @@ export SequenceMapping export species, uniprotX export try_download_alphafold, download_alphafolds, getchain, findall_subseq export filter_species!, filter_long!, sortperm_msa, chimerax_script -export project_sequences, columnwise_entropy, align, residue_centroid, residue_centroid_matrix, mapclosest, chargelocations, positive_locations, negative_locations +export project_sequences, columnwise_entropy, align, residue_centroid, residue_centroid_matrix, mapclosest +export chargelocations, positive_locations, negative_locations, charge_magnitudes, fields export StructAlign, residueindex, ismapped export BWScheme, lookupbw export aa_properties, aa_properties_zscored diff --git a/src/analyze.jl b/src/analyze.jl index d832859..bd79e4e 100644 --- a/src/analyze.jl +++ b/src/analyze.jl @@ -146,3 +146,81 @@ end positive_locations(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] for loc in chargelocs if loc[3] ∈ ("ARG", "LYS", "HIS")] negative_locations(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] for loc in chargelocs if loc[3] ∈ ("ASP", "GLU")] + +_charge_magnitude(aa) = aa == "LYS" ? 1.0 : + aa ∈ ("ARG", "HIS") ? 0.5 : + aa ∈ ("ASP", "GLU") ? -0.5 : 0.0 +charge_magnitudes(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] => _charge_magnitude(loc[3]) for loc in chargelocs] + +struct Fields{T} <: FieldVector{2, T} + electrostatic::T + steric::T +end + +gvwr(resname, a::PDBAtom) = get(vanderwaalsradius, + (resname, a.atom), + get(vanderwaalsradius, (resname, a.element), nothing)) +steric(x::Coordinates, y::Coordinates, σ) = exp(-distance(x, y)^2 / (2 * σ^2)) +function steric(x::Coordinates, pdb::AbstractVector{PDBResidue}, σfac=1) + s = 0.0 + for aa in pdb + resname = aa.id.name + for a in aa.atoms + vdwr = gvwr(resname, a) + vdwr === nothing && continue + s += steric(x, a.coordinates, vdwr * σfac) + end + end + return s +end +function steric!(ss::AbstractArray, xs::AbstractArray{Coordinates}, pdb::AbstractVector{PDBResidue}, σfac=1) + # Performance optimization of `steric` + axes(ss) == axes(xs) || throw(DimensionMismatch("ss and xs must have commensurate axes, got $(axes(ss)) and $(axes(xs))")) + fill!(ss, 0.0) + for aa in pdb + resname = aa.id.name + for a in aa.atoms + vdwr = gvwr(resname, a) + vdwr === nothing && continue + for (i, x) in pairs(xs) + ss[i] += steric(x, a.coordinates, vdwr * σfac) + end + end + end + return ss +end + +const membrane_width = 34 # Å, hydrophobic thickness (see https://pubs.rsc.org/en/content/articlehtml/2016/sm/c6sm01777k) +const λB1 = 568 # Bjerrum length for ϵr = 1, in Å + +function dielectric(density::Real, location::Coordinates) + extramembrane = 1 / (1 + exp(abs(2*location[3]/membrane_width)^2)) + ϵext = 25 * extramembrane + 2.2 * (1-extramembrane) + return 6.5 * density + ϵext * (1 - density) +end + +function fields(pdb::AbstractVector{PDBResidue}, locations, r0=1.0) + sterics = similar(locations, Float64) + steric!(sterics, locations, pdb) + # Get the local dielectric constant: see + # On the Dielectric "Constant" of Proteins: Smooth Dielectric Function for Macromolecular Modeling and Its Implementation in DelPhi, + # Lin Li, Chuan Li, Zhe Zhang, Emil Alexov + # J Chem Theory Comput. 2013 Apr 9;9(4):2126-2136. doi: 10.1021/ct400065j + # https://pubmed.ncbi.nlm.nih.gov/23585741/ + # We use the "density" as an indicator of whether we are interior or exterior + density = similar(locations, Float64) + steric!(density, locations, pdb, 5.0) + dmin, dmax = extrema(density) + dielec = dielectric.((density .- dmin) ./ (dmax - dmin), locations) + cmags = charge_magnitudes(chargelocations(pdb)) + electrostatic = similar(locations, Float64) + fill!(electrostatic, 0.0) + for (i, x) in pairs(locations) + λ = λB1 / dielec[i] + for (y, v) in cmags + r = distance(x, y) + r0 + electrostatic[i] += v/r * exp(-r/λ) + end + end + return StructArray{Fields{Float64}}((; electrostatic, steric=sterics)) +end From 215e4172477da27bca56cc5cf35894af3b21bc0e Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Wed, 7 Dec 2022 10:44:39 -0600 Subject: [PATCH 2/3] Switch to Lennard-Jones Also tweaked a couple of constants via visual inspection --- src/analyze.jl | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/analyze.jl b/src/analyze.jl index bd79e4e..24fc32d 100644 --- a/src/analyze.jl +++ b/src/analyze.jl @@ -160,20 +160,31 @@ end gvwr(resname, a::PDBAtom) = get(vanderwaalsradius, (resname, a.atom), get(vanderwaalsradius, (resname, a.element), nothing)) -steric(x::Coordinates, y::Coordinates, σ) = exp(-distance(x, y)^2 / (2 * σ^2)) -function steric(x::Coordinates, pdb::AbstractVector{PDBResidue}, σfac=1) + +# Repulsive "meta-potentials" to estimate steric interactions +gaussian(x::Coordinates, y::Coordinates, σ) = exp(-distance(x, y)^2 / (2 * σ^2)) + +function lennardjones(x::Coordinates, y::Coordinates, σ) + # Lennard-Jones potential + r = distance(x, y) + k = σ / r + k6 = k^6 + return iszero(r) ? Inf : k6^2 - k6 +end + +function steric(f::F, x::Coordinates, pdb::AbstractVector{PDBResidue}, σfac=1) where F s = 0.0 for aa in pdb resname = aa.id.name for a in aa.atoms vdwr = gvwr(resname, a) vdwr === nothing && continue - s += steric(x, a.coordinates, vdwr * σfac) + s += f(x, a.coordinates, vdwr * σfac) end end return s end -function steric!(ss::AbstractArray, xs::AbstractArray{Coordinates}, pdb::AbstractVector{PDBResidue}, σfac=1) +function steric!(f::F, ss::AbstractArray, xs::AbstractArray{Coordinates}, pdb::AbstractVector{PDBResidue}, σfac=1) where F # Performance optimization of `steric` axes(ss) == axes(xs) || throw(DimensionMismatch("ss and xs must have commensurate axes, got $(axes(ss)) and $(axes(xs))")) fill!(ss, 0.0) @@ -183,7 +194,7 @@ function steric!(ss::AbstractArray, xs::AbstractArray{Coordinates}, pdb::Abstrac vdwr = gvwr(resname, a) vdwr === nothing && continue for (i, x) in pairs(xs) - ss[i] += steric(x, a.coordinates, vdwr * σfac) + ss[i] += f(x, a.coordinates, vdwr * σfac) end end end @@ -194,14 +205,15 @@ const membrane_width = 34 # Å, hydrophobic thickness (see https://pubs.rsc.or const λB1 = 568 # Bjerrum length for ϵr = 1, in Å function dielectric(density::Real, location::Coordinates) - extramembrane = 1 / (1 + exp(abs(2*location[3]/membrane_width)^2)) - ϵext = 25 * extramembrane + 2.2 * (1-extramembrane) + extramembrane = 1 / (1 + exp(abs(2*location[3]/membrane_width)^8)) + # @assert 0.0 <= extramembrane <= 1.0 + ϵext = 30 * extramembrane + 2.2 * (1-extramembrane) return 6.5 * density + ϵext * (1 - density) end -function fields(pdb::AbstractVector{PDBResidue}, locations, r0=1.0) +function fields(pdb::AbstractVector{PDBResidue}, locations, r0=0.0) sterics = similar(locations, Float64) - steric!(sterics, locations, pdb) + steric!(lennardjones, sterics, locations, pdb) # Get the local dielectric constant: see # On the Dielectric "Constant" of Proteins: Smooth Dielectric Function for Macromolecular Modeling and Its Implementation in DelPhi, # Lin Li, Chuan Li, Zhe Zhang, Emil Alexov @@ -209,7 +221,7 @@ function fields(pdb::AbstractVector{PDBResidue}, locations, r0=1.0) # https://pubmed.ncbi.nlm.nih.gov/23585741/ # We use the "density" as an indicator of whether we are interior or exterior density = similar(locations, Float64) - steric!(density, locations, pdb, 5.0) + steric!(gaussian, density, locations, pdb) dmin, dmax = extrema(density) dielec = dielectric.((density .- dmin) ./ (dmax - dmin), locations) cmags = charge_magnitudes(chargelocations(pdb)) From 766d056376f81ba15d4f448e501ef79f5bdd51b4 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Sat, 10 Dec 2022 08:32:12 -0600 Subject: [PATCH 3/3] Add [compat] for StructArrays --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 2b8b598..e03176e 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ MIToS = "2" MultivariateStats = "0.9, 0.10" Rotations = "1" StaticArrays = "1" +StructArrays = "0.6" TravelingSalesmanHeuristics = "0.3" julia = "1.6"