diff --git a/R/colocboost_assemble.R b/R/colocboost_assemble.R index fddc13c..9ccf7cd 100644 --- a/R/colocboost_assemble.R +++ b/R/colocboost_assemble.R @@ -47,7 +47,7 @@ colocboost_assemble <- function(cb_obj, if (data_info$n_outcomes == 1 & output_level == 1) { output_level <- 2 } - if (cb_obj$cb_model_para$num_updates == 1) { + if (cb_obj$cb_model_para$num_updates == 0) { cb_output <- list( "cos_summary" = NULL, "vcp" = NULL, @@ -59,10 +59,10 @@ colocboost_assemble <- function(cb_obj, if (output_level != 1) { tmp <- get_full_output(cb_obj = cb_obj, past_out = NULL, variables = NULL) if (output_level == 2) { - cb_output$ucos_details <- tmp$ucos_detials + cb_output <- c(cb_output, list("ucos_details" = tmp$ucos_detials)) cb_output <- cb_output[c("cos_summary", "vcp", "cos_details", "data_info", "model_info", "ucos_details")] } else { - cb_output$ucos_details <- tmp$ucos_detials + cb_output <- c(cb_output, list("ucos_details" = tmp$ucos_detials)) cb_output$diagnostic_details <- tmp[-1] cb_output <- cb_output[c("cos_summary", "vcp", "cos_details", "data_info", "model_info", "ucos_details", "diagnostic_details")] } diff --git a/R/colocboost_init.R b/R/colocboost_init.R index 767e7bf..dc1b21d 100644 --- a/R/colocboost_init.R +++ b/R/colocboost_init.R @@ -294,7 +294,7 @@ colocboost_init_para <- function(cb_data, cb_model, tau = 0.01, "variables" = cb_data$variable.names, "focal_outcome_idx" = focal_outcome_idx, "coveraged" = TRUE, - "num_updates" = 1, + "num_updates" = 0, "coveraged_outcome" = coveraged_outcome, "num_updates_outcome" = num_updates_outcome ) diff --git a/R/colocboost_one_causal.R b/R/colocboost_one_causal.R index e680fc9..0d70f1d 100644 --- a/R/colocboost_one_causal.R +++ b/R/colocboost_one_causal.R @@ -67,7 +67,7 @@ colocboost_one_iteration <- function(cb_model, cb_model_para, cb_data) { } } # -- remove redundant parameters - cb_model_para$num_updates <- 0 + cb_model_para$num_updates <- 1 rm_elements <- c("update_temp", "update_y") if (!is.null(cb_model_para$need_more)) { rm_elements <- c(rm_elements, "need_more") @@ -295,7 +295,7 @@ colocboost_diagLD <- function(cb_model, cb_model_para, cb_data) { } } # -- remove redundant parameters - cb_model_para$num_updates <- 0 + cb_model_para$num_updates <- 1 rm_elements <- c("update_temp", "update_y") if (!is.null(cb_model_para$need_more)) { rm_elements <- c(rm_elements, "need_more") diff --git a/R/colocboost_workhorse.R b/R/colocboost_workhorse.R index 7167ee4..da0fa0f 100644 --- a/R/colocboost_workhorse.R +++ b/R/colocboost_workhorse.R @@ -102,7 +102,6 @@ colocboost_workhorse <- function(cb_data, } if (sum(cb_model_para$update_y == 1) == 0) { - cb_model_para$num_updates <- 1 cb_obj <- list("cb_data" = cb_data, "cb_model" = cb_model, "cb_model_para" = cb_model_para) class(cb_obj) <- "colocboost" return(cb_obj) diff --git a/tests/testthat/test_colocboost.R b/tests/testthat/test_colocboost.R index c31d914..58db554 100644 --- a/tests/testthat/test_colocboost.R +++ b/tests/testthat/test_colocboost.R @@ -294,4 +294,6 @@ test_that("colocboost handles missing/invalid inputs appropriately", { ), "Please provide the sample index of X and Y, since they do not have the same samples!" ) -}) \ No newline at end of file +}) + + diff --git a/tests/testthat/test_one_causal.R b/tests/testthat/test_one_causal.R new file mode 100644 index 0000000..9128c06 --- /dev/null +++ b/tests/testthat/test_one_causal.R @@ -0,0 +1,293 @@ +library(testthat) + +# Helper function for test data generation +generate_one_causal_test_data <- function(n = 200, p = 30, L = 2, seed = 123) { + set.seed(seed) + + # Generate X with strong LD structure to test LD handling + sigma <- matrix(0, p, p) + for (i in 1:p) { + for (j in 1:p) { + sigma[i, j] <- 0.95^abs(i - j) # Higher LD than standard test + } + } + X <- MASS::mvrnorm(n, rep(0, p), sigma) + colnames(X) <- paste0("SNP", 1:p) + + # Generate true effects - clear single causal variant per trait + true_beta <- matrix(0, p, L) + # Different causal variants for each trait + true_beta[5, 1] <- 0.8 # Strong effect for trait 1 + true_beta[15, 2] <- 0.9 # Strong effect for trait 2 + + # Generate Y with some noise + Y <- matrix(0, n, L) + for (l in 1:L) { + Y[, l] <- X %*% true_beta[, l] + rnorm(n, 0, 0.7) # Less noise for clearer signals + } + + # Calculate beta, se, and z-scores + beta <- matrix(0, ncol(X), ncol(Y)) + se <- matrix(0, ncol(X), ncol(Y)) + z <- matrix(0, ncol(X), ncol(Y)) + + for (i in 1:ncol(Y)) { + output = matrix(0, ncol(X), 2) + for (j in 1:ncol(X)) { + fit = summary(lm(Y[,i] ~ X[,j]))$coef + if (nrow(fit) == 2) + output[j,] = as.vector(fit[2,1:2]) + else + output[j,] = c(0,0) + } + beta[,i] <- output[,1] + se[,i] <- output[,2] + z[,i] <- beta[,i]/se[,i] + } + + # Create summary statistics data frames + sumstat_list <- list() + for (i in 1:ncol(Y)) { + sumstat_list[[i]] <- data.frame( + beta = beta[,i], + sebeta = se[,i], + z = z[,i], + n = nrow(X), + variant = colnames(X) + ) + } + + # Add effect matrices for HyPrColoc format + rownames(beta) <- rownames(se) <- colnames(X) + colnames(beta) <- colnames(se) <- paste0("Y", 1:ncol(Y)) + + # Return all formats + list( + X = X, + Y = Y, + LD = cor(X), + sumstat = sumstat_list, + effect_est = beta, + effect_se = se, + effect_n = rep(n, ncol(Y)) + ) +} + +# Create common test data - do this outside of test_that blocks +test_data <- generate_one_causal_test_data() +three_trait_data <- generate_one_causal_test_data(L = 3) +non_coloc_data <- generate_one_causal_test_data() # Different seed would be better for true non-coloc + +# Test 1: LD-free mode (diagonal LD matrix) +test_that("colocboost runs in LD-free mode correctly", { + # Run colocboost without LD matrix + warnings <- capture_warnings({ + result_ld_free <- colocboost( + sumstat = test_data$sumstat, + M = 1 # One iteration for LD-free mode + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_ld_free, "colocboost") + + # Check that it used the LD-free mode + expect_true(result_ld_free$model_info$model_coveraged == "one_causal") + + # Check that there's exactly one iteration + expect_equal(result_ld_free$model_info$n_updates, 1) + + # Check dimensions + expect_equal(length(result_ld_free$data_info$variables), ncol(test_data$X)) + expect_equal(result_ld_free$data_info$n_outcomes, 2) +}) + +# Test 2: One iteration mode with LD +test_that("colocboost runs in one iteration mode correctly", { + # Run colocboost with LD matrix but M=1 + warnings <- capture_warnings({ + result_one_iter <- colocboost( + sumstat = test_data$sumstat, + LD = test_data$LD, + M = 1 # One iteration + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_one_iter, "colocboost") + + # Check that it used the one_causal mode + expect_true(result_one_iter$model_info$model_coveraged == "one_causal") + + # Check that the number of updates is 1 + expect_equal(result_one_iter$model_info$n_updates, 1) + + # It should have the update_status matrix + expect_true(!is.null(result_one_iter$model_info$jk_star)) +}) + +# Test 3: Focal outcome functionality with one causal +test_that("colocboost one causal mode handles focal outcome correctly", { + # Run colocboost with focal_outcome_idx = 1 + warnings <- capture_warnings({ + result_focal <- colocboost( + sumstat = test_data$sumstat, + LD = test_data$LD, + M = 1, + focal_outcome_idx = 1 + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_focal, "colocboost") + + # Check focal outcome is correctly set + expect_equal(result_focal$data_info$outcome_info$is_focal[1], TRUE) + expect_equal(result_focal$data_info$outcome_info$is_focal[2], FALSE) + + # Check that it used the one_causal mode + expect_true(result_focal$model_info$model_coveraged == "one_causal") +}) + +# Test 4: Different lambda values for focal outcome +test_that("colocboost one causal mode uses different lambda for focal outcome", { + # Run with different lambda for focal outcome + warnings <- capture_warnings({ + result_lambda <- colocboost( + sumstat = test_data$sumstat, + LD = test_data$LD, + M = 1, + focal_outcome_idx = 1, + lambda = 0.3, # Base lambda + lambda_focal_outcome = 0.8 # Different lambda for focal outcome + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_lambda, "colocboost") + + # Result should use one_causal approach + expect_true(result_lambda$model_info$model_coveraged == "one_causal") +}) + +# Test 5: Multiple traits with no true colocalization +test_that("colocboost one causal mode correctly handles non-colocalized traits", { + # Run colocboost + warnings <- capture_warnings({ + result_non_coloc <- colocboost( + sumstat = non_coloc_data$sumstat, + LD = non_coloc_data$LD, + M = 1 + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_non_coloc, "colocboost") + + # Result should use one_causal approach + expect_true(result_non_coloc$model_info$model_coveraged == "one_causal") +}) + +# Test 6: With three or more traits +test_that("colocboost one causal mode handles three or more traits", { + # Run colocboost + warnings <- capture_warnings({ + result_three_traits <- colocboost( + sumstat = three_trait_data$sumstat, + LD = three_trait_data$LD, + M = 1 + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_three_traits, "colocboost") + + # Check dimensions + expect_equal(result_three_traits$data_info$n_outcomes, 3) + + # Result should use one_causal approach + expect_true(result_three_traits$model_info$model_coveraged == "one_causal") +}) + +# Test 7: Effect of different jk_equiv parameters +test_that("colocboost one causal mode respects jk_equiv parameters", { + # Run with stricter jk_equiv settings + warnings <- capture_warnings({ + result_strict <- colocboost( + sumstat = test_data$sumstat, + LD = test_data$LD, + M = 1, + jk_equiv_corr = 0.95, # Very high correlation required + jk_equiv_loglik = 0.5 # More stringent loglik difference + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_strict, "colocboost") + + # Should still be one_causal mode + expect_true(result_strict$model_info$model_coveraged == "one_causal") + + # Run with lenient jk_equiv settings + warnings <- capture_warnings({ + result_lenient <- colocboost( + sumstat = test_data$sumstat, + LD = test_data$LD, + M = 1, + jk_equiv_corr = 0.5, # Lower correlation required + jk_equiv_loglik = 2.0 # Larger loglik difference allowed + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_lenient, "colocboost") +}) + +# Test 8: Interaction with output level parameter +test_that("colocboost one causal mode works with different output_level values", { + for (level in 1:3) { + warnings <- capture_warnings({ + result <- colocboost( + sumstat = test_data$sumstat, + LD = test_data$LD, + M = 1, + output_level = level + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result, "colocboost") + + # Check output level-specific fields + if (level >= 2) { + expect_true("ucos_details" %in% names(result)) + } + if (level >= 3) { + expect_true("diagnostic_details" %in% names(result)) + } + } +}) + +# Test 10: Test with HyPrColoc compatible format +test_that("colocboost one causal mode works with HyPrColoc format", { + # Run colocboost with effect matrices + warnings <- capture_warnings({ + result_hypr <- colocboost( + effect_est = test_data$effect_est, + effect_se = test_data$effect_se, + effect_n = test_data$effect_n, + LD = test_data$LD, + M = 1 + ) + }) + + # Test that we get a colocboost object + expect_s3_class(result_hypr, "colocboost") + + # Check that it used the one_causal mode + expect_true(result_hypr$model_info$model_coveraged == "one_causal") + + # Check dimensions + expect_equal(length(result_hypr$data_info$variables), ncol(test_data$X)) + expect_equal(result_hypr$data_info$n_outcomes, 2) +}) \ No newline at end of file diff --git a/tests/testthat/test_plot.R b/tests/testthat/test_plot.R index 7c10cdf..775d3fe 100644 --- a/tests/testthat/test_plot.R +++ b/tests/testthat/test_plot.R @@ -75,12 +75,12 @@ generate_test_result <- function(n = 500, p = 60, L = 4, seed = 42, output_level return(result) } +# Generate a test colocboost result +cb_res <- generate_test_result() + # Test colocboost_plot function with basic options test_that("colocboost_plot basic functionality works", { - # Generate a test colocboost result - cb_res <- generate_test_result() - # Basic plot call expect_error(suppressWarnings(colocboost_plot(cb_res)), NA) @@ -92,9 +92,6 @@ test_that("colocboost_plot basic functionality works", { # Test colocboost_plot with different y-axis options test_that("colocboost_plot handles different y-axis options", { - # Generate a test colocboost result - cb_res <- generate_test_result() - # Test with different y-axis values expect_error(suppressWarnings(colocboost_plot(cb_res, y = "log10p")), NA) expect_error(suppressWarnings(colocboost_plot(cb_res, y = "z_original")), NA) @@ -110,9 +107,6 @@ test_that("colocboost_plot handles different y-axis options", { # Test colocboost_plot with plot filtering options test_that("colocboost_plot handles filtering options", { - # Generate a test colocboost result - cb_res <- generate_test_result() - # Test with outcome index filtering expect_error(suppressWarnings(colocboost_plot(cb_res, outcome_idx = 1)), NA) expect_error(suppressWarnings(colocboost_plot(cb_res, outcome_idx = 2)), NA) @@ -140,9 +134,6 @@ test_that("colocboost_plot handles filtering options", { # Test colocboost_plot with visual customization options test_that("colocboost_plot handles visual customization options", { - # Generate a test colocboost result - cb_res <- generate_test_result() - # Test with custom colors expect_error(suppressWarnings(colocboost_plot(cb_res, points_color = "red")), NA) expect_error(suppressWarnings(colocboost_plot(cb_res, cos_color = c("blue", "green", "orange"))), NA) @@ -165,9 +156,6 @@ test_that("colocboost_plot handles visual customization options", { # Test colocboost_plot with layout options test_that("colocboost_plot handles layout options", { - # Generate a test colocboost result - cb_res <- generate_test_result() - # Test with different plot_cols values expect_error(suppressWarnings(colocboost_plot(cb_res, plot_cols = 1)), NA) expect_error(suppressWarnings(colocboost_plot(cb_res, plot_cols = 3)), NA) @@ -187,9 +175,6 @@ test_that("colocboost_plot handles layout options", { # Test colocboost_plot with additional visualization options test_that("colocboost_plot handles additional visualization options", { - # Generate a test colocboost result - cb_res <- generate_test_result() - # Test with vertical line options expect_error(suppressWarnings(colocboost_plot(cb_res, add_vertical = TRUE, add_vertical_idx = c(5, 10))), NA) @@ -201,6 +186,51 @@ test_that("colocboost_plot handles additional visualization options", { expect_error(suppressWarnings(colocboost_plot(cb_res, show_variable = TRUE)), NA) }) +# Test colocboost_plot with custom outcome names +test_that("colocboost_plot handles custom outcome names", { + + # Test with custom outcome names + expect_error(suppressWarnings(colocboost_plot(cb_res, outcome_names = c("Trait1", "Trait2"))), NA) +}) + +# Test colocboost_plot with a specific range +test_that("colocboost_plot handles zoom-in with grange", { + + # Test with grange option to zoom in + expect_error(suppressWarnings(colocboost_plot(cb_res, grange = 5:15)), NA) +}) + +# Test colocboost_plot with focal outcome in L=4 case +test_that("colocboost_plot handles focal outcome in complex cases", { + + # Generate a test colocboost result with 4 traits and focal outcome set + cb_res_focal <- generate_test_result(L = 4, output_level = 3) + + # Basic plot call with focal outcome + expect_error(suppressWarnings(colocboost_plot(cb_res_focal)), NA) + + # Test plot_focal_only option + expect_error(suppressWarnings(colocboost_plot(cb_res_focal, plot_focal_only = TRUE)), NA) + + # Test plot_focal_cos_outcome_only option + expect_error(suppressWarnings(colocboost_plot(cb_res_focal, plot_focal_cos_outcome_only = TRUE)), NA) + + # Combine focal outcome filtering with other options + expect_error(suppressWarnings(colocboost_plot(cb_res_focal, + plot_focal_only = TRUE, + y = "cos_vcp")), NA) + + expect_error(suppressWarnings(colocboost_plot(cb_res_focal, + plot_focal_cos_outcome_only = TRUE, + plot_ucos = TRUE)), NA) + + # Test focusing only on outcomes colocalized with focal outcome + expect_error(suppressWarnings(colocboost_plot(cb_res_focal, + plot_focal_cos_outcome_only = TRUE, + outcome_idx = 1:3)), NA) +}) + + # Test colocboost_plot with single trait (finemapping) results test_that("colocboost_plot handles single trait results", { @@ -215,25 +245,6 @@ test_that("colocboost_plot handles single trait results", { expect_error(suppressWarnings(colocboost_plot(cb_res_single, plot_cols = 1)), NA) }) -# Test colocboost_plot with uncolocalized visualization options -test_that("colocboost_plot handles uncolocalized visualization options", { - - # Generate a test colocboost result with high output level to include ucos details - cb_res <- generate_test_result(output_level = 3) - - # Test with plot_ucos options - expect_error(suppressWarnings(colocboost_plot(cb_res, plot_ucos = TRUE)), NA) - - # Test with show_cos_to_uncoloc options - expect_error(suppressWarnings(colocboost_plot(cb_res, show_cos_to_uncoloc = TRUE)), NA) - - # Generate a different colocboost result to test the warning for plot_ucos - cb_res_low <- generate_test_result(output_level = 1) - # This should give a warning but not an error - expect_warning(colocboost_plot(cb_res_low, plot_ucos = TRUE), - "Since you want to plot trait-specific \\(uncolocalized\\) sets with plot_ucos = TRUE") -}) - # Test colocboost_plot with L=4 case for complex colocalization and trait-specific effects test_that("colocboost_plot handles L=4 case with complex colocalization patterns", { @@ -273,54 +284,25 @@ test_that("colocboost_plot handles L=4 case with complex colocalization patterns show_cos_to_uncoloc = TRUE)), NA) }) -# Test colocboost_plot with custom outcome names -test_that("colocboost_plot handles custom outcome names", { - - # Generate a test colocboost result - cb_res <- generate_test_result() - - # Test with custom outcome names - expect_error(suppressWarnings(colocboost_plot(cb_res, outcome_names = c("Trait1", "Trait2"))), NA) -}) -# Test colocboost_plot with a specific range -test_that("colocboost_plot handles zoom-in with grange", { - - # Generate a test colocboost result - cb_res <- generate_test_result() - - # Test with grange option to zoom in - expect_error(suppressWarnings(colocboost_plot(cb_res, grange = 5:15)), NA) -}) +# Generate a test colocboost result with high output level to include ucos details +cb_res <- generate_test_result(output_level = 3) -# Test colocboost_plot with focal outcome in L=4 case -test_that("colocboost_plot handles focal outcome in complex cases", { - - # Generate a test colocboost result with 4 traits and focal outcome set - cb_res_focal <- generate_test_result(L = 4, output_level = 3) - - # Basic plot call with focal outcome - expect_error(suppressWarnings(colocboost_plot(cb_res_focal)), NA) - - # Test plot_focal_only option - expect_error(suppressWarnings(colocboost_plot(cb_res_focal, plot_focal_only = TRUE)), NA) +# Test colocboost_plot with uncolocalized visualization options +test_that("colocboost_plot handles uncolocalized visualization options", { - # Test plot_focal_cos_outcome_only option - expect_error(suppressWarnings(colocboost_plot(cb_res_focal, plot_focal_cos_outcome_only = TRUE)), NA) - # Combine focal outcome filtering with other options - expect_error(suppressWarnings(colocboost_plot(cb_res_focal, - plot_focal_only = TRUE, - y = "cos_vcp")), NA) + # Test with plot_ucos options + expect_error(suppressWarnings(colocboost_plot(cb_res, plot_ucos = TRUE)), NA) - expect_error(suppressWarnings(colocboost_plot(cb_res_focal, - plot_focal_cos_outcome_only = TRUE, - plot_ucos = TRUE)), NA) + # Test with show_cos_to_uncoloc options + expect_error(suppressWarnings(colocboost_plot(cb_res, show_cos_to_uncoloc = TRUE)), NA) - # Test focusing only on outcomes colocalized with focal outcome - expect_error(suppressWarnings(colocboost_plot(cb_res_focal, - plot_focal_cos_outcome_only = TRUE, - outcome_idx = 1:3)), NA) + # Generate a different colocboost result to test the warning for plot_ucos + cb_res_low <- generate_test_result(output_level = 1) + # This should give a warning but not an error + expect_warning(colocboost_plot(cb_res_low, plot_ucos = TRUE), + "Since you want to plot trait-specific \\(uncolocalized\\) sets with plot_ucos = TRUE") }) # Test colocboost_plot with varying cutoff settings from get_robust_colocalization @@ -369,4 +351,6 @@ test_that("colocboost_plot handles varying cutoff settings", { plot_ucos = TRUE, y = "cos_vcp")), NA) } -}) \ No newline at end of file +}) + +