Add aa_properties_matrix (#37)

This is useful for covariance analyses.

Author: timholy

Project.toml

@@ -1,13 +1,14 @@
 name = "GPCRAnalysis"
 uuid = "c1d73f9e-d42a-418a-8d5b-c7b00ec0358f"
-authors = ["Tim Holy <[email protected]> and contributors"]
 version = "0.5.1"
+authors = ["Tim Holy <[email protected]> and contributors"]
 
 [deps]
 BioStructures = "de9282ab-8554-53be-b2d6-f6c222edabfc"
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+FASTX = "c2308a5c-f048-11e8-3e8a-31650f418d12"
 FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 GaussianMixtureAlignment = "f2431ed1-b9c2-4fdb-af1b-a74d6c93b3b3"
@@ -38,6 +39,7 @@ BioStructures = "4.2"
 ColorTypes = "0.11, 0.12"
 Distances = "0.10"
 Downloads = "1"
+FASTX = "2"
 FixedPointNumbers = "0.8"
 GZip = "0.6"
 GaussianMixtureAlignment = "0.2.2, 0.3"

src/GPCRAnalysis.jl

@@ -5,6 +5,7 @@ using Statistics
 using LinearAlgebra
 
 using BioStructures
+using FASTX
 
 using MIToS: MIToS, Pfam, MSA
 using MIToS.MSA: AbstractMultipleSequenceAlignment, AnnotatedAlignedSequence, AnnotatedMultipleSequenceAlignment,
@@ -51,7 +52,7 @@ export project_sequences, columnwise_entropy, align, residue_centroid, residue_c
        alphacarbon_coordinates_matrix, chargelocations, positive_locations, negative_locations
 export StructAlign, residueindex, ismapped, skipnothing
 export BWScheme, lookupbw
-export aa_properties, aa_properties_zscored
+export aa_properties, aa_properties_zscored, aa_properties_matrix
 export sidechaincentroid, scvector, inward_tm_residues, inward_ecl_residues
 export features_from_structure
 export forcecomponents, optimize_weights, forcedict

src/properties.jl

@@ -1,12 +1,23 @@
+# Biophysical properties of amino acids
+
 struct AAProperties <: AbstractVector{Float64}
     charge::Float64
     hydropathy::Float64   # from https://www.sciencedirect.com/science/article/pii/0022283682905150?via%3Dihub & https://en.wikipedia.org/wiki/Amino_acid
     volume::Float64       # from http://proteinsandproteomics.org/content/free/tables_1/table08.pdf (van der Waals volume in ų)
 end
-Base.size(v::AAProperties) = (3,)
+Base.size(::AAProperties) = (3,)
 Base.getindex(v::AAProperties, i::Int) = i == 1 ? v.charge :
                                          i == 2 ? v.hydropathy :
                                          i == 3 ? v.volume : Base.throw_boundserror(v, i)
+
+StaticArrays.SVector(v::AAProperties) = SVector{3}(v.charge, v.hydropathy, v.volume)
+
+Base.:(+)(v::AAProperties, w::AAProperties) = AAProperties(v.charge + w.charge, v.hydropathy + w.hydropathy, v.volume + w.volume)
+Base.:(-)(v::AAProperties, w::AAProperties) = AAProperties(v.charge - w.charge, v.hydropathy - w.hydropathy, v.volume - w.volume)
+Base.:(/)(v::AAProperties, r::Real) = AAProperties(v.charge / r, v.hydropathy / r, v.volume / r)
+Base.:(*)(u::AAProperties, v::Adjoint{Float64, GPCRAnalysis.AAProperties}) = SVector(u) * SVector(v')'
+Base.:(\)(L::LowerTriangular{T, StaticArraysCore.SMatrix{3, 3, T, 9}}, v::AAProperties)  where T<:Real = L \ SVector(v)
+
 const aa_properties = Dict(
     MSA.Residue('A') => AAProperties(0, 1.8, 67),
     MSA.Residue('R') => AAProperties(1, -4.5, 148),
@@ -35,4 +46,18 @@ const featμ = sum(p for (r, p) in aa_properties)/length(aa_properties)
 const featC = sum((Δp = p - featμ; Δp*Δp') for (r, p) in aa_properties)/(length(aa_properties)-1)
 const featChol = cholesky(featC)
 
-const aa_properties_zscored = Dict(r => SVector{3}(featChol.L \ (p - featμ)) for (r, p) in aa_properties)
+const aa_properties_zscored = Dict(r => (featChol.L \ (p - featμ)) for (r, p) in aa_properties)
+
+function aa_properties_matrix(msa::AbstractMultipleSequenceAlignment)
+    props = copy(aa_properties_zscored)
+    props[GAP] = zero(valtype(props))
+    props[MSA.Residue('X')] = zero(valtype(props))
+    return [props[residue] for residue in permutedims(msa)]
+end
+
+function aa_properties_matrix(seqs::AbstractVector{FASTX.FASTA.Record})
+    props = Dict(Char(r) => v for (r, v) in aa_properties_zscored)
+    props['-'] = zero(valtype(props))
+    props['X'] = zero(valtype(props))
+    return reduce(hcat, [[props[r] for r in sequence(rec)] for rec in seqs])
+end

test/runtests.jl

@@ -1,8 +1,9 @@
 using GPCRAnalysis
 using MIToS: MSA, Pfam
 using MIToS.MSA: @res_str, three2residue, coverage, GappedAlphabet, filtersequences!, percentsimilarity,
-                 nsequences, sequencenames, getsequencemapping
+                 nsequences, sequencenames, getsequencemapping, getcolumnmapping
 using BioStructures
+using FASTX
 using GaussianMixtureAlignment
 using InvertedIndices
 using Statistics
@@ -85,6 +86,22 @@ using Test
         @test MSACode(msa, AccessionCode("Q3V4T1")) == MSACode("Y070_ATV/2-70")
         @test msa[MSACode("Y070_ATV/2-70")][8] == msa[AccessionCode("Q3V4T1")][8] == MSA.Residue('V')
     end
+    @testset "Properties" begin
+        pf09645_sto = "PF09645_full.stockholm"
+        msa = MSA.read_file(pf09645_sto, Pfam.Stockholm)
+        X = aa_properties_matrix(msa)
+        ΔX = X .- mean(X, dims=2)
+        i = findfirst(==(14), getcolumnmapping(msa))
+        @test all(iszero, ΔX[i, :])
+        @test !all(iszero, ΔX[i-1, :])
+        seqs = FASTAReader(open("test.fasta")) do io
+            collect(io)
+        end
+        X = aa_properties_matrix(seqs)
+        ΔX = X .- mean(X, dims=2)
+        @test all(iszero, ΔX[1, :])
+        @test !all(iszero, ΔX[3, :])
+    end
     @testset "ChimeraX" begin
         tmpfile = tempname() * ".cxc"
         chimerax_script(tmpfile, ["ABCD.pdb", "EFGH.pdb"], [[5, 10, 15], [7, 11]])

test/test.fasta

@@ -0,0 +1,4 @@
+> Seq:1
+ACA
+> Seq:2
+ACC