diff --git a/DESCRIPTION b/DESCRIPTION index 4531e2c..53a08a0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,6 +25,7 @@ Imports: checkmate, dplyr, coconat (> 0.1.2), + Matrix, stringr, magrittr, bit64, diff --git a/R/cosine.R b/R/cosine.R index 3acab23..acf43df 100644 --- a/R/cosine.R +++ b/R/cosine.R @@ -18,7 +18,8 @@ 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) } @@ -26,7 +27,8 @@ multi_cosine_matrix <- function(x, partners, nas, group='type') { 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 diff --git a/R/multihop.R b/R/multihop.R new file mode 100644 index 0000000..d3fcf41 --- /dev/null +++ b/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) +} diff --git a/man/cf_cosine_plot.Rd b/man/cf_cosine_plot.Rd index d386b90..f716a4a 100644 --- a/man/cf_cosine_plot.Rd +++ b/man/cf_cosine_plot.Rd @@ -10,6 +10,9 @@ cf_cosine_plot( ..., 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, @@ -28,6 +31,9 @@ multi_connection_table( 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, @@ -48,6 +54,20 @@ option since neuprint sometimes returns 0 weight edges.} \item{partners}{Whether to return inputs or outputs} +\item{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}.} + +\item{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}.} + +\item{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}.} + \item{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 @@ -137,6 +157,21 @@ examples. \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 diff --git a/tests/testthat/test-multihop.R b/tests/testthat/test-multihop.R new file mode 100644 index 0000000..b4b2be4 --- /dev/null +++ b/tests/testthat/test-multihop.R @@ -0,0 +1,120 @@ +# Offline tests of the multihop effective-connectivity core. cf_partners is +# mocked with a small synthetic 3-layer network so no live service is needed. +# +# Network (outputs / downstream): +# query q1 (fw:1), q2 (fw:2) +# layer1 i11 (fw:11, type A), i12 (fw:12, type B) +# layer2 t21 (fw:21, type X), t22 (fw:22, type Y) + +synthetic_edges <- function() { + data.frame( + pre_key = c("fw:1","fw:2","fw:2", "fw:11","fw:11","fw:12"), + post_key = c("fw:11","fw:11","fw:12", "fw:21","fw:22","fw:22"), + weight = c(10, 5, 5, 8, 2, 4), + stringsAsFactors = FALSE) +} + +# type of each neuron (partner side) +synthetic_type <- function(key) { + lu <- c("fw:11"="A","fw:12"="B","fw:21"="X","fw:22"="Y", + "fw:1"="Q","fw:2"="Q") + unname(lu[key]) +} + +mock_cf_partners <- function(ids, threshold=1L, partners="outputs", + MoreArgs=list(), ...) { + e <- synthetic_edges() + keycol <- if(partners=="outputs") "pre_key" else "post_key" + partcol <- if(partners=="outputs") "post_key" else "pre_key" + e <- e[e[[keycol]] %in% as.character(ids) & e$weight >= threshold, , drop=FALSE] + if(nrow(e)==0) return(NULL) + e$type <- synthetic_type(e[[partcol]]) + e$dataset <- "flywire" + e +} + +test_that("multihop_effective_matrix matches hand-computed one-hop", { + testthat::local_mocked_bindings(cf_partners = mock_cf_partners) + eff <- coconatfly:::multihop_effective_matrix( + c("fw:1","fw:2"), partners="outputs", nhops=1L, + threshold=1L, min_frac=0, group="type") + m <- as.matrix(eff) + # expected (see header): q1 X=2/3 Y=2/9 ; q2 X=1/3 Y=7/9 + expect_setequal(rownames(m), c("fw:1","fw:2")) + expect_setequal(colnames(m), c("X","Y")) + m <- m[c("fw:1","fw:2"), c("X","Y")] + expect_equal(unname(m), + matrix(c(2/3, 2/9, 1/3, 7/9), nrow=2, byrow=TRUE)) +}) + +test_that("frontier pruning by min_frac drops weak interneuron types", { + testthat::local_mocked_bindings(cf_partners = mock_cf_partners) + # interneuron type A max effective input-fraction = 2/3, type B = 1. + # min_frac 0.7 prunes A, keeping only i12 (type B) -> only target Y via q2. + eff <- coconatfly:::multihop_effective_matrix( + c("fw:1","fw:2"), partners="outputs", nhops=1L, + threshold=1L, min_frac=0.7, group="type") + m <- as.matrix(eff) + expect_equal(colnames(m), "Y") + # q2 -> i12 -> t22 fully: 1 ; q1 has no path through i12 + expect_equal(m["fw:2","Y"], 1) + expect_true(!"fw:1" %in% rownames(m) || m["fw:1","Y"]==0) +}) + +test_that("remove_query drops query neurons from the target set", { + # add a recurrent path i11 -> q1 (fw:1, type Q) so a query neuron is a target + edges2 <- rbind(synthetic_edges(), + data.frame(pre_key="fw:11", post_key="fw:1", weight=6, + stringsAsFactors=FALSE)) + typ2 <- c("fw:11"="A","fw:12"="B","fw:21"="X","fw:22"="Y","fw:1"="Q","fw:2"="Q") + mock2 <- function(ids, threshold=1L, partners="outputs", MoreArgs=list(), ...) { + kc <- if(partners=="outputs")"pre_key" else "post_key" + pc <- if(partners=="outputs")"post_key" else "pre_key" + e <- edges2[edges2[[kc]] %in% as.character(ids) & edges2$weight>=threshold,,drop=FALSE] + if(nrow(e)==0) return(NULL) + e$type <- unname(typ2[e[[pc]]]); e$dataset <- "flywire"; e + } + testthat::local_mocked_bindings(cf_partners = mock2) + # without removal, the query neuron fw:1 (type Q) shows up as a feature + eff0 <- coconatfly:::multihop_effective_matrix( + c("fw:1","fw:2"), "outputs", 1L, 1L, 0, "type", remove_query=FALSE) + expect_true("Q" %in% colnames(eff0)) + # with removal it is dropped (no other Q target exists) + eff1 <- coconatfly:::multihop_effective_matrix( + c("fw:1","fw:2"), "outputs", 1L, 1L, 0, "type", remove_query=TRUE) + expect_false("Q" %in% colnames(eff1)) +}) + +test_that("nhops=2 handles an interneuron dead end (row alignment)", { + # q1 -> {i1(A), i2(B)} ; i1 -> j1(X) ; i2 has NO outputs (dead end) ; + # j1 -> t1(Z). With alignment, i2's missing onward row must become zero. + edges3 <- data.frame( + pre_key = c("fw:1","fw:1", "fw:11", "fw:21"), + post_key = c("fw:11","fw:12","fw:21", "fw:31"), + weight = c(10, 5, 8, 6), stringsAsFactors=FALSE) + typ3 <- c("fw:11"="A","fw:12"="B","fw:21"="X","fw:31"="Z") + mock3 <- function(ids, threshold=1L, partners="outputs", MoreArgs=list(), ...) { + e <- edges3[edges3$pre_key %in% as.character(ids) & edges3$weight>=threshold,,drop=FALSE] + if(nrow(e)==0) return(NULL) + e$type <- unname(typ3[e$post_key]); e$dataset <- "flywire"; e + } + testthat::local_mocked_bindings(cf_partners = mock3) + eff <- coconatfly:::multihop_effective_matrix( + "fw:1", "outputs", nhops=2L, threshold=1L, min_frac=0, group="type") + m <- as.matrix(eff) + # only surviving path q1 -> i1 -> j1 -> t1(Z); each step fully normalised -> 1 + expect_equal(colnames(m), "Z") + expect_equal(unname(m["fw:1","Z"]), 1) +}) + +test_that(".mh_matrix2df melts to a cf_partners-like table", { + m <- Matrix::Matrix(c(2/3, 2/9, 1/3, 7/9), nrow=2, byrow=TRUE, + dimnames=list(c("fw:1","fw:2"), c("X","Y"))) + df <- coconatfly:::.mh_matrix2df(m, dataset="flywire", partners="outputs") + expect_true(all(c("pre_key","post_key","type","weight","dataset","partners", + "pre_id","post_id") %in% colnames(df))) + expect_setequal(df$pre_key, c("fw:1","fw:2")) + expect_setequal(df$type, c("X","Y")) + expect_equal(nrow(df), 4L) + expect_true(is.mct(transform(df)) || all(c("pre_id","post_id") %in% colnames(df))) +})