---
title: "NeighbourNet: Workflow Example"
author: "Yidi Deng"
date: "2025-11-10"
output: html_document
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(Seurat)
library(NeighbourNet)
library(ggplot2)
library(igraph)
library(patchwork)

# Data of murine hematopoietic progenitors generated by Nestorowa et al. (2016)
# We acquired the data from the Seurat vignette that demonstrate cell-cycle scoring 
# https://satijalab.org/seurat/articles/cell_cycle_vignette.html
exp.mat <- read.table(
  file = "./nestorawa_forcellcycle_expressionMatrix.txt",
  header = TRUE,
  row.names = 1,
  check.names = FALSE
)

# Cell-cycle markers from Tirosh et al. (2015) are available via Seurat.
# These are separated into S-phase and G2/M-phase markers.
s.genes  <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

# Create a Seurat object and run basic preprocessing
obj <- CreateSeuratObject(
  counts = Matrix::Matrix(as.matrix(exp.mat), sparse = TRUE)
)
obj <- NormalizeData(obj)
obj <- CellCycleScoring(
  obj,
  s.features   = s.genes,
  g2m.features = g2m.genes,
  set.ident    = TRUE
)
```

## Preprocessing

We first perform quality control and gene selection, then prepare the object for NeighbourNet.

```{r gene-selection}
# Filter out lowly expressed genes and annotate TFs/targets/receptors/ligands
genes <- NeighbourNet::select.gene(obj, min.cells = 10)
```

Run PCA on the selected genes. Co-expression will be measured within this gene set.

```{r prepare-seurat}
obj <- NeighbourNet::prepare.seurat(obj, genes = genes$genes)
```

Construct the KNN graph used by NeighbourNet.

```{r prepare-graph}
obj <- NeighbourNet::prepare.graph(obj)
```

Optionally, for data sets with many cells (> 5000), you can run the regression on a subset of representative cells.

```{r select-cell, eval = FALSE}
# Not run
obj <- NeighbourNet::select.cell(obj)
```

Finally, set up the regression problem and compute local variance.

```{r prepare-reg}
# Local variance is computed for predictor and response genes.
# Response genes are low-rank approximated.
# Here we measure co-expression between TFs (predictors) and cell-cycle targets (responses).
targets <- c(s.genes, g2m.genes)

obj <- NeighbourNet::prepare.reg(
  obj,
  predictors = genes$tfs,
  responses  = targets
)
```

The NeighbourNet regression settings are stored in the `misc` slot:

```{r inspect-setting}
Seurat::Misc(obj, "NNet.setting") %>% names()
```

## NeighbourNet regression

Run NeighbourNet regression to obtain per-cell co-expression networks.

```{r run-nn-reg}
obj <- NeighbourNet::run.nn.reg(
  obj,
  responses    = targets,
  return.p.val = TRUE
)
```

The resulting network ensemble is stored in `NNet.mod` within the `misc` slot.

```{r inspect-mod}
# effect: (response × predictor × cell) tensor of co-expression effect sizes
# p.val : corresponding significance tensor
Seurat::Misc(obj, "NNet.mod") %>% names()
```

Build meta-networks summarising the per-cell networks.

```{r build-meta-network}
obj <- NeighbourNet::build.meta.network(obj)
```

Meta-networks are stored as a (response × predictor × meta-cell) tensor:

```{r meta-dim, eval = FALSE}
Seurat::Misc(obj, "NNet.mod")$meta.network$meta.network %>% dim()
```

## Visualisation

We next visualise key TF–target subnetworks and upstream signalling structure.

### Select central genes

First, identify central genes based on eigenvector centrality across meta-networks.

```{r central-genes}
# All responses will be visualised when keep.responses = TRUE.
# n.net controls the number of meta-networks to examine.
# For each meta-network, k singular vectors are used, and
# n.per.component top-loading genes per component are selected.
central.genes <- NeighbourNet::select.central.genes(
  obj,
  keep.responses   = FALSE,
  n.per.component  = 4,
  n.net            = 5,
  k                = 2
)
```

### Prepare visualisation settings

```{r prepare-visualise}
# TFs (or the chosen g2 layer) are clustered into n.clu groups (colour-coded)
obj <- NeighbourNet::prepare.visualise(
  obj,
  central.genes = central.genes,
  n.clu         = 4
)
```

### Visualise a per-cell network

```{r plot-cell, fig.height = 10, fig.width = 10}
i <- 1
NeighbourNet::visualise.network(
  obj,
  i,
  cutoff      = 0.95,
  radius      = c(0.4, 0.7, 0.85, 1),
  pie.radius  = 0.04,
  text.size   = 5
)
```

### Visualise a meta-network component

```{r plot-meta, fig.height = 10, fig.width = 10}
i <- 1
meta.p <- Seurat::Misc(obj, "NNet.mod")$meta.network$p.val[,, i]
cutoff <- apply(meta.p, 1, max) %>%
  sort(decreasing = TRUE) %>%
  mean()

NeighbourNet::visualise.network(
  obj,
  i,
  meta.network = TRUE,
  cutoff       = cutoff,
  radius       = c(0.4, 0.7, 0.85, 1),
  pie.radius   = 0.04,
  text.size    = 5
)
```

### Visualise meta-network components in PC space

We can also inspect which cells contribute to a given meta-network component in the
meta-network PC space.

```{r plot-meta-pcs, fig.height = 5, fig.width = 6}
i <- 1
pcs     <- Seurat::Misc(obj, "NNet.mod")$meta.network$pcs
weights <- Seurat::Misc(obj, "NNet.mod")$meta.network$npca.loadings[, i]

ggplot(pcs) +
  geom_point(aes(component_1, component_2, colour = weights)) +
  scale_color_continuous("Weights")+
  theme_classic() +
  theme(text = element_text(size = 15))+
  ggtitle("Cell Weights on Meta-network 1")

ggplot(pcs) +
  geom_point(aes(component_1, component_2, colour = obj$Phase)) +
  scale_color_discrete("Phase")+
  theme_classic() +
  theme(text = element_text(size = 15)) +
  ggtitle("Cell by Predicted Phase")
```

## Receptor activity

We now compute receptor activity scores using NeighbourNet's receptor-TF-target
integration.

```{r receptor-activity}
act <- NeighbourNet::receptor.activity(obj)
```

### Visualise NOTCH2 activity and identify mediating TFs

```{r notch2-activity, fig.height = 5, fig.width = 6}
notch2.act <- act$receptor.act["NOTCH2", ] %>% as.numeric()

ggplot(pcs) +
  geom_point(aes(component_1, component_2, colour = notch2.act)) +
  scale_color_continuous("Act")+
  theme_classic() +
  theme(text = element_text(size = 15))+
  ggtitle("Notch2")

tf.cor <- cor(notch2.act, t(act$tf.act)) %>%
  drop() %>%
  sort(decreasing = TRUE)
head(tf.cor)
```

We can inspect the inferred signalling path from NOTCH2 to a candidate TF
(e.g. E2F1) using the prior signalling graph:

```{r notch2-e2f1-path}
igraph::shortest_paths(
  NeighbourNet::sig.graph,
  from    = "NOTCH2",
  to      = "E2F1",
  weights = 1 / igraph::E(NeighbourNet::sig.graph),
  output  = "both"
)$vpath
```
