Merge pull request #58 from m3g/fix_residue_ticks_map2D

Improved 2D density maps: fixed residue ticks for repeated residues, support heatmap and contour

Project.toml

@@ -1,7 +1,7 @@
 name = "ComplexMixtures"
 uuid = "6f35c628-ac57-5bae-8ea9-703a8964f6e9"
 authors = ["Leandro Martinez <[email protected]>"]
-version = "2.10.2"
+version = "2.11.0-DEV"
 
 [deps]
 CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"

docs/src/assets/scripts/example2/script3.jl

@@ -12,11 +12,11 @@ using PDBTools
 
 # Chemical groups of the polymer monomers, defined by the atom types:
 groups = (
-    "CH_3" => ["CF", "HF1", "HF2", "HF3"], # methyles
+    "CH_{3L}" => ["CF", "HF1", "HF2", "HF3"], # left-terminal methyl
     "CO" => ["OE1", "CD"], # carbonyl
     "NH_2" => ["NE2", "HE22", "HE21"], # amine 
     "CHCH_2" => ["C", "H2", "H1", "CA", "HA"], # backbone
-    "CH_3" => ["CL", "HL1", "HL2", "HL3"], # terminal methyles
+    "CH_{3R}" => ["CL", "HL1", "HL2", "HL3"], # right-terminal methyles
 )
 
 system = readPDB("./equilibrated.pdb")
@@ -27,14 +27,12 @@ results = load("./mddf.json")
 # each chemical group of each polymer mer independently
 group_contribs = Vector{Float64}[]
 labels = LaTeXString[]
-for (imer, mer) in enumerate(eachresidue(acr))
+for mer in eachresidue(acr)
     for (group_label, group_atoms) in groups
-        # only first residue has a terminal CH3
-        if imer != 1 && group_label == "CH_3"
+        if resnum(mer) != 1 && group_label == "CH_{3L}" # first residue only
             continue
         end
-        # only last residue has a terminal CH3
-        if imer != 5 && group_label == "CH_3"
+        if resnum(mer) != 5 && group_label == "CH_{3R}" # last residue only
             continue
         end
         # Filter the atoms of this mer that belong to the group

ext/Plotting.jl

@@ -16,9 +16,11 @@ using PDBTools: Residue, residue_ticks, Atom, eachresidue, resnum
         xlabel="Residue",
         ylabel="r / Å",
     )
+    contour(rc::ResidueContributions, args...; kargs...)
+    heatmap(rc::ResidueContributions, args...; kargs...)
 
-Plot the contribution of each residue to the solute-solvent pair distribution function as a contour plot.
-This function requires loading the `Plots` package.
+Plot the contribution of each residue to the solute-solvent pair distribution function as a 2D density map.
+This function requires loading the `Plots` package. The calling syntax for `contour` and `heatmap` is the same as for `contourf`.
 
 # Arguments
 
@@ -27,7 +29,7 @@ This function requires loading the `Plots` package.
 
 # Optional arguments
 
-- `step`: The step of the residue ticks in the x-axis of the plot. Default is 1.
+- `step`: The step of the residue ticks in the x-axis of the plot. Default is 1 or will be set to show at most 20 ticks labels.
 - `oneletter::Bool`: Use one-letter residue codes. Default is `false`. One-letter codes are only available for the 20 standard amino acids.
 - `xlabel` and `ylabel`: Labels for the x and y axes. Default is `"Residue"` and `"r / Å"`.
 
@@ -59,15 +61,22 @@ julia> plt = contourf(rc; step=5, size=(800,400), title="Title", clims=(-0.1, 0.
 ```
 
 !!! compat
-    This function requires loading the `Plots` package and is available in
-    ComplexMixtures v2.5.0 or greater.
+    This function requires loading the `Plots` package.
 
-    Support for all `Plots.contourf` parameters was introduced in ComplexMixtures v2.6.0.
+    Support for all `Plots.contourf` parameters was introduced in ComplexMixtures v2.6.0, and support for
+    `contour` and `heatmap` were introduced in ComplexMixtures v2.11.0. 
 
 """
-function Plots.contourf(
+Plots.contourf, Plots.contour, Plots.heatmap
+
+Plots.contourf(rc::ResidueContributions, args...; kargs...) = _density2D(Plots.contourf, rc, args...; kargs...)
+Plots.contour(rc::ResidueContributions, args...; kargs...) = _density2D(Plots.contour, rc, args...; kargs...)
+Plots.heatmap(rc::ResidueContributions, args...; kargs...) = _density2D(Plots.heatmap, rc, args...; kargs...)
+
+function _density2D(
+    plot_type::Function,
     rc::ResidueContributions;
-    step::Int=1,
+    step::Union{Nothing,Integer}=nothing,
     oneletter=false,
     xlabel="Residue",
     ylabel="r / Å",
@@ -76,9 +85,20 @@ function Plots.contourf(
     kargs...
 )
 
-    # Plot a contour courves with the density at each distance from each residue
+    # Plot a contour curves with the density at each distance from each residue
     # colors, linewidths, etc. are defined here and can be tuned
-    tick_range = 1:step:length(rc.xticks[1])
+    input_step = isnothing(step) ? 1 : step
+    tick_range = firstindex(rc.xticks[1]):input_step:lastindex(rc.xticks[1])
+    nticks = 50
+    if isnothing(step) && length(tick_range) > nticks
+        step = length(rc.resnums) ÷ nticks
+        @warn """\n
+            Consider using a step to set the number of residue ticks in the plot. 
+            - step will be set to $step to display $nticks residue ticks.
+
+        """ _line=nothing _file=nothing
+        tick_range = first(tick_range):step:last(tick_range)
+    end
     tick_marks = rc.xticks[1][tick_range]
     tick_labels = rc.xticks[2][tick_range]
     tick_resnums = rc.resnums[tick_range]
@@ -110,10 +130,20 @@ function Plots.contourf(
         levels = 12
     end
 
+    # density to plot
+    rc_range = 1:length(rc)
+    if length(rc_range) > 2000 
+        rc_step = length(rc_range) ÷ 2000
+        rc_range = 1:rc_step:length(rc)
+        @warn """\n
+            The number of residues to plot is too large. Will plot every $rc_step residues.
+
+        """ _line=nothing _file=nothing
+    end
     Plots.default(fontfamily="Computer Modern")
-    plt = Plots.contourf(
-        rc.resnums,
-        rc.d, hcat(rc.residue_contributions...);
+    plt = plot_type(
+        rc.xticks[1][rc_range],
+        rc.d, hcat(rc[rc_range].residue_contributions...);
         color=Plots.cgrad(colorscale),
         linewidth=1,
         linecolor=:black,

src/tools/residue_contributions.jl

@@ -201,8 +201,6 @@ function ResidueContributions(
     # Obtain pretty labels for the residues in the x-axis (using PDBTools)
     resnums = PDBTools.resnum.(residues)
     xticks = PDBTools.residue_ticks(atoms;
-        first=first(resnums),
-        last=last(resnums),
         oneletter=false,
         serial=false,
     )