diff --git a/.Rbuildignore b/.Rbuildignore index f73ddc0..bdf2b2d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -20,3 +20,4 @@ ^docs$ ^_pkgdown\.yml$ ^CRAN-SUBMISSION$ +^\.superpowers$ diff --git a/DESCRIPTION b/DESCRIPTION index ca7aa6c..4e16d32 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: colocboost Type: Package -Date: 2026-06-07 +Date: 2026-09-08 Title: Multi-Context Colocalization Analysis for QTL and GWAS Studies -Version: 1.0.9 +Version: 1.0.10 Authors@R: c( person(given = "Xuewei", family = "Cao", email = "xc2270@cumc.columbia.edu", role = c("cre", "aut", "cph")), person(given = "Haochen", family = "Sun", email = "hs3393@cumc.columbia.edu", role = c("aut", "cph")), diff --git a/R/colocboost_assemble_cos.R b/R/colocboost_assemble_cos.R index 5912070..254c652 100644 --- a/R/colocboost_assemble_cos.R +++ b/R/colocboost_assemble_cos.R @@ -1,4 +1,10 @@ #' @importFrom stats as.dist cutree hclust +.group_coloc_candidates <- function(update, pos.coloc) { + signatures <- apply(update[, pos.coloc, drop = FALSE], 2, paste0, collapse = ",") + group_ids <- match(signatures, unique(signatures)) + split(seq_along(pos.coloc), group_ids) +} + colocboost_assemble_cos <- function(cb_obj, coverage = 0.95, weight_fudge_factor = 1.5, @@ -22,6 +28,7 @@ colocboost_assemble_cos <- function(cb_obj, cb_model <- cb_obj$cb_model cb_model_para <- cb_obj$cb_model_para cb_data <- cb_obj$cb_data + purity_outcomes <- .cb_unique_purity_outcomes(cb_data, seq_len(cb_model_para$L)) # define the confident sets for colocalization update <- cb_model_para$update_status @@ -124,22 +131,15 @@ colocboost_assemble_cos <- function(cb_obj, } } } else { - coloc_candidate <- update[, pos.coloc] - coloc_candidate <- apply(coloc_candidate, 2, paste0, collapse = ",") - coloc_temp <- table(coloc_candidate) - # iterations for each colocalization sets - pos_coloc_sets <- lapply(1:length(coloc_temp), function(x) { - which(coloc_candidate == names(coloc_temp)[x]) - }) - names(pos_coloc_sets) <- names(coloc_temp) + pos_coloc_sets <- .group_coloc_candidates(update, pos.coloc) # - define coloc_sets coloc_sets <- avWeight_coloc_sets <- total_change_Loglik_coloc <- evidence_strength_coloc <- cs_change_coloc <- coloc_outcomes_sets <- list() flag <- 0 - for (i in 1:length(coloc_temp)) { + for (i in seq_along(pos_coloc_sets)) { pos_temp_coloc_each <- pos_coloc_sets[[i]] - coloc_outcomes <- which(unlist(strsplit(names(coloc_temp)[i], ",")) == 1) + coloc_outcomes <- which(update[, pos.coloc[pos_temp_coloc_each[1]]] == 1) # - if only one iteration for this coloc_set if (length(pos_temp_coloc_each) == 1) { @@ -349,14 +349,17 @@ colocboost_assemble_cos <- function(cb_obj, # calculate between purity ncsets <- length(coloc_sets) min_between <- max_between <- ave_between <- matrix(0, nrow = ncsets, ncol = ncsets) - for (i.between in 1:(ncsets - 1)) { - for (j.between in (i.between + 1):ncsets) { + overlap_pairs <- .merge_ucos_overlap_pairs(coloc_sets) + if (nrow(overlap_pairs) > 0L) { + for (pair_idx in seq_len(nrow(overlap_pairs))) { + i.between <- overlap_pairs[pair_idx, 1L] + j.between <- overlap_pairs[pair_idx, 2L] cset1 <- coloc_sets[[i.between]] cset2 <- coloc_sets[[j.between]] res <- list() - for (i in 1:cb_model_para$L) { + for (i in purity_outcomes) { X_dict <- cb_data$dict[i] - res[[i]] <- get_between_purity(cset1, cset2, + res[[length(res) + 1L]] <- get_between_purity(cset1, cset2, X = cb_data$data[[X_dict]]$X, Xcorr = cb_data$data[[X_dict]]$XtX, miss_idx = cb_data$data[[i]]$variable_miss, diff --git a/R/colocboost_check_update_jk.R b/R/colocboost_check_update_jk.R index aef6b0b..8639114 100644 --- a/R/colocboost_check_update_jk.R +++ b/R/colocboost_check_update_jk.R @@ -5,6 +5,62 @@ #' #' @return update_status and real_update_jk for each trait #' @noRd +.cb_append_update_history <- function(cb_model_para, update_jk, update_status, real_update_jk) { + used <- attr(cb_model_para, "update_history_n") + used <- if (is.null(used)) 0L else used + capacity <- attr(cb_model_para, "update_history_capacity") + required <- used + 1L + + if (is.null(capacity) || capacity < required) { + new_capacity <- max(128L, required, if (is.null(capacity)) 0L else capacity * 2L) + L <- cb_model_para$L + + update_status_new <- matrix(0, nrow = L, ncol = new_capacity) + real_update_jk_new <- matrix(NA_real_, nrow = new_capacity, ncol = L) + jk_new <- matrix(NA_real_, nrow = new_capacity, ncol = L + 1L) + + if (used > 0L) { + used_idx <- seq_len(used) + update_status_new[, used_idx] <- cb_model_para$update_status[, used_idx, drop = FALSE] + real_update_jk_new[used_idx, ] <- cb_model_para$real_update_jk[used_idx, , drop = FALSE] + jk_new[used_idx, ] <- cb_model_para$jk[used_idx, , drop = FALSE] + } + + cb_model_para$update_status <- update_status_new + cb_model_para$real_update_jk <- real_update_jk_new + cb_model_para$jk <- jk_new + attr(cb_model_para, "update_history_capacity") <- new_capacity + } + + cb_model_para$update_status[, required] <- update_status + cb_model_para$real_update_jk[required, ] <- real_update_jk + cb_model_para$jk[required, ] <- update_jk + attr(cb_model_para, "update_history_n") <- required + cb_model_para +} + +.cb_trim_update_history <- function(cb_model_para) { + used <- attr(cb_model_para, "update_history_n") + if (is.null(used)) { + return(cb_model_para) + } + + if (used == 0L) { + cb_model_para$update_status <- c() + cb_model_para$real_update_jk <- c() + cb_model_para$jk <- c() + } else { + used_idx <- seq_len(used) + cb_model_para$update_status <- cb_model_para$update_status[, used_idx, drop = FALSE] + cb_model_para$real_update_jk <- cb_model_para$real_update_jk[used_idx, , drop = FALSE] + cb_model_para$jk <- cb_model_para$jk[used_idx, , drop = FALSE] + } + + attr(cb_model_para, "update_history_n") <- NULL + attr(cb_model_para, "update_history_capacity") <- NULL + cb_model_para +} + colocboost_check_update_jk <- function(cb_model, cb_model_para, cb_data) { pos.update <- which(cb_model_para$update_y == 1) @@ -299,9 +355,12 @@ boost_check_update_jk_nofocal <- function(cb_model, cb_model_para, cb_data) { } # - update cb_model and report results - cb_model_para$jk <- rbind(cb_model_para$jk, update_jk) - cb_model_para$update_status <- cbind(cb_model_para$update_status, as.matrix(update_status)) - cb_model_para$real_update_jk <- rbind(cb_model_para$real_update_jk, real_update_jk) + cb_model_para <- .cb_append_update_history( + cb_model_para, + update_jk = update_jk, + update_status = update_status, + real_update_jk = real_update_jk + ) update_temp <- list( "update_status" = update_status, @@ -478,9 +537,12 @@ boost_check_update_jk_focal <- function(cb_model, cb_model_para, cb_data, } # - update cb_model and report results - cb_model_para$jk <- rbind(cb_model_para$jk, update_jk) - cb_model_para$update_status <- cbind(cb_model_para$update_status, as.matrix(update_status)) - cb_model_para$real_update_jk <- rbind(cb_model_para$real_update_jk, real_update_jk) + cb_model_para <- .cb_append_update_history( + cb_model_para, + update_jk = update_jk, + update_status = update_status, + real_update_jk = real_update_jk + ) update_temp <- list( "update_status" = update_status, @@ -554,70 +616,80 @@ check_pair_jkeach <- function(jk_each, jk_equiv_corr = 0.8, jk_equiv_loglik = 0.001) { + n_pair <- length(jk_each) + if (n_pair <= 1) { + return(matrix(0, nrow = n_pair, ncol = n_pair)) + } #' @importFrom stats cor get_LD_jk_each <- function(jk_each, X = NULL, XtX = NULL, N = NULL, remain_jk = NULL, ref_label = "LD") { + jk_unique <- unique(jk_each) + unique_idx <- match(jk_each, jk_unique) + if (!is.null(X)) { - LD_temp <- suppressWarnings({ - get_cormat(X[, jk_each]) + LD_unique <- suppressWarnings({ + get_cormat(X[, jk_unique, drop = FALSE]) }) - LD_temp[which(is.na(LD_temp))] <- 0 - # LD_temp <- LD_temp[1, 2] + LD_temp <- LD_unique[unique_idx, unique_idx, drop = FALSE] } else if (!is.null(XtX)) { if (identical(ref_label, "No_ref")) { LD_temp <- matrix(0, length(jk_each), length(jk_each)) } else { - jk.remain <- match(jk_each, remain_jk) + jk.remain <- match(jk_unique, remain_jk) if (identical(ref_label, "X_ref")) { - LD_temp <- suppressWarnings({ get_cormat(XtX[, jk.remain]) }) + LD_unique <- suppressWarnings({ + get_cormat(XtX[, jk.remain, drop = FALSE]) + }) } else { - LD_temp <- XtX[jk.remain, jk.remain] + LD_unique <- XtX[jk.remain, jk.remain, drop = FALSE] } - LD_temp[which(is.na(LD_temp))] <- 0 + LD_temp <- LD_unique[unique_idx, unique_idx, drop = FALSE] } } return(LD_temp) } - detect_func <- function(idx, LD_all, jk_i, jk_j, i, j){ - change_log_jk_i <- model_update[[idx]]$change_loglike[jk_i] - change_log_jk_j <- model_update[[idx]]$change_loglike[jk_j] - change_each <- abs(change_log_jk_i - change_log_jk_j) - LD_temp <- LD_all[[idx]][i, j] - return((change_each <= jk_equiv_loglik) & (abs(LD_temp) >= jk_equiv_corr)) + get_ld_key <- function(idx) { + ref_idx <- X_dict[idx] + ref_label <- cb_data$data[[ref_idx]]$ref_label + missing_key <- paste(data_update[[idx]]$variable_miss, collapse = ",") + paste(ref_idx, ref_label, length(model_update[[idx]]$res), missing_key, sep = "|") } data_update <- cb_data$data[pos.update] - LD_all <- lapply(1:length(jk_each), function(idx){ - get_LD_jk_each(jk_each, - X = cb_data$data[[X_dict[idx]]]$X, - XtX = cb_data$data[[X_dict[idx]]]$XtX, - N = data_update[[idx]]$N, - remain_jk = setdiff(1:length(model_update[[idx]]$res), data_update[[idx]]$variable_miss), - ref_label = cb_data$data[[X_dict[idx]]]$ref_label - ) - }) + ld_keys <- vapply(seq_along(jk_each), get_ld_key, character(1)) + ld_cache <- list() + for (idx in seq_along(jk_each)) { + ld_key <- ld_keys[idx] + if (is.null(ld_cache[[ld_key]])) { + ref_idx <- X_dict[idx] + ld_cache[[ld_key]] <- get_LD_jk_each(jk_each, + X = cb_data$data[[ref_idx]]$X, + XtX = cb_data$data[[ref_idx]]$XtX, + N = data_update[[idx]]$N, + remain_jk = setdiff(seq_along(model_update[[idx]]$res), data_update[[idx]]$variable_miss), + ref_label = cb_data$data[[ref_idx]]$ref_label + ) + ld_cache[[ld_key]][which(is.na(ld_cache[[ld_key]]))] <- 0 + } + } # -- check if jk_i ~ jk_j - change_each_pair <- matrix(FALSE, nrow = length(jk_each), ncol = length(jk_each)) - for (i in 1:length(jk_each)) { - jk_i <- jk_each[i] - for (j in i:length(jk_each)) { - if (j != i) { - jk_j <- jk_each[j] - change_each_pair[i, j] <- detect_func(idx = i, LD_all, jk_i, jk_j, i, j) - # if jk_i and jk_j are equivalent on dataset i, then we don't need to check dataset j - if ( !change_each_pair[i, j] ){ - change_each_pair[j, i] <- detect_func(idx = j, LD_all, jk_i, jk_j, i, j) - } - } else { - change_each_pair[i, j] <- FALSE - } - } + change_loglike <- t(vapply(seq_along(jk_each), function(idx) { + model_update[[idx]]$change_loglike[jk_each] + }, numeric(length(jk_each)))) + change_ok <- abs(sweep(change_loglike, 1, diag(change_loglike), "-")) <= jk_equiv_loglik + + detected <- matrix(FALSE, nrow = n_pair, ncol = n_pair) + for (ld_key in unique(ld_keys)) { + rows <- which(ld_keys == ld_key) + ld_ok <- abs(ld_cache[[ld_key]][rows, , drop = FALSE]) >= jk_equiv_corr + detected[rows, ] <- change_ok[rows, , drop = FALSE] & ld_ok } - change_each_pair <- change_each_pair + t(change_each_pair) + diag(detected) <- FALSE + change_each_pair <- (detected | t(detected)) * 1 return(change_each_pair) } diff --git a/R/colocboost_init.R b/R/colocboost_init.R index 403e114..a164243 100644 --- a/R/colocboost_init.R +++ b/R/colocboost_init.R @@ -160,7 +160,7 @@ colocboost_init_model <- function(cb_data, "z" = NULL, "learning_rate_init" = learning_rate_init, "stop_thresh" = stop_thresh, - "ld_jk" = c(), + "ld_jk" = list(), "jk" = c(), "scaling_factor" = if (!is.null(cb_data$data[[i]]$N)) (cb_data$data[[i]]$N - 1) else 1, "beta_scaling" = if (!is.null(cb_data$data[[i]]$N)) 1 else 100, diff --git a/R/colocboost_output.R b/R/colocboost_output.R index f98953b..49a9fea 100644 --- a/R/colocboost_output.R +++ b/R/colocboost_output.R @@ -238,7 +238,7 @@ get_robust_colocalization <- function(cb_output, cos_details$cos_purity$min_abs_cor <- as.matrix(cos_details$cos_purity$min_abs_cor)[-remove_idx, -remove_idx, drop = FALSE] cos_details$cos_purity$median_abs_cor <- as.matrix(cos_details$cos_purity$median_abs_cor)[-remove_idx, -remove_idx, drop = FALSE] cos_details$cos_purity$max_abs_cor <- as.matrix(cos_details$cos_purity$max_abs_cor)[-remove_idx, -remove_idx, drop = FALSE] - vcp <- as.vector(1 - apply(1 - do.call(cbind, cos_details$cos_vcp), 1, prod)) + vcp <- as.vector(1 - apply(1 - do.call(cbind, unname(cos_details$cos_vcp)), 1, prod)) names(vcp) <- cb_output$data_info$variables cb_output$vcp <- vcp cb_output$cos_details <- cos_details @@ -323,7 +323,7 @@ get_robust_colocalization <- function(cb_output, use_entropy = use_entropy, residual_correlation = residual_correlation) names(int_weight) <- names(cos_weights) <- colocset_names cos_details$cos_weights <- cos_weights - vcp <- as.vector(1 - apply(1 - do.call(cbind, int_weight), 1, prod)) + vcp <- as.vector(1 - apply(1 - do.call(cbind, unname(int_weight)), 1, prod)) names(vcp) <- cb_output$data_info$variables cb_output$vcp <- vcp cos_details$cos_vcp <- int_weight diff --git a/R/colocboost_plot.R b/R/colocboost_plot.R index ebe1644..8702cf9 100644 --- a/R/colocboost_plot.R +++ b/R/colocboost_plot.R @@ -404,7 +404,7 @@ get_input_plot <- function(cb_output, plot_cos_idx = NULL, coloc_index <- cb_output$cos_details$cos_outcomes$outcome_index # top_variables coloc_hits <- lapply(names(coloc_cos), function(cn) { - p <- grep(cn, rownames(cb_output$cos_details$cos_top_variables)) + p <- which(startsWith(rownames(cb_output$cos_details$cos_top_variables), cn)) cb_output$cos_details$cos_top_variables$top_index[p] }) names(coloc_hits) <- names(coloc_cos) @@ -413,7 +413,7 @@ get_input_plot <- function(cb_output, plot_cos_idx = NULL, cos_vcp <- lapply(1:length(analysis_outcome), function(iy) { pos <- which(sapply(coloc_index, function(idx) iy %in% idx)) if (length(pos) != 0) { - w <- do.call(cbind, cb_output$cos_details$cos_vcp[pos]) + w <- do.call(cbind, unname(cb_output$cos_details$cos_vcp[pos])) return(1 - apply(1 - w, 1, prod)) } else { return(rep(0, length(variables))) @@ -810,4 +810,3 @@ plot_initial <- function(cb_plot_input, y = "log10p", return(args) } - diff --git a/R/colocboost_update.R b/R/colocboost_update.R index 4037214..29ed8b9 100644 --- a/R/colocboost_update.R +++ b/R/colocboost_update.R @@ -12,6 +12,16 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { pos.update <- which(cb_model_para$update_temp$update_status != 0) focal_outcome_idx <- cb_model_para$focal_outcome_idx tau = cb_model_para$tau + ld_jk_cache <- list() + make_ld_jk_cache_key <- function(update_jk, outcome_idx, ref_idx) { + paste( + update_jk, + ref_idx, + cb_data$data[[ref_idx]]$ref_label, + paste(cb_data$data[[outcome_idx]]$variable_miss, collapse = ","), + sep = "|" + ) + } for (i in pos.update) { update_jk <- cb_model_para$update_temp$real_update_jk[i] @@ -20,20 +30,27 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { ########## BEGIN: MAIN CALCULATION ################### # - calucalate LD between update_jk and other variables - if (update_jk %in% unlist(cb_model[[i]]$jk)) { - pos <- which(unlist(cb_model[[i]]$jk) == update_jk) - ld_jk <- cb_model[[i]]$ld_jk[pos, ] - } else { + ld_jk_key <- as.character(update_jk) + ld_jk_cache_key <- make_ld_jk_cache_key(update_jk, i, X_dict) + ld_jk <- cb_model[[i]]$ld_jk[[ld_jk_key]] + if (!is.null(ld_jk)) { + ld_jk_cache[[ld_jk_cache_key]] <- ld_jk + } + if (is.null(ld_jk)) { cb_model[[i]]$jk <- c(cb_model[[i]]$jk, update_jk) - ld_jk <- get_LD_jk(update_jk, - X = cb_data$data[[X_dict]]$X, - XtX = cb_data$data[[X_dict]]$XtX, - N = cb_data$data[[i]]$N, - remain_idx = setdiff(1:cb_model_para$P, cb_data$data[[i]]$variable_miss), - P = cb_model_para$P, - ref_label = cb_data$data[[X_dict]]$ref_label - ) - cb_model[[i]]$ld_jk <- rbind(cb_model[[i]]$ld_jk, ld_jk) + ld_jk <- ld_jk_cache[[ld_jk_cache_key]] + if (is.null(ld_jk)) { + ld_jk <- get_LD_jk(update_jk, + X = cb_data$data[[X_dict]]$X, + XtX = cb_data$data[[X_dict]]$XtX, + N = cb_data$data[[i]]$N, + remain_idx = setdiff(seq_len(cb_model_para$P), cb_data$data[[i]]$variable_miss), + P = cb_model_para$P, + ref_label = cb_data$data[[X_dict]]$ref_label + ) + ld_jk_cache[[ld_jk_cache_key]] <- ld_jk + } + cb_model[[i]]$ld_jk[[ld_jk_key]] <- ld_jk } ld_feature <- sqrt(abs(ld_jk)) @@ -57,7 +74,7 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { x_tmp <- cb_data$data[[X_dict]]$X scaling_factor <- cb_model[[i]]$scaling_factor cov_Xtr <- if (!is.null(x_tmp)) { - t(x_tmp) %*% as.matrix(cb_model[[i]]$res) / scaling_factor + crossprod(x_tmp, cb_model[[i]]$res) / scaling_factor } else { cb_model[[i]]$res / scaling_factor } @@ -99,10 +116,7 @@ colocboost_update <- function(cb_model, cb_model_para, cb_data) { prediction_beta <- cb_data$data[[X_dict]]$X %*% (beta_grad) cb_model[[i]]$res <- cb_model[[i]]$res - step1 * prediction_beta # - profile-loglikelihood - x <- cb_data$data[[X_dict]]$X - y <- cb_data$data[[i]]$Y - beta <- cb_model[[i]]$beta - profile_log <- mean((y - x %*% beta)^2) + profile_log <- mean(cb_model[[i]]$res^2) } else if (!is.null(cb_data$data[[X_dict]]$XtX)) { beta_scaling <- cb_model[[i]]$beta_scaling # - summary statistics @@ -316,7 +330,7 @@ boost_obj_last <- function(cb_data, cb_model, cb_model_para) { x_tmp <- cb_data$data[[X_dict]]$X scaling_factor <- cb_model[[i]]$scaling_factor cov_Xtr <- if (!is.null(x_tmp)) { - t(x_tmp) %*% as.matrix(cb_model[[i]]$res) / scaling_factor + crossprod(x_tmp, cb_model[[i]]$res) / scaling_factor } else { cb_model[[i]]$res / scaling_factor } diff --git a/R/colocboost_utils.R b/R/colocboost_utils.R index 7e956bc..eeab354 100644 --- a/R/colocboost_utils.R +++ b/R/colocboost_utils.R @@ -13,12 +13,33 @@ #' @rdname colocboost_refine_cos #' @keywords cb_refine_cos #' @noRd +.cb_unique_purity_outcomes <- function(cb_data, outcomes) { + if (length(outcomes) == 0L) { + return(integer(0)) + } + + keys <- vapply(outcomes, function(i) { + ref_idx <- cb_data$dict[i] + ref_label <- cb_data$data[[ref_idx]]$ref_label + ref_label <- if (is.null(ref_label)) "" else ref_label + missing_idx <- cb_data$data[[i]]$variable_miss + missing_idx <- if (is.null(missing_idx)) integer(0) else missing_idx + paste(ref_idx, ref_label, paste(missing_idx, collapse = ","), sep = "|") + }, character(1)) + + outcomes[!duplicated(keys)] +} + merge_cos_ucos <- function(cb_obj, out_cos, out_ucos, coverage = 0.95, min_abs_corr = 0.5, tol = 1e-9, median_cos_abs_corr = 0.8) { change_obj_each <- out_ucos$change_obj_each coloc_sets <- out_cos$cos$cos ucos_each <- out_ucos$ucos_each + purity_outcomes <- .cb_unique_purity_outcomes( + cb_obj$cb_data, + seq_len(cb_obj$cb_model_para$L) + ) # - remove overlap between coloc_sets and single_sets is_overlap <- is_highLD <- c() @@ -60,10 +81,13 @@ merge_cos_ucos <- function(cb_obj, out_cos, out_ucos, coverage = 0.95, # - if fine_y not in coloc_y, we check overlap and also min_purity change_obj_coloc <- out_cos$cos$cs_change cset1 <- coloc_sets[[i]] + if (length(intersect(cset1, cset2)) == 0L) { + next + } res <- list() - for (ii in 1:cb_obj$cb_model_para$L) { + for (ii in purity_outcomes) { X_dict <- cb_obj$cb_data$dict[ii] - res[[ii]] <- get_between_purity(cset1, cset2, + res[[length(res) + 1L]] <- get_between_purity(cset1, cset2, X = cb_obj$cb_data$data[[X_dict]]$X, Xcorr = cb_obj$cb_data$data[[X_dict]]$XtX, miss_idx = cb_obj$cb_data$data[[ii]]$variable_miss, @@ -120,6 +144,40 @@ merge_cos_ucos <- function(cb_obj, out_cos, out_ucos, coverage = 0.95, return(ll) } +# Build only uCoS pairs that can merge because they share at least one variant. +.merge_ucos_overlap_pairs <- function(ucos_each) { + empty_pairs <- matrix(integer(0), ncol = 2L, dimnames = list(NULL, c("i", "j"))) + if (length(ucos_each) < 2L) { + return(empty_pairs) + } + + variants <- unlist(ucos_each, use.names = FALSE) + if (length(variants) == 0L) { + return(empty_pairs) + } + + ucos_idx <- rep(seq_along(ucos_each), lengths(ucos_each)) + pairs <- lapply(split(ucos_idx, variants), function(idx) { + idx <- sort(unique(idx)) + if (length(idx) < 2L) { + return(NULL) + } + do.call(rbind, lapply(seq_len(length(idx) - 1L), function(i) { + cbind(idx[i], idx[(i + 1L):length(idx)]) + })) + }) + pairs <- pairs[!vapply(pairs, is.null, logical(1))] + if (length(pairs) == 0L) { + return(empty_pairs) + } + + pairs <- do.call(rbind, pairs) + pairs <- pairs[!duplicated(paste(pairs[, 1L], pairs[, 2L], sep = "\r")), , drop = FALSE] + pairs <- pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE] + colnames(pairs) <- c("i", "j") + pairs +} + #' @importFrom stats na.omit merge_ucos <- function(cb_obj, past_out, min_abs_corr = 0.5, @@ -131,12 +189,37 @@ merge_ucos <- function(cb_obj, past_out, out_cos <- past_out$cos ucos_each <- out_ucos$ucos_each change_obj_each <- out_ucos$change_obj_each + get_top_abs_corr <- function(pos1, pos2, X = NULL, Xcorr = NULL, + miss_idx = NULL, P = NULL, ref_label = "LD") { + if (is.null(Xcorr)) { + value <- suppressWarnings(stats::cor(X[, pos1], X[, pos2])) + } else if (identical(ref_label, "No_ref") || (length(Xcorr) == 1 && Xcorr == 1)) { + value <- 0 + } else { + if (length(miss_idx) != 0) { + pos1 <- match(pos1, setdiff(seq_len(P), miss_idx)) + pos2 <- match(pos2, setdiff(seq_len(P), miss_idx)) + } + if (is.na(pos1) || is.na(pos2)) { + value <- 0 + } else if (identical(ref_label, "X_ref")) { + value <- suppressWarnings(stats::cor(Xcorr[, pos1], Xcorr[, pos2])) + } else { + value <- Xcorr[pos1, pos2] + } + } + if (is.na(value)) value <- 0 + abs(value) + } # calculate between purity ncsets <- length(ucos_each) min_between <- max_between <- ave_between <- matrix(0, nrow = ncsets, ncol = ncsets) - for (i.between in 1:(ncsets - 1)) { - for (j.between in (i.between + 1):ncsets) { + overlap_pairs <- .merge_ucos_overlap_pairs(ucos_each) + if (nrow(overlap_pairs) > 0L) { + for (pair_idx in seq_len(nrow(overlap_pairs))) { + i.between <- overlap_pairs[pair_idx, 1L] + j.between <- overlap_pairs[pair_idx, 2L] cset1 <- ucos_each[[i.between]] cset2 <- ucos_each[[j.between]] y.i <- out_ucos$ucos_outcome[i.between] @@ -145,6 +228,22 @@ merge_ucos <- function(cb_obj, past_out, next } yy <- c(y.i, y.j) + # Top-top LD is one element of the full between-set LD matrix; if it + # cannot pass the merge cutoff, the full min-between check cannot pass. + top_abs_corr <- vapply(yy, function(ii) { + X_dict <- cb_obj$cb_data$dict[ii] + get_top_abs_corr(cset1[1], cset2[1], + X = cb_obj$cb_data$data[[X_dict]]$X, + Xcorr = cb_obj$cb_data$data[[X_dict]]$XtX, + miss_idx = cb_obj$cb_data$data[[ii]]$variable_miss, + P = cb_obj$cb_model_para$P, + ref_label = cb_obj$cb_data$data[[X_dict]]$ref_label + ) + }, numeric(1)) + top_abs_corr <- if (min_abs_corr == 0) min(top_abs_corr) else max(top_abs_corr) + if (top_abs_corr <= median_cos_abs_corr) { + next + } res <- list() flag <- 1 for (ii in yy) { @@ -662,7 +761,7 @@ get_cos_details <- function(cb_obj, coloc_out, data_info = NULL) { coloc_out$purity <- purity_all[is_pure,,drop = FALSE] colocset_names <- colocset_names[is_pure] } - vcp <- as.vector(1 - apply(1 - do.call(cbind, int_weight), 1, prod)) + vcp <- as.vector(1 - apply(1 - do.call(cbind, unname(int_weight)), 1, prod)) names(vcp) <- data_info$variables @@ -851,7 +950,7 @@ get_full_output <- function(cb_obj, past_out = NULL, variables = NULL, cb_output cb$change_loglike <- cb$change_loglike[ordered] cb$correlation <- as.numeric(cb$correlation[ordered]) cb$z <- as.numeric(cb$z[ordered]) - cb$ld_jk <- cb$ld_jk[, ordered] + cb$ld_jk <- lapply(cb$ld_jk, function(x) x[ordered]) cb$z_univariate <- as.numeric(cb$z_univariate[ordered]) cb$beta_hat <- as.numeric(cb$beta_hat[ordered]) cb$multi_correction <- as.numeric(cb$multi_correction[ordered]) @@ -1094,4 +1193,3 @@ pseudo_inverse <- function(mat) { diag(1 / eig$values[1:keep], keep, keep) %*% t(eig$vectors[, 1:keep, drop = FALSE]) } - diff --git a/R/colocboost_workhorse.R b/R/colocboost_workhorse.R index 2f0a4f7..a709b79 100644 --- a/R/colocboost_workhorse.R +++ b/R/colocboost_workhorse.R @@ -254,6 +254,7 @@ colocboost_workhorse <- function(cb_data, cb_model_para$coveraged <- FALSE } + cb_model_para <- .cb_trim_update_history(cb_model_para) # -- remove redundant parameters rm_elements <- c("update_temp", "update_y") diff --git a/_pkgdown.yml b/_pkgdown.yml index b180b8c..b050857 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -27,16 +27,22 @@ articles: - Disease_Prioritized_Colocalization - title: Interpretation and Visualization - desc: "Tutorials on how to interpret and visualize the output from ColocBoost." + desc: "Tutorials on how to interpret and visualize ColocBoost results." contents: - Interpret_ColocBoost_Output - Visualization_ColocBoost_Output - - title: Advanced Topics - desc: "Advanced topics and special cases in colocalization and fine-mapping analysis." + - title: Colocalization Concepts and Methodology + desc: "Educational tutorials on the statistical concepts, assumptions, and algorithms underlying ColocBoost." contents: - - Partial_Overlap_Variants + - Conceptual_Multi_Trait_Colocalization + - Advanced_Colocalization_Scenarios - ColocBoost_Update + + - title: Advanced Workflows and Special Cases + desc: "Practical guidance for specialized data structures, diagnostic scenarios, and extended ColocBoost workflows." + contents: + - Partial_Overlap_Variants - ColocBoost_Wrapper_Pipeline - LD_Free_Colocalization - Ambiguous_Colocalization diff --git a/cran-comments.md b/cran-comments.md index fae0132..ba0f099 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,47 +1,30 @@ -## colocboost 1.0.9 release comments +## colocboost 1.0.10 release comments -This is a CRAN-requested patch update to colocboost 1.0.8. +This is a CRAN-requested patch update to colocboost 1.0.10. This patch includes: -* Fixed the CRAN-reported macOS arm64 test issue in the uCoS robustness tests. +* Fixed the CRAN-reported MKL test issue in the CoS robustness tests. The tests now use a stronger simulation setting and aligned inputs for `get_robust_ucos()` and `get_ucos_evidence()`, so they no longer depend on weak or platform-sensitive simulated signals. -No R source code or package dependency changes were made for this CRAN-requested -patch. +* R source code changes were limited to targeted computational optimizations and robustness fixes. No package dependencies were changed. -## CRAN-requested macOS arm64 fix -CRAN reported test failures for colocboost 1.0.8 on macOS arm64 and requested a -correction before 2026-06-21. The failing test checked robust trait-specific -uncolocalized event filtering and evidence calculation. The issue has been -addressed by strengthening the simulation used in the test suite and by using -matched `cb_obj` and `cb_res` inputs for the evidence check. +## Advanced documentation and tutorials -CRAN also pointed to the M1mac check service for arm64 issues: -https://www.stats.ox.ac.uk/pub/bdr/M1mac/README.txt. CRAN noted that this -service runs a much older OS/toolchain and that toolchain differences often -matter. This submission therefore fixes the test design itself, rather than -relying on a platform-specific workaround. +* Added conceptual and advanced-scenario vignettes to help users understand multi-trait colocalization events and interpret CoS, VCP, and NPC. +* Added practical illustrations of multiple causal variants and weaker disease GWAS signals, together with supporting figures and guidance. ## R CMD check results There is one NOTE about installed package size: * checking installed package size ... NOTE - installed size is 5.0 MB - sub-directories of 1Mb or more: - data 2.0 MB - doc 1.9 MB + installed size is 7.1Mb + sub-directories of 1Mb or more: + data 2.0Mb + doc 4.1Mb -This NOTE is expected. The installed size is mainly due to reduced example datasets and rendered vignettes with figures. These files are kept to make the tutorials reproducible and self-contained for multi-trait colocalization workflows. No external data are downloaded during examples or vignette rendering. - -## Previous comments - -* This package implements methods described in our paper "ColocBoost" - (Cao et al., 2025), now cited in DESCRIPTION. -* Previous CRAN-requested fixes addressed tarball size, user option handling, - LICENSE metadata, and accepted package/domain terms in inst/WORDLIST. -* The examples and vignettes use small datasets to avoid long check times. +This NOTE is expected. The increase from the previous release is primarily due to two substantive documentation vignettes that introduce the conceptual framework for multi-trait colocalization and provide tutorials on advanced analysis scenarios in ColocBoost. These vignettes include six high-level summary figures: three explain the key conceptual definitions and their relationships to existing methods, and three provide empirical illustrations based on representative simulation results. The figures are embedded in the rendered, self-contained vignettes in the `doc` directory to help users understand the methodology and interpret ColocBoost results. The `data` directory contains reduced example datasets used in reproducible tutorials. No external data are downloaded during examples or vignette rendering. diff --git a/inst/WORDLIST b/inst/WORDLIST index a11345d..b8e2940 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -8,8 +8,11 @@ SuSiE # Sum of Single Effects regression model # Statistical and Genetic Terms bim # PLINK BIM variant information file +combinatorial # Describing combinations in a hypothesis space +combinatorially # Growing according to combinations eQTL # Expression Quantitative Trait Loci GWAS # Genome-Wide Association Study +heritability # Proportion of phenotypic variance attributable to genetic variation INDELs # Insertions and Deletions LD # Linkage Disequilibrium ld # Linkage Disequilibrium @@ -51,6 +54,10 @@ YI # Yang I. Li et # and more al # and more +# Journal Abbreviations +JRSS # Journal of the Royal Statistical Society +PLOS # Public Library of Science + # Technical Terms changelog # Release change log cis # Referring to nearby location of a regulatory element @@ -70,11 +77,13 @@ pos # Position in genome precomputed # Computed in advance precomputes # Computes in advance probabilistically # Based on probability theory +priori # Used in the phrase a priori qc # Quality Control rcond # Reciprocal condition number reconciliate # Process of resolving discrepancies repo # Repository rss # Residual Sum of Squares +scalability # Ability to handle increasing problem size subsampled # Analyzed using data subsets uncolocalized # Not showing colocalization uCoS # Trait-specific (uncolocalization) confidence set @@ -91,6 +100,8 @@ colocalize # Having undergone colocalization analysis colocalized # Having undergone colocalization analysis CoS # Colocalization confidence set in our proposed ColocBoost method NPC # Normalization probability of colocalization in our proposed ColocBoost method +NPUC # Normalized probability of uncolocalization +PPFC # Posterior probability of full colocalization PIPs # Posterior Inclusion Probabilities SEL # Single-effect learner in our proposed ColocBoost method uS # The number of uCoS @@ -99,3 +110,12 @@ Pre # Before pre # Before jk # Index used in ColocBoost nd # Second + +# LaTeX Commands +ge +geq +ldots +mathrm +neq +substack +underbrace diff --git a/tests/testthat/test_Xref.R b/tests/testthat/test_Xref.R index b5cb8a6..c3af41c 100644 --- a/tests/testthat/test_Xref.R +++ b/tests/testthat/test_Xref.R @@ -25,7 +25,7 @@ generate_xref_test_data <- function(n = 200, n_ref = 50, p = 30, L = 2, seed = 4 true_beta <- matrix(0, p, L) true_beta[5, 1] <- 0.7 # SNP5 affects trait 1 true_beta[5, 2] <- 0.6 # SNP5 also affects trait 2 (colocalized) - true_beta[20, 2] <- 0.5 # SNP20 only affects trait 2 + true_beta[20, 2] <- 1.0 # SNP20 only affects trait 2 # Generate Y with some noise Y <- matrix(0, n, L) @@ -498,8 +498,9 @@ test_that("X_ref model has XtX_beta_cache in diagnostic output", { # ============================================================================ test_that("purity functions dispatch correctly for X_ref", { - get_purity <- get("get_purity", envir = asNamespace("colocboost")) - get_between_purity <- get("get_between_purity", envir = asNamespace("colocboost")) + source_env <- environment(colocboost) + get_purity <- get("get_purity", envir = source_env) + get_between_purity <- get("get_between_purity", envir = source_env) set.seed(42) n_ref <- 50 @@ -586,11 +587,14 @@ test_that("get_robust_ucos works with X_ref results", { sumstat = test_data$sumstat, X_ref = test_data$X_ref, M = 10, - output_level = 2 + output_level = 2, + pvalue_cutoff = NULL, + cos_npc_cutoff = 0, + npc_outcome_cutoff = 0 ) })) - skip_if(is.null(result$ucos_details), "No ucos detected") + expect_false(is.null(result$ucos_details)) expect_error( suppressMessages( @@ -728,6 +732,44 @@ test_that("get_LD_jk and get_LD_jk1_jk2 dispatch correctly for X_ref", { expect_equal(ld_pair_noref, 0) }) +test_that("colocboost stores update-jk LD cache as a named list", { + set.seed(42) + n <- 80 + p <- 25 + sigma <- 0.8^abs(outer(1:p, 1:p, "-")) + X <- MASS::mvrnorm(n, rep(0, p), sigma) + colnames(X) <- paste0("SNP", seq_len(p)) + Y <- cbind(X[, 5] * 0.8 + rnorm(n), X[, 5] * 0.7 + rnorm(n)) + + sumstat <- lapply(seq_len(ncol(Y)), function(i) { + beta <- se <- z <- numeric(p) + for (j in seq_len(p)) { + fit <- summary(lm(Y[, i] ~ X[, j]))$coef + beta[j] <- fit[2, 1] + se[j] <- fit[2, 2] + z[j] <- beta[j] / se[j] + } + data.frame(beta = beta, sebeta = se, z = z, n = n, variant = colnames(X)) + }) + + suppressWarnings(suppressMessages({ + result <- colocboost( + sumstat = sumstat, + LD = cor(X), + M = 15, + output_level = 3, + stop_thresh = 0 + ) + })) + + cb_model <- result$diagnostic_details$cb_model + for (model in cb_model) { + expect_type(model$ld_jk, "list") + expect_true(all(nzchar(names(model$ld_jk)))) + expect_equal(length(model$ld_jk), length(unique(model$jk))) + } +}) + # ============================================================================ # Test 22: ref_label is never NULL in internal cb_data after processing @@ -776,4 +818,4 @@ test_that("ref_label is always explicitly set in cb_data, never NULL", { for (i in seq_along(cb_data_ind$data)) { expect_equal(cb_data_ind$data[[i]]$ref_label, "individual") } -}) \ No newline at end of file +}) diff --git a/tests/testthat/test_large_outcome_inference.R b/tests/testthat/test_large_outcome_inference.R new file mode 100644 index 0000000..f511e8b --- /dev/null +++ b/tests/testthat/test_large_outcome_inference.R @@ -0,0 +1,156 @@ +library(testthat) + +test_that("colocalization candidate grouping supports more than 10000-byte signatures", { + n_outcomes <- 4589L + update <- matrix(-1L, nrow = n_outcomes, ncol = 3L) + update[seq_len(100L), 1:2] <- 1L + update[101:200, 3] <- 1L + + old_signature <- paste0(update[, 1], collapse = ",") + expect_gt(nchar(old_signature, type = "bytes"), 10000L) + + result <- .group_coloc_candidates(update, seq_len(ncol(update))) + + expect_equal(unname(result), list(c(1L, 2L), 3L)) + expect_true(all(nchar(names(result), type = "bytes") < 10000L)) +}) + +test_that("CoS details support more than 10000-byte outcome-set names", { + n_outcomes <- 4589L + cos_name <- paste0("cos1:", paste0("y", seq_len(n_outcomes), collapse = "_")) + expect_gt(nchar(cos_name, type = "bytes"), 10000L) + + local_mocked_bindings( + get_cos_evidence = function(...) { + list( + normalization_evidence = list(data.frame(npc_outcome = 1)), + npc = 1 + ) + }, + get_integrated_weight = function(...) c(1, 0), + get_in_cos = function(...) list(1L), + get_purity = function(...) c(1, 1, 1), + .package = "colocboost" + ) + + data_entry <- list( + X = matrix(0, nrow = 2, ncol = 2), + XtX = NULL, + N = 2L, + variable_miss = integer(0), + ref_label = "individual" + ) + cb_obj <- list( + cb_model_para = list( + outcome_names = paste0("trait", seq_len(n_outcomes)), + variables = c("v1", "v2"), + P = 2L, + coverage = 0.95, + min_abs_corr = 0.5, + median_abs_corr = NULL, + n_purity = 100L, + weight_fudge_factor = 1.5, + use_entropy = FALSE, + residual_correlation = NULL + ), + cb_data = list( + dict = rep(1L, n_outcomes), + data = rep(list(data_entry), n_outcomes) + ) + ) + coloc_out <- list( + cos = list(cos1 = 1L), + coloc_outcomes = list(seq_len(n_outcomes)), + avWeight = list(matrix(1, nrow = 2L, ncol = n_outcomes)), + cs_change = matrix(1, nrow = 1L, ncol = n_outcomes) + ) + data_info <- list(variables = c("v1", "v2")) + + expect_no_error(result <- get_cos_details(cb_obj, coloc_out, data_info)) + expect_equal(names(result$cos_results$cos$cos_index), cos_name) +}) + +test_that("robust CoS filtering supports more than 10000-byte outcome-set names", { + n_outcomes <- 4589L + outcome_idx <- seq_len(n_outcomes) + outcome_names <- paste0("trait", outcome_idx) + weights <- matrix(rep(c(0.9, 0.1), n_outcomes), nrow = 2L) + colnames(weights) <- paste0("outcome", outcome_idx) + outcome_npc <- data.frame( + relative_logLR = rep(1, n_outcomes), + npc_outcome = rep(1, n_outcomes), + outcomes_index = outcome_idx, + row.names = outcome_names + ) + purity <- matrix( + 1, + 2L, + 2L, + dimnames = list(c("cos1", "cos2"), c("cos1", "cos2")) + ) + + cb_output <- list( + data_info = list( + variables = c("v1", "v2"), + z = rep(list(c(10, 0)), n_outcomes), + coef = rep(list(c(1, 0)), n_outcomes), + n_outcomes = n_outcomes, + n_variables = 2L, + outcome_info = data.frame( + outcome_names = outcome_names, + is_focal = rep(FALSE, n_outcomes) + ) + ), + vcp = c(v1 = 0.9, v2 = 0.1), + cos_details = list( + cos_top_variables = data.frame( + top_index = c(1L, 2L), + top_variables = c("v1", "v2"), + row.names = c("cos1", "cos2") + ), + cos = list( + cos_index = list(cos1 = 1L, cos2 = 2L), + cos_variables = list(cos1 = "v1", cos2 = "v2") + ), + cos_outcomes = list( + outcome_index = list(cos1 = outcome_idx, cos2 = 1L), + outcome_name = list(cos1 = outcome_names, cos2 = outcome_names[1]) + ), + cos_outcomes_npc = list( + cos1 = outcome_npc, + cos2 = data.frame( + relative_logLR = 0, + npc_outcome = 0, + outcomes_index = 1L, + row.names = outcome_names[1] + ) + ), + cos_vcp = list(cos1 = c(0.9, 0.1), cos2 = c(0.1, 0.9)), + cos_weights = list( + cos1 = weights, + cos2 = matrix(c(0.1, 0.9), ncol = 1L, dimnames = list(NULL, "outcome1")) + ), + cos_npc = c(cos1 = 1, cos2 = 1), + cos_min_npc_outcome = c(cos1 = 1, cos2 = 1), + cos_purity = list( + min_abs_cor = purity, + median_abs_cor = purity, + max_abs_cor = purity + ) + ) + ) + class(cb_output) <- "colocboost" + + expect_no_error( + result <- get_robust_colocalization( + cb_output, + cos_npc_cutoff = 0, + npc_outcome_cutoff = 0, + pvalue_cutoff = 1 + ) + ) + expect_equal(length(result$cos_details$cos_outcomes$outcome_index[[1]]), n_outcomes) + expect_gt(nchar(names(result$cos_details$cos$cos_index)[1], type = "bytes"), 10000L) + expect_no_error(plot_input <- get_input_plot(result)) + expect_equal(length(plot_input$cos_vcp), n_outcomes) +}) diff --git a/tests/testthat/test_optimization_phase1.R b/tests/testthat/test_optimization_phase1.R index a82716a..006aadd 100644 --- a/tests/testthat/test_optimization_phase1.R +++ b/tests/testthat/test_optimization_phase1.R @@ -1,5 +1,7 @@ library(testthat) +colocboost_test_env <- environment(colocboost) + # ---- Shared test data generators ---- generate_test_data_opt <- function(n = 200, p = 30, L = 2, seed = 42) { @@ -253,6 +255,94 @@ test_that("optimization works with 3 outcomes", { expect_equal(result$data_info$n_outcomes, 3) }) +test_that("post-assembly purity optimization handles non-contiguous outcomes", { + X <- matrix(rep(seq_len(8), 4), nrow = 8, ncol = 4) + colnames(X) <- paste0("SNP", 1:4) + cb_obj <- list( + cb_data = list( + dict = c(1, 1, 3), + data = list( + list(X = X, XtX = NULL, variable_miss = integer(0), ref_label = "individual", N = nrow(X)), + list(X = X, XtX = NULL, variable_miss = integer(0), ref_label = "individual", N = nrow(X)), + list(X = X, XtX = NULL, variable_miss = integer(0), ref_label = "individual", N = nrow(X)) + ) + ), + cb_model_para = list(L = 3, P = ncol(X)) + ) + out_cos <- list(cos = list( + cos = list(c(1, 2)), + avWeight = list(matrix(c(0.7, 0.3, 0, 0), ncol = 1, + dimnames = list(colnames(X), "outcome1"))), + coloc_outcomes = list(1), + cs_change = matrix(0, nrow = 1, ncol = 3) + )) + out_ucos <- list( + ucos_each = list(c(2, 3)), + ucos_outcome = 3, + avW_ucos_each = matrix(c(0, 0.7, 0.3, 0), ncol = 1, + dimnames = list(colnames(X), "Y3")), + change_obj_each = matrix(0, nrow = 1, ncol = 3), + purity_each = matrix(1, nrow = 1, ncol = 3) + ) + + expect_error({ + merged <- merge_cos_ucos(cb_obj, out_cos, out_ucos, + coverage = 0.8, median_cos_abs_corr = 0.8) + expect_equal(merged$cos$cos$coloc_outcomes[[1]], c(1, 3)) + }, NA) +}) + +test_that("colocboost_assemble_cos handles non-contiguous purity outcomes", { + p <- 4 + L <- 3 + LD <- matrix(0.95, nrow = p, ncol = p) + diag(LD) <- 1 + cb_obj <- list( + cb_data = list( + dict = c(1L, 1L, 3L), + data = lapply(seq_len(L), function(i) { + list(X = NULL, XtX = LD, variable_miss = integer(0), + ref_label = "LD", N = 20) + }) + ), + cb_model = list( + list(weights_path = rbind(c(0.6, 0.4, 0, 0), c(0.6, 0, 0.4, 0))), + list(weights_path = rbind(c(0.5, 0.5, 0, 0))), + list(weights_path = rbind(c(0.5, 0, 0.5, 0))) + ), + cb_model_para = list( + L = L, + P = p, + update_status = rbind(c(1, 1), c(1, 0), c(0, 1)) + ) + ) + class(cb_obj) <- "colocboost" + + local_mocked_bindings( + check_null_post = function(cb_obj, coloc_sets_temp, ...) { + list( + cs_change = matrix(1, nrow = length(coloc_sets_temp), ncol = L), + is_non_null = seq_along(coloc_sets_temp) + ) + }, + get_between_purity = function(...) { + c(min_abs_cor = 0.9, max_abs_cor = 1, median_abs_cor = 0.9) + }, + .package = "colocboost" + ) + + expect_error({ + assembled <- colocboost_assemble_cos( + cb_obj, + coverage = 0.8, + min_abs_corr = 0, + median_cos_abs_corr = 0.8 + ) + expect_equal(length(assembled$cos$cos), 1L) + expect_equal(assembled$cos$coloc_outcomes[[1]], c(1L, 2L, 3L)) + }, NA) +}) + # ---- Test: Missing variants ---- test_that("optimization handles missing variants correctly", { @@ -416,3 +506,544 @@ test_that("optimization does not break convergence behavior", { # Model should have converged (model_info should exist) expect_type(result$model_info, "list") }) + +test_that("pairwise jk checks reuse LD calculation for shared reference data", { + set.seed(202) + n <- 80 + p <- 20 + n_outcomes <- 6 + X <- matrix(rnorm(n * p), nrow = n) + colnames(X) <- paste0("SNP", seq_len(p)) + jk_each <- c(3, 6, 9, 12, 15, 18) + pos.update <- seq_len(n_outcomes) + 1L + cb_data <- list( + data = c( + list(list(X = X, XtX = NULL, ref_label = "individual")), + lapply(seq_len(n_outcomes), function(i) { + list(N = n, variable_miss = integer(0)) + }) + ) + ) + X_dict <- rep(1L, n_outcomes) + model_update <- lapply(seq_len(n_outcomes), function(i) { + change_loglike <- seq_len(p) / p + i * 1e-4 + change_loglike[jk_each] <- change_loglike[jk_each] + rnorm(n_outcomes, sd = 1e-5) + list(change_loglike = change_loglike, res = numeric(p)) + }) + + ns <- colocboost_test_env + original_get_cormat <- get("get_cormat", envir = ns) + call_count <- 0L + local_mocked_bindings( + get_cormat = function(...) { + call_count <<- call_count + 1L + original_get_cormat(...) + }, + .package = "colocboost" + ) + + pair_check <- get("check_pair_jkeach", envir = ns) + res <- pair_check(jk_each, pos.update, model_update, cb_data, X_dict) + + expect_equal(call_count, 1L) + expect_equal(dim(res), c(n_outcomes, n_outcomes)) + expect_true(all(res == t(res))) + expect_true(all(diag(res) == 0)) +}) + +make_shared_update_fixture <- function(n = 60, p = 25, n_outcomes = 5, update_jk = 4, seed = 303) { + set.seed(seed) + X <- matrix(rnorm(n * p), nrow = n) + X <- scale(X) + colnames(X) <- paste0("SNP", seq_len(p)) + Y <- matrix(rnorm(n * n_outcomes), nrow = n) + Y <- scale(Y) + + cb_data <- list( + data = lapply(seq_len(n_outcomes), function(i) { + list( + X = if (i == 1L) X else NULL, + Y = as.matrix(Y[, i]), + N = n, + variable_miss = integer(0), + ref_label = "individual" + ) + }), + dict = rep(1L, n_outcomes), + variable.names = colnames(X) + ) + class(cb_data) <- "colocboost" + + cb_model <- lapply(seq_len(n_outcomes), function(i) { + correlation <- rnorm(p) + list( + res = as.matrix(Y[, i]), + beta = rep(0, p), + weights_path = list(), + profile_loglike_each = mean(Y[, i]^2), + obj_path = 999999, + obj_single = 999999, + change_loglike = abs(rnorm(p)), + correlation = correlation, + z = correlation, + learning_rate_init = 0.01, + stop_thresh = 1e-6, + ld_jk = list(), + jk = integer(0), + scaling_factor = n - 1, + beta_scaling = 1, + XtX_beta_cache = NULL + ) + }) + class(cb_model) <- "colocboost" + + cb_model_para <- list( + update_temp = list( + update_status = rep(1, n_outcomes), + real_update_jk = rep(update_jk, n_outcomes) + ), + focal_outcome_idx = NULL, + tau = 0.01, + lambda = 0.5, + lambda_focal_outcome = 1, + func_simplex = "LD_z2z", + LD_free = FALSE, + dynamic_learning_rate = FALSE, + learning_rate_decay = 1, + P = p + ) + class(cb_model_para) <- "colocboost" + + list(cb_model = cb_model, cb_model_para = cb_model_para, cb_data = cb_data) +} + +test_that("colocboost_update reuses LD_jk for outcomes sharing reference data", { + fixture <- make_shared_update_fixture(n_outcomes = 5, update_jk = 6) + + ns <- colocboost_test_env + original_get_LD_jk <- get("get_LD_jk", envir = ns) + call_count <- 0L + local_mocked_bindings( + get_LD_jk = function(...) { + call_count <<- call_count + 1L + original_get_LD_jk(...) + }, + .package = "colocboost" + ) + + updated <- colocboost_update(fixture$cb_model, fixture$cb_model_para, fixture$cb_data) + + expect_equal(call_count, 1L) + expect_true(all(vapply(updated, function(model) "6" %in% names(model$ld_jk), logical(1)))) +}) + +test_that("individual profile log equals explicit and residual calculations", { + fixture <- make_shared_update_fixture(n_outcomes = 1, update_jk = 5) + fixture$cb_model[[1]]$ld_jk[["5"]] <- rep(1, fixture$cb_model_para$P) + matmul_count <- 0L + fixture$cb_data$data[[1]]$X <- structure( + fixture$cb_data$data[[1]]$X, + class = c("counted_matrix", "matrix") + ) + assign("%*%.counted_matrix", function(x, y) { + matmul_count <<- matmul_count + 1L + NextMethod() + }, envir = .GlobalEnv) + on.exit(rm("%*%.counted_matrix", envir = .GlobalEnv), add = TRUE) + + updated <- colocboost_update(fixture$cb_model, fixture$cb_model_para, fixture$cb_data) + profile_log <- tail(updated[[1]]$profile_loglike_each, n = 1) + X_plain <- matrix(fixture$cb_data$data[[1]]$X, nrow = nrow(fixture$cb_data$data[[1]]$X)) + explicit_profile <- mean((fixture$cb_data$data[[1]]$Y - X_plain %*% updated[[1]]$beta)^2) + residual_profile <- mean(updated[[1]]$res^2) + + expect_equal(matmul_count, 1L) + expect_equal( + as.numeric(profile_log), + explicit_profile, + tolerance = 1e-12 + ) + expect_equal( + as.numeric(profile_log), + residual_profile, + tolerance = 1e-12 + ) + expect_equal( + explicit_profile, + residual_profile, + tolerance = 1e-12 + ) +}) + +test_that("merge_ucos skips between-purity checks for disjoint uCoS pairs", { + set.seed(404) + p <- 30 + n <- 20 + n_outcomes <- 4 + LD <- diag(p) + LD[1, 2] <- LD[2, 1] <- 0.9 + cb_obj <- list( + cb_data = list( + data = lapply(seq_len(n_outcomes), function(i) { + list(X = NULL, XtX = LD, variable_miss = integer(0), ref_label = "LD", N = n) + }), + dict = seq_len(n_outcomes) + ), + cb_model_para = list(P = p, L = n_outcomes) + ) + class(cb_obj) <- "colocboost" + + ucos_each <- list(c(1, 2), c(10, 11), c(2, 3), c(20, 21)) + names(ucos_each) <- paste0("sets:Y", seq_len(n_outcomes), ":ucos1") + avW <- matrix(runif(p * length(ucos_each)), nrow = p) + colnames(avW) <- names(ucos_each) + past_out <- list( + ucos = list( + ucos_each = ucos_each, + avW_ucos_each = avW, + change_obj_each = matrix(0.1, nrow = length(ucos_each), ncol = n_outcomes), + purity_each = matrix(1, nrow = length(ucos_each), ncol = 3), + ucos_outcome = seq_len(n_outcomes) + ), + cos = list(cos = list()) + ) + + ns <- colocboost_test_env + between_calls <- 0L + local_mocked_bindings( + get_between_purity = function(pos1, pos2, ...) { + if (length(intersect(pos1, pos2)) == 0) { + stop("disjoint uCoS pair should not require between-purity") + } + between_calls <<- between_calls + 1L + c(min_abs_cor = 0.9, max_abs_cor = 1, median_abs_cor = 0.9) + }, + get_purity = function(...) c(1, 1, 1), + .package = "colocboost" + ) + + result <- get("merge_ucos", envir = ns)( + cb_obj, past_out, + min_abs_corr = 0.5, + median_cos_abs_corr = 0.8 + ) + + expect_equal(between_calls, 2L) + expect_equal(length(result$cos$cos$cos), 1L) + expect_equal(length(result$ucos$ucos_each), 2L) +}) + +test_that("merge_ucos skips full between-purity when top variants cannot pass merge cutoff", { + p <- 10 + n <- 20 + LD <- diag(p) + LD[1, 3] <- LD[3, 1] <- 0.1 + cb_obj <- list( + cb_data = list( + data = lapply(seq_len(2), function(i) { + list(X = NULL, XtX = LD, variable_miss = integer(0), ref_label = "LD", N = n) + }), + dict = seq_len(2) + ), + cb_model_para = list(P = p, L = 2) + ) + class(cb_obj) <- "colocboost" + + ucos_each <- list(c(1, 2), c(3, 2)) + names(ucos_each) <- c("sets:Y1:ucos1", "sets:Y2:ucos1") + avW <- matrix(runif(p * length(ucos_each)), nrow = p) + colnames(avW) <- names(ucos_each) + past_out <- list( + ucos = list( + ucos_each = ucos_each, + avW_ucos_each = avW, + change_obj_each = matrix(0.1, nrow = length(ucos_each), ncol = 2), + purity_each = matrix(1, nrow = length(ucos_each), ncol = 3), + ucos_outcome = seq_len(2) + ), + cos = list(cos = list()) + ) + + ns <- colocboost_test_env + local_mocked_bindings( + get_between_purity = function(...) { + stop("top-variant prefilter should skip full between-purity") + }, + .package = "colocboost" + ) + + result <- get("merge_ucos", envir = ns)( + cb_obj, past_out, + min_abs_corr = 0.5, + median_cos_abs_corr = 0.8 + ) + + expect_equal(length(result$cos$cos$cos), 0L) + expect_equal(length(result$ucos$ucos_each), 2L) +}) + +test_that("overlap candidate pairs match pairwise intersect scan", { + ucos_each <- list( + c(1, 2), + c(3, 4), + c(2, 5), + c(5, 6), + 7, + c(2, 6) + ) + reference <- do.call(rbind, lapply(seq_len(length(ucos_each) - 1L), function(i) { + do.call(rbind, lapply((i + 1L):length(ucos_each), function(j) { + if (length(intersect(ucos_each[[i]], ucos_each[[j]])) == 0) { + return(NULL) + } + c(i, j) + })) + })) + colnames(reference) <- c("i", "j") + + candidate_pairs <- get(".merge_ucos_overlap_pairs", envir = colocboost_test_env)(ucos_each) + + expect_equal(candidate_pairs, reference) +}) + +test_that("merge_ucos candidate pairs preserve pairwise intersect output", { + set.seed(505) + p <- 30 + n_outcomes <- 5 + LD <- diag(p) + LD[1, 2] <- LD[2, 1] <- 0.9 + LD[1, 20] <- LD[20, 1] <- 0.95 + cb_obj <- list( + cb_data = list( + data = lapply(seq_len(n_outcomes), function(i) { + list(X = NULL, XtX = LD, variable_miss = integer(0), ref_label = "LD", N = 20) + }), + dict = seq_len(n_outcomes) + ), + cb_model_para = list(P = p, L = n_outcomes) + ) + class(cb_obj) <- "colocboost" + + ucos_each <- list(c(1, 2), c(10, 11), c(2, 3), c(20, 21), c(3, 4)) + names(ucos_each) <- paste0("sets:Y", seq_len(n_outcomes), ":ucos1") + avW <- matrix(runif(p * length(ucos_each)), nrow = p) + colnames(avW) <- names(ucos_each) + past_out <- list( + ucos = list( + ucos_each = ucos_each, + avW_ucos_each = avW, + change_obj_each = matrix(0.1, nrow = length(ucos_each), ncol = n_outcomes), + purity_each = matrix(1, nrow = length(ucos_each), ncol = 3), + ucos_outcome = seq_len(n_outcomes) + ), + cos = list(cos = list()) + ) + + ns <- colocboost_test_env + local_mocked_bindings( + get_between_purity = function(...) { + c(min_abs_cor = 0.9, max_abs_cor = 1, median_abs_cor = 0.9) + }, + get_purity = function(...) c(1, 1, 1), + .package = "colocboost" + ) + + pairwise_pairs <- function(ucos_each) { + if (length(ucos_each) < 2L) { + return(matrix(integer(0), ncol = 2L, dimnames = list(NULL, c("i", "j")))) + } + pairs <- do.call(rbind, lapply(seq_len(length(ucos_each) - 1L), function(i) { + do.call(rbind, lapply((i + 1L):length(ucos_each), function(j) { + if (length(intersect(ucos_each[[i]], ucos_each[[j]])) == 0) { + return(NULL) + } + c(i, j) + })) + })) + if (is.null(pairs)) { + return(matrix(integer(0), ncol = 2L, dimnames = list(NULL, c("i", "j")))) + } + colnames(pairs) <- c("i", "j") + pairs + } + + candidate_result <- get("merge_ucos", envir = ns)( + cb_obj, past_out, + min_abs_corr = 0.5, + median_cos_abs_corr = 0.8 + ) + + local_mocked_bindings( + .merge_ucos_overlap_pairs = pairwise_pairs, + .package = "colocboost" + ) + pairwise_result <- get("merge_ucos", envir = ns)( + cb_obj, past_out, + min_abs_corr = 0.5, + median_cos_abs_corr = 0.8 + ) + + expect_equal(candidate_result, pairwise_result) + expect_equal(length(candidate_result$cos$cos$cos), 1L) + expect_equal(length(candidate_result$ucos$ucos_each), 3L) +}) + +test_that("duplicate purity contexts are collapsed for post-assembly checks", { + cb_data <- list( + data = lapply(seq_len(4), function(i) { + list( + X = NULL, + XtX = diag(5), + variable_miss = if (i <= 2) integer(0) else 5L, + ref_label = "LD" + ) + }), + dict = rep(1L, 4) + ) + + unique_outcomes <- get(".cb_unique_purity_outcomes", envir = colocboost_test_env)( + cb_data, + seq_len(4) + ) + + expect_equal(unique_outcomes, c(1L, 3L)) +}) + +test_that("merge_cos_ucos reuses duplicate purity contexts", { + p <- 6 + L <- 3 + cb_obj <- list( + cb_data = list( + data = lapply(seq_len(L), function(i) { + list( + X = NULL, + XtX = diag(p), + variable_miss = integer(0), + ref_label = "LD" + ) + }), + dict = rep(1L, L) + ), + cb_model_para = list(P = p, L = L) + ) + class(cb_obj) <- "colocboost" + + out_cos <- list(cos = list( + cos = list(cos1 = c(1, 2)), + coloc_outcomes = list(1L), + avWeight = list(matrix(runif(p), nrow = p, dimnames = list(NULL, "outcome1"))), + cs_change = matrix(0.1, nrow = 1, ncol = L) + )) + out_ucos <- list( + ucos_each = list(ucos1 = c(2, 3)), + avW_ucos_each = matrix(runif(p), nrow = p, dimnames = list(NULL, "ucos1")), + change_obj_each = matrix(0.2, nrow = 1, ncol = L), + purity_each = matrix(1, nrow = 1, ncol = 3), + ucos_outcome = 2L + ) + + ns <- colocboost_test_env + between_calls <- 0L + local_mocked_bindings( + get_between_purity = function(...) { + between_calls <<- between_calls + 1L + c(min_abs_cor = 0.9, max_abs_cor = 1, median_abs_cor = 0.9) + }, + .package = "colocboost" + ) + + result <- get("merge_cos_ucos", envir = ns)( + cb_obj, + out_cos, + out_ucos, + median_cos_abs_corr = 0.8 + ) + + expect_equal(between_calls, 1L) + expect_null(result$ucos$ucos_each) + expect_equal(result$cos$cos$coloc_outcomes[[1]], c(1L, 2L)) +}) + +test_that("merge_cos_ucos skips purity checks for disjoint different-outcome sets", { + p <- 6 + L <- 2 + cb_obj <- list( + cb_data = list( + data = lapply(seq_len(L), function(i) { + list( + X = NULL, + XtX = diag(p), + variable_miss = integer(0), + ref_label = "LD" + ) + }), + dict = rep(1L, L) + ), + cb_model_para = list(P = p, L = L) + ) + class(cb_obj) <- "colocboost" + + out_cos <- list(cos = list( + cos = list(cos1 = c(1, 2)), + coloc_outcomes = list(1L), + avWeight = list(matrix(runif(p), nrow = p, dimnames = list(NULL, "outcome1"))), + cs_change = matrix(0.1, nrow = 1, ncol = L) + )) + out_ucos <- list( + ucos_each = list(ucos1 = c(4, 5)), + avW_ucos_each = matrix(runif(p), nrow = p, dimnames = list(NULL, "ucos1")), + change_obj_each = matrix(0.2, nrow = 1, ncol = L), + purity_each = matrix(1, nrow = 1, ncol = 3), + ucos_outcome = 2L + ) + + ns <- colocboost_test_env + local_mocked_bindings( + get_between_purity = function(...) { + stop("disjoint different-outcome sets should not require between-purity") + }, + .package = "colocboost" + ) + + result <- get("merge_cos_ucos", envir = ns)( + cb_obj, + out_cos, + out_ucos, + median_cos_abs_corr = 0.8 + ) + + expect_equal(length(result$ucos$ucos_each), 1L) + expect_equal(result$cos$cos$coloc_outcomes[[1]], 1L) +}) + +test_that("chunked update history preserves legacy matrix output", { + updates <- list( + list(update_status = c(1, 0, -1), real_update_jk = c(5, NA, 7), jk = c(5, 5, NA, 7)), + list(update_status = c(0, 1, 1), real_update_jk = c(NA, 4, 4), jk = c(4, NA, 4, 4)), + list(update_status = c(-1, 0, 0), real_update_jk = c(2, NA, NA), jk = c(2, 2, NA, NA)) + ) + + legacy <- list(update_status = c(), real_update_jk = c(), jk = c()) + chunked <- list(L = 3, update_status = c(), real_update_jk = c(), jk = c()) + + for (update in updates) { + legacy$update_status <- cbind(legacy$update_status, as.matrix(update$update_status)) + legacy$real_update_jk <- rbind(legacy$real_update_jk, update$real_update_jk) + legacy$jk <- rbind(legacy$jk, update$jk) + chunked <- get(".cb_append_update_history", envir = colocboost_test_env)( + chunked, + update_jk = update$jk, + update_status = update$update_status, + real_update_jk = update$real_update_jk + ) + } + + expect_gt(attr(chunked, "update_history_capacity"), length(updates)) + chunked <- get(".cb_trim_update_history", envir = colocboost_test_env)(chunked) + + expect_equal(chunked$update_status, legacy$update_status) + expect_equal(chunked$real_update_jk, legacy$real_update_jk) + expect_equal(chunked$jk, legacy$jk) + expect_null(attr(chunked, "update_history_capacity")) + expect_null(attr(chunked, "update_history_n")) +}) diff --git a/tests/testthat/test_utils.R b/tests/testthat/test_utils.R index e62fcf9..320e5d4 100644 --- a/tests/testthat/test_utils.R +++ b/tests/testthat/test_utils.R @@ -23,9 +23,9 @@ generate_test_result <- function(n = 100, p = 20, L = 2, seed = 42) { true_beta[10, 1] <- 1 # SNP10 also affects the trait } else { # Multi-trait case - true_beta[5, 1] <- 1 # SNP5 affects trait 1 - true_beta[5, 2] <- 1 # SNP5 also affects trait 2 (colocalized) - true_beta[10, 2] <- 0.5 # SNP10 only affects trait 2 + true_beta[5, 1] <- 2 # SNP5 affects trait 1 + true_beta[5, 2] <- 2 # SNP5 also affects trait 2 (colocalized) + true_beta[10, 2] <- 1.2 # SNP10 only affects trait 2 } # Generate Y with some noise @@ -401,8 +401,7 @@ test_that("get_cos extracts CoS correctly with generated test results", { expect_named(result_median_purity, c("cos", "cos_purity")) # Test empty colocalization results - empty_cb_output <- cb_output - empty_cb_output$cos_details$cos <- NULL + empty_cb_output <- list(cos_details = list(cos = NULL)) expect_warning( result_empty <- get_cos(empty_cb_output, coverage = 0.95), @@ -895,4 +894,3 @@ test_that("get_hierarchical_clusters handles extreme correlation structures", { expect_named(result_mixed, c("cluster", "Q_modularity")) expect_equal(nrow(result_mixed$cluster), P_mixed) }) - diff --git a/vignettes/Advanced_Colocalization_Scenarios.Rmd b/vignettes/Advanced_Colocalization_Scenarios.Rmd new file mode 100644 index 0000000..359b5c4 --- /dev/null +++ b/vignettes/Advanced_Colocalization_Scenarios.Rmd @@ -0,0 +1,84 @@ +--- +title: "Advanced Colocalization Scenarios with ColocBoost" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Advanced Colocalization Scenarios with ColocBoost} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + dpi = 70 +) +``` + +This vignette uses representative simulation studies to illustrate two advanced colocalization scenarios addressed by ColocBoost: + +- **Multiple causal variants per trait**: Regions containing multiple causal variants, beyond the one-causal-variant-per-trait assumption. +- **Weaker effects in disease GWAS**: Shared causal signals for which the disease trait contributes weaker association evidence than the accompanying molecular traits. + +# 1. Multiple causal variants within a genomic region + +To reduce the combinatorial hypothesis space, Bayesian multi-trait colocalization methods commonly assume that *each trait has at most one causal variant within a genomic region* (**one-causal-variant-per-trait assumption**). +This assumption becomes increasingly restrictive as the number of phenotypes increases and more distinct signals and trait-sharing patterns must be resolved. +Collapsing these signals into a single-signal representation can obscure event-specific sharing patterns, leading to missed or incorrectly localized colocalization events. +This concern has also been emphasized and evaluated for pairwise colocalization using COLOC (V5) (Wallace, 2021, *PLOS Genetics*). + + +## Scenario 1: Heterogeneous effects across traits + +A common multi-signal scenario arises when multiple causal variants are shared across traits but have **heterogeneous** effects. +Consider two traits influenced by two causal variants. Under the *one-causal-variant-per-trait* assumption, each trait is represented only by its strongest signal (Figure 2b(i)): + +- **Trait 1** is represented by causal variant 1, which has the strongest association with Trait 1. +- **Trait 2** is represented by causal variant 2, which has the strongest association with Trait 2. + +The resulting single-signal representations appear as two distinct trait-specific signals, +leading to a false conclusion of no colocalization even though both causal variants are shared across the two traits. +ColocBoost instead resolves the two shared signals as distinct colocalization events. + +Heterogeneous effects of two causal variants across traits. + + +## Scenario 2: Non-causal strongest marginal effect + +Another multi-signal scenario occurs when a non-causal variant tags multiple causal variants through LD and consequently has the strongest marginal association. +Distinguishing marginal association from causal attribution motivates multi-effect fine-mapping methods such as SuSiE (Wang et al., 2020, *JRSS B*). + +Consider two traits sharing the same two causal variants. Under the *one-causal-variant-per-trait* assumption (Figure 2b(ii)): + +- **Trait 1 and Trait 2** are represented by the non-causal marginal lead (green dot), which has a stronger marginal association than either true causal variant (red dots). + +The resulting single-signal representation incorrectly localizes the colocalized signal to a non-causal variant, whereas ColocBoost resolves the two shared causal signals as distinct colocalization events. + +A non-causal variant has the strongest marginal association. + +# 2. Colocalization with weaker effects in GWAS + +In practice, it is often of interest to colocalize a disease GWAS with multiple molecular QTL traits to elucidate the functional basis of disease associations. +An important technical aspect of GWAS-xQTL colocalization is that GWAS traits often have lower per-variant contributions to heritability than molecular xQTL traits. + +## Scenario 3: Weaker effects in disease GWAS + +Consider a disease GWAS and an xQTL sharing the same two causal variants (Figure 2b(iii)): + +- **xQTL** shows strong association evidence for both causal variants. +- **Disease GWAS** shows strong evidence for causal variant 1 but weaker evidence for causal variant 2. + +COLOC (V5) identifies the event supported by the stronger GWAS signal but misses the second event with weaker GWAS evidence. +As a two-stage approach that performs fine-mapping before colocalization, +it may have reduced sensitivity to weaker signals with limited support in the initial single-trait analysis. +ColocBoost identifies both shared signals using its disease-prioritized colocalization approach. + +Colocalization with a weaker causal effect in the disease GWAS. + +See [Mixed Data-type and Disease Prioritized Colocalization](https://statfungen.github.io/colocboost/articles/Disease_Prioritized_Colocalization.html) for practical guidance on GWAS-xQTL analysis with the ColocBoost disease-prioritized mode. diff --git a/vignettes/ColocBoost_Update.Rmd b/vignettes/ColocBoost_Update.Rmd index a7de526..408ad3d 100644 --- a/vignettes/ColocBoost_Update.Rmd +++ b/vignettes/ColocBoost_Update.Rmd @@ -26,8 +26,8 @@ The animation below demonstrates how ColocBoost iteratively updates across multi Observe how the blue points (proxies) distribute around the selected best update, showing the algorithm's ability to capture LD structure during optimization. ```{r, out.width="80%"} -knitr::include_graphics("../man/figures/ColocBoost_update.gif") +knitr::include_graphics("figures/ColocBoost_update.gif") ``` -See more details in the [ColocBoost Manuscript](https://www.medrxiv.org/content/10.1101/2025.04.17.25326042v1) and the [GitHub repository](https://github.com/StatFunGen/colocboost). \ No newline at end of file +See more details in the [ColocBoost Manuscript](https://www.medrxiv.org/content/10.1101/2025.04.17.25326042v1) and the [GitHub repository](https://github.com/StatFunGen/colocboost). diff --git a/vignettes/Conceptual_Multi_Trait_Colocalization.Rmd b/vignettes/Conceptual_Multi_Trait_Colocalization.Rmd new file mode 100644 index 0000000..6560f04 --- /dev/null +++ b/vignettes/Conceptual_Multi_Trait_Colocalization.Rmd @@ -0,0 +1,259 @@ +--- +title: "Conceptual Framework for Multi-trait Colocalization and ColocBoost" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Conceptual Framework for Multi-trait Colocalization and ColocBoost} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + dpi = 70 +) +``` + +This vignette introduces key conceptual definitions for multi-trait colocalization, describes their practical implementation in ColocBoost, +and relates them to established concepts in fine-mapping and existing colocalization methods. + +# 1. Multi-trait colocalization problem as a colocalization event decomposition + +The multi-trait colocalization problem concerns identifying *shared causal signals* across \(L \geq 2\) traits within a genomic region of interest. +To formalize this problem, ColocBoost introduces a new conceptual definition of a multi-trait ***Colocalization Event*** as a triplet \( \{ g(s), T(s), CoS_{\alpha}(s) \} \) for each detected event \(s\), comprising + +- \(g(s)\): ***Genomic Region*** underlying the colocalization analysis. +- \(T(s)\): ***Trait Configuration*** specifying which subset of \(L\) traits share the same causal variant for event \(s\). +- \(CoS_{\alpha}(s)\): ***Colocalization Confidence Set*** representing the smallest set of variants containing the shared causal variant with probability at least \(\alpha\) (default \(\alpha = 0.95\)). + +A *central challenge* in multi-trait colocalization is that multiple events can arise within the same genomic region, involving overlapping or distinct sets of traits and multiple causal variants. + +- In Bayesian formulations that specify all possible trait configurations a *priori*, the hypothesis space grows combinatorially and rapidly becomes computationally prohibitive. +- ColocBoost avoids this enumeration by reformulating colocalization as a multi-task learning problem optimized through gradient boosting. This formulation jointly resolves multiple colocalization events and scales to hundreds of traits. + + +Illustration of colocalization events in ColocBoost. + + +# 2. Conceptual ColocBoost summaries and their analogies to existing methods + +ColocBoost characterizes each event using two complementary summaries: variant-level localization and event-level evidence for sharing across the corresponding trait configuration. +These summaries have direct structural analogies to established quantities in fine-mapping and colocalization, while being defined specifically for multi-trait colocalization events. + +- **Variant-level evidence:** The ***colocalization confidence set*** (CoS) and ***variant colocalization probability*** (VCP) localize the variants underlying each colocalization event. +They are structurally analogous to the credible set (CS) and posterior inclusion probability (PIP), respectively, in statistical fine-mapping methods such as SuSiE. + +- **Colocalization evidence:** The ***normalized probability of colocalization*** (NPC) quantifies support for colocalization. +It is structurally analogous to PP.H4 in COLOC for pairwise colocalization and PPFC in HyPrColoc for multi-trait colocalization. + + +## 2.1. Variant-level evidence + +For each event \(s\), representing a single causal signal shared by the subset of traits \(T(s)\), + +- \(\xi^s = (\xi_1^s, \ldots, \xi_P^s)\) denotes the vector of *single-effect colocalization probabilities* across \(P\) variants in the region, quantifying the probability that each variant is the shared causal variant underlying traits \(T(s)\). +- Structural analogy: single-effect posterior \( \alpha_l \) in SuSiE for single-effect \(l\) in a single-trait model. + +ColocBoost then defines an \(\alpha\)-level *Colocalization Confidence Set*, \(CoS_{\alpha}(s)\) (default \(\alpha = 0.95\)), including the candidate causal variant and its high-LD proxies, +with construction based on \(\xi^s\): + +\[ + CoS_{\alpha}(s) = \left\{ v_1, v_2, \ldots, v_{p_0} : p_0 = \min (p: \sum_{j=1}^p \xi_{(j)}^s \geq \alpha ) \right\}, +\] + +where \(\xi_{(1)}^s \geq \ldots \geq \xi_{(P)}^s\) are the sorted single-effect colocalization probabilities. +ColocBoost discards the CoS with low *purity* (minimum absolute correlation between all pairs of variants within CoS, default threshold \(purity < 0.5\)). + +- Structural analogy: \(\alpha\)-level credible set (CS) in SuSiE in a single-trait model. + +ColocBoost also defines the *variant colocalization probability* (VCP) for each variant \(j\) in the region as +\[ + VCP_j = 1 - \prod_{s=1}^{S} (1 - \xi_j^s). +\] +The construction of VCP is based on the assumption that each event \(s\) is *conditionally* independent. + +- Structural analogy: posterior inclusion probability (PIP) in SuSiE in a single-trait model. +- Structural analogy: variant-level posterior probability for H4 (SNP.PP.H4) in COLOC in a pairwise colocalization analysis. + +Illustration of variant-level analogues in ColocBoost. + + +## 2.2. Colocalization evidence + +### Narrative definition of normalized probability of colocalization (NPC) + +For each event \(s\), ColocBoost defines the *normalized probability of colocalization* (NPC) as an empirical, event-level measure comparing shared colocalization with trait-specific alternatives. + +- Higher NPC values indicate stronger evidence that the event represents a causal signal shared by at least two traits. +- Lower NPC values indicate weaker evidence for sharing, while greater consistency with the single-trait specific causal signals. + +NPC is structurally analogous, but not probabilistically equivalent, to + +- PP.H4: posterior probability of H4, two traits shared the same causal variant, in COLOC for pairwise colocalization. +- PPFC: posterior probability of full colocalization in HyPrColoc for multi-trait colocalization. + +NPC is assigned to each detected colocalization event \(s\), allowing multiple distinct events within the same genomic region and a potentially different trait configuration \(T(s)\). + +::: {.note-box} +**Note:** PP.H4 is defined **only** for two traits and denotes posterior probability for \(H_4\), the hypothesis that both traits are associated and share the same causal variant, +relative to hypotheses \(H_0\), \(H_1\), \(H_2\), and \(H_3\). For \(L > 2\), the colocalization hypothesis space must encompass all possible trait configurations across shared and distinct causal variants (Foley et al., 2021, *Nature Communications*). +For example, \(H_{(L-2,1,1)}\) represents a configuration in which \(L-2\) traits share one causal variant, while the remaining two traits have distinct causal variants. +The full hypothesis space grows combinatorially, comprising \(\mathrm{Bell}(L+1)\) hypotheses. +::: + +ColocBoost avoids this combinatorial explosion by discovering supported configurations \(T(s)\) and their corresponding \(CoS_{\alpha}(s)\) in a data-driven manner, +without enumerating all possible trait configurations (details in ColocBoost paper). + + +### Mathematical definition of normalized probability of colocalization (NPC) + +For a detected colocalization event \(s\) with a triplet \( \{ g(s), T(s), CoS_{\alpha}(s) \} \), +ColocBoost first defines the trait-level normalized evidence score \(NP_l^s\) for each trait \(l \in T(s)\): + +\[ + NP_l^s = 1 - \exp(-\lambda_l LRT_l^s), +\] + +where \(\lambda_l \) is a trait-specific rate that adjusts the scale of normalization based on a baseline log-likelihood ratio for each trait \(l\) (details in Supplementary Note of ColocBoost paper). + +Here, \(LRT_l^s\) is a log-likelihood ratio test statistic between two models, \(M^s_{0,l}\) and \(M^s_{1,l}\). + +- \(M^s_{0,l}\): the null model where all variants have zero effects on trait \(l\) (\( \beta_l=0 \)). +- \(M^s_{1,l}\): the alternative model that variants in \(CoS_{\alpha}(s)\) have non-zero effects (\( \beta^s_l \neq 0 \)). + +::: {.note-box} +***Rationale:*** NPC is evaluated only after event \(s\) has been detected; therefore, at least one trait is expected to provide strong evidence for the event. +NPC provides event-level evidence that the detected event is jointly supported by at least two traits, rather than being driven by evidence from only one trait. +::: + +**Two-trait example** + +For the two-trait case, without loss of generality, let \(NP_1^s \ge NP_2^s\), so that trait 1 is the leading trait for event \(s\). +ColocBoost approximates the evidence for the single-trait, non-colocalized configuration, in which trait 1 contributes but trait 2 does not, as + +\[ +NPUC_s += +\underbrace{NP_1^s}_{ +\substack{\text{evidence supporting}\\ +\text{trait 1 contribution}} +} +\times +\underbrace{(1-NP_2^s)}_{ +\substack{\text{lack of evidence supporting}\\ +\text{trait 2 contribution}} +}. +\] + +Accordingly, a large NPUC indicates that event \(s\) is supported **primarily** by the leading trait. +ColocBoost subsequently defines the event-level colocalization evidence NPC as + +\[ + NPC_s = 1 - NPUC_s. +\] + +NPUC and NPC are normalized colocalization evidence scores rather than posterior probabilities. + +**Multi-trait generalization** + +For \(L>2\), ColocBoost generalizes NPUC as the *leading-trait-only explanation* across all traits: + +\[ +NPUC_s = NP^s_{l_{\max}} \prod_{l\neq l_{\max}} \left(1-NP_{l}^{s}\right), \, \mathrm{and} \, NPC_s = 1 - NPUC_s. +\] + +Here, \(NP^s_{l_{\max}}\) represents the normalized evidence from the leading trait, whereas the product term captures the lack of support from the remaining traits. +Together, these terms quantify the extent to which event \(s\) is supported primarily by the leading trait. +Accordingly, a high NPC indicates support from more than one trait. + + +::: {.note-box} +**Highlight:** ColocBoost provides complementary evidence at two levels: \(NPC_s\) evaluates the overall colocalization event, +whereas \(NP_l^s\) quantifies each trait's support for that event. +In our numerical studies, the lenient thresholds \(NPC_s \geq 0.5\) and \(NP_l^s \geq 0.2\) maintained well-controlled FDR while retaining reasonable detection power. +More stringent thresholds may be applied to prioritize the strongest colocalization signals. +::: + +See more details about filtering colocalization events by relative strength of evidence +using ['get_robust_colocalization'](https://statfungen.github.io/colocboost/articles/Interpret_ColocBoost_Output.html#filter-colocalization-events-by-relative-strength-of-evidence) function. + + +Illustration of event-level analogues in ColocBoost. + + +Concordance between NPC and PP.H4 or PPFC was assessed only for detected CoS with at least 95% overlap variants between ColocBoost and COLOC or HyPrColoc, respectively. + + +# 3. Practical interpretation for `colocboost` output + +This section maps the concepts above to the corresponding `colocboost` output fields. +See [Interpret ColocBoost Output](https://statfungen.github.io/colocboost/articles/Interpret_ColocBoost_Output.html) for detailed guidance. + +After running `res = colocboost()`, + +- `res$cos_summary`: a summary of all colocalization events. Each row corresponds to one colocalization event \(s\) and includes columns with + - `colocalized_outcomes`: *Trait Configuration* -- \(T(s)\); + - `colocalized_variables`: *Colocalization Confidence Set* -- \(CoS_{\alpha}(s)\); + - `colocalized_variables_vcp`: *Variant Colocalization Probability* for shared variants in \(CoS_{\alpha}(s)\); + - `cos_npc`: *Normalized Probability of Colocalization* -- \(NPC_s\); + - `purity`: minimum absolute correlation between all pairs of variants within \(CoS_{\alpha}(s)\). +- `res$cos_details`: detailed information about all colocalization events, including sublists with + - `vcp`: a length-\(P\) vector of *Variant Colocalization Probability* for all \(P\) variants; + - `cos_vcp`: a list of single-effect colocalization probabilities for each event \(s\); + - `cos_outcomes_npc`: a list of trait-level evidence (\(NP^s_l\)) for each event \(s\). + +## Example: Causal variant structure +The dataset features two causal variants with indices 194 and 589. + +- Causal variant 194 is associated with traits 1, 2, 3, and 4. +- Causal variant 589 is associated with traits 2, 3, and 5. + +```{r run-colocboost} +library(colocboost) +# Loading the Dataset +data(Ind_5traits) +# Run colocboost +res <- colocboost(X = Ind_5traits$X, Y = Ind_5traits$Y) +``` + +Colocalization events summary: +```{r events_summary} +cos_summary <- res$cos_summary +cos_summary[, + c( + "colocalized_outcomes", + "colocalized_variables", + "colocalized_variables_vcp", + "cos_npc", + "purity" + ) +] +``` + +Trait-level evidence: +```{r trait_level} +res$cos_details$cos_outcomes_npc +``` diff --git a/vignettes/Interpret_ColocBoost_Output.Rmd b/vignettes/Interpret_ColocBoost_Output.Rmd index b927391..b650031 100644 --- a/vignettes/Interpret_ColocBoost_Output.Rmd +++ b/vignettes/Interpret_ColocBoost_Output.Rmd @@ -7,6 +7,22 @@ vignette: > %\VignetteEncoding{UTF-8} --- + + ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, @@ -73,20 +89,30 @@ cos_interest_outcome <- get_cos_summary(res, interest_outcome = c("Y1", "Y2")) ## 2. Filter colocalization events by relative strength of evidence -In `cos_summary`, for each 95% CoS, the `cos_npc` column provides a normalized probability of colocalization and -`min_npc_outcome` column provides the minimum normalized probability among colocalized traits. -Those two metrics are measured as an empirical evidence of colocalization both in CoS-level and in trait-level. -To obtain the best minimal colocalization configuration can be defined by using both `cos_npc` and `npc_outcome`. -See the detailed usage of this function in [link](https://statfungen.github.io/colocboost/reference/get_robust_colocalization.html). +For each 95% colocalization confidence set (CoS) reported in `cos_summary`, two complementary quantities summarize the strength of colocalization evidence: +- `cos_npc` reports the event-level normalized probability of colocalization (NPC), which quantifies support for sharing beyond a single-trait explanation. +- `min_npc_outcome` reports the minimum trait-level normalized evidence among the traits assigned to the event. + +These quantities assess whether an event is supported by multiple traits and whether each trait contributes sufficient evidence. +The `get_robust_colocalization()` function applies thresholds at both levels to obtain a parsimonious, well-supported trait configuration. + +See [Conceptual Framework for Multi-trait Colocalization and ColocBoost](https://statfungen.github.io/colocboost/articles/Conceptual_Multi_Trait_Colocalization.html) for detailed definitions and interpretation of these quantities. ```{r run-strong-colocalization} -filter_res <- get_robust_colocalization(res, cos_npc_cutoff = 0.5, npc_outcome_cutoff = 0.2) +filter_res <- get_robust_colocalization( + res, + cos_npc_cutoff = 0.5, + npc_outcome_cutoff = 0.2 +) ``` -- The output from `get_robust_colocalization` is the same as output from `colocboost`, which can be directly used in any post inference and visualization. -- `npc=0.5` or `npc_outcome = 0.2` maintains robust colocalization signals for cases when many traits are evaluated. -Higher thresholds can be specified if users want to focus only on strong colocalization events. +::: {.note-box} +**Highlight:** + +- The output of `get_robust_colocalization()` retains the same structure as the output of `colocboost()` and can be used directly in downstream inference and visualization workflows. +- In our numerical studies, `cos_npc_cutoff = 0.5` and `npc_outcome_cutoff = 0.2` retained robust colocalization signals when many traits were evaluated. Higher thresholds may be used to prioritize events with stronger colocalization evidence. +::: ## 3. More details on ColocBoost output diff --git a/vignettes/Partial_Overlap_Variants.Rmd b/vignettes/Partial_Overlap_Variants.Rmd index b74d7f3..a1055cd 100644 --- a/vignettes/Partial_Overlap_Variants.Rmd +++ b/vignettes/Partial_Overlap_Variants.Rmd @@ -1,8 +1,8 @@ --- -title: "Handling partial overlapping variants across traits in ColocBoost" +title: "Handling Partial Overlapping Variants across Traits in ColocBoost" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{Handling partial overlapping variants across traits in ColocBoost} + %\VignetteIndexEntry{Handling Partial Overlapping Variants across Traits in ColocBoost} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -22,7 +22,7 @@ This vignette demonstrates how ColocBoost handles partial overlapping variants a library(colocboost) ``` -![Illustration of partial overlapping variants across traits](../man/figures/missing_representation.png) +![Illustration of partial overlapping variants across traits](figures/missing_representation.png) @@ -107,4 +107,3 @@ res$data_info$n_variables # Plotting the results colocboost_plot(res) ``` - diff --git a/vignettes/announcements.Rmd b/vignettes/announcements.Rmd index 2b1fe29..498aa2b 100644 --- a/vignettes/announcements.Rmd +++ b/vignettes/announcements.Rmd @@ -10,10 +10,13 @@ vignette: > ## ColocBoost on the media -- *April 20, 2025*: Manuscript describing ColocBoost methods with applications is posted on [medRxiv](https://www.medrxiv.org/content/10.1101/2025.04.17.25326042v1). +- *April 20, 2025*: Manuscript describing ColocBoost methods with applications is posted on [medRxiv](https://doi.org/10.1101/2025.04.17.25326042). - *May 2, 2025*: `colocboost` R package is available on [CRAN](https://CRAN.R-project.org/package=colocboost). ## Software updates +- `v1.0.10` [**Important**] Major enhancements to conceptual interpretation and computational scalability. + - Expanded the event-based framework for multi-trait colocalization, supported by new tutorials and practical guidance for interpreting CoS, VCP, and NPC. + - Improved computational efficiency relative to v1.0.9 by reducing redundant LD calculations, purity evaluations, matrix operations, and memory reallocation, particularly for analyses involving large numbers of phenotypes. - `v1.0.9` Improvements to summary-statistics workflows, trait-specific result filtering, and computational efficiency. - Added `X_ref` support as a memory-efficient alternative to precomputed LD matrices for large summary-statistics analyses. - Added `get_robust_ucos` to recalibrate and summarize robust trait-specific, uncolocalized events. diff --git a/man/figures/ColocBoost_update.gif b/vignettes/figures/ColocBoost_update.gif similarity index 100% rename from man/figures/ColocBoost_update.gif rename to vignettes/figures/ColocBoost_update.gif diff --git a/vignettes/figures/Colocalization_Events.png b/vignettes/figures/Colocalization_Events.png new file mode 100644 index 0000000..dacf69a Binary files /dev/null and b/vignettes/figures/Colocalization_Events.png differ diff --git a/vignettes/figures/Event_Level_Analogues.png b/vignettes/figures/Event_Level_Analogues.png new file mode 100644 index 0000000..8246e82 Binary files /dev/null and b/vignettes/figures/Event_Level_Analogues.png differ diff --git a/vignettes/figures/Figure2b_i.png b/vignettes/figures/Figure2b_i.png new file mode 100644 index 0000000..bde00b0 Binary files /dev/null and b/vignettes/figures/Figure2b_i.png differ diff --git a/vignettes/figures/Figure2b_ii.png b/vignettes/figures/Figure2b_ii.png new file mode 100644 index 0000000..294ad5e Binary files /dev/null and b/vignettes/figures/Figure2b_ii.png differ diff --git a/vignettes/figures/Figure2b_iii.png b/vignettes/figures/Figure2b_iii.png new file mode 100644 index 0000000..2ddff72 Binary files /dev/null and b/vignettes/figures/Figure2b_iii.png differ diff --git a/vignettes/figures/Variant_Level_Analogues.png b/vignettes/figures/Variant_Level_Analogues.png new file mode 100644 index 0000000..7c12a8d Binary files /dev/null and b/vignettes/figures/Variant_Level_Analogues.png differ diff --git a/man/figures/missing_representation.png b/vignettes/figures/missing_representation.png similarity index 100% rename from man/figures/missing_representation.png rename to vignettes/figures/missing_representation.png