最近这个方向很火，前来学习报道一下，虚拟细胞敲除数据分析，已经有很多博主也出过教程了，我这边也只是参考官网学习一下，大家感兴趣的也可以去看原文。scTenifoldKnk，一个R/MATLAB/Python软件包，用于对单细胞基因调控网络进行虚拟敲除实验。scTenifoldKnk 是一种机器学习工作流程，利用野生型（WT）对照样本中的单细胞 RNA 测序（scRNAseq）数据作为输入，进行虚拟敲除实验。构建单细胞基因调控网络（scGRN），并通过将WT scGRN的邻接基质中敲除目标基因，将该基因的出度边设为零。scTenifoldKnk 随后将敲除的 scGRN 与 WT scGRN 进行比较，以识别差异调控的基因，称为虚拟敲除扰动基因，用于评估基因敲除的影响并揭示该基因在分析细胞中的功能。 scTenifoldKnk是一个基于机器学习的计算工具，能够仅利用野生型（Wild-Type, WT）样本的单细胞转录组数据，在计算机上模拟基因敲除实验，从而预测特定基因的功能。它的核心优势在于，无需进行真实、昂贵且耗时的动物敲除实验，即可系统性、高通量地评估基因功能。
它包含以下三个主要步骤：
构建野生型（WT）单细胞基因调控网络
方法首先利用你提供的WT样本单细胞表达矩阵，通过一个包含随机采样、主成分回归和基于CANDECOMP/PARAFAC (CP)张量分解去噪的流程，构建一个高质量、有方向、有权重的基因调控网络。这个网络不仅记录了基因之间的调控关系，还量化了相互作用的强度。
执行虚拟基因敲除（Virtual Knockout）
这是该工具的核心操作。在获得WT的调控网络矩阵（Wd）后，scTenifoldKnk会复制一份该矩阵，然后将目标基因所在的那一行所有元素全部置零（Wd˜）。这模拟了该基因无法再去调控其他基因的状态，即“被敲除”了。
比较网络并识别差异调控基因
现在有了Wd（WT网络）和Wd˜（伪KO网络）两个矩阵。scTenifoldKnk采用一种称为拟流形对齐（quasi-manifold alignment）的非线性学习方法，将两个网络中的所有基因投影到同一个低维空间。
在这个新空间中，每个基因都会有两个坐标，一个来自WT网络，一个来自伪KO网络。如果一个基因在两个网络中扮演的角色基本不变，它的两个坐标点就会非常接近；反之，如果它的调控关系因目标基因被敲除而发生了显著变化，它的两个点就会分得很开。
这个距离（distance）被计算出来并经过统计检验，最终得到差异调控基因（differentially regulated genes）列表。这些基因就是受虚拟敲除影响最显著的基因，通过对它们进行功能富集分析，就可以推断出被敲除的目标基因的可能功能。
scTenifoldKnk 所需的输入是一个表达矩阵，行中包含基因，列中包含单元格（条形码）。scTenifoldKnk的运行时间很大程度上取决于从子采样表达式矩阵构造scGRN所需的时间。时间的增加与输入数据中细胞和基因的数量成正比。scTenifoldKnk的Python版本可在以下网站获取：https://github.com/qwerty239qwe/scTenifoldpy
MATLAB版本可下载：https://github.com/jamesjcai/scGEAToolbox
scTenifoldKnk 的输出是一个包含 3 个槽位的列表，具体如下：张量网络：在CANDECOMP/PARAFAC（CP）张量分解后计算出的权重平均去噪声基因调控网络。它包含两个名额，包括：流形对齐：非线性流形对齐生成的低维特征。它是一个数据框架，行中基因数量是2倍，列中为d维（默认=2）diffRegulation：差分调控分析的结果。它是一个包含6列的数据框架，具体如下：gene：一个特征向量，基因ID从manifoldAlignment输出中识别。距离：在同一基因坐标之间计算出的欧几里得距离数值向量，适用于这两种条件。Z：Box-Cox幂变换后计算出的Z分数数值向量。p.value：与折叠变化相关的p值数值向量，概率符号为P[X>x]，采用卡方分布，自由度为一。p.adj：使用 Benjamini 和 Hochberg（1995）FDR 校正的调整后 p 值数值向量。# 1.安装包library(remotes)install_github('cailab-tamu/scTenifoldKnk')library(scTenifoldKnk)# 2.准备数据# 这是一个测试用的单细胞数据集scRNAseq <- system.file("single-cell/example.csv", package = "scTenifoldKnk")scRNAseq <- read.csv(scRNAseq, row.names = 1)# 查看数据的基本信息dim(scRNAseq) # 查看维度（基因数 x 细胞数）scRNAseq[1:5, 1:5] # 查看前5行5列# 假设你有Seurat对象 'seurat_obj'# 提取表达矩阵（counts或data均可）count_matrix <- GetAssayData(seurat_obj, layer = "counts") # Seurat v5语法# 或对于旧版本: count_matrix <- GetAssayData(seurat_obj, slot = "counts")# 确保格式正确：基因为行名，细胞为列名# 此时count_matrix已经可以直接用于scTenifoldKnk# 3.运行虚拟敲除分析# 运行scTenifoldKnk，这里我们敲除基因 "G100"# 注意：示例数据中的基因名是G100、G101等格式results <- scTenifoldKnk( countMatrix = scRNAseq, gKO = "G100", # 要敲除的目标基因 qc = TRUE, # 进行质量控制 qc_minLSize = 0, # 最小文库大小（示例数据中设为0） nc_nNet = 10, # 构建网络的数量 nc_nCells = 500, # 每次采样的细胞数 td_K = 3, # CP张量分解的秩 ma_nDim = 2, # 降维后的维度 nCores = 1 # 使用的CPU核心数)
# 4. 查看结果的整体结构str(results, max.level = 1)# 结果包含三个主要部分：# 1. tensorNetworks - 基因调控网络# 2. manifoldAlignment - 流形对齐结果# 3. diffRegulation - 差异调控分析结果（最重要）# 5. 提取差异调控结果diff_reg <- results$diffRegulation# 查看前10行head(diff_reg, 10)# 按p值排序，查看最显著受影响的基因diff_reg_sorted <- diff_reg[order(diff_reg$p.value), ]head(diff_reg_sorted, 10)# 6. 筛选FDR < 0.05的显著差异基因sig_genes <- diff_reg[diff_reg$p.adj < 0.05, ]nrow(sig_genes) # 查看有多少基因显著受敲除影响# 按距离排序，找出受影响最大的基因top_affected <- sig_genes[order(sig_genes$distance, decreasing = TRUE), ]head(top_affected, 20)
自动选择显著差异调控的基因；从WT网络中提取它们的相互作用；按权重分位数筛选重要边；可选地进行功能富集注释，在节点上显示分类信息# 7.使用plotKO函数绘制以敲除基因为中心的调控网络# 注意：这个函数在包的文档中有说明，可能需要从GitHub获取[citation:1]plotKO( X = results, # scTenifoldKnk的输出结果 gKO = "G100", # 被敲除的基因 q = 0.99, # 边权重的分位数阈值（只显示最强的1%边） annotate = TRUE, # 进行功能注释 fdrThreshold = 0.05 # FDR阈值)
# 提取top 20受影响基因top20 <- head(diff_reg_sorted, 20)# 绘制条形图library(ggplot2)ggplot(top20, aes(x = reorder(gene, -distance), y = distance, fill = -log10(p.adj))) + geom_bar(stat = "identity") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + labs( title = paste("Top 20 genes affected by", "G100", "knockout"), x = "Gene", y = "Distance (Regulatory Change)", fill = "-log10(FDR)" ) + scale_fill_gradient(low = "blue", high = "red")
# 创建火山图所需的列diff_reg$log10p <- -log10(diff_reg$p.adj)diff_reg$Significant <- diff_reg$p.adj < 0.05ggplot(diff_reg, aes(x = distance, y = log10p, color = Significant)) + geom_point(alpha = 0.6, size = 1) + scale_color_manual(values = c("grey", "red")) + theme_minimal() + labs( title = paste("Volcano plot for", "G100", "virtual knockout"), x = "Distance (Regulatory Change)", y = "-log10(Adjusted p-value)" ) + geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "blue") + theme(legend.position = "bottom")
# 提取显著差异基因列表sig_genes_list <- as.character(sig_genes$gene)# 使用clusterProfiler进行GO富集分析（需要安装）if (!require("clusterProfiler", quietly = TRUE)) { BiocManager::install("clusterProfiler")}library(clusterProfiler)library(org.Hs.eg.db) # 人类基因，请根据物种选择对应数据库# 转换基因名为ENTREZ ID（如果需要）# GO富集分析ego <- enrichGO( gene = sig_genes_list, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP", # Biological Process pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1)# 查看结果head(ego, 10)# 绘制气泡图dotplot(ego, showCategory = 15) + ggtitle(paste("GO Enrichment for", "G100", "knockout"))
文末有话说（你的关注、点赞、在看、评论非常重要！！！）欢迎大家留言讨论，各种合作欢迎咨询，联系管理员：kriswcyYQ。平台宗旨：提供（空间）单细胞组学、常规组学个性化解决方案；提供数据分析、R及Python绘图优化、代码复现等科研服务；提供数据分析课程服务（包括基于AI的深度学习课程、影像组学等），欢迎大家前来咨询。目前冲5000粉，免费数据分析咨询一次！！！
10个月宝宝每天需要喝多少奶粉？