make targene specific outputs (#194)

* make targene specific outputs

* add missing test dep

* fix outputs to remove SAMPLE_IDS

* add svp into TargeneCore

* add MKL

* fix missing test dep

* add missing cmd function map

* fix test

* add new pvalue section per component

Author: olivierlabayle

Manifest.toml

@@ -15,6 +15,9 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -34,10 +37,12 @@ CairoMakie = "0.12"
 CategoricalArrays = "0.10"
 Combinatorics = "1.0"
 DataFrames = "1.2"
+OrderedCollections = "1.6.3"
 PackageCompiler = "2.1.17"
 SnpArrays = "0.3"
 StableRNGs = "1.0.1"
 Statistics = "1.10"
 TMLE = "0.17"
+MKL = "0.7"
 YAML = "0.4.9"
 julia = "1.10, 1"

Project.toml

@@ -1,5 +1,10 @@
 module TargeneCore
 
+
+if occursin("Intel", Sys.cpu_info()[1].model)
+    using MKL
+end
+
 using ArgParse
 using DataFrames
 using CSV
@@ -18,6 +23,8 @@ using JLD2
 using CairoMakie
 using TMLECLI
 using HypothesisTests
+using OrderedCollections
+using Mmap
 
 ###############################################################################
 ###                               INCLUDES                                  ###
@@ -26,10 +33,11 @@ using HypothesisTests
 include("utils.jl")
 include("confounders.jl")
 include("dataset.jl")
-include("plots.jl")
+include("outputs.jl")
 include("inputs_from_estimands.jl")
 include("inputs_from_config.jl")
 include("estimation_inputs.jl")
+include("sieve_variance.jl")
 include("cli.jl")
 
 ###############################################################################
@@ -38,7 +46,7 @@ include("cli.jl")
 
 export estimation_inputs
 export make_dataset
-export summary_plots
+export make_outputs
 export filter_chromosome, merge_beds, adapt_flashpca
 
 export get_outcome, get_treatments, get_outcome_extra_covariates, get_confounders, get_all_confounders

src/TargeneCore.jl

@@ -1,7 +1,7 @@
 function cli_settings()
     s = ArgParseSettings(description="TargeneCore CLI.")
 
-    @add_arg_table s begin
+    @add_arg_table! s begin
         "estimation-inputs"
             action = :command
             help = "Generate estimation inputs"
@@ -10,7 +10,7 @@ function cli_settings()
             action = :command
             help = "Generate a dataset."
 
-        "summary-plots"
+        "make-outputs"
             action = :command
             help = "Generate summary plots."
 
@@ -25,9 +25,13 @@ function cli_settings()
         "adapt-flashpca"
             action = :command
             help = "Adapts FlashPCA output."
+        
+        "svp"
+            action = :command
+            help = "Run Sieve Variance Plateau."
     end
 
-    @add_arg_table s["estimation-inputs"] begin
+    @add_arg_table! s["estimation-inputs"] begin
         "config-file"
             arg_type = String
             help = "Configuration file (.yaml)."
@@ -72,7 +76,7 @@ function cli_settings()
             default = 0
     end
 
-    @add_arg_table s["make-dataset"] begin
+    @add_arg_table! s["make-dataset"] begin
         "--genotypes-prefix"
             help = "Prefix path to BGEN/BED chromosomes."
             required = true
@@ -108,12 +112,12 @@ function cli_settings()
             default = 0
     end
 
-    @add_arg_table s["summary-plots"] begin
-        "results-file"
-            help = "Path to the results file."
+    @add_arg_table! s["make-outputs"] begin
+        "input-prefix"
+            help = "Input prefix to the results file."
             required = true
 
-        "--outprefix"
+        "--output-prefix"
             help = "Prefix to output plots."
             default = "."
 
@@ -123,7 +127,7 @@ function cli_settings()
             default = 0
     end
 
-    @add_arg_table s["filter-chromosome"] begin
+    @add_arg_table! s["filter-chromosome"] begin
         "input"
             help = "Prefix to bed files."
 
@@ -154,22 +158,58 @@ function cli_settings()
             required = false
     end
 
-    @add_arg_table s["merge-beds"] begin
+    @add_arg_table! s["merge-beds"] begin
         "input-prefix"
             help = "Prefix to bed files."
 
         "output"
             help = "Output file."
     end
 
-    @add_arg_table s["adapt-flashpca"] begin
+    @add_arg_table! s["adapt-flashpca"] begin
         "input"
             help = "Input file."
 
         "output"
             help = "Output file."
     end
 
+    @add_arg_table! s["svp"] begin
+        "input-prefix"
+            arg_type = String
+            help = "Input prefix to HDF5 files generated by the tmle CLI."
+
+        "--out"
+            arg_type = String
+            help = "Output filename."
+            default = "svp.hdf5"
+
+        "--grm-prefix"
+            arg_type = String
+            help = "Prefix to the aggregated GRM."
+            default = "GRM"
+        
+        "--verbosity"
+            arg_type = Int
+            default = 0
+            help = "Verbosity level"
+        
+        "--n-estimators"
+            arg_type = Int
+            default = 10
+            help = "Number of variance estimators to build for each estimate."
+
+        "--max-tau"
+            arg_type = Float64
+            default = 0.8
+            help = "Maximum distance between any two individuals."
+
+        "--estimator-key"
+            arg_type = String
+            help = "Estimator to use to proceed with sieve variance correction."
+            default = "1"
+    end
+
     return s
 end
 
@@ -199,10 +239,10 @@ function julia_main()::Cint
             call_threshold=cmd_settings["call-threshold"],
             verbosity=cmd_settings["verbosity"]
         )
-    elseif cmd == "summary-plots"
-        summary_plots(
-            cmd_settings["results-file"],
-            outprefix=cmd_settings["outprefix"],
+    elseif cmd == "make-outputs"
+        make_outputs(
+            cmd_settings["input-prefix"];
+            output_prefix=cmd_settings["output-prefix"],
             verbosity=cmd_settings["verbosity"],
         )
     elseif cmd == "filter-chromosome"
@@ -218,6 +258,15 @@ function julia_main()::Cint
         adapt_flashpca(
             cmd_settings["input"], cmd_settings["output"]
         )
+    elseif cmd =="svp"
+        sieve_variance_plateau(cmd_settings["input-prefix"];
+            out=cmd_settings["out"],
+            grm_prefix=cmd_settings["grm-prefix"],
+            verbosity=cmd_settings["verbosity"], 
+            n_estimators=cmd_settings["n-estimators"], 
+            max_tau=cmd_settings["max-tau"],
+            estimator_key=cmd_settings["estimator-key"]
+        )
     else
         throw(ArgumentError(string("No function matching command:", cmd)))
     end

src/cli.jl

@@ -0,0 +1,108 @@
+function add_pvalue_columns!(results)
+    for colname ∈ names(results)
+        results[!, Symbol(colname, :_PVALUE)] = [pvalue_or_nan(Ψ̂) for Ψ̂ ∈ results[!, colname]]
+    end
+end
+
+TMLE.estimate(Ψ̂::TMLECLI.FailedEstimate) = "Failed"
+
+get_treatment_changes(Ψ) = Dict(treatment_name => join(treatment_values, " => ") for (treatment_name, treatment_values) ∈ Ψ.treatment_values)
+get_treatment_changes(Ψ::TMLE.JointEstimand) = [get_treatment_changes(Ψᵢ) for Ψᵢ ∈ Ψ.args]
+
+get_estimand_type(Ψ) = replace(string(typeof(Ψ)), "TMLE.Statistical" => "")
+get_estimand_type(Ψ::TMLE.JointEstimand) = get_estimand_type(first(Ψ.args))
+
+log10_uniform_quantiles(n) = -log10.(collect(LinRange(0., 1., n + 1))[2:end])
+
+log10_beta_quantiles(n, alpha) = -log10.([quantile(Beta(k, n + 1 − k), alpha) for k in 1:n])
+
+function qqplot(output, results)
+    fig = Figure()
+    ax = Axis(fig[1, 1], 
+        title="Q-Q Plot",
+        xlabel = "-log₁₀(Uniform)",
+        ylabel = "-log₁₀(Observed)"
+    )
+    # QQ plots
+    pvalue_cols = [colname for colname ∈ names(results) if endswith(colname, "PVALUE")]
+    for pvalue_col ∈ pvalue_cols
+        pvalues = -log10.(filter(x -> !isnan(x), results[!, pvalue_col]))
+        n = length(pvalues)
+        unif_quantiles = log10_uniform_quantiles(n)
+        qqplot!(ax, unif_quantiles, pvalues, qqline=:identity, label=replace(string(pvalue_col), "_PVALUE" => ""))
+    end
+    # Confidence bands (from: https://lbelzile.github.io/lineaRmodels/qqplot.html)
+    n = size(results, 1)
+    unif_quantiles = log10_uniform_quantiles(n)
+    ub = log10_beta_quantiles(n, 0.025)
+    lb = log10_beta_quantiles(n, 0.975)
+    band!(ax, unif_quantiles, lb, ub, color=(:grey, 0.6))
+    # Legend
+    axislegend(position=:lt)
+    # Save figure
+    CairoMakie.save(output, fig)
+    return fig
+end
+
+summary_statistics(Ψ̂) = OrderedDict(
+        :PVALUE => pvalue_or_nan(Ψ̂),
+        :EFFECT_SIZE => TMLE.estimate(Ψ̂)
+)
+
+function summary_statistics(Ψ̂::TMLE.JointEstimate)
+    return OrderedDict(
+        :PVALUE => pvalue_or_nan(Ψ̂),
+        :COMPONENTS => [summary_statistics(Ψ̂ᵢ) for Ψ̂ᵢ ∈ Ψ̂.estimates]
+    )
+end
+
+function save_summary_yaml(output, results)
+    estimator_names = Symbol.(filter(x -> !endswith(x, "PVALUE"), names(results)))
+    estimator_1 = first(estimator_names)
+    summary_dicts = []
+    for row in eachrow(results)
+        Ψ = row[estimator_1].estimand
+        estimand_results = OrderedDict(
+            :EFFECT_TYPE => get_estimand_type(Ψ),
+            :OUTCOME => get_outcome(Ψ),
+            :TREATMENTS => get_treatment_changes(Ψ),            
+        )
+        for estimator_name in estimator_names
+            estimand_results[estimator_name] = summary_statistics(row[estimator_name])
+        end
+        push!(summary_dicts, estimand_results)
+    end
+    YAML.write_file(output, summary_dicts)
+end
+
+function load_results_as_df(input_prefix)
+    # Load from all files
+    input_files = files_matching_prefix(input_prefix)
+    results = DataFrame(TMLECLI.read_results_from_files(input_files))
+    # Drop potential SAMPLE_IDS columns
+    hasproperty(results, :SAMPLE_IDS) && select!(results, Not(:SAMPLE_IDS))
+    # Add p-values columns for each estimator
+    add_pvalue_columns!(results)
+    return results
+end
+
+function make_outputs(input_prefix;
+    output_prefix="final",
+    verbosity=1
+    )
+    # Load results
+    verbosity > 0 && @info("Loading results from files.")
+    results = load_results_as_df(input_prefix)
+    # Write full results
+    verbosity > 0 && @info("Saving aggregated results.")
+    jldsave(make_filepath_from_prefix(output_prefix; filename="results.hdf5"); results)
+    # Write Summary YAML file
+    verbosity > 0 && @info("Saving summary YAML file.")
+    save_summary_yaml(make_filepath_from_prefix(output_prefix; filename="results.summary.yaml"), results)
+    # Save QQ plot
+    verbosity > 0 && @info("Making QQ plot.")
+    qqplot(make_filepath_from_prefix(output_prefix; filename="QQ.png"), results)
+    
+    verbosity > 0 && @info("Done.")
+    return 0
+end