scTenifoldNet is an R package for constructing and comparing single-cell gene regulatory networks (scGRNs) from single-cell RNA-seq (scRNA-seq) data collected under different conditions. It employs a machine-learning workflow combining principal component regression, tensor decomposition, and manifold alignment to accurately identify even subtly shifted gene expression programs.
Implementations in other languages are also available:
scTenifoldNet is available on CRAN:
install.packages("scTenifoldNet")
To install the development version from GitHub:
# install.packages("remotes")
# remotes::install_github("cailab-tamu/scTenifoldNet")
install.packages("scTenifoldNet", repos = c('https://cailab-tamu.r-universe.dev', 'https://cloud.r-project.org'))
The package is written in R only, so installing from source needs no compiler.
The scTenifoldNet() function orchestrates a six-step pipeline. Each step reports progress to the console via the cli package, with sample-level labels ([X] / [Y]) so that operations on each condition are clearly distinguishable.
| Step | Function | Description |
|---|---|---|
| 1 | scQC |
Quality control — filters cells by library size, outlier detection, minimum gene expression fraction, and mitochondrial read ratio |
| 2 | cpmNormalization |
Counts-per-million (CPM) normalization |
| 3 | makeNetworks |
Constructs gene regulatory networks from subsampled cells using principal component regression (pcNet). All the per-gene regressions come from one eigendecomposition, which gives the same networks as fitting each gene separately |
| 4 | tensorDecomposition |
CANDECOMP/PARAFAC (CP) tensor decomposition for network denoising |
| 5 | manifoldAlignment |
Non-linear manifold alignment of the two denoised networks |
| 6 | dRegulation |
Differential regulation testing via Box-Cox transformation and chi-square statistics |
Individual functions are exported and fully documented, allowing users to run or modify any step independently.
The required input is a raw counts matrix with genes as rows and cells (barcodes) as columns. Data should be unnormalized when using the main scTenifoldNet() function. The modular design of the package allows users to substitute custom preprocessing at any step.
scTenifoldNet() returns a list with three elements:
tensorNetworks — Weight-averaged denoised gene regulatory networks after CP tensor decomposition, containing:
X: The network for sample X (sparse matrix of class dgCMatrix).Y: The network for sample Y (sparse matrix of class dgCMatrix).The diagonal (self-loops) of both networks is set to zero. The CP decomposition can produce non-zero values there, but they are not used by the manifold alignment.
manifoldAlignment — A data frame of low-dimensional features from the non-linear manifold alignment, with 2 × n genes rows and d columns (default d = 30).
diffRegulation — A data frame with six columns:
gene: Gene identifier.distance: Euclidean distance between the gene's coordinates in the two conditions.Z: Z-score after Box-Cox power transformation.FC: Fold change with respect to the expectation.p.value: P-value from the chi-square distribution with one degree of freedom.p.adj: Adjusted p-value (Benjamini & Hochberg FDR correction).Network construction subsamples cells at random, so results depend on the random seed. The seed argument (default 1) is set before each random step. The same input and parameters therefore give the same result on every run, whatever the caller's RNG state, and that state is restored when the function returns.
seed = 1, seed = 2, ...) to see how much results vary between runs.seed = NULL to let a set.seed() call made before the function control the result.tensorDecomposition() accepts the same seed argument.
Running time grows mainly with the number of genes. The number of cells matters little, because each network is built from a fixed-size subsample of cells (nc_nCells). Benchmarks with the default parameters (10 networks of 500 cells) on simulated counts, measured on an Apple M2 Pro (16 GB RAM) with R 4.5 and its reference BLAS:
| Cells | Genes | Time |
|---|---|---|
| 300 | 1,000 | 35 s |
| 1,000 | 1,000 | 33 s |
| 1,000 | 5,000 | 8.4 min |
| 2,500 | 5,000 | 7.9 min |
| 5,000 | 5,000 | 8.0 min |
Before version 1.4.1, networks were built by fitting one SVD per gene, and the earlier benchmarks in this README (measured on a different machine) reported about 3 hours for 5,000 genes.
Peak memory grows with the square of the number of genes, because the networks are stacked into a tensor of genes x genes x nc_nNet entries. Estimated peak memory with 10 networks:
| Genes | Peak memory |
|---|---|
| 1,000 | 0.8 GB |
| 2,000 | 2.2 GB |
| 5,000 | 11.8 GB |
| 10,000 | 46 GB |
| 15,000 | 103 GB |
scTenifoldNet() compares this estimate with the memory available after quality control and warns, before building the networks, when it does not fit, reporting the largest number of genes that does. The estimate can also be checked beforehand with scTenifoldNet::checkMemory(nGenes, nNet = 10, nConditions = 2).
We create a sparse count matrix of 2,000 cells and 100 genes drawn from a negative binomial distribution (~67 % zeros). The last ten genes are prefixed with mt- to simulate mitochondrial genes.
library(scTenifoldNet)
nCells <- 2000
nGenes <- 100
set.seed(1)
X <- matrix(
rnbinom(n = nGenes * nCells, size = 20, prob = 0.98),
ncol = nCells
)
rownames(X) <- c(paste0("ng", 1:90), paste0("mt-", 1:10))
Genes 10, 2, and 3 are replaced with the expression profiles of genes 50, 11, and 5, respectively.
Y <- X
Y[10, ] <- Y[50, ]
Y[2, ] <- Y[11, ]
Y[3, ] <- Y[5, ]
Under H₀ (no differential regulation) we supply the same matrix for both conditions. Under H_A we supply the control and perturbed matrices.
outputH0 <- scTenifoldNet(X = X, Y = X,
nc_nNet = 10, nc_nCells = 500,
td_K = 3, qc_minLibSize = 30)
outputHA <- scTenifoldNet(X = X, Y = Y,
nc_nNet = 10, nc_nCells = 500,
td_K = 3, qc_minLibSize = 30)
The denoised networks are stored as sparse matrices that can be converted to igraph objects for downstream topological analysis:
igraph::graph_from_adjacency_matrix(outputHA$tensorNetworks$X, weighted = TRUE)
igraph::graph_from_adjacency_matrix(outputHA$tensorNetworks$Y, weighted = TRUE)
Under H₀, no genes reach significance (FDR < 0.1). Under H_A, the six genes involved in the perturbation are correctly identified:
head(outputH0$diffRegulation, n = 6)
# gene distance Z FC p.value p.adj
# 59 ng59 1.138659e-15 2.497859 4.338493 0.03725989 0.6411916
# 32 ng32 9.970451e-16 2.081719 3.326448 0.06817394 0.6411916
# 72 ng72 9.452170e-16 1.917581 2.989608 0.08380044 0.6411916
# 89 ng89 8.924023e-16 1.742758 2.664849 0.10258758 0.6411916
# 12 ng12 8.851226e-16 1.718018 2.621550 0.10542145 0.6411916
# 17 ng17 8.784453e-16 1.695183 2.582145 0.10807510 0.6411916
head(outputHA$diffRegulation, n = 6)
# gene distance Z FC p.value p.adj
# 2 ng2 0.023526702 2.762449 12.193413 0.0004795855 0.02414332
# 50 ng50 0.023514429 2.761550 12.180695 0.0004828665 0.02414332
# 11 ng11 0.022443941 2.681598 11.096894 0.0008647241 0.02882414
# 3 ng3 0.020263415 2.508478 9.045415 0.0026335445 0.06583861
# 10 ng10 0.019194561 2.417929 8.116328 0.0043868326 0.07711821
# 5 ng5 0.019079975 2.407977 8.019712 0.0046270923 0.07711821
Quantile-quantile plots provide a clear summary. Genes with FDR < 0.1 are highlighted in red.
par(mfrow = c(2, 1), mar = c(3, 3, 1, 1), mgp = c(1.5, 0.5, 0))
set.seed(1)
qChisq <- rchisq(100, 1)
geneColor <- rev(ifelse(outputH0$diffRegulation$p.adj < 0.1, "red", "black"))
qqplot(qChisq, outputH0$diffRegulation$FC, pch = 16, main = expression(H[0]),
col = geneColor, xlab = expression(chi^2 ~ Quantiles), ylab = "FC",
xlim = c(0, 8), ylim = c(0, 13))
qqline(qChisq)
legend("bottomright", legend = "FDR < 0.1", pch = 16, col = "red", bty = "n", cex = 0.7)
geneColor <- rev(ifelse(outputHA$diffRegulation$p.adj < 0.1, "red", "black"))
qqplot(qChisq, outputHA$diffRegulation$FC, pch = 16, main = expression(H[A]),
col = geneColor, xlab = expression(chi^2 ~ Quantiles), ylab = "FC",
xlim = c(0, 8), ylim = c(0, 13))
qqline(qChisq)
legend("bottomright", legend = "FDR < 0.1", pch = 16, col = "red", bty = "n", cex = 0.7)
If you use scTenifoldNet in your research, please cite:
Osorio, D., Zhong, Y., Li, G., Huang, J. Z., & Cai, J. J. (2020). scTenifoldNet: A Machine Learning Workflow for Constructing and Comparing Transcriptome-wide Gene Regulatory Networks from Single-Cell Data. Patterns, 1(9), 100139. doi:10.1016/j.patter.2020.100139
BibTeX:
@Article{osorio2020sctenifoldnet,
title = {scTenifoldNet: A Machine Learning Workflow for Constructing and
Comparing Transcriptome-wide Gene Regulatory Networks from
Single-Cell Data},
author = {Daniel Osorio and Yan Zhong and Guanxun Li and Jianhua Z. Huang
and James J. Cai},
journal = {Patterns},
year = {2020},
volume = {1},
number = {9},
pages = {100139},
issn = {2666-3899},
doi = {10.1016/j.patter.2020.100139},
}
© The Texas A&M University System. All rights reserved.