Prepare ColocBoost 1.0.10: scalability, robustness, and tutorials

.Rbuildignore

@@ -20,3 +20,4 @@
 ^docs$
 ^_pkgdown\.yml$
 ^CRAN-SUBMISSION$
+^\.superpowers$

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")),

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,

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)
 }
 

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,

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

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)
 }
-