add docs for plotqq and twas

docs/src/examples/gwas.md

@@ -69,4 +69,30 @@ rowgap!(f.layout, 10)
 resize_to_layout!(f)
 f
 ```
-![](../figs/manhattan-chromosomes.png)
+![](../figs/manhattan-chromosomes.png)
+
+We can then use `GeneticsMakie.plotqq!` to draw QQ plots.
+```julia
+f = Figure(resolution = (408, 792))
+axs = [Axis(f[2, i]) for i in 1:length(titles)]
+for i in eachindex(titles)
+    GeneticsMakie.plotqq!(axs[i], dfs[i]; ystep = 5)
+    axs[i].xlabel = ""
+    axs[i].ylabel = ""
+    ylims!(axs[i], 0, 40)
+    i > 1 ? hideydecorations!(axs[i]) : nothing
+end
+for (i, title) in enumerate(titles)
+    Box(f[1, i], color = :gray90)
+    Label(f[1, i], title, tellwidth = false, textsize = 8, padding = (0, 0, 3, 3))
+end
+Label(f[3, 1:length(titles)], text = "Expected -log[p]", textsize = 8)
+Label(f[2, 0], text = "Observed -log[p]", textsize = 8, rotation = pi / 2, tellheight = false)
+rowsize!(f.layout, 2, Aspect(2, 1))
+colgap!(f.layout, 5)
+rowgap!(f.layout, 1, 0)
+rowgap!(f.layout, 2, 5)
+resize_to_layout!(f)
+f
+```
+![](../figs/qq.png)

docs/src/examples/twas.md

@@ -0,0 +1,48 @@
+# Plotting TWAS
+
+One can also visualize gene-level association results using `GeneticsMakie.plotgwas!`.
+Here we focus on the results of meta analysis of gene-level burden test from case and control studies
+and gene-level Poisson test from family trio studies for schizophrenia.
+
+```julia
+using Pkg
+Pkg.add(["GeneticsMakie", "CairoMakie", "CSV", "DataFrames"])
+```
+
+```julia
+using GeneticsMakie, CairoMakie, CSV, DataFrames
+url = "https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_39/GRCh37_mapping/gencode.v39lift37.annotation.gtf.gz"
+gencode = Arrow.Table("data/gencode/$(splitext(basename(url))[1]).arrow")|> DataFrame
+Downloads.download("https://atgu-exome-browser-data.s3.amazonaws.com/SCHEMA/SCHEMA_gene_results.tsv.bgz", "data/gwas/SCHEMA_gene_results.tsv.bgz")
+df = CSV.read("data/gwas/SCHEMA_gene_results.tsv.bgz", DataFrame)
+```
+
+The key is to calculate coordinates for each gene in the summary statistics.
+
+```julia
+df.CHR .= "1"
+df.BP .= -1.5
+storage = filter(x -> x.feature == "gene", gencode)
+for i in 1:nrow(df)
+    ind = findfirst(isequal(df.gene_id[i]), storage.gene_id)
+    isnothing(ind) ? continue : nothing
+    df.CHR[i] = storage.seqnames[ind]
+    df.BP[i] = (storage.start[ind] + storage.end[ind]) / 2
+end
+rename!(df, "P meta" => "P")
+dropmissing!(df, :P)
+filter!(x -> x.P != "NA", df)
+df.P = parse.(Float64, df.P)
+```
+
+```julia
+f = Figure(resolution = (408, 792))
+ax = Axis(f[1, 1])
+GeneticsMakie.plotgwas!(ax, df; ymax = 13, p = 2.2e-6)
+hidespines!(ax, :t, :r)
+Label(f[1, 1, Top()], text = "SCZ (2022): SCHEMA", textsize = 8)
+rowsize!(f.layout, 1, 50)
+resize_to_layout!(f)
+f
+```
+![](../figs/schema.png)

src/plotgwas.jl

@@ -48,7 +48,7 @@ function plotgwas!(
         scatter!(ax, dfsig.x, dfsig.P, markersize = 1.5, color = scattercolor)
     end
     if !isnothing(linecolor)
-        hlines!(ax, -log(10, p); xmin = 0.0, xmax = xmax, linewidth = 0.75, color = linecolor)
+        lines!(ax, [0.0, xmax], fill(-log(10, p), 2), color = linecolor, linewidth = 0.75)
     end
     xlims!(ax, 0, xmax)
     ylims!(ax, 0, ymax)

src/plotqq.jl

@@ -1,3 +1,11 @@
+"""
+    plotqq!(ax::Axis, gwas::DataFrame; kwargs)
+    plotqq!(ax::Axis, P::AbstractVector; kwargs)
+
+Plot QQ plot of `P` values where the expected distribution is the uniform distribution.
+
+Keyword arguments include `xstep::Real` and `ystep::Real` for x and y axes ticks step sizes.
+"""
 function plotqq!(
     ax::Axis,
     P::AbstractVector; 
@@ -45,12 +53,4 @@ function plotqq!(
     ax.yticks = yticks
 end
 
-"""
-    plotqq!(ax::Axis, gwas::DataFrame; kwargs)
-    plotqq!(ax::Axis, P::AbstractVector; kwargs)
-
-Plot QQ plot of `P` values where the expected distribution is the uniform distribution.
-
-Keyword arguments include `xstep::Real` and `ystep::Real` for x and y axes ticks step sizes.
-"""
 plotqq!(ax::Axis, gwas::DataFrame; kwargs...) = plotqq!(ax, gwas.P; kwargs...)