fix and add tests (#173)

Author: olivierlabayle

src/sieve_plateau.jl

@@ -57,7 +57,8 @@ function build_work_list(prefix, grm_ids; pval=0.05)
             end
         end
     end
-    return reduce(vcat, transpose(influence_curves)), n_obs, file_tmlereport_pairs
+    influence_curves = length(influence_curves) > 0 ? reduce(vcat, transpose(influence_curves)) : influence_curves 
+    return influence_curves, n_obs, file_tmlereport_pairs
 end
 
 
@@ -152,6 +153,15 @@ function save_results(outfilename, τs, variances, std_errors, file_queryreport_
     end
 end
 
+function save_empty_results(outfilename, τs, file_queryreport_pairs)
+    jldopen(outfilename, "w") do io
+        io["TAUS"] = τs
+        io["VARIANCES"] = []
+        io["SOURCEFILE_REPORTID_PAIRS"] = file_queryreport_pairs
+        io["STDERRORS"] = []
+    end
+end
+
 
 corrected_stderrors(variances, n_obs) =
     sqrt.(view(maximum(variances, dims=1), 1, :) ./ n_obs)
@@ -166,10 +176,13 @@ function sieve_variance_plateau(parsed_args)
     grm, grm_ids = readGRM(parsed_args["grm-prefix"])
 
     influence_curves, n_obs, file_queryreport_pairs = build_work_list(prefix, grm_ids; pval=pval)
-    variances = compute_variances(influence_curves, grm, τs, n_obs)
-    std_errors = corrected_stderrors(variances, n_obs)
-
-    save_results(outfilename, τs, variances, std_errors, file_queryreport_pairs)
+    # If the work list is not empty
+    if length(influence_curves) > 0
+        variances = compute_variances(influence_curves, grm, τs, n_obs)
+        std_errors = corrected_stderrors(variances, n_obs)
+        save_results(outfilename, τs, variances, std_errors, file_queryreport_pairs)
+    # Otherwise save empty structure
+    else
+        save_empty_results(outfilename, τs, file_queryreport_pairs)
+    end
 end
-
-

test/sieve_plateau.jl

@@ -211,4 +211,34 @@ end
     clean_hdf5estimates_files(prefix)
 end
 
+@testset "Test sieve_variance_plateau empty work list" begin
+    # In this case no p-value will make the threshold
+    grm_ids = TMLEEpistasis.GRMIDs("data/grm/test.grm.id")
+    prefix = "rs12_rs45"
+    nb_estimators = 10
+    build_results_files(grm_ids, prefix)
+    outfilename = "rs12_rs45_sieve_variance.hdf5"
+    parsed_args = Dict(
+        "prefix" => prefix,
+        "pval" => -1,
+        "grm-prefix" => "data/grm/test.grm",
+        "out" => outfilename,
+        "nb-estimators" => nb_estimators,
+        "max-tau" => 0.75
+    )
+
+    sieve_variance_plateau(parsed_args)
+
+    results_file = jldopen(outfilename)
+    @test size(results_file["TAUS"]) == (nb_estimators,)
+    @test size(results_file["VARIANCES"]) == (0,)
+    @test results_file["SOURCEFILE_REPORTID_PAIRS"] == []
+    @test size(results_file["STDERRORS"]) == (0,)
+
+    close(results_file)
+
+    rm(outfilename)
+    clean_hdf5estimates_files(prefix)
+end
+
 end

test/summary.jl

@@ -84,6 +84,48 @@ end
 end
 
 
+@testset "Test build_summary with empty sieve" begin
+    grm_ids = TMLEEpistasis.GRMIDs("data/grm/test.grm.id")
+    prefix = "rs12_rs45"
+    summaryfilename = summary_filename(prefix)
+    build_results_files(grm_ids, prefix)
+
+    TMLEEpistasis.save_empty_results(sieve_filename(prefix), [], [])
+    parsed_args = Dict(
+        "prefix" => prefix,
+        "out" => summary_filename(prefix),
+        "sieve" => true
+    )
+
+    build_summary(parsed_args)
+    
+    summary = CSV.read(summaryfilename, DataFrame)
+
+    @test summary.PHENOTYPE == ["cancer", "cancer", "height", "height", "bmi", "bmi"]
+    @test summary.PHENOTYPE_TYPE == ["BINARY", "BINARY", "CONTINUOUS", "CONTINUOUS", "CONTINUOUS", "CONTINUOUS"]
+    @test summary.QUERYNAME == ["QUERY_1", "QUERY_2", "QUERY_1", "QUERY_2", "QUERY_1", "QUERY_2"]
+    @test summary.FILENAME_ORIGIN == ["rs12_rs45_batch_1_Bool.hdf5", 
+                                      "rs12_rs45_batch_1_Bool.hdf5", 
+                                      "rs12_rs45_batch_1_Real.hdf5", 
+                                      "rs12_rs45_batch_1_Real.hdf5", 
+                                      "rs12_rs45_batch_1_Real.hdf5", 
+                                      "rs12_rs45_batch_1_Real.hdf5"]
+    @test summary.REPORT_KEY == ["1_1", "1_2", "1_1", "1_2", "2_1", "2_2"]
+    @test eltype(summary.INITIAL_ESTIMATE) == Float64
+    @test eltype(summary.ESTIMATE) == Float64
+    @test eltype(summary.PVAL) == Float64
+    @test eltype(summary.LWB) == Float64
+    @test eltype(summary.UPB) == Float64
+    for col in (:SIEVE_STDERR, :SIEVE_PVAL, :SIEVE_LWB, :SIEVE_UPB)
+        @test eltype(summary[!, col]) == Missing
+    end
+
+    clean_hdf5estimates_files(prefix)
+    rm(sieve_filename(prefix))
+    rm(summaryfilename)
+
+end
+
 end
 
 true