Multihop (effective) connectivity clustering in cf_cosine_plot

DESCRIPTION

@@ -25,6 +25,7 @@ Imports:
     checkmate,
     dplyr,
     coconat (> 0.1.2),
+    Matrix,
     stringr,
     magrittr,
     bit64,

R/cosine.R

@@ -18,15 +18,17 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') {
       x[['outputs']],
       inputcol = 'pre_key',
       outputcol = groupcol,
-      inputids = ids)
+      inputids = ids,
+      standardise_input = FALSE)
     cm[['cout']] = coconat::cosine_sim(oam, transpose = T)
     cm[['wout']]=sum(x[['outputs']]$weight)
   }
   if('inputs' %in% partners) {
     groupcol <- if(isFALSE(group)) "pre_key" else group
     iam <- coconat::partner_summary2adjacency_matrix(
       x[['inputs']],
-      inputcol = groupcol, outputcol = 'post_key', outputids = ids)
+      inputcol = groupcol, outputcol = 'post_key', outputids = ids,
+      standardise_input = FALSE)
     cm[['cin']] = coconat::cosine_sim(iam, transpose = F)
     cm[['win']]=sum(x[['inputs']]$weight)
   }
@@ -53,6 +55,21 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') {
 #'   \code{hemilineage}, \code{serial} (serially homologous cell group) if
 #'   available. This can reveal other interesting features of organisation.
 #'
+#'   \bold{Multihop (effective) connectivity}. By default (\code{nhops=0})
+#'   neurons are clustered by their direct partners. Setting \code{nhops>0}
+#'   instead clusters them by their \emph{effective} connectivity onto partners
+#'   reached through one or more intermediate interneuron layers, following the
+#'   method of Schlegel et al. (2021). At each step the connectivity is
+#'   input-normalised (so the inputs to every postsynaptic cell sum to 1) and
+#'   the successive matrices are multiplied, computed separately within each
+#'   dataset. \code{nhops=1} corresponds to their \dQuote{one-hop} (2nd-order)
+#'   pathways, \code{nhops=2} to 3rd-order and so on. Because the intermediate
+#'   partner sets grow quickly, \code{min_frac} prunes each layer to types that
+#'   receive at least that fraction of effective input (per-type, applied at
+#'   every hop including the final one); it defaults to a small non-zero value
+#'   when \code{nhops>0}. A grouping column (e.g. the default \code{group="type"})
+#'   is required for \code{nhops>0}.
+#'
 #'   The \code{labRow} argument is most conveniently specified as a length 1
 #'   string to be interpolated by \code{\link[glue]{glue}}; this will happen in
 #'   the context of a data frame generated by \code{\link{cf_meta}}. One reason
@@ -79,6 +96,17 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') {
 #' @param group The name or the grouping column for partner connectivity
 #'   (defaults to \code{"type"}) or a logical where \code{group=FALSE} means no
 #'   grouping (see details).
+#' @param nhops Number of intermediate interneuron layers to traverse when
+#'   computing effective connectivity. \code{0} (the default) uses direct
+#'   partners; \code{1} is the \dQuote{one-hop} (2nd-order) case, etc. See
+#'   \bold{details}.
+#' @param min_frac Per-type fractional threshold (default \code{0.005}) used
+#'   when \code{nhops>0} to prune each layer (including the final one) to partner
+#'   types receiving at least this fraction of effective input. A scalar or a
+#'   vector with one value per hop. Ignored when \code{nhops=0}.
+#' @param remove_query Whether to exclude the query neurons from the partners at
+#'   every hop, both as intermediate interneurons and as final targets (default
+#'   \code{FALSE}). Only relevant when \code{nhops>0}.
 #' @param labRow Optionally, either string that can be interpolated by
 #'   \code{\link[glue]{glue}} \emph{or} a character vector matching the number
 #'   of neurons specified by \code{ids}. See \bold{details} for an important
@@ -236,6 +264,8 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') {
 #' }
 cf_cosine_plot <- function(ids=NULL, ..., threshold=5,
                            partners = c("outputs", "inputs"),
+                           nhops=0L, min_frac=0.005,
+                           remove_query=FALSE,
                            labRow='{type}_{coconatfly::abbreviate_datasets(dataset)}{side}',
                            group='type',
                            heatmap=TRUE,
@@ -253,7 +283,10 @@ cf_cosine_plot <- function(ids=NULL, ..., threshold=5,
     x=ids
     partners=unique(x$partners)
   } else
-    x=multi_connection_table(ids, partners = partners, threshold = threshold, group=group, min_datasets = min_datasets)
+    x=multi_connection_table(ids, partners = partners, threshold = threshold,
+                             group=group, min_datasets = min_datasets,
+                             nhops=nhops, min_frac=min_frac,
+                             remove_query=remove_query)
 
   cm <- multi_cosine_matrix(x, partners = partners, group=group, nas=nas)
 
@@ -342,6 +375,8 @@ cf_cosine_plot <- function(ids=NULL, ..., threshold=5,
 #'  output neurons.
 multi_connection_table <- function(ids, partners=c("inputs", "outputs"),
                                    threshold=1L, group='type',
+                                   nhops=0L, min_frac=0.005,
+                                   remove_query=FALSE,
                                    check_missing=TRUE,
                                    min_datasets=Inf,
                                    prefer.foreign=NA,
@@ -352,10 +387,15 @@ multi_connection_table <- function(ids, partners=c("inputs", "outputs"),
   if(isTRUE(group))
     group='type'
   partners=match.arg(partners, several.ok = T)
+  if(nhops>0 && !is.character(group))
+    stop("Multihop connectivity (nhops>0) requires a grouping column ",
+         "(e.g. group='type') to define partner layers across hops.")
   kk=keys(ids)
   if(length(partners)>1) {
     l=sapply(partners, simplify = F, function(p)
       multi_connection_table(kk, partners=p, threshold = threshold, group=group,
+                             nhops=nhops, min_frac=min_frac,
+                             remove_query=remove_query,
                              check_missing=F, min_datasets = min_datasets,
                              prefer.foreign=prefer.foreign, MoreArgs=MoreArgs,
                              keep.all=keep.all, ...))
@@ -383,8 +423,14 @@ multi_connection_table <- function(ids, partners=c("inputs", "outputs"),
        ((length(datasets)>1 && "malecns" %in% datasets) && is.na(prefer.foreign)))
       MoreArgs=list(malecns=list(prefer.foreign=TRUE))
   }
-  x <- cf_partners(kk, threshold = threshold, partners = partners,
-                   MoreArgs = MoreArgs, keep.all=keep.all, ...)
+  x <- if(nhops>0)
+    multihop_partner_summary(kk, partners=partners, nhops=nhops,
+                             threshold=threshold, min_frac=min_frac,
+                             group=group, remove_query=remove_query,
+                             MoreArgs=MoreArgs)
+  else
+    cf_partners(kk, threshold = threshold, partners = partners,
+                MoreArgs = MoreArgs, keep.all=keep.all, ...)
   if(is.character(group))
     x <- match_types(x, group, partners=partners, min_datasets = min_datasets)
   # mark which column was used for the query

R/multihop.R

@@ -0,0 +1,134 @@
+# Multihop (effective) connectivity
+#
+# Compute effective connectivity from a set of query neurons to their n-th order
+# partners, passing through intermediate interneuron layers, following the
+# input-normalise-then-multiply approach of Schlegel et al. 2021
+# (doi:10.7554/eLife.62576). See multihop-plan.md / cf_cosine_plot docs.
+
+# Build a raw pre x post sparse adjacency matrix from a cf_partners table.
+# Data is already standardised by cf_partners so we skip standardise_input.
+.mh_adjacency <- function(tbl) {
+  coconat::partner_summary2adjacency_matrix(
+    tbl, inputcol = "pre_key", outputcol = "post_key",
+    standardise_input = FALSE, sparse = TRUE)
+}
+
+# Effective query x target-type matrix for a SINGLE dataset and direction.
+# Returns a sparse Matrix (query neurons as rows, target types as columns) or
+# NULL if the walk dies out.
+multihop_effective_matrix <- function(dskeys, partners, nhops, threshold,
+                                      min_frac, group = "type",
+                                      remove_query = FALSE, MoreArgs = list()) {
+  stopifnot(nhops >= 1L)
+  # one min_frac per stage: nhops intermediate prunes + 1 terminal cut
+  mf <- rep_len(min_frac, nhops + 1L)
+  query_keys <- unique(as.character(dskeys))
+  # which side of a hop table is the advancing frontier / query anchor
+  newfrontier_key <- if (partners == "outputs") "post_key" else "pre_key"
+
+  frontier <- query_keys
+  running <- NULL          # query x current-frontier (neuron resolution)
+  far_types <- NULL        # named vector: terminal neuron key -> type
+
+  for (h in seq_len(nhops + 1L)) {
+    if (length(frontier) == 0) return(NULL)
+    tbl <- cf_partners(frontier, threshold = threshold, partners = partners,
+                       MoreArgs = MoreArgs)
+    if (is.null(tbl) || nrow(tbl) == 0) return(NULL)
+    if (!group %in% colnames(tbl))
+      stop("Grouping column `", group, "` not present in cf_partners result!")
+
+    Mn <- coconat::colScaleM(.mh_adjacency(tbl))
+    # advance the query x frontier running matrix (see multihop-plan.md):
+    # outputs -> multiply by Mn (old x new); inputs -> by t(Mn) (old x new)
+    # advance the query x frontier running matrix by one hop. effective_
+    # connectivity aligns step's rows to colnames(running) by name (dead-end
+    # paths -> zero). normalise=FALSE: running is already a product of
+    # normalised matrices and step was normalised above, so neither is rescaled.
+    step <- if (partners == "outputs") Mn else Matrix::t(Mn)
+    running <- if (is.null(running)) step
+               else coconat::effective_connectivity(list(running, step),
+                                                     normalise = FALSE)
+
+    # type of each new-frontier neuron (partner side of this hop)
+    ntypes <- tbl[[group]][!duplicated(tbl[[newfrontier_key]])]
+    names(ntypes) <- tbl[[newfrontier_key]][!duplicated(tbl[[newfrontier_key]])]
+    ntypes <- ntypes[colnames(running)]
+
+    if (h <= nhops) {
+      # intermediate layer: prune frontier by type (selection only)
+      G <- coconat::grouping_matrix(colnames(running), ntypes)
+      grp <- running %*% G
+      keep_types <- colnames(grp)[apply(as.matrix(grp), 2, max) >= mf[h]]
+      surviving <- colnames(running)[ntypes %in% keep_types]
+      if (remove_query)
+        surviving <- setdiff(surviving, query_keys)
+      running <- running[, surviving, drop = FALSE]
+      frontier <- surviving
+    } else {
+      far_types <- ntypes
+    }
+  }
+
+  # terminal: optionally drop query neurons from the target set (a query type
+  # can still appear as a feature if non-query neurons of that type are targets)
+  if (remove_query) {
+    keepcols <- setdiff(colnames(running), query_keys)
+    running <- running[, keepcols, drop = FALSE]
+    far_types <- far_types[keepcols]
+  }
+  # group far neurons to type (per-neuron normalisation already done)
+  eff <- running %*% coconat::grouping_matrix(colnames(running), far_types)
+  # final per-type cut
+  eff <- Matrix::drop0(eff * (eff >= mf[nhops + 1L]))
+  eff
+}
+
+# Melt an effective query x target-type matrix into a cf_partners-like long
+# table for one direction, with the query on the appropriate key column and the
+# target type in `group`. Feeds straight into multi_cosine_matrix.
+.mh_matrix2df <- function(eff, dataset, partners, group = "type") {
+  if (is.null(eff) || length(eff) == 0 || sum(eff != 0) == 0)
+    return(NULL)
+  s <- Matrix::summary(Matrix::drop0(eff))
+  qk <- rownames(eff)[s$i]
+  tp <- colnames(eff)[s$j]
+  df <- data.frame(stringsAsFactors = FALSE, weight = s$x)
+  df[[group]] <- tp
+  # query lives in pre_key for outputs, post_key for inputs
+  if (partners == "outputs") {
+    df$pre_key <- qk
+    df$post_key <- tp
+  } else {
+    df$post_key <- qk
+    df$pre_key <- tp
+  }
+  # pre_id/post_id kept so the table is recognised by is.mct(); they are not
+  # coerced downstream (multi_cosine_matrix uses standardise_input=FALSE)
+  df$pre_id <- sub("^[a-z]+:", "", df$pre_key)
+  df$post_id <- df$post_key
+  df$dataset <- dataset
+  df$partners <- partners
+  df
+}
+
+# Effective connectivity long table for ONE direction across all datasets.
+# Mirrors the per-direction output of cf_partners (before match_types) so it can
+# be dropped into multi_connection_table.
+multihop_partner_summary <- function(kk, partners, nhops, threshold, min_frac,
+                                     group = "type", remove_query = FALSE,
+                                     MoreArgs = list()) {
+  kk <- keys(kk)
+  kdf <- keys2df(kk)
+  datasets <- unique(kdf$dataset)
+  res <- lapply(datasets, function(ds) {
+    dskeys <- kk[kdf$dataset == ds]
+    ma <- if (!is.null(MoreArgs[[ds]])) MoreArgs[ds] else list()
+    eff <- multihop_effective_matrix(
+      dskeys, partners = partners, nhops = nhops, threshold = threshold,
+      min_frac = min_frac, group = group, remove_query = remove_query,
+      MoreArgs = ma)
+    .mh_matrix2df(eff, dataset = ds, partners = partners, group = group)
+  })
+  dplyr::bind_rows(res)
+}