还来得及！虚拟敲除(1)-scTenifoldKnk包scRNA数据单基因虚拟敲除【详细解释】

KS科研分享与服务-TS的美梦

发布于 2026-06-12 14:29:05

370

偷偷问一下，关注了吗

从最简单的开始虚拟敲除这个热点，还来得及，来得及！

1scTenifoldKnk简介

现在（2025年中-2026年5月）铺天盖地的热度都是虚拟敲除（Virtual Knockout），相信对于scTenifoldKnk也不陌生了。总体而言，在众多虚拟敲除工具中，scTenifoldKnk是一款比较简单的虚拟敲除工具，所以介绍最多，使用范围最广，引用很多，但并不是说它最优秀，最值得信赖，这一点需要清楚。scTenifoldKnk文章发的很早，近两年热度上来之后频繁使用。scTenifoldKnk是一个基于单细胞RNA测序数据执行“虚拟”基因敲除实验的机器学习工作流，它以野生型（WT）对照样本的单细胞RNA 测序数据为输入，执行虚拟基因敲除实验。该方法首先构建单细胞基因调控网络（scGRN），然后通过将目标基因在WT scGRN邻接矩阵中的出度边权（outdegree edges）设为0，来实现对该基因的“敲除”。接着，scTenifoldKnk 将敲除后的scGRN与WT的scGRN进行比较，识别出差异调控的基因，这些基因被称为虚拟敲除扰动基因（perturbed genes）。这些基因可用于评估目标基因敲除的影响，并揭示该基因在所分析细胞中的功能。

github官网：https://github.com/cailab-tamu/scTenifoldKnk

paper链接：https://www.sciencedirect.com/science/article/pii/S2666389922000010

安装包，可能需要先安装依赖包scTenifoldNet：

install.packages("scTenifoldNet")
library(remotes)
install_github('cailab-tamu/scTenifoldKnk')
library(scTenifoldKnk)

1.1scTenifoldKnk虚拟敲除原理简述

ScTenifoldKnk的大致流程是以来自野生型（WT）样本的基因-细胞表达矩阵作为输入，不需要真实的敲除样本数据。该工作流程首先从输入表达矩阵构建野生型单细胞基因调控网络（scGRN），然后通过从野生型scGRN中敲除某个基因来生成一个伪敲除（pseudo-KO）scGRN。最后，它采用网络比较方法比较伪敲除scGRN和野生型scGRN，以识别差异调控（ifferentially regulated，DR）基因。这些DR基因即为虚拟敲除扰动基因。通过对这些虚拟敲除扰动基因进行功能富集分析，可以推断被敲除基因（即被虚拟敲除的基因）的功能。

具体而言，其核心原理分为以下三个主要步骤：

构建野生型基因调控网络 (Constructing WT scGRN)

输入野生型样本的单细胞RNA测序基因表达矩阵。通过随机下采样产生多个细胞子集，并对每个子集运行主成分回归，构建出多个初始的基因调控网络，储存为一个signed, weighted, directional graph，用pxp的邻接矩阵Wi来表示。最后利用张量分解（Tensor Decomposition）技术对这些网络进行去噪和融合，最终生成一个稳定、加权的野生型基因调控网络。

生成“伪敲除”网络 (Generating Pseudo-KO scGRN)

在得到的野生型邻接矩阵中，找到目标敲除基因所在的行，将该基因对应行中的所有权值全部设为0。这一步在计算上模拟了该基因功能的缺失，即切断了该基因对网络中所有其他基因的输出调控作用，从而产生一个“伪敲除”网络。

流形对齐与差异评估 (Manifold Alignment)

使用一种称为准流形对齐（Quasi-manifold alignment)的方法，将野生型网络和伪敲除网络进行对齐，两个scGRN中包含的所有基因都被投影到k维空间中，投影后，每个基因都有两个低维表示：一个对应于Wd（WT scGRN），另一个对应于~Wd（对应于Pseudo-KO scGRN）。然后计算每个基因在两个网络投影之间的欧几里得距离，距离越大的基因，说明其在敲除目标基因后受到的调控影响越显著。根据距离值对基因进行排序，生成一个排序基因列表，作为基因集富集分析（GSEA）的输入，最后，应用卡方检验来识别显著的DR基因，即虚拟敲除（KO）扰动基因。

1.2关于运行时间

scTenifoldKnk因为需要进行网络构建和重复的统计模拟，所以比较耗费时间。至于影响运行时间的因素，作者已经做了比较，参照下面的表格，不同的基因数和细胞数在默认参数下的运行时间，时间的增加与用作输入的数据集中的细胞和基因的数量成正比。总体而言，基因数的影响更大，因为基因数的多少决定构建网络的大小，相同的细胞数，基因数量增加，运行时间显著延长。

2单细胞数据准备

演示的数据来源于GEO公开数据库GSE293902（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE293902）,目前没有看到文章信息。数据基本信息是来源于小鼠脾脏中分离的Germinal centers (GC) B细胞，并通过流式细胞荧光分选技术进行分选，最后对4个对照（WT）sample，4个实验组（Dtx1基因敲除）样本获取的单细胞转录组测序数据。这里WT和实验组各选了2个sample进行演示。接下来就是正常的seurat流程，按照一般的质控进行了分析，并去除了双细胞，最后选择两个cluster演示敲除。

library(Seurat)
library(dplyr)
library(patchwork)
#===============================================================================
setwd("/home/tq_ziv/data_analysis/虚拟敲除/scRNA_knock/")
#读入数据创建seurat obj
paths <- list(
  WT1  = "WT1",
  WT2  = "WT2",
  DTX1 = "DTX1",
  DTX2 = "DTX2"
)
objs <- lapply(names(paths), function(sample){
  dat <- Read10X(data.dir = paths[[sample]])
  seu <- CreateSeuratObject(
    counts = dat,
    project = sample,
    min.cells = 3,
    min.features = 200
  )
  seu$sample   <- sample
  seu$genotype <- ifelse(grepl("^WT", sample), "WT", "DTX")
  return(seu)
})
names(objs) <- names(paths)
#QC
for(x in names(objs)){
  objs[[x]][["percent.mt"]] <- PercentageFeatureSet(objs[[x]], pattern = "^mt-")
  objs[[x]] <- subset(objs[[x]], subset = nCount_RNA < 50000 & nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 20)
}
#SCT+整合
sce_ko <- merge(x = objs[[1]],y = objs[2: length(objs)],add.cell.ids = names(objs))
sce_ko <- SCTransform(sce_ko)
sce_ko <- RunPCA(sce_ko, verbose=F)
sce_ko <- IntegrateLayers(object = sce_ko, method = CCAIntegration, normalization.method = "SCT", verbose = F)
sce_ko[["RNA"]] <- JoinLayers(sce_ko[["RNA"]])
#降维聚类
sce_ko <- RunUMAP(sce_ko, dims = 1:20)
sce_ko <- FindNeighbors(sce_ko, reduction = "integrated.dr", dims = 1:20)
sce_ko <- FindClusters(sce_ko, resolution = 0.4)
#===============================================================================
#去除双细胞

library(DoubletFinder)
seurat_object <- SplitObject(sce_ko,split.by ='sample')
seurat_object[[1]] <- ks_detectDoublet(seurat_object[[1]],dims = 1:30,estDubRate=0.065,
                                       ncores = 2, SCTransform=T, Homotypic=F, annotation="seurat_clusters")
seurat_object[[2]] <- ks_detectDoublet(seurat_object[[2]],dims = 1:30,estDubRate=0.065,
                                       ncores = 2, SCTransform=T, Homotypic=F, annotation="seurat_clusters")
seurat_object[[3]] <- ks_detectDoublet(seurat_object[[3]],dims = 1:30,estDubRate=0.065,
                                       ncores = 2, SCTransform=T, Homotypic=F, annotation="seurat_clusters")
seurat_object[[4]] <- ks_detectDoublet(seurat_object[[4]],dims = 1:30,estDubRate=0.065,
                                       ncores = 2, SCTransform=T, Homotypic=F, annotation="seurat_clusters")
for (i in seq_along(seurat_object)) {

  seurat_object[[i]] <- subset(seurat_object[[i]], subset = (DF.classify == "Singlet"))

}
#===============================================================================
#合并重新分析
#SCT+整合
sce_konck <- merge(x = seurat_object[[1]],y = seurat_object[2: length(seurat_object)])
sce_konck <- SCTransform(sce_konck)
sce_konck <- RunPCA(sce_konck, verbose=F)
sce_konck <- IntegrateLayers(object = sce_konck, method = CCAIntegration, normalization.method = "SCT", verbose = F)
sce_konck[["RNA"]] <- JoinLayers(sce_konck[["RNA"]])
#降维聚类
sce_konck <- RunUMAP(sce_konck, dims = 1:25,n.neighbors = 60,min.dist = 0.3)
sce_konck <- FindNeighbors(sce_konck, reduction = "integrated.dr", dims = 1:25)
sce_konck <- FindClusters(sce_konck, resolution = 0.2)
# DimPlot(sce_konck,label = T)
# FeaturePlot(sce_konck,features = c("Dtx1"))
# VlnPlot(sce_konck,features = "Dtx1",group.by = 'sample',pt.size = 0.1)
# DotPlot(sce_konck,features = "Dtx1",group.by = 'sample')
#提取比较聚集的细胞
seu <- subset(sce_konck, idents = c("0", "1"))
DimPlot(seu,label = T)
selected_cells <- CellSelector(plot = DimPlot(seu, label = TRUE))
sce_KO <- subset(seu, cells = selected_cells)
sce_KO@meta.data$seurat_clusters = droplevels(sce_KO@meta.data$seurat_clusters, 
                                       exclude = setdiff(levels(sce_KO@meta.data$seurat_clusters),
                                                         unique(sce_KO@meta.data$seurat_clusters)))
save(sce_KO,file = 'sce_KO.RData')

3、scTenifoldKnk虚拟敲除

3.1、虚拟敲除函数参数解释

函数相关参数解析。大多数一般使用默认参数，关于网络构建的某些参数，值越大可能对于结果更精确，但是耗费的时间也会大大增加。

#核心函数参数解释
scTenifoldKnk(
  countMatrix,#counts表达矩阵，行是基因，列是sample/cellid
  qc = TRUE,#是否使用函数内置的QC流程，QC指标有线粒体比例 (qc_mtThreshold)；library size (qc_minLSize)；基因在多少细胞中表达 (qc_minCells)，这些参数都可以设置，如果是Seurat流程已经完成了相关质控流程，这一步可以不需要。也可以自行调整。
  gKO = NULL,#需要敲除的基因，需要保证名称正确，包含在countMatrix中

  #QC相关参数

  qc_mtThreshold = 0.1,#如QC为TRUE，单个sample/cell线粒体基因比例。
  qc_minLSize = 1000,#如QC为TRUE，library size的最小值，即每个细胞总UMI数，低于该值的细胞被滤掉。
  qc_minCells = 25,#如QC为TRUE，某个基因检测到有表达的最小细胞数。

  #网络相关参数

  nc_lambda = 0,#0–1 之间的小数，用来增强网络的“方向性”。默认0（不调整方向性）。
  nc_nNet = 10,#基于主成分回归生成的网络数量（即要构建多少张网络，最后整合），越大估计更稳健。
  nc_nCells = 500,#每次子采样用于建一张网络的细胞数。
  nc_nComp = 3,#整数值。PCA中生成网络的主成分的数量。应该大于2，并且低于基因总数。
  nc_scaleScores = TRUE,#是否把网络权重缩放到最大绝对值为 1。
  nc_symmetric = FALSE,#是否强制输出对称的网络（无向）。
  nc_q = 0.9,#0-1之间的小数，保留前百分之多少最强关系，越大，网络越密。

  #张量分解相关参数

  td_K = 3,#CP分解中的秩（rank），一般3-5即可。
  td_maxIter = 1000,#CP分解迭代的最大迭代次数。
  td_maxError = 1e-05,
  td_nDecimal = 3,#内部结果保留的小数位数
  ma_nDim = 2,
  nCores = parallel::detectCores()#并行计算线程数
)

3.2运行虚拟敲除

提取WT数据表达用于虚拟分析。

#提取WT数据进行虚拟敲除
sce_WT <- subset(sce_KO,genotype=='WT')

#获取表达矩阵
WT_countMatrix <- GetAssayData(sce_WT,layer = 'counts',assay = 'SCT')

过滤一下基因，我这里并没有过滤太多，保证演示比较快速，但是实际情况中，不应该对基因进行过滤，如果过滤了，那么构建的网络就不准确了，结果也会有点问题。

#过滤基因，为了演示减少基因数目，那些在很少细胞中检测到的基因过滤
#减少矩阵
minCells <- rowSums(WT_countMatrix)
countMatrix <- WT_countMatrix[minCells >= 100,]
dim(countMatrix)

## [1] 9265 3588

#精确匹配检测目的基因是否存在于矩阵中
gene_name <- "Dtx1"
is_present <- gene_name %in% rownames(countMatrix)
if(is_present) {
  cat(gene_name, "is detcted in matrix\n")
} else {
  cat(gene_name, "is not detcted in matrix\n")
}

运行敲除:

KE <- scTenifoldKnk(countMatrix = countMatrix, gKO = 'Dtx1',qc = F, nCores = 8)
#save(KE,file = 'Dtx1_Vknock.RData')
#运行日志
── scTenifoldKnk ────────────────────────────────────────────────────────────────────────────
ℹ Simulating Dtx1 gene knockout
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=08h 32m 30s
✔ Network construction complete (nNet = 10, nCells per net = 500)
  |===================================================================================| 100%
✔ Tensor decomposition completed (K = 3)
✔ Prepared WT adjacency matrix for KO simulation
✔ Simulated knockout for Dtx1
✔ Manifold alignment completed (d = 2)
✔ Differential regulation computed for Dtx1
✔ Finished scTenifoldKnk for Dtx1

3.3查看虚拟敲除结果

可以看到，我演示的9265个基因3588个细胞的矩阵运行时间是08h 32m 30s。得到的分析结果是一个list，tensorNetworks是经过ANDECOMP/PARAFAC (CP) tensor decomposition计算分析后的WT和KO网络。manifoldAlignment是非线性流形对齐生成的、用于低维特征的降维结果。它是一个数据框，其行数为基因数的2倍，列数为d维（默认 d = 2）。diffRegulation是本次虚拟敲除后差异调控分析的结果。

names(KE)

## [1] "tensorNetworks"    "manifoldAlignment" "diffRegulation"

请注意，在diffRegulation文件中，这个差异结果不是传统转录组意义上的差异分析结果，比如FC不是logFC，而是看哪些基因显著受到虚拟敲除的“扰动”，所以FC没有符号。本次数据集的演示结果不是很好，显著性结果很少。

head(KE[['diffRegulation']])

##       gene     distance        Z           FC      p.value        p.adj
## 2556  Dtx1 6.489131e-06 5.952583 1.026183e+05 0.000000e+00 0.000000e+00
## 8701   Fau 1.904783e-06 5.188502 8.841842e+03 0.000000e+00 0.000000e+00
## 7244 Rpl37 3.343029e-07 4.157130 2.723529e+02 3.484358e-61 1.076086e-57
## 6960 Rps24 7.088095e-08 3.287898 1.224364e+01 4.668456e-04 9.999996e-01
## 8302 Rpl36 6.941641e-08 3.276507 1.174292e+01 6.107517e-04 9.999996e-01
## 2050  Cd52 4.010504e-08 2.980083 3.919672e+00 4.772420e-02 9.999996e-01

scTenifoldKnk也提供了可视化函数plotKO，从scTenifoldKnk的输出结果中生成并绘制以敲除基因为中心的子网络。该函数首先筛选出差异调控显著的基因，然后从重建的野生型网络中提取这些基因之间的相互作用，再根据权重的分位数对边进行过滤，最后展示该网络。当annotate =TRUE时，该函数会查询富集分析数据库，并在网络图的节点上叠加表示功能类别的饼图，同时显示显著功能条目的图例。annotate =False，则只绘图，不进行富集分析。

scTenifoldKnk::plotKO(KE, gKO = 'Dtx1',annotate = F)

scTenifoldKnk虚拟敲除的原理和演示就到此结束，因为结果不好，所以不再做对比，实验敲除和虚拟敲除之间的差别还是非常大的。对于好的结果，可以对虚拟扰动比较大的基因进行排序，进行功能富集分析，筛选靶点，为自己研究增色。觉得我们分享有用的，点个赞再走呗！

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-06-09，如有侵权请联系 cloudcommunity@tencent.com 删除

排序

本文分享自 KS科研分享与服务微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度