Modernize chimera script generation (#18)

Old hacky behavior is deprecated in favor of something saner.
Also leverages the new latest-version alphafold file machinery.

src/GPCRAnalysis.jl

@@ -54,5 +54,6 @@ include("properties.jl")
 include("chimerax.jl")
 include("pocket.jl")
 include("features.jl")
+include("deprecated.jl")
 
 end

src/chimerax.jl

@@ -1,5 +1,5 @@
 function chimerax_script(scriptfilename, struct_filenames, ridxs::AbstractVector{<:AbstractVector{<:Integer}};
-                         chain_transparency=80, styles=Dict{Tuple{Int,Int},String}(), extras=String[])
+                         align::Bool=true, chain_transparency=80, styles=Dict{Tuple{Int,Int},String}(), extras=String[])
     open(scriptfilename, "w") do io
         for (i, afname) in enumerate(struct_filenames)
             println(io, "open ", afname)
@@ -13,29 +13,45 @@ function chimerax_script(scriptfilename, struct_filenames, ridxs::AbstractVector
             end
         end
         println(io, "transparency #1-", length(struct_filenames), ' ', chain_transparency, " target c")
-        if (length(struct_filenames) > 1 &&
-            !any(str -> startswith(str, "align"), struct_filenames) &&
-            !any(str -> startswith(str, "align"), extras))
-            println(io, "matchmaker #2-", length(struct_filenames), " to #1")
+        if (length(struct_filenames) > 1 && align)
+            if (any(str -> startswith(str, "align"), struct_filenames) ||
+                any(str -> startswith(str, "align"), extras))
+                @warn "skipping alignment based on file names is deprecated, use `align=false` instead" maxlog=1
+            else
+                println(io, "matchmaker #2-", length(struct_filenames), " to #1")
+            end
         end
-        for ex in extras
-            println(io, ex)
+        if isa(extras, AbstractString)
+            println(io, extras) 
+        else
+            for ex in extras
+                println(io, ex)
+            end
         end
     end
 end
 
 """
     chimerax_script(scriptfilename, uprot_list, msa::AnnotatedMultipleSequenceAlignment, colidxs;
-                    dir=pwd(), chain_transparency=80, styles=Dict{Int,String}(), extras=String[])
+                    dir=pwd(), align=true, chain_transparency=80, styles=Dict{Int,String}(), extras=String[])
 
-Create a [chimerax]() visualization script with name `scriptfilename`. `uprot_list` is a list of UniProtX names that you
-want to visualize.  The AlphaFold PDB files for these proteins should be `dir`. `msa` is a Multiple Sequence alignment
-and `colidxs` specifies the column indices in `msa` corresponding to amino acid side chains that you'd like to visualize.
+Create a [chimerax](https://www.cgl.ucsf.edu/chimerax/) visualization script
+with name `scriptfilename`. `uprot_list` is a list of UniProtX names that you
+want to visualize. `msa` is a Multiple Sequence alignment and `colidxs`
+specifies the column indices in `msa` corresponding to amino acid side chains
+that you'd like to visualize.
 
-`chain_transparency` sets the transparency on the ribbon diagrams (0 = not transparent).
-`styles` can be used to affect the display, e.g., `Dict(k => "@SD sphere")` would cause methionines at column index `k`
-to be displayed with the sulfur in sphere mode. `extras` can be used to hand-specify a number of additional commands;
-this can be useful if, for example, the `msa` has occasional misalignments.
+Keyword arguments:
+- `dir` is the directory with the protein structure files
+- `align` determines whether to align the structures to the first one (uses the
+  `matchmaker` tool)
+- `chain_transparency` sets the transparency on the ribbon diagrams (0 = not
+transparent)
+- `styles` can be used to affect the display, e.g., `Dict(k => "@SD sphere")`
+would cause methionines at column index `k` to be displayed with the sulfur in
+sphere mode.
+- `extras` can be used to hand-specify a number of additional commands; this can
+be useful if, for example, the `msa` has occasional misalignments.
 
 # Examples
 
@@ -52,9 +68,11 @@ function chimerax_script(scriptfilename, uprot_list, msa::AnnotatedMultipleSeque
     ridxs = [Int[] for _ in 1:length(uprot_list)]
     struct_filenames = Vector{String}(undef, length(uprot_list))
     rcstyles = Dict{Tuple{Int,Int},String}()
+    afs = alphafoldfile(msa, dir; join=true)
+    uprot2msaidx = Dict{AccessionCode,Int}(AccessionCode(msa, name) => i for (i, name) in enumerate(sequencenames(msa)))
     for (i, p) in enumerate(uprot_list)
-        struct_filenames[i] = joinpath(dir, alphafoldfilename(p))
-        seqidx = findfirst(str -> startswith(str, p), sequencenames(msa))
+        j = uprot2msaidx[AccessionCode(p)]
+        struct_filenames[i] = afs[MSACode(sequencenames(msa, j))]
         sm = getsequencemapping(msa, seqidx)
         for (j, c) in enumerate(colidxs)
             ridx = sm[c]
@@ -72,6 +90,12 @@ function chimerax_script(scriptfilename, uprot_list, msa::AnnotatedMultipleSeque
     return chimerax_script(scriptfilename, struct_filenames, ridxs; styles=rcstyles, kwargs...)
 end
 
-function markers(modelnum::Integer, positions, radius::Real, color)
-    return ["marker #$modelnum position $(pos[1]),$(pos[2]),$(pos[3]) radius $radius color $color" for pos in positions]
+"""
+    str = marker(modelnum, pos, radius, color)
+
+Return a string that represents a marker in ChimeraX. `modelnum` is the model number, `pos` is a 3D position,
+`radius` is the radius of the marker, and `color` is the color of the marker.
+"""
+@noinline function marker(modelnum::Integer, pos::AbstractVector{<:Real}, radius::Real, color)
+    return "marker #$modelnum position $(pos[1]),$(pos[2]),$(pos[3]) radius $radius color $color"
 end

src/deprecated.jl

@@ -0,0 +1,4 @@
+function markers(modelnum::Integer, positions::AbstractVector{<:AbstractVector{<:Real}}, radius::Real, color)
+    Base.depwarn("`markers` is deprecated, use `marker.(modelnum, positions, radius, color)` instead", :markers)
+    return [marker(modelnum, pos, radius, color) for pos in positions]
+end

test/deprecated.jl

@@ -0,0 +1,11 @@
+@testset "Deprecated" begin
+    @info "Testing deprecated functions (warnings are expected)"
+    @test only(GPCRAnalysis.markers(2, [Coordinates(5, 4, 3)], 0.5, "blue")) == "marker #2 position 5.0,4.0,3.0 radius 0.5 color blue"
+    dot = GPCRAnalysis.marker(2, Coordinates(5, 4, 3), 0.5, "blue")
+    tmpfile = tempname() * ".cxc"
+    chimerax_script(tmpfile, ["align_ABCD.pdb", "align_EFGH.pdb"], [[5, 10, 15], [7, 11]]; extras=[dot])
+    script = read(tmpfile, String)
+    @test !occursin("matchmaker #2-2 to #1", script)
+    @test occursin("marker #2 position 5.0,4.0,3.0 radius 0.5 color blue", script)
+    rm(tmpfile)
+end

test/runtests.jl

@@ -86,11 +86,12 @@ using Test
         @test occursin("show #2 :11", script)
         @test occursin("transparency #1-2 80 target c", script)
         @test occursin("matchmaker #2-2 to #1", script)
-        dots = GPCRAnalysis.markers(2, [Coordinates(5, 4, 3)], 0.5, "blue")
-        chimerax_script(tmpfile, ["align_ABCD.pdb", "align_EFGH.pdb"], [[5, 10, 15], [7, 11]]; extras=dots)
+        dot = GPCRAnalysis.marker(2, Coordinates(5, 4, 3), 0.5, "blue")
+        chimerax_script(tmpfile, ["ABCD.pdb", "EFGH.pdb"], [[5, 10, 15], [7, 11]]; align=false, extras=[dot])
         script = read(tmpfile, String)
         @test !occursin("matchmaker #2-2 to #1", script)
         @test occursin("marker #2 position 5.0,4.0,3.0 radius 0.5 color blue", script)
+        rm(tmpfile)
     end
     @testset "Needleman-Wunsch" begin
         function testscore(S, P, D, gapcosts)
@@ -348,3 +349,5 @@ using Test
         end
     end
 end
+
+include("deprecated.jl")