WIP: Add BioStockholm extension

Author: timholy

Project.toml

@@ -1,7 +1,7 @@
 name = "GPCRAnalysis"
 uuid = "c1d73f9e-d42a-418a-8d5b-c7b00ec0358f"
-authors = ["Tim Holy <[email protected]> and contributors"]
 version = "0.6.0"
+authors = ["Tim Holy <[email protected]> and contributors"]
 
 [deps]
 BioStructures = "de9282ab-8554-53be-b2d6-f6c222edabfc"
@@ -27,15 +27,18 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TravelingSalesmanHeuristics = "8c8f4381-2cdd-507c-846c-be2bcff6f45f"
 
 [weakdeps]
+BioStockholm = "eeb925a3-6f9d-43e6-829e-e0ea03b76ecf"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MIToS = "51bafb47-8a16-5ded-8b04-24ef4eede0b5"
 
 [extensions]
+GPCRAnalysisBioStockholmExt = "BioStockholm"
 GPCRAnalysisJuMPExt = ["JuMP", "HiGHS"]
 GPCRAnalysisMIToSExt = "MIToS"
 
 [compat]
+BioStockholm = "0.2.1"
 BioStructures = "4.2"
 ColorTypes = "0.11, 0.12"
 Distances = "0.10"
@@ -63,11 +66,12 @@ TravelingSalesmanHeuristics = "0.3"
 julia = "1.10"
 
 [extras]
+BioStockholm = "eeb925a3-6f9d-43e6-829e-e0ea03b76ecf"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 InvertedIndices = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MIToS = "51bafb47-8a16-5ded-8b04-24ef4eede0b5"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["HiGHS", "InvertedIndices", "JuMP", "MIToS", "Test"]
+test = ["BioStockholm", "HiGHS", "InvertedIndices", "JuMP", "MIToS", "Test"]

ext/GPCRAnalysisBioStockholmExt.jl

@@ -0,0 +1,96 @@
+module GPCRAnalysisBioStockholmExt
+
+using GPCRAnalysis
+using BioStockholm
+using BioStockholm: OrderedDict   # from OrderedCollections.jl
+
+function conscols(msa::MSA)
+    ss = msa.GC["SS_cons"]
+    return findfirst(!=( '.'), ss):findlast(!=( '.'), ss)
+end
+
+# Low-level API implementation
+# GPCRAnalysis.sequenceindexes(msaseq::AnnotatedAlignedSequence) = getsequencemapping(msaseq)
+# GPCRAnalysis.sequenceindexes(msaseq::MSA, i::Int) = getsequencemapping(msaseq, i)
+# GPCRAnalysis.isgap(res::MSA.Residue) = res == GAP
+# GPCRAnalysis.isunknown(res::MSA.Residue) = res == XAA
+GPCRAnalysis.sequencekeys(msa::MSA) = collect(keys(msa.seq))
+# GPCRAnalysis.msasequence(msa::AbstractMultipleSequenceAlignment, key) = getsequence(msa, key)
+function GPCRAnalysis.residuematrix(msa::MSA)
+    keepcols = conscols(msa)
+    reduce(vcat, [permutedims(seq[keepcols]) for (_, seq) in msa.seq])
+end
+GPCRAnalysis.subseqs(msa::MSA{T}, rowmask::AbstractVector{Bool}) where T = MSA{T}(OrderedDict(pr for (pr, keep) in zip(msa.seq, rowmask) if keep), msa.GF, OrderedDict(pr for (pr, keep) in zip(msa.GS, rowmask) if keep), msa.GC, msa.GR)
+function GPCRAnalysis.subseqs!(msa::MSA, rowmask::AbstractVector{Bool})
+    for ((key, _), keep) in zip(msa.seq, rowmask)
+        if !keep
+            delete!(msa.seq, key)
+            delete!(msa.GS, key)
+        end
+    end
+    return msa
+end
+
+# These are piracy
+function Base.eachrow(msa::MSA)
+    keepcols = conscols(msa)
+    return Iterators.map(values(msa.seq)) do s
+        return s[keepcols]
+    end
+end
+Base.size(msa::MSA) = (length(msa.seq), length(conscols(msa)))
+Base.size(msa::MSA, dim::Int) = dim == 1 ? length(msa.seq) : length(conscols(msa))
+
+Base.getindex(msa::MSA, seqname::MSACode) = msa.seq[seqname.name]
+Base.getindex(msa::MSA, seqname::AccessionCode) = msa[MSACode(msa, seqname)]
+
+
+function GPCRAnalysis.AccessionCode(msa::MSA, seqname::AbstractString)
+    AccessionCode(split(msa.GS[seqname]["AC"], '.')[1])
+end
+GPCRAnalysis.AccessionCode(msa::MSA, seqname::MSACode) = AccessionCode(msa, seqname.name)
+GPCRAnalysis.AccessionCode(::MSA, seqname::AccessionCode) = seqname
+
+function GPCRAnalysis.MSACode(msa::MSA, accession::AbstractString)
+    acs = [split(ac["AC"], '.')[1] for (_, ac) in msa.GS]
+    i = findfirst(==(accession), acs)
+    return MSACode(GPCRAnalysis.sequencekeys(msa)[i])
+end
+GPCRAnalysis.MSACode(msa::MSA, accession::AccessionCode) = MSACode(msa, accession.name)
+GPCRAnalysis.MSACode(::MSA, accession::MSACode) = accession
+
+
+function reduced_alphabet(r::Char)
+    if r in ('A','I','L','M','V')
+        return 1  # hydrophobic
+    elseif r in ('N','Q','S','T')
+        return 2  # polar
+    elseif r in ('R','H','K')
+        return 3  # charged
+    elseif r in ('D','E')
+        return 4  # charged
+    elseif r in ('F','W','Y')
+        return 5  # aromatic
+    end
+    offset = findfirst(==(r), ('C','G','P','-'))
+    offset === nothing && throw(ArgumentError("Unknown residue '$r'"))
+    return 5 + offset  # special or unknown
+end
+
+GPCRAnalysis.columnwise_entropy(msa) = columnwise_entropy(reduced_alphabet, msa)
+
+function GPCRAnalysis.percent_similarity(msa)
+    M = reduced_alphabet.(GPCRAnalysis.residuematrix(msa))
+    n = size(M, 1)
+    S = zeros(Float64, n, n)
+    for i in 1:n
+        for j in i:n
+            S[i, j] = 100 * sum(M[i, :] .== M[j, :]) / size(M, 2)
+            S[j, i] = S[i, j]
+        end
+    end
+    return S
+end
+
+
+end

test/PF09645_full.stockholm

@@ -11,4 +11,4 @@ F112_SSV1/3-112                .....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRA
 #=GR F112_SSV1/3-112     SS    .....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX.....
 #=GC SS_cons                   .....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX.....
 #=GC seq_cons                  ........NshphAclhaKILppKtElolEDIlAQFEISsosAYsI.+sL+hICEpH.-ECpsppKsRKTlhh.hKpEphppptpEp..ppItKIhsAp................h....
-//
+//

test/runtests.jl

@@ -5,7 +5,7 @@ using GPCRAnalysis: sequenceindexes, isgap, isunknown, sequencekeys, msasequence
                      subseqs, subseqs!, percent_similarity
 using MIToS: MSA, Pfam
 using MIToS.MSA: coverage, GappedAlphabet, nsequences
-                 nsequences, sequencenames, getsequencemapping, getcolumnmapping
+using BioStockholm: BioStockholm
 using BioStructures
 using FASTX
 using GaussianMixtureAlignment
@@ -17,6 +17,7 @@ using ColorTypes
 using Test
 
 columnindexes(msa::MSA.AbstractMultipleSequenceAlignment) = MSA.getcolumnmapping(msa)
+columnindexes(msa::BioStockholm.MSA) = eachindex(first(msa.seq).second)
 
 # skip the network-hitting components by setting `skip_download = true` in the global namespace
 
@@ -66,6 +67,7 @@ columnindexes(msa::MSA.AbstractMultipleSequenceAlignment) = MSA.getcolumnmapping
     @testset "MSA" begin
         # The test file is copied from MIToS/test/data, with gratitude
         pf09645_sto = "PF09645_full.stockholm"
+        ## First in MIToS format
         msa = MSA.read_file(pf09645_sto, Pfam.Stockholm)
         @test MSA.nsequences(filter_species!(deepcopy(msa), "ATV")) == 1
         @test MSA.nsequences(filter_long!(deepcopy(msa), 70)) == 3
@@ -91,15 +93,43 @@ columnindexes(msa::MSA.AbstractMultipleSequenceAlignment) = MSA.getcolumnmapping
         @test AccessionCode(msa, MSACode("Y070_ATV/2-70")) == AccessionCode("Q3V4T1")
         @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')
+
+        ## Now in BioStockholm format
+        msa = read(pf09645_sto, BioStockholm.MSA)
+        @test length(sequencekeys(filter_species!(deepcopy(msa), "ATV"))) == 1
+        @test length(sequencekeys(filter_long!(deepcopy(msa), 70))) == 3
+
+        idx = SequenceMapping([0, 4, 5, 0])
+        seqvals = fill(NaN, 9)
+        seqvals[idx] = [0.1, 0.2, 0.3, 0.4]
+        @test seqvals[4] == 0.2
+        @test seqvals[5] == 0.3
+        @test all(isnan, seqvals[1:3])
+        @test all(isnan, seqvals[6:end])
+
+        # analyze
+        e = columnwise_entropy(msa)
+        @test length(e) == size(msa, 2) && e[9] == 0
+        e2 = columnwise_entropy(identity, msa)
+        @test all(e2 .>= e)
+        @test !all(e2 .== e)
+
+        @test size(project_sequences(msa)) == (3, 4)
+        @test size(project_sequences(msa; fracvar=0.5)) == (1, 4)
+
+        @test AccessionCode(msa, MSACode("Y070_ATV/2-70")) == AccessionCode("Q3V4T1")
+        @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), columnindexes(msa))
-        @test all(iszero, ΔX[i, :])
-        @test !all(iszero, ΔX[i-1, :])
+        for msa in (MSA.read_file(pf09645_sto, Pfam.Stockholm), read(pf09645_sto, BioStockholm.MSA))
+            X = aa_properties_matrix(msa)
+            ΔX = X .- mean(X, dims=2)
+            i = findfirst(==(14), columnindexes(msa))
+            @test all(iszero, ΔX[i, :])
+            @test !all(iszero, ΔX[i-1, :])
+        end
         seqs = FASTAReader(open("test.fasta")) do io
             collect(io)
         end