Package {scTenifoldKnk}

Type:	Package
Title:	In-Silico Knockout Experiments from Single-Cell Gene Regulatory Networks
Version:	2.0.0
Description:	A workflow based on 'scTenifoldNet' to perform in-silico knockout experiments using single-cell RNA sequencing (scRNA-seq) data from wild-type (WT) control samples as input. First, the package constructs a single-cell gene regulatory network (scGRN) and knocks out a target gene from the adjacency matrix of the WT scGRN by setting the gene’s outdegree edges to zero. Then, it compares the knocked out scGRN with the WT scGRN to identify differentially regulated genes, called virtual-knockout perturbed genes, which are used to assess the impact of the gene knockout and reveal the gene’s function in the analyzed cells. It also predicts the direction (up or down) of the response of each gene from the WT expression, and reads all knockouts of a network from a single heat kernel, which makes transcriptome-wide knockout screens practical.
URL:	https://github.com/cailab-tamu/scTenifoldKnk
BugReports:	https://github.com/cailab-tamu/scTenifoldKnk/issues
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
RoxygenNote:	7.3.3
Imports:	Matrix, methods, stats, MASS, scTenifoldNet (≥ 1.4.3), cli, enrichR, igraph, reshape2, grDevices, graphics
Suggests:	testthat (≥ 2.1.0), locfdr
NeedsCompilation:	no
Packaged:	2026-10-08 13:13:03 UTC; runner
Author:	Daniel Osorio [aut, cre], Yan Zhong [aut, ctb], Guanxun Li [aut, ctb], Qian Xu [aut, ctb], Yongjian Yang [aut, ctb], Yanan Tian [aut, ctb], Robert Chapkin [aut, ctb], Jianhua Huang [aut, ctb], James J. Cai [aut, ctb, ths]
Maintainer:	Daniel Osorio <dcosorioh@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-09 13:10:02 UTC

Evaluates gene differential regulation based on manifold alignment distances.

Description

Using the output of the non-linear manifold alignment, this function computes the Euclidean distance between the coordinates for the same gene in both conditions. Calculated distances are then transformed using Box-Cox power transformation, and standardized to ensure normality. P-values are assigned following the chi-square distribution over the fold-change of the squared distance computed with respect to the expectation, the mean squared distance of the genes that were not knocked out (gKO), or, when empiricalNull = TRUE, using Efron's empirical null estimated from the Z-scores with locfdr. Genes whose distance is at the level of floating-point noise (at most sqrt(.Machine$double.eps) times the largest absolute coordinate) did not move between conditions and get a p-value of 1; if this applies to every gene, for example when the knockout leaves the network unchanged, a warning is raised.

Usage

dRegulation(
  manifoldOutput,
  gKO = NULL,
  empiricalNull = FALSE,
  direction = NULL
)

Arguments

Value

A data frame with 6 columns (8 when direction is given) as follows:

• gene A character vector with the gene id identified from the manifoldAlignment output.

• distance A numeric vector of the Euclidean distance computed between the coordinates of the same gene in both conditions.

• Z A numeric vector of the Z-scores computed after Box-Cox power transformation.

• FC A numeric vector of the FC computed with respect to the expectation, the mean squared distance of the genes not in gKO.

• p.value A numeric vector of the p-values associated to the fold-changes, probabilities are assigned as P[X > x] using the Chi-square distribution with one degree of freedom.

• p.adj A numeric vector of adjusted p-values using Benjamini & Hochberg (1995) FDR correction.

• direction Only when direction is given: "up" or "down", the predicted direction of the change of the gene; the knocked-out genes are "down" by construction.

• directionScore Only when direction is given: the signed direction score (NA for the knocked-out genes).

References

• Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics, 29, 1165-1188. doi: 10.1214/aos/1013699998.

Examples

library(scTenifoldKnk)

# Simulating a dataset following a negative binomial distribution with high sparsity (~67%)
nCells = 2000
nGenes = 100
set.seed(1)
X <- rnbinom(n = nGenes * nCells, size = 20, prob = 0.98)
X <- round(X)
X <- matrix(X, ncol = nCells)
rownames(X) <- c(paste0('ng', 1:90), paste0('mt-', 1:10))

# Performing single-cell quality control
qcOutput <- scQC(
  X = X,
  minLibSize = 30,
  removeOutlierCells = TRUE,
  minPCT = 0.05,
  maxMTratio = 0.1
)

# Computing 3 gene regulatory networks from subsamples of 500 cells
xNetworks <- scTenifoldNet::makeNetworks(
  X = qcOutput,
  nNet = 3,
  nCells = 500,
  nComp = 3,
  scaleScores = TRUE,
  symmetric = FALSE,
  q = 0.95
)

# Computing a K = 3 CANDECOMP/PARAFAC (CP) Tensor Decomposition
tdOutput <- scTenifoldNet::tensorDecomposition(xNetworks, K = 3, maxError = 1e5, maxIter = 1e3)

## Not run: 
# Computing the manifold alignment
maOutput <- scTenifoldNet::manifoldAlignment(tdOutput$X, tdOutput$X)

# Evaluating differential regulation
drOutput <- dRegulation(maOutput)
head(drOutput)

# Plotting — genes with FDR < 0.05 colored in red
geneColor <- ifelse(drOutput$p.adj < 0.05, 'red', 'black')
qqnorm(drOutput$Z, main = 'Standardized Distance', pch = 16, col = geneColor)
qqline(drOutput$Z)

## End(Not run)

Spectral heat kernel of a gene-gene matrix

Description

Computes the heat kernel

H = \sum_k \exp(t (\lambda_k / \lambda_{max} - 1)) v_k v_k^T

from the eigendecomposition X = \sum_k \lambda_k v_k v_k^T of a symmetric gene-gene matrix (a gene regulatory network or a gene-gene correlation matrix). Row i of H describes how a change in gene i diffuses over the whole network, summing paths of every length; larger t concentrates the kernel on the dominant network modules.

Usage

heatKernel(X, t = 10, symmetric = TRUE)

Arguments

Value

A square numeric matrix with the same dimension names as X.

Examples

set.seed(1)
A <- matrix(rnorm(25), 5, 5, dimnames = list(letters[1:5], letters[1:5]))
H <- heatKernel(A, t = 10)
isSymmetric(H)

Heat-kernel (hk) manifold alignment for in-silico knockouts

Description

Kernel counterpart of the non-linear manifold alignment used by scTenifoldKnk. The heat kernel of the WT network is computed once (heatKernel), and the knockout of each gene (or gene set) is read from its rows: removing gene x lowers it by its mean expression and the change diffuses over the network,

\Delta_g = -\frac{mean(x)}{sd(x)} H[x, g] \, sd(g)

on log1p(CPM) expression (summed over the genes of a multi-gene knockout). The absolute value |\Delta_g| plays the role of the manifold alignment distance: it ranks the perturbed genes like the manifold alignment does, without recomputing an alignment for every knockout, which makes transcriptome-wide perturbation efficient.

Usage

hkManifoldAlignment(WT, X, gKO = NULL, t = 10, H = NULL)

Arguments

Value

A numeric matrix of perturbation distances |\Delta_g| (for the knocked-out genes themselves, their mean log1p(CPM) expression, the amount they lose) with one row per knockout (named by the knocked-out genes, joined with "+" for multi-gene knockouts) and one column per gene of the network.

See Also

heatKernel, knockoutDirection

Direction of the response to an in-silico knockout

Description

Predicts whether each gene goes up or down after knocking out gKO, using only the WT expression data: the heat kernel (heatKernel) of the gene-gene Pearson correlation matrix of log1p(CPM) expression is diffused from the knocked-out gene(s),

s_g = -\sum_{x \in gKO} \frac{mean(x)}{sd(x)} H[x, g] \, sd(g)

, and the sign of s_g is the predicted direction.

Evaluated against bulk knockdown/knockout profiles of the same cell lines, this direction is correct more often than chance, but it mostly reflects the response shared by most knockdowns along the dominant WT expression program rather than regulation specific to gKO, and its accuracy varies between cell types and WT data sets.

In single-cell data, log1p(CPM) expression often keeps a dependence on sequencing depth that makes nearly all genes correlate positively; the heat kernel then follows that axis and predicts almost every gene to go down. With regressLibSize = TRUE, log library size is regressed out of each gene before computing the correlations, which removes the depth axis. It is off by default: it corrected tissue data where nearly all genes were predicted down, but in cell lines it left the direction unchanged or made it slightly worse, because the expression that follows library size can be biological.

Usage

knockoutDirection(X, gKO, genes = rownames(X), t = 5, regressLibSize = FALSE)

Arguments

Value

A numeric matrix of direction scores s_g (positive: predicted up, negative: predicted down) with one row per knockout and one column per gene.

See Also

heatKernel, dRegulation

Map of transcriptome-wide perturbation profiles

Description

Places every virtual knockout from a transcriptome-wide run of scTenifoldKnk on a two-dimensional map together with a reference signature, such as a disease-vs-control expression signature, and highlights the genes whose knockout reproduces it. Each knockout is described by its signed perturbation profile: the rank of the perturbation distance of every gene, scaled to (0, 1], multiplied by its predicted direction. The similarity of a profile to the signature is the cosine between both over the genes they share; the knocked-out gene itself is left out of its own profile.

Usage

perturbationMap(
  X,
  signature,
  genes = NULL,
  layout = c("signature", "pca"),
  nLabels = 10,
  plot = TRUE
)

Arguments

Value

Invisibly, a data frame with one row per knockout, sorted by decreasing similarity: gene, similarity (cosine to the signature), rank (1 = most similar), percentile (fraction of knockouts at least as similar), and the map coordinates x and y. With layout = "pca", the coordinates of the signature are stored in the "signature" attribute.

Examples

library(scTenifoldKnk)

# Counts with three gene modules
set.seed(1)
modules <- matrix(rgamma(3 * 600, shape = 2), 3, 600)
loadings <- matrix(0, 60, 3)
loadings[cbind(1:60, rep(1:3, length.out = 60))] <- runif(60, 0.5, 2)
X <- matrix(rpois(60 * 600, lambda = 5 * loadings %*% modules), 60, 600)
rownames(X) <- paste0('g', 1:60)
colnames(X) <- paste0('c', 1:600)

# Knock out every gene of the network
O <- scTenifoldKnk(X, transcriptomeWide = TRUE, qc = FALSE,
                   nc_nNet = 3, nc_nCells = 300, nCores = 1)

# A signature that resembles the knockout of g1
sig <- O$perturbationDistances['g1', ] * O$perturbationDirections['g1', ]
sig <- sig + rnorm(length(sig), sd = sd(sig))

M <- perturbationMap(O, signature = sig, genes = 'g1')
head(M)

Plot KO network

Description

Generate and plot a KO-centered subnetwork from the output of scTenifoldKnk. The function selects genes with significant differential regulation (FDR < 0.05), extracts their interactions from the reconstructed WT network, filters edges by weight quantile, and displays the network using igraph. When annotate = TRUE the function queries enrichment databases via enrichR and overlays category pies on nodes with a legend of significant terms.

Usage

plotKO(
  X,
  gKO,
  q = 0.99,
  annotate = TRUE,
  nCategories = 20,
  fdrThreshold = 0.05
)

Arguments

Value

Invisibly returns NULL. Called for the side effect of plotting the network.

Examples

## Not run: 
library(scTenifoldKnk)

# Load example data
scRNAseq <- system.file("single-cell/example.csv", package = "scTenifoldKnk")
scRNAseq <- read.csv(scRNAseq, row.names = 1)

# Run scTenifoldKnk
output <- scTenifoldKnk(countMatrix = scRNAseq, gKO = "G100", qc_minLibSize = 0)

# Plot the KO-centered subnetwork with enrichment annotation
plotKO(output, gKO = "G100")

# Plot without enrichment annotation
plotKO(output, gKO = "G100", annotate = FALSE)

## End(Not run)

Performs single-cell data quality control

Description

This function performs quality control filters over the provided input matrix. It checks for minimum cell library size, mitochondrial ratio, outlier cells, and the fraction of cells where a gene is expressed.

Usage

scQC(
  X,
  minLibSize = 1000,
  removeOutlierCells = TRUE,
  minPCT = 0.05,
  maxMTratio = 0.1,
  label = NULL
)

Arguments

Value

A dgCMatrix object with the cells and the genes that pass the quality control filters.

References

Ilicic, Tomislav, et al. "Classification of low quality cells from single-cell RNA-seq data." Genome biology 17.1 (2016): 29.

Examples

library(scTenifoldKnk)

# Simulating a dataset following a negative binomial distribution with high sparsity (~67%)
nCells = 2000
nGenes = 100
set.seed(1)
X <- rnbinom(n = nGenes * nCells, size = 20, prob = 0.98)
X <- round(X)
X <- matrix(X, ncol = nCells)
rownames(X) <- c(paste0('ng', 1:90), paste0('mt-', 1:10))

# Performing single-cell quality control
qcOutput <- scQC(
  X = X,
  minLibSize = 30,
  removeOutlierCells = TRUE,
  minPCT = 0.05,
  maxMTratio = 0.1
)

# Comparing dimensions before and after QC
dim(X)
dim(qcOutput)

scTenifoldKNK

Description

Predict gene perturbations using in-silico knockout experiments from single-cell gene regulatory networks.

Usage

scTenifoldKnk(
  countMatrix,
  gKO = NULL,
  transcriptomeWide = FALSE,
  qc = TRUE,
  qc_minLibSize = 1000,
  qc_removeOutlierCells = TRUE,
  qc_minPCT = 0.05,
  qc_maxMTratio = 0.1,
  nc_lambda = 0,
  nc_nNet = 10,
  nc_nCells = 500,
  nc_nComp = 3,
  nc_scaleScores = TRUE,
  nc_symmetric = FALSE,
  nc_q = 0.9,
  nc_priorNetwork = NULL,
  td_K = 3,
  td_maxIter = 1000,
  td_maxError = 1e-05,
  td_nDecimal = 3,
  ma_nDim = 2,
  ma_method = NULL,
  ma_heatT = 10,
  dr_empiricalNull = FALSE,
  dr_direction = TRUE,
  dr_directionT = 5,
  dr_directionRegressLibSize = FALSE,
  nCores = parallel::detectCores(),
  seed = 1
)

Arguments

Value

In single knockout mode (transcriptomeWide = FALSE), a list with 3 slots as follows:

• tensorNetworks: A list with the WT weight-averaged denoised gene regulatory network. The KO network is the WT network with the rows of gKO set to 0; it is not returned, which roughly halves the size of the output.

• manifoldAlignment: The generated low-dimensional features result of the non-linear manifold alignment (only when ma_method = "manifold").

• diffRegulation: The results of the differential regulation analysis (see dRegulation), including the direction and directionScore columns when dr_direction = TRUE.

In transcriptome-wide mode (transcriptomeWide = TRUE), a list with 2 slots as follows:

• tensorNetworks: A list with the WT weight-averaged denoised gene regulatory network.

• perturbationDistances: A numeric matrix of perturbation distances (heat manifold alignment by default, manifold alignment when ma_method = "manifold") with the perturbed genes as rows and all genes in the WT network as columns.

• perturbationDirections: Only when dr_direction = TRUE: a numeric matrix of the same dimensions with the predicted direction of each gene (1 up, -1 down, 0 undetermined).

Author(s)

Daniel Osorio <dcosorioh@gmail.com>

Examples

library(scTenifoldKnk)

# Simulating a dataset following a negative binomial distribution with high sparsity (~67%)
nCells = 2000
nGenes = 100
set.seed(1)
X <- rnbinom(n = nGenes * nCells, size = 20, prob = 0.98)
X <- round(X)
X <- matrix(X, ncol = nCells)
rownames(X) <- c(paste0('ng', 1:90), paste0('mt-', 1:10))

## Not run: 
# Running scTenifoldKnk — simulating knockout of gene ng10
output <- scTenifoldKnk(
  countMatrix = X,
  gKO = "ng10",
  nc_nNet = 10,
  nc_nCells = 500,
  td_K = 3,
  qc_minLibSize = 30
)

# Structure of the output
str(output)

# Accessing the WT gene regulatory network, and rebuilding the KO network
dim(output$tensorNetworks$WT)
KO <- output$tensorNetworks$WT
KO['ng10', ] <- 0

# Accessing the manifold alignment result
head(output$manifoldAlignment)

# Differential regulation results — top perturbed genes
head(output$diffRegulation, n = 10)

# Plotting the KO-centered subnetwork
plotKO(output, gKO = "ng10")

# Multi-gene knockout — ng10 and ng20 knocked out together
dkoOutput <- scTenifoldKnk(
  countMatrix = X,
  gKO = c("ng10", "ng20"),
  nc_nNet = 10,
  nc_nCells = 500,
  td_K = 3,
  qc_minLibSize = 30
)
head(dkoOutput$diffRegulation, n = 10)
plotKO(dkoOutput, gKO = c("ng10", "ng20"))

# Transcriptome-wide perturbation — knock out every gene in the WT network
twOutput <- scTenifoldKnk(
  countMatrix = X,
  transcriptomeWide = TRUE,
  nc_nNet = 10,
  nc_nCells = 500,
  td_K = 3,
  qc_minLibSize = 30
)

# Matrix of distances: perturbed genes (rows) by all genes (columns)
dim(twOutput$perturbationDistances)
twOutput$perturbationDistances[1:5, 1:5]

## End(Not run)