diff --git a/src/GeneticsMakie.jl b/src/GeneticsMakie.jl
index 2b9335c..0291138 100644
--- a/src/GeneticsMakie.jl
+++ b/src/GeneticsMakie.jl
@@ -10,14 +10,14 @@ using Distributions
 include("parsegtf.jl")
 include("plotgenes.jl")
 include("plotisoforms.jl")
-include("plotld.jl")
 include("mungesumstats.jl")
 include("plotlocus.jl")
 include("plotloops.jl")
 include("plotqq.jl")
 include("plotgwas.jl")
-include("plotrg.jl")
 include("gwas.jl")
 include("genome.jl")
+include("plotld.jl")
+include("plotrg.jl")
 
-end
+end
\ No newline at end of file
diff --git a/src/plotld.jl b/src/plotld.jl
index a7e9db7..d643912 100644
--- a/src/plotld.jl
+++ b/src/plotld.jl
@@ -1,5 +1,6 @@
 """
     plotld(LD::AbstractMatrix; kwargs)
+    plotld!(ax::Axis, LD::AbstractMatrix; kwargs)
 
 Heatmap of symmetric correlation matrix `LD` with the diagonal elements on the x-axis.
 
diff --git a/src/plotrg.jl b/src/plotrg.jl
index 4f7c0da..b70f280 100644
--- a/src/plotrg.jl
+++ b/src/plotrg.jl
@@ -1,18 +1,39 @@
-function plotrg!(
-    ax::Axis,
-    r::AbstractMatrix,
-    colnames::AbstractVector; 
-    p::Union{AbstractMatrix, Nothing} = nothing,
-    circle::Bool = true,
-    diag::Bool = false,
+"""
+    plotrg(r::AbstractMatrix; kwargs)
+    plotrg!(ax::Axis, r::AbstractMatrix; kwargs)
+
+Correlation plot of matrix `r`.
+
+# Keyword arguments
+- `circle`      : whether to draw cicles instead of rectangles (default to `true`)
+- `diagonal`    : whether to visualize diagonal elements (default to `false`)
+- `colormap`    : colormap of values (default to `:RdBu_10`)
+- `colorrange`  : start and end points of `colormap` (default to `(-1, 1)`)
+- `strokewidth` : width of outline around surrounding boxes (default to `0.5`)
+"""
+@recipe(Plotrg, r) do scene
+    Attributes(
+        circle      = true,
+        diagonal    = false,
+        colormap    = :RdBu_10,
+        colorrange  = (-1, 1),
+        strokewidth = 0.5
     )
+end
 
+function Makie.plot!(plot::Plotrg{<:Tuple{<:AbstractMatrix}})
+    r           = plot[:r][]
+    circle      = plot[:circle][]
+    diagonal    = plot[:diagonal][]
+    colormap    = plot[:colormap][]
+    colorrange  = plot[:colorrange][]
+    strokewidth = plot[:strokewidth][]
     n = size(r, 1)
-    ps = Vector{Polygon}(undef, n^2)
+    polys = Vector{Polygon}(undef, n^2)
     counter = 1
     for i in 1:n
         for j in 1:n
-            ps[counter] = Polygon([
+            polys[counter] = Polygon([
                 Point2f(i - 1, 1 - j),
                 Point2f(i - 1, -j),
                 Point2f(i, -j),
@@ -21,9 +42,9 @@ function plotrg!(
             counter += 1
         end
     end
-    psd = Vector{Polygon}(undef, n)
+    polysd = Vector{Polygon}(undef, n)
     for i in 1:n
-        psd[i] = Polygon([
+        polysd[i] = Polygon([
             Point2f(i - 1, 1 - i),
             Point2f(i - 1, -i),
             Point2f(i, -i),
@@ -45,7 +66,7 @@ function plotrg!(
                 counter += 1
             end
         end
-        if diag
+        if diagonal
             csd = Vector{Circle}(undef, n)
             colord = Vector{Float32}(undef, n)
             for i in 1:n
@@ -70,47 +91,21 @@ function plotrg!(
                 counter += 1
             end
         end
-        if diag
-            csd = Vector{Circle}(undef, n)
+        if diagonal
+            csd = Vector{Rect}(undef, n)
             colord = Vector{Float32}(undef, n)
             for i in 1:n
                 widthd = r[i, i] * 0.95
-                csd[i] = Rect(i - 1 + (1 - withd) / 2, 1 - i - (1 - withd) / 2, withd, -withd)
+                csd[i] = Rect(i - 1 + (1 - widthd) / 2, 1 - i - (1 - widthd) / 2, widthd, -widthd)
                 colord[i] = r[i, i]
             end
         end
     end
-    poly!(ax, ps, color = :white, strokewidth = 0.5, strokecolor = "#C3C3C3")
-    # poly!(ax, psd, color = "#C3C3C3", strokewidth = 0.5, strokecolor = "#C3C3C3")
-    poly!(ax, csl, color = colorl, colorrange = (-1, 1), colormap = :RdBu_10)
-    poly!(ax, csu, color = coloru, colorrange = (-1, 1), colormap = :RdBu_10)
-    if !isnothing(p)
-        pvalues = Matrix{Float64}(undef, n^2, 3)
-        counter = 1
-        for i in 1:n
-            for j in 1:n
-                j == i ? pvalues[counter, 1] = 1 : pvalues[counter, 1] = p[j, i]
-                pvalues[counter, 2] = i - 0.5
-                pvalues[counter, 3] = -j + 0.5
-                counter += 1
-            end
-        end
-        ind = findall(x -> x < 0.05, pvalues[:, 1]) 
-        scatter!(ax, pvalues[ind, 2], pvalues[ind, 3], marker = '✳', markersize = 7, color = :black)
-    end
-    if diag
-        poly!(ax, csd, color = colord, colorrange = (-1, 1), colormap = :RdBu_10)
+    poly!(plot, polys, color = :white, strokewidth = strokewidth, strokecolor = "#C3C3C3")
+    # poly!(plot, polysd, color = "#C3C3C3", strokewidth = strokewidth, strokecolor = "#C3C3C3")
+    poly!(plot, csl, color = colorl, colorrange = colorrange, colormap = colormap)
+    poly!(plot, csu, color = coloru, colorrange = colorrange, colormap = colormap)
+    if diagonal
+        poly!(plot, csd, color = colord, colorrange = colorrange, colormap = colormap)
     end
-    ax.xticks = (0.5:(n - 0.5), colnames)
-    ax.yticks = (-0.5:-1:-(n - 0.5), colnames)
-    ax.xticklabelsize = 6
-    ax.yticklabelsize = 6
-    ax.xticklabelpad = 0
-    ax.yticklabelpad = 1
-    ax.xaxisposition = :top
-    ax.yaxisposition = :left
-    hidedecorations!(ax, ticklabels = false)
-    xlims!(ax, 0, n)
-    ylims!(ax, -n, 0)
-    ax.aspect = DataAspect()
 end
\ No newline at end of file
diff --git a/src/plotwas.jl b/src/plotwas.jl
deleted file mode 100644
index 4b7bfbc..0000000
--- a/src/plotwas.jl
+++ /dev/null
@@ -1,295 +0,0 @@
-function plotwas!(
-    ax::Axis,
-    twas::DataFrame,
-    gencode::DataFrame; 
-    ymax::Real = 0,
-    sigline::Bool = false,
-    sigcolor::Bool = true,
-    zscore::Bool = false,
-    gene_id::Bool = false
-)
-    gencodeₛ = filter(x -> x.feature == "gene", gencode)
-    gencodeₛ.x = (gencodeₛ.start + gencodeₛ.end) / 2
-    if gene_id
-        select!(gencodeₛ, :gene_id, :x, :CHR)
-        df = innerjoin(twas, gencodeₛ; on = :gene_id)
-    else
-        select!(gencodeₛ, :gene_name, :x, :CHR)
-        df = innerjoin(twas, gencodeₛ; on = :gene_name)
-    end
-    df = select(gwas, :CHR, :BP, :P)
-    df.P = -log.(10, df.P)
-    if ymax == 0
-        ymax = maximum(df.P) / 4 * 5
-        ymax <= 10 ? ymax = 10 : nothing
-    end
-    storage = DataFrame(CHR = vcat(string.(1:22), ["X", "Y"]))
-    storage.maxpos = [GRCh37_totlength[chr] for chr in storage.CHR]
-    storage.add = cumsum(storage.maxpos) - storage.maxpos
-    df = leftjoin(df, storage; on = :CHR)
-    df.x = df.BP + df.add
-    xmax = sum(unique(df.maxpos))
-    indeven = findall(in(vcat(string.(1:2:22), "X")), df.CHR)
-    indodd = findall(in(vcat(string.(2:2:23), "Y")), df.CHR)
-    scatter!(ax, view(df, indeven, :x), view(df, indeven, :P), markersize = 1.5, color = "#0D0D66")
-    scatter!(ax, view(df, indodd, :x), view(df, indodd, :P), markersize = 1.5, color = "#7592C8")
-    if sigcolor
-        ind = df.P .> -log(10, 5e-8)
-        dfsig = view(df, ind, :)
-        scatter!(ax, dfsig.x, dfsig.P,  markersize = 1.5, color = "#4DB069")
-    end
-    if sigline
-        hlines!(ax, -log(10, 5e-8); xmin = 0.0, xmax = xmax, linewidth = 0.75, color = :red2)
-    end
-    xlims!(ax, 0, xmax)
-    ylims!(ax, 0, ymax)
-    hideydecorations!(ax, label = false, ticklabels = false, ticks = false)
-    hidexdecorations!(ax, label = false, ticklabels = false)
-    ax.xlabel = "Chromosome"
-    ax.ylabel = "-log[p]"
-    ax.xticks = ((cumsum(storage.maxpos) + storage.add) / 2, storage.CHR)
-    ax.yticks = setticks(ymax)
-    ax.xticklabelsize = 6
-    ax.yticklabelsize = 6
-    ax.xlabelsize = 8
-    ax.ylabelsize = 8
-    ax.xticksize = 3
-    ax.yticksize = 3
-    ax.xtickwidth = 0.75
-    ax.ytickwidth = 0.75
-    ax.spinewidth = 0.75
-    ax.xticklabelpad = 0
-    ax.yticklabelpad = 2.5
-end
-
-function plotwas!(
-    ax::Axis,
-    chromosome::AbstractString,
-    range1::Real,
-    range2::Real,
-    twas::DataFrame,
-    gencode::DataFrame
-)
-
-end
-
-plotwas!(ax::Axis, chromosome::AbstractString, bp::Real, twas::DataFrame, gencode::DataFrame; window::Real = 1e6) =
-    plotwas!(ax, chromosome, bp - window, bp + window, twas, gencode)
-
-function plotwas!(ax::Axis, gene::AbstractString, twas::DataFrame, gwas::DataFrame; window::Real = 1e6)
-    chr, start, stop = findgene(gene, gencode)
-    plotwas!(ax, chr, start - window, stop + window, twas, gencode)
-end
-
-
-# g <- ggplot() + geom_point(data=df, aes(x=pos, y=TWAS.Z, col=(CHR %% 2 == 0))) + 
-#   geom_point(data=df[df$FOCUS==TRUE,], aes(x=pos, y=TWAS.Z, col="3")) + 
-#   scale_color_manual(values=c("FALSE"="gray60","TRUE"="gray70","3"="#E41A1C")) +
-#   scale_x_continuous("Chromosome", breaks=floor(0.5*(ticks+ticks2)), labels=1:22) +
-#   geom_text_repel(data=df[df$FOCUS==TRUE,], aes(x=pos, y=TWAS.Z, label=ID), 
-#                   color="black", size=3, segment.size=.1) + 
-#   theme_bw() + theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank(), 
-#                      panel.grid.major.y=element_blank(), panel.grid.minor.y=element_blank(), legend.position="none") + 
-#   ylim(-10,10) + ylab("TWAS Z-score") + geom_hline(yintercept=4.63, linetype="dashed", size=0.375, color="#377EB8") + 
-#   geom_hline(yintercept=-4.63, linetype="dashed", size=0.375, color="#377EB8") + 
-#   geom_vline(xintercept=ticks2-1, linetype="dotted", size=0.375, color="black")
-
-function coordinateisforms(
-    gene::AbstractString, 
-    gencode::DataFrame,
-    orderby::Union{Nothing, AbstractVector{<:AbstractString}},
-    height::Real,
-    text::Union{Bool, Symbol},
-    shrink::Real
-    )
-
-    df = filter(x -> x.gene_name == gene, gencode)
-    nrow(df) == 0 ? error("Cannot find $(gene) in the annotation.") : nothing
-    dfi = view(df, df.feature .== "transcript", :)
-    dfe = view(df, df.feature .== "exon", :)
-    chromosome = df.seqnames[1]
-    isoforms = unique(dfi.transcript_id)
-    n = length(isoforms)
-    ps = Vector{Vector{Polygon}}(undef, n) 
-    bs = Matrix{Float64}(undef, n, 2)
-    rows = collect(1:n)
-    if shrink != 1
-        for (exon1, exon2) in zip(storage.start[2:end], storage.end[1:end])
-            (exon1 - exon2) * shrink
-        end
-    end
-    if !isnothing(orderby)
-        prior = []
-        for k in 1:length(orderby)
-            if any(isoforms .== orderby[k])
-                push!(prior, orderby[k])
-            end
-        end
-        if length(prior) > 0
-            dforder = DataFrame(transcript_id = [prior; filter(!in(orderby), isoforms)], rank = 1:n)
-            dfi = leftjoin(dfi, dforder; on = :transcript_id)
-            sort!(dfi, [order(:rank)])
-            isoforms = unique(dfi.transcript_id)
-        end
-    end
-    for j in eachindex(isoforms)
-        ranges = view(dfe, findall(isequal(isoforms[j]), dfe.transcript_id), [:start, :end])
-        m = size(ranges, 1)
-        p = Vector{Polygon}(undef, m)
-        if text == true || text == :top || text == :t || text == :bottom || text == :b
-            for i = 1:m
-                p[i] = Polygon(
-                    [Point2f(ranges[i, 1], 1 - height - (rows[j] - 1) * (0.25 + height)),
-                    Point2f(ranges[i, 1], 1 - (rows[j] - 1) * (0.25 + height)),
-                    Point2f(ranges[i, 2], 1 - (rows[j] - 1) * (0.25 + height)),
-                    Point2f(ranges[i, 2], 1 - height - (rows[j] - 1) * (0.25 + height))]
-                )
-            end
-        else
-            for i = 1:m
-                p[i] = Polygon(
-                    [Point2f(ranges[i, 1], 1 - height - (rows[j] - 1) * (0.025 + height)),
-                    Point2f(ranges[i, 1], 1 - (rows[j] - 1) * (0.025 + height)),
-                    Point2f(ranges[i, 2], 1 - (rows[j] - 1) * (0.025 + height)),
-                    Point2f(ranges[i, 2], 1 - height - (rows[j] - 1) * (0.025 + height))]
-                )
-            end
-        end
-        ps[j] = p
-        ind = findfirst(isequal(isoforms[j]), dfi.transcript_id)
-        start = dfi.start[ind]
-        stop = dfi[ind, :end]
-        bs[j, 1] = start
-        bs[j, 2] = stop
-    end
-    return isoforms, ps, bs, rows, chromosome
-end
-
-"""
-    plotisoforms!(ax::Axis, gene::AbstractString, gencode::DataFrame; kwargs)
-
-Plot each isoform of a given `gene` on a separate row.
-
-# Arguments
-- `orderby::Union{Nothing, AbstractVector{<:AbstractString}} = nothing`: the order of isoforms.
-- `highlight::Union{Nothing, Tuple{AbstractVector, AbstractVector}} = nothing`: isoforms to be highlighted and their colors.
-- `height::Real = 0.25`: the height of exons.
-- `isoformcolor = :royalblue`: the color of isoforms.
-- `textcolor = :black`: the color of isoform labels.
-- `text::Union{Bool, Symbol} = :top`: the position of isoform labels. 
-"""
-function plotisoforms!(
-    ax::Axis,
-    gene::AbstractString,
-    gencode::DataFrame;
-    orderby::Union{Nothing, AbstractVector{<:AbstractString}} = nothing,
-    highlight::Union{Nothing, Tuple{AbstractVector, AbstractVector}} = nothing, 
-    height::Real = 0.25,
-    isoformcolor = :royalblue,
-    textcolor = :black,
-    text::Union{Bool, Symbol} = :top,
-    shrink::Real = 1
-    )
-
-    isoforms, ps, bs, rows, chromosome = coordinateisforms(gene, gencode, orderby, height, text, shrink)
-    isnothing(highlight) ? highlight = ([nothing], [nothing]) : nothing
-    if text == true || text == :top || text == :t
-        for j in 1:size(ps, 1)
-            ind = findfirst(isequal(isoforms[j]), highlight[1])
-            if isnothing(ind)
-                poly!(ax, ps[j], color = isoformcolor, strokewidth = 0)
-                lines!(ax, [bs[j, 1], bs[j, 2]], 
-                    [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
-                    color = isoformcolor, linewidth = 0.5)
-                text!(ax, "$(isoforms[j])", 
-                    position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - (rows[j] - 1) * (0.25 + height)), 
-                    align = (:center, :bottom), fontsize = 6, color = textcolor)
-            else
-                poly!(ax, ps[j], color = highlight[2][ind], strokewidth = 0)
-                lines!(ax, [bs[j, 1], bs[j, 2]], 
-                    [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
-                    color = highlight[2][ind], linewidth = 0.5)
-                text!(ax, "$(isoforms[j])", 
-                    position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - (rows[j] - 1) * (0.25 + height)), 
-                    align = (:center, :bottom), fontsize = 6, color = textcolor)
-            end
-        end
-    elseif text == :bottom || text == :b
-        for j in 1:size(ps, 1)
-            ind = findfirst(isequal(isoforms[j]), highlight[1])
-            if isnothing(ind)
-                poly!(ax, ps[j], color = isoformcolor, strokewidth = 0)
-                lines!(ax, [bs[j, 1], bs[j, 2]], 
-                    [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
-                    color = isoformcolor, linewidth = 0.5)
-                text!(ax, "$(isoforms[j])", 
-                    position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - height - (rows[j] - 1) * (0.25 + height)), 
-                    align = (:center, :top), fontsize = 6, color = textcolor)
-            else
-                poly!(ax, ps[j], color = highlight[2][ind], strokewidth = 0)
-                lines!(ax, [bs[j, 1], bs[j, 2]], 
-                    [1 - height / 2 - (rows[j] - 1) * (0.25 + height), 1 - height / 2 - (rows[j] - 1) * (0.25 + height)],
-                    color = highlight[2][ind], linewidth = 0.5)
-                text!(ax, "$(isoforms[j])", 
-                    position = ((bs[j, 1] + bs[j, 2]) / 2, 1 - height - (rows[j] - 1) * (0.25 + height)), 
-                    align = (:center, :top), fontsize = 6, color = textcolor)
-            end
-        end
-    else
-        for j in 1:size(ps, 1)
-            ind = findfirst(isequal(isoforms[j]), highlight[1])
-            if isnothing(ind)
-                poly!(ax, ps[j], color = isoformcolor, strokewidth = 0)
-                lines!(ax, [bs[j, 1], bs[j, 2]], 
-                    [1 - height / 2 - (rows[j] - 1) * (0.025 + height), 1 - height / 2 - (rows[j] - 1) * (0.025 + height)],
-                    color = isoformcolor, linewidth = 0.5)
-            else
-                poly!(ax, ps[j], color = highlight[2][ind], strokewidth = 0)
-                lines!(ax, [bs[j, 1], bs[j, 2]], 
-                    [1 - height / 2 - (rows[j] - 1) * (0.025 + height), 1 - height / 2 - (rows[j] - 1) * (0.025 + height)],
-                    color = highlight[2][ind], linewidth = 0.5)
-            end
-        end
-    end
-    range1 = minimum(bs[:, 1])
-    range2 = maximum(bs[:, 2])
-    diff = range2 - range1
-    prop = 23000 * diff / 2.2e6
-    if text == true || text == :top || text == :t || text == :bottom || text == :b
-        storage = bs[:, 1] + bs[:, 2]
-        ind = argmin(storage)
-        low = storage[ind] / 2 - prop * 12
-        low < range1 ? range1 = low : range1 -= prop * 4.5
-        ind = argmax(storage)
-        high = storage[ind] / 2 + prop * 12
-        high > range2 ? range2 = high : range2 += prop * 4.5
-    else 
-        range1 = range1 - diff / 50
-        range2 = range2 + diff / 50
-    end
-    ax.spinewidth = 0.75
-    hidexdecorations!(ax)
-    if text == :left || text == :l || text == :right || text == :r
-        ax.yticks = ([1 - height / 2 - (rows[j] - 1) * (0.025 + height) for j in 1:length(isoforms)], isoforms)
-        hideydecorations!(ax, ticklabels = false)
-        ax.yticklabelsize = 4
-        ax.yticklabelpad = 3
-        (text == :right) || (text == :r) ? ax.yaxisposition = :right : ax.yaxisposition = :left
-        rs = 18 * maximum(rows) * (0.025 + height) / 0.5
-    elseif text == true || text == :top || text == :t || text == :bottom || text == :b
-        hideydecorations!(ax)
-        rs = 18 * maximum(rows) * (0.25 + height) / 0.5
-    else
-        hideydecorations!(ax)
-        rs = 18 * maximum(rows) * (0.025 + height) / 0.5
-    end
-    xlims!(ax, range1, range2)
-    if text == true || text == :top || text == :t
-        ylims!(ax, 0.875 - height  - (maximum(rows) - 1) * (0.25 + height), 1.375)
-    elseif text == :bottom || text == :b
-        ylims!(ax, 0.65 - height  - (maximum(rows) - 1) * (0.25 + height), 1.05)
-    else
-        ylims!(ax, 0.975 - height  - (maximum(rows) - 1) * (0.025 + height), 1.05)
-    end
-    return rs, chromosome, range1, range2
-end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 98498d4..479f27f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ using CSV
 using DataFrames
 using SnpArrays
 using Statistics
+using Distributions
 
 @testset "Plotting genes/isoforms" begin
     gencode = CSV.read("data/gencode.gtf", DataFrame)
@@ -162,18 +163,6 @@ end
     rm("manhattan.png")
 end
 
-@testset "Plotting LD" begin
-    LD = rand(10, 10)
-    LD = LD + transpose(LD)
-    [LD[i, i] = 1 for i in 1:size(LD, 1)]
-    f = Figure()
-    ax = Axis(f[1, 1])
-    GeneticsMakie.plotld!(ax, LD)
-    save("ld.png", f)
-    @test isfile("ld.png")
-    rm("ld.png")
-end
-
 @testset "Plotting LocusZoom" begin
     kgp = SnpData("data/kgp")
     gwas = CSV.read("data/locus.csv", DataFrame)
@@ -233,21 +222,6 @@ end
     rm("qq.png")
 end
 
-@testset "Plotting correlation" begin
-    f = Figure()
-    ax = Axis(f[1, 1])
-    n = 10
-    GeneticsMakie.plotrg!(ax, (rand(n, n) .* 2 .- 1), string.(1:n), circle = true)
-    colsize!(f.layout, 1, Aspect(1, 1))
-    rowsize!(f.layout, 1, 18 * n)
-    save("cor.png", f)
-    @test isfile("cor.png")
-    rm("cor.png")
-end
-
-@testset "Plotting TWAS" begin
-end
-
 @testset "Plotting loops" begin
     loopdf = CSV.read("data/loops.csv", DataFrame)
     loopdf.chr1 = string.(loopdf.chr1)
@@ -306,3 +280,31 @@ end
     @test isfile("loops.png")
     rm("loops.png")
 end
+
+@testset "Plotting LD" begin
+    LD = rand(50, 50)
+    LD = (LD + transpose(LD)) / 2
+    [LD[i, i] = 1 for i in 1:size(LD, 1)]
+    f, ax, p = plotld(LD)
+    ax.aspect = DataAspect()
+    xlims!(0, 10)
+    ylims!(0, 5)
+    hidedecorations!(ax)
+    hidespines!(ax, :r, :l, :b)
+    Colorbar(f[1, 2], p, label = "LD", height = Relative(0.2))
+    save("ld.png", f)
+    @test isfile("ld.png")
+    rm("ld.png")
+end
+
+@testset "Plotting correlation" begin
+    r = rand(Uniform(-1, 1), 10, 10)
+    f, ax, p = plotrg(r)
+    ax.aspect = DataAspect()
+    xlims!(0, size(r, 1))
+    ylims!(-size(r, 1), 0)
+    hidedecorations!(ax)
+    save("corr.png", f)
+    @test isfile("corr.png")
+    rm("corr.png")
+end
\ No newline at end of file