R/blisa.R

In blisa: Infer Cell-Cell Communication from Spatial Transcriptomics

#' Run BLISA spatial cell-cell communication analysis
#'
#' Generic function for running BLISA (Bivariate Local Indicator of Spatial
#' Association). Dispatches on the class of \code{x}:
#' \itemize{
#'   \item \code{blisa.default} accepts a pre-binned gene-by-bin count
#'     matrix and a matching \code{bins} polygon object.
#'   \item \code{blisa.SpatialExperiment} accepts a cell-level
#'     \code{SpatialExperiment} object and bins cells into hexagonal tiles
#'     internally via \code{\link{hexBinCells}} before running the analysis.
#' }
#'
#' @param x A gene-by-bin count matrix (for \code{blisa.default}) or a
#'   cell-level \code{SpatialExperiment} object (for
#'   \code{blisa.SpatialExperiment}).
#' @param ... Additional arguments passed to the relevant method.
#'
#' @return A list; see individual method documentation for details.
#' @examples
#' \dontrun{
#' # Download and cache the example Xenium breast cancer dataset (~21 MB)
#' data_url  <- paste0(
#'   "https://github.com/ChenLaboratory/example_data/releases/",
#'   "download/v1.0.0/spe_xenium_bc_s1rep1.rds"
#' )
#' cache_dir <- tools::R_user_dir("blisa", "cache")
#' data_file <- file.path(cache_dir, "spe_xenium_bc_s1rep1.rds")
#' if (!file.exists(data_file)) {
#'   dir.create(cache_dir, recursive = TRUE, showWarnings = FALSE)
#'   download.file(data_url, data_file, mode = "wb")
#' }
#' spe <- readRDS(data_file)
#'
#' # blisa.SpatialExperiment method: bins cells, runs LISA, scores CCI
#' result <- blisa(spe, bin_size = 50, group = "cell_type")
#' result
#' }
#' @export
blisa <- function(x, ...) UseMethod("blisa")


#' @describeIn blisa Method for a gene-by-bin count matrix.
#'
#' @param bins An \code{sf} object of bin polygons. Row order must match the
#'   column order of \code{x}.
#' @param LR_df Data frame of ligand-receptor pairs with columns
#'   \code{ligand.symbol} and \code{receptor.symbol}. When \code{NULL},
#'   CellChatDB for the chosen \code{species} is downloaded automatically.
#' @param bin_size Numeric. Bin spacing used to define queen adjacency for
#'   nearby-mode interactions. Default \code{50}.
#' @param dmax Numeric. Maximum distance for diffuse-mode neighbours.
#'   Default \code{250}.
#' @param nsim Integer. Number of permutations for Moran's I significance.
#'   Default \code{999}.
#' @param p_cutoff Numeric. P-value threshold for High-High hotspots.
#'   Default \code{0.05}.
#' @param min_ligand Numeric. Minimum ligand count threshold. Default \code{10}.
#' @param min_receptor Numeric. Minimum receptor count threshold.
#'   Default \code{10}.
#' @param min_cells Integer. Minimum number of cells per bin. For a
#'   \code{SpatialExperiment}, bins below this threshold are dropped during
#'   binning by \code{\link{hexBinCells}}. For pre-binned input supplied
#'   directly to \code{blisa.default}, they are dropped from \code{bins},
#'   \code{x} and \code{counts_by_group} before the spatial weights are built;
#'   this requires \code{n_cells_col} and is ignored when
#'   \code{n_cells_col = NA}. Default \code{1}.
#' @param n_cells_col Character or \code{NA}. Column in \code{bins} holding
#'   per-bin cell counts used for \code{min_cells} filtering.
#'   Set to \code{NA} to skip (default).
#' @param annotation_col Character. Column in \code{LR_df} specifying
#'   interaction category used for communication-mode assignment. Default
#'   \code{"annotation"}.
#' @param default_mode Character. CCC mode assigned to LR pairs whose
#'   annotation does not match \code{diffuse_category}. Default
#'   \code{"diffuse"}.
#' @param diffuse_category Character vector of annotation categories treated
#'   as diffuse signalling.
#' @param species Character. Which CellChatDB to download when
#'   \code{LR_df = NULL}. One of \code{"human"} (default) or \code{"mouse"}.
#' @param genes Character vector of gene names to consider when matching
#'   ligand-receptor pairs. Only LR pairs whose ligand and receptor genes are
#'   all present in \code{genes} and pass the expression thresholds will be
#'   tested. Defaults to \code{rownames(x)} (all genes in the count matrix).
#'   Genes specified here but absent from \code{rownames(x)} are silently
#'   ignored.
#' @param counts_by_group Named list of gene-by-bin count matrices, one per
#'   group level (e.g. cell type), as returned by \code{\link{hexBinCells}}
#'   when \code{group} is supplied. When provided, \code{\link{runCCI}} is
#'   called automatically after the BLISA loop and its output is included in
#'   the result as \code{CCI_scores}. Default \code{NULL}.
#' @param fast Logical. When \code{TRUE} (default), uses
#'   \code{fastLISA::local_moran_bv} (a fast C/OpenMP backend) for the
#'   bivariate local Moran's I computation. When \code{FALSE}, uses the
#'   original \code{spdep::localmoran_bv} + \code{spdep::hotspot} pipeline.
#' @param cpu_threads Integer. Number of OpenMP threads used by
#'   \code{fastLISA::local_moran_bv}. Only used when \code{fast = TRUE}.
#'   Ignored on platforms without OpenMP. Default \code{4L}.
#' @param verbose Logical. If \code{TRUE}, print progress messages (per-LR-pair
#'   progress, bin filtering, CCI status). Default \code{FALSE} (silent).
#'   A progress bar over the ligand-receptor pairs is shown automatically in
#'   interactive sessions (and includes the current LR pair label when
#'   \code{verbose = TRUE}).
#'
#' @return An object of class \code{blisa} with four components:
#' \describe{
#'   \item{LR_results}{Data frame of BLISA results for each LR pair, including
#'     \code{ccc_mode}, \code{sig_numbers}, \code{sig_index}, \code{sig_pval},
#'     \code{all_pval}, \code{all_lisa}, \code{all_quadrant}, and original
#'     columns from \code{LR_df}. \code{all_quadrant} is a character vector of
#'     hotspot quadrant labels (\code{"High-High"}, \code{"Low-Low"}, etc.)
#'     for every bin; non-tested bins are \code{NA}.}
#'   \item{bins}{Bin-level \code{sf} object of hexagonal polygons.}
#'   \item{spatial_weights}{Spatial weights list from \code{\link{computeSpatialWeights}}.}
#'   \item{CCI_scores}{\code{NULL} unless \code{counts_by_group} is supplied,
#'     in which case a wide data frame of interaction scores from
#'     \code{\link{runCCI}}: rows are \code{"Sender->Receiver"} group pairs,
#'     columns are LR pairs.}
#' }
#' @export
blisa.default <- function(
    x,
    bins,
    LR_df             = NULL,
    bin_size         = 50,
    dmax             = 250,
    nsim             = 999,
    p_cutoff         = 0.05,
    min_ligand       = 10,
    min_receptor     = 10,
    min_cells        = 1,
    n_cells_col      = NA,
    annotation_col   = "annotation",
    default_mode      = "diffuse",
    diffuse_category  = c("Secreted Signaling", "Non-protein Signaling"),
    species           = c("human", "mouse"),
    genes             = NULL,
    counts_by_group   = NULL,
    fast              = TRUE,
    cpu_threads       = 4L,
    verbose           = FALSE,
    ...
) {
  msg   <- function(...) if (verbose) message(...)
  quiet <- function(expr)
    if (verbose) expr else suppressWarnings(suppressMessages(expr))

  ## ---------------------------
  ## Drop low-cell bins
  ## For users who supply their own binned data and want cell-number filtering
  ## without going through the blisa binning pipeline, where hexBinCells() has
  ## already applied this filter. 
  ## ---------------------------
  if (!is.na(n_cells_col)) {
    if (!n_cells_col %in% colnames(bins))
      stop("Column '", n_cells_col, "' not found in bins.")

    keep_bins <- bins[[n_cells_col]] >= min_cells

    if (all(keep_bins)) {
      msg("No bins below ", min_cells, " cells, subsetting skipped.")
    } else {
      if (ncol(x) != nrow(bins))
        stop("ncol(x) (", ncol(x), ") must match nrow(bins) (", nrow(bins),
             ") to subset bins.")
      msg(sum(!keep_bins), " bins dropped (< ", min_cells, " cells, column '",
          n_cells_col, "').")

      bins <- bins[keep_bins, , drop = FALSE]
      x    <- x[, keep_bins, drop = FALSE]
      if (!is.null(counts_by_group))
        counts_by_group <- lapply(counts_by_group,
                                  function(m) m[, keep_bins, drop = FALSE])
    }
  }

  sw             <- quiet(computeSpatialWeights(bins, bin_size, dmax))
  queen_wt       <- sw$queen_wt
  dist_wt        <- sw$dist_wt
  keep_idx_queen <- sw$keep_idx_queen
  keep_idx_dist  <- sw$keep_idx_dist

  # Restrict gene panel for LR pair matching; default to all genes in x
  genes <- if (is.null(genes)) rownames(x) else intersect(genes, rownames(x))

  LR_df_filtered <- filterLRpairs(
    counts       = x[genes, , drop = FALSE],
    min_ligand   = min_ligand,
    min_receptor = min_receptor,
    LR_df        = LR_df,
    species      = species
  )

  LR_results <- quiet(LR_df_add_mode(LR_df_filtered, annotation_col,
                                     default_mode, diffuse_category))

  # Replace the CellChatDB interaction_name row names with informative IDs that
  # spell out ligand/receptor subunits (e.g. "TGFB1_TGFBR2|TGFBR1"). The original
  # name is preserved in a column for cross-referencing.
  LR_results$interaction_name <- rownames(LR_results)
  rownames(LR_results)        <- make_lr_id(LR_results$ligand.symbol,
                                            LR_results$receptor.symbol)

  LR_results$sig_numbers   <- integer(nrow(LR_results))
  LR_results$sig_index     <- vector("list", nrow(LR_results))
  LR_results$sig_pval      <- vector("list", nrow(LR_results))
  LR_results$all_pval      <- vector("list", nrow(LR_results))
  LR_results$all_lisa      <- vector("list", nrow(LR_results))
  LR_results$all_quadrant  <- vector("list", nrow(LR_results))

  n_pairs <- nrow(LR_results)
  show_pb <- interactive() && n_pairs > 0
  if (verbose || show_pb) message("Testing ", n_pairs, " LR pairs...")

  for (i in seq_len(n_pairs)) {
    ligand   <- LR_results$ligand.symbol[i]
    receptor <- LR_results$receptor.symbol[i]
    mode     <- LR_results$ccc_mode[i]

    if (mode == "nearby") {
      wt       <- queen_wt
      keep_idx <- keep_idx_queen
    } else {
      wt       <- dist_wt
      keep_idx <- keep_idx_dist
    }

    if (show_pb) {
      # Under knitr, carriage-return updates are each captured as separate
      # lines, so only the final (100%) line is drawn there; in an interactive
      # console the bar updates in place every iteration.
      in_knit <- isTRUE(getOption("knitr.in.progress"))
      if (!in_knit || i == n_pairs) {
        frac   <- i / n_pairs
        width  <- 40L
        filled <- round(frac * width)
        lab    <- if (verbose)
          sprintf("%-42s", paste0(rownames(LR_results)[i], " (", mode, ")")) else ""
        bar <- sprintf("  %s|%s%s| %3.0f%%",
                       lab, strrep("=", filled), strrep(" ", width - filled),
                       frac * 100)
        if (in_knit) {
          cat(bar, "\n", sep = "")
        } else {
          cat("\r", bar, sep = "")
          utils::flush.console()
          if (i == n_pairs) cat("\n")
        }
      }
    } else {
      msg(rownames(LR_results)[i])
      msg("ccc mode is ", mode)
    }

    x_full <- get_min_expr(receptor, x)
    y_full <- get_min_expr(ligand, x)

    if (fast) {
      # fastLISA backend: C/OpenMP
      res_bv <- fastLISA::local_moran_bv(
        x_full[keep_idx], y_full[keep_idx], wt,
        nsim    = nsim,
        p.value = p_cutoff,
        n.cores = cpu_threads
      )
      hs <- attr(res_bv, "cluster")
    } else {
      # spdep backend (original implementation)
      res_bv <- spdep::localmoran_bv(x_full[keep_idx], y_full[keep_idx], wt,
                                     nsim = nsim)
      hs <- spdep::hotspot(
        res_bv,
        Prname        = "Pr(folded) Sim",
        cutoff        = p_cutoff,
        quadrant.type = "pysal",
        p.adjust      = "none"
      )
    }

    full_pval <- rep(1, ncol(x))
    full_lisa <- rep(0, ncol(x))
    full_pval[keep_idx] <- res_bv[, "Pr(folded) Sim"]
    full_lisa[keep_idx] <- res_bv[, "Ibvi"]

    hs_idx_hh <- !is.na(hs) & hs == "High-High"
    HH_idx    <- keep_idx[which(hs_idx_hh)]
    HH_pval   <- full_pval[HH_idx]

    full_quadrant <- rep(NA_character_, ncol(x))
    full_quadrant[keep_idx] <- as.character(hs)

    LR_results$sig_numbers[i]   <- length(HH_idx)
    LR_results$sig_index[[i]]   <- HH_idx
    LR_results$sig_pval[[i]]    <- HH_pval
    LR_results$all_pval[[i]]    <- full_pval
    LR_results$all_lisa[[i]]    <- full_lisa
    LR_results$all_quadrant[[i]] <- full_quadrant
    LR_results$ccc_mode[i]      <- mode
  }

  LR_results <- LR_results[order(-LR_results$sig_numbers), , drop = FALSE]
  front_cols <- c("ccc_mode", "sig_numbers", "sig_index", "sig_pval")
  LR_results <- LR_results[, c(front_cols, setdiff(colnames(LR_results), front_cols))]

  CCI_scores <- NULL
  if (!is.null(counts_by_group)) {
    msg("Running CCI analysis...")
    tmp <- new_blisa(LR_results, bins, sw)
    CCI_scores <- runCCI.default(tmp, counts_by_group)
  }

  new_blisa(LR_results, bins, sw, CCI_scores)
}


#' @describeIn blisa Method for a cell-level SpatialExperiment object.
#'   For cell-level platforms (e.g. Xenium) it bins cells into hexagonal tiles
#'   via \code{\link{hexBinCells}}; for Visium it treats each spot as a bin via
#'   \code{\link{visiumSpotBins}}. Either way it then delegates to
#'   \code{blisa.default}. The path is chosen by \code{platform}.
#'
#' @param platform Character. \code{"auto"} (default) inspects
#'   \code{colData(x)}: if \code{array_col}/\code{array_row} are present the
#'   object is treated as \code{"visium"}, otherwise \code{"xenium"}.
#'   Set explicitly to override the detection.
#' @param bin_size Numeric. Width of each hexagonal bin in coordinate units
#'   (e.g. microns). Passed to \code{\link{hexBinCells}} and
#'   \code{\link{computeSpatialWeights}}. Default \code{50}; for the Visium
#'   path it defaults to \code{100} (spot center-to-center spacing) when left
#'   unset.
#' @param group Character. Column name in \code{colData(x)} to use as the
#'   grouping variable (e.g. cell type) for per-group bin aggregation and
#'   downstream CCI analysis via \code{\link{runCCI}}. If the column is not
#'   found in \code{colData(x)}, a message is issued and CCI is skipped.
#'   Default \code{"cell_type"}.
#' @param genes Character vector of gene names to consider when matching
#'   ligand-receptor pairs. Defaults to \code{rownames(x)} (all genes in the
#'   \code{SpatialExperiment} object).
#' @param min_total_counts Numeric. Bins whose total counts (summed over all
#'   genes) fall below this threshold are dropped during binning by
#'   \code{\link{hexBinCells}}. Set to \code{0} to disable. Default \code{10}.
#' @param verbose Logical. If \code{TRUE}, print progress messages. Default
#'   \code{FALSE} (silent).
#'
#' @export
blisa.SpatialExperiment <- function(x, bin_size = 50, LR_df = NULL,
                                       group = "cell_type", genes = NULL,
                                       min_cells = 1, min_total_counts = 10,
                                       platform = c("auto", "xenium", "visium"),
                                       verbose = FALSE, ...) {
  platform  <- match.arg(platform)
  cd_cols   <- colnames(SummarizedExperiment::colData(x))
  has_array <- all(c("array_col", "array_row") %in% cd_cols)

  if (platform == "auto") {
    platform <- if (has_array) "visium" else "xenium"
    if (verbose) message("platform auto-detected as '", platform, "'.")
  }
  if (platform == "visium" && !has_array && verbose)
    message("platform = 'visium' but no array_col/array_row in colData; ",
            "spatialCoords will be used directly (no de-tilt) \u2014 ensure ",
            "bin_size matches the coordinate units.")

  ## ---- Visium: each spot is already a bin ------------------------------
  if (platform == "visium") {
    if (!is.null(group) && group %in% cd_cols && verbose)
      message("Visium spots have no per-cell labels \u2014 group '", group,
              "' ignored; CCI analysis skipped.")

    # bin_size unset -> let visiumSpotBins pick it (100 um for the array-index
    # lattice, or measured from the coordinates in the fallback).
    sp <- if (missing(bin_size)) NULL else bin_size
    binned <- visiumSpotBins(x,
                             spot_pitch       = sp,
                             min_total_counts = min_total_counts,
                             verbose          = verbose)

    # Scale dmax to the coordinate units unless the caller supplied one.
    dots <- list(...)
    if (is.null(dots$dmax)) dots$dmax <- 2.5 * binned$pitch

    return(do.call(blisa.default, c(list(
      x        = binned$counts_matrix,
      bins     = binned$bins,
      LR_df    = LR_df,
      bin_size = binned$pitch,
      genes    = if (is.null(genes)) rownames(x) else genes,
      verbose  = verbose
    ), dots)))
  }

  ## ---- Xenium (cell-level): bin cells into hexes -----------------------
  coords <- as.data.frame(SpatialExperiment::spatialCoords(x))

  # Resolve group vector from colData
  if (!is.null(group) && group %in% cd_cols) {
    group_vec <- SummarizedExperiment::colData(x)[[group]]
  } else {
    if (!is.null(group) && verbose)
      message("Column '", group, "' not found in colData(x) \u2014 CCI analysis will be skipped.")
    group_vec <- NULL
  }

  if (is.null(genes)) genes <- rownames(x)

  binned <- hexBinCells(coords, SummarizedExperiment::assay(x, "counts"),
                        bin_size = bin_size, min_cells = min_cells,
                        min_total_counts = min_total_counts, group = group_vec,
                        verbose = verbose)

  blisa.default(
    x               = binned$counts_matrix,
    bins          = binned$bins,
    LR_df           = LR_df,
    bin_size        = bin_size,
    n_cells_col     = "n_cells",
    genes           = genes,
    counts_by_group = binned$counts_by_group,
    verbose         = verbose,
    ...
  )
}