Add BioStockholm extension

This supports MSAs loaded by BioStockholm.

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,99 @@
+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.sequencekeys(msa::MSA) = collect(keys(msa.seq))
+GPCRAnalysis.msasequence(msa::MSA, key) = msa.seq[key][conscols(msa)]
+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
+GPCRAnalysis.columnindexes(msa::BioStockholm.MSA) = conscols(msa)
+
+Base.getindex(msa::MSA, seqname::MSACode) = msa.seq[seqname.name][conscols(msa)]
+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 == '-'
+        return 0
+    elseif 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(f, msa::MSA)
+    # This mimics MIToS's implementation
+    function pctsim(v1, v2)
+        same = l = 0
+        for (a, b) in zip(v1, v2)
+            a == b == 0 && continue  # skip gaps
+            same += a == b
+            l += 1
+        end
+        return 100 * same / l
+    end
+
+    M = f.(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] = pctsim(M[i, :], M[j, :])
+            S[j, i] = S[i, j]
+        end
+    end
+    return S
+end
+GPCRAnalysis.percent_similarity(msa::MSA) = GPCRAnalysis.percent_similarity(reduced_alphabet, msa)
+
+
+end

ext/GPCRAnalysisMIToSExt.jl

@@ -11,7 +11,7 @@ using MIToS: MIToS, Pfam, MSA
 using MIToS.MSA: AbstractMultipleSequenceAlignment, AnnotatedAlignedSequence, AnnotatedMultipleSequenceAlignment,
                  ReducedAlphabet, ResidueAlphabet, GAP, XAA
 using MIToS.MSA: getsequence, getannotsequence, getsequencemapping, getresidues, three2residue, sequencenames,
-                 filtersequences, filtersequences!, percentsimilarity
+                 filtersequences, filtersequences!, percentsimilarity, getcolumnmapping
 
 
 # Low-level API implementation
@@ -25,6 +25,7 @@ GPCRAnalysis.residuematrix(msa::AbstractMultipleSequenceAlignment) = getresidues
 GPCRAnalysis.subseqs(msa::AbstractMultipleSequenceAlignment, rowmask)  = filtersequences(msa, rowmask)
 GPCRAnalysis.subseqs!(msa::AbstractMultipleSequenceAlignment, rowmask) = filtersequences!(msa, rowmask)
 GPCRAnalysis.percent_similarity(msa::AbstractMultipleSequenceAlignment) = percentsimilarity(msa)
+GPCRAnalysis.columnindexes(msa::MSA.AbstractMultipleSequenceAlignment) = getcolumnmapping(msa)
 
 Base.getindex(msa::AbstractMultipleSequenceAlignment, seqname::MSACode) = getsequence(msa, seqname.name)
 Base.getindex(msa::AbstractMultipleSequenceAlignment, seqname::AccessionCode) = getsequence(msa, MSACode(msa, seqname).name)

src/msa.jl

@@ -14,6 +14,13 @@ The two-argument form retrieves the sequenceindexes for the `i`th sequence in `m
 """
 function sequenceindexes end
 
+"""
+    idxs = columnindexes(msa)
+
+Return the indices (within the reference sequence) covered by the conserved columns of the MSA.
+"""
+function columnindexes end
+
 """
     isgap(res)
 
@@ -108,11 +115,11 @@ Remove all sequences from `msa` with fewer than `minres` matching residues.
 """
 function filter_long!(msa, minres::Real)
     # Get rid of short sequences
-    nresidues = map(eachrow(msa)) do v
-        sum(!isgap, v)
+    nresidues = map(sequencekeys(msa)) do key
+        sum(!isgap, msasequence(msa, key))
     end
     rowmask = nresidues .> minres
-    subseqs(msa, rowmask)
+    subseqs!(msa, rowmask)
 end
 
 struct SequenceMapping <: AbstractVector{Int}

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

@@ -1,11 +1,11 @@
 using GPCRAnalysis
 using GPCRAnalysis: three2char
 # MSA interface functions
-using GPCRAnalysis: sequenceindexes, isgap, isunknown, sequencekeys, msasequence, residuematrix,
+using GPCRAnalysis: sequenceindexes, columnindexes, isgap, isunknown, sequencekeys, msasequence, residuematrix,
                      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
@@ -16,8 +16,6 @@ using JuMP, HiGHS
 using ColorTypes
 using Test
 
-columnindexes(msa::MSA.AbstractMultipleSequenceAlignment) = MSA.getcolumnmapping(msa)
-
 # skip the network-hitting components by setting `skip_download = true` in the global namespace
 
 @testset "GPCRAnalysis.jl" begin
@@ -66,6 +64,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 +90,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(residuematrix(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] == '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