本节概览：
- 虚拟基因敲除简介：了解什么是虚拟敲除及其与实验敲除的本质区别。
- scTenifoldKnk原理：通过比较“野生型”和“敲除型”基因共表达网络，识别受扰动影响的基因。
- 数据准备：从Seurat对象出发，提取表达矩阵并筛选目标基因。
- 核心分析：使用scTenifoldKnk进行虚拟敲除，计算差异调控基因。
- 结果解读：识别受目标基因调控的下游响应基因。
- 结果可视化：条形图、火山图、MA图、点图全面展示扰动效应。
在前几节中，我们完成了细胞类型注释、差异表达分析、功能富集和转录因子调控网络分析。这些分析帮助我们回答了“细胞是什么”、“它们表达了什么差异基因”、“谁在驱动这些变化”等问题。
但在实际研究中，我们经常会遇到一个更深层次的问题：“如果我把某个基因删掉，会有什么后果？”
在传统实验中，回答这个问题需要CRISPR-Cas9基因敲除、RNAi干扰或条件性基因敲除小鼠——这些实验耗时耗力、成本高昂。但在单细胞数据分析中，我们是否可以通过计算模拟来预测基因敲除的效应呢？
scTenifoldKnk 正是为此而生的。它通过比较“野生型”和“虚拟敲除型”的基因共表达网络，预测某个基因被扰动后对整个调控网络的影响。简单来说，它在计算机上“删除”一个基因，然后观察哪些基因的调控关系发生了改变。
本文将带领大家从已注释好的Seurat对象出发，使用scTenifoldKnk完成虚拟基因敲除分析的完整流程。
一、虚拟基因敲除的核心概念
1.1 什么是虚拟敲除（in silico knockout）？
虚拟敲除是指利用计算模拟预测基因功能缺失的后果，而不进行实际的生物学实验。其核心逻辑是：
- 基于单细胞表达数据构建野生型基因共表达网络
- 从表达矩阵中移除目标基因
- 基于剩余基因重新构建敲除型基因共表达网络
- 比较两个网络的差异，识别受扰动影响的基因
1.2 scTenifoldKnk的工作原理
scTenifoldKnk的算法流程分为三个核心步骤：
第一步：构建野生型基因共表达网络
从单细胞表达矩阵中，通过计算基因间的相关性（或使用更复杂的网络推断方法），构建一个基因共表达网络。
第二步：构建虚拟敲除型基因共表达网络
从表达矩阵中移除目标基因，基于剩余基因重新构建共表达网络。
第三步：比较两个网络，识别差异调控基因
通过张量分解和流形对齐等技术，比较两个网络的差异，识别受目标基因调控的下游响应基因。
二、准备工作：安装与加载R包
2.1 安装scTenifoldKnk
# scTenifoldKnk安装
install.packages("scTenifoldKnk")
install.packages(c("openxlsx", "ggrepel"))
2.2 加载R包与数据
# 1. 加载R包
rm(list = ls())
options(stringsAsFactors = FALSE)
library(scTenifoldKnk) # 核心包
library(Seurat)
library(dplyr)
library(ggplot2)
library(openxlsx) # 导出Excel
library(ggrepel) # 标签防重叠
# 2. 设置工作目录
dir.create("cko", showWarnings = FALSE)
setwd("cko")
# 3. 读取数据
scRNA <- readRDS("../4_scRNA_celltype.rds")
DefaultAssay(scRNA) <- "RNA"
table(scRNA$cell_type)
三、数据准备
3.1 提取高变基因与目标基因
虚拟敲除分析的计算量较大，我们需要聚焦于高变基因（信息最丰富的基因），同时确保目标基因（待敲除的基因）被包含在分析中。
# 4. 提取高变基因与目标基因
set.seed(123)
# 获取高变基因
variable_genes <- VariableFeatures(scRNA)
# 定义目标基因（待敲除的基因）
target_gene <- "Ptn" # 可以替换为你感兴趣的基因
# 确保目标基因在数据中
if (!target_gene %in% rownames(scRNA)) {
stop(paste("目标基因", target_gene, "不在表达矩阵中！"))
}
# 合并高变基因与目标基因
com_genes <- union(variable_genes, target_gene)
# 子集化Seurat对象
scRNA_sub <- subset(scRNA, features = com_genes)
cat("高变基因数:", length(variable_genes), "\n")
cat("目标基因:", target_gene, "\n")
cat("分析基因总数:", nrow(scRNA_sub), "\n")
3.2 随机抽样（控制计算量）
如果细胞数较多（>10000），建议随机抽样以减少计算时间：
# 5. 随机抽样（可选）
# 设置抽样细胞数（可根据机器性能调整）
n_cells <- 1000
if (ncol(scRNA_sub) > n_cells) {
set.seed(456)
sampled_cells <- sample(1:ncol(scRNA_sub), n_cells, replace = FALSE)
scRNA_sub <- scRNA_sub[, sampled_cells]
cat("抽样后细胞数:", ncol(scRNA_sub), "\n")
} else {
cat("使用全部细胞:", ncol(scRNA_sub), "\n")
}
3.3 提取表达矩阵
scTenifoldKnk的输入是一个基因为行、细胞为列的counts矩阵：
# 6. 提取表达矩阵
expression_matrix <- LayerData(
scRNA_sub,
assay = "RNA",
layer = "counts"
)
# 确认目标基因存在
gene_presence <- target_gene %in% rownames(expression_matrix)
print(paste("目标基因", target_gene, "存在于表达矩阵中:", gene_presence))
# 查看矩阵维度
dim(expression_matrix)
四、核心分析：scTenifoldKnk虚拟敲除
4.1 运行虚拟敲除
这是整个流程的核心步骤——scTenifoldKnk函数会完成网络构建、网络比较和差异基因识别：
# 7. 运行虚拟敲除
set.seed(456)
# 创建目标基因文件夹
dir.create(target_gene, showWarnings = FALSE)
setwd(target_gene)
# 运行scTenifoldKnk
perturbation_results <- scTenifoldKnk(
countMatrix = expression_matrix,
gKO = target_gene,
qc = TRUE,
qc_mtThreshold = 0.1,
qc_minLSize = 1000,
nc_lambda = 0,
nc_nNet = 10, # 构建10个子网络
nc_nCells = 500, # 每个子网络的细胞数
nc_nComp = 3,
nc_scaleScores = TRUE,
nc_symmetric = FALSE,
nc_q = 0.9,
td_K = 3,
td_maxIter = 1000,
td_maxError = 1e-05,
td_nDecimal = 3,
ma_nDim = 2,
nCores = parallel::detectCores()
)
参数解读：
| 参数 | 说明 |
|---|---|
countMatrix |
基因×细胞的表达矩阵（counts） |
gKO |
要敲除的目标基因名称 |
nc_nNet |
构建的子网络数量（越大越稳定，计算越慢） |
nc_nCells |
每个子网络使用的细胞数 |
td_K |
张量分解的秩 |
ma_nDim |
流形对齐的维度 |
⏱️ 运行时间：取决于细胞数和网络数量。
4.2 结果提取与导出
# 8. 结果处理与导出
# 提取差异调控基因
differential_genes <- perturbation_results$diffRegulation
differential_genes$log_fold_change <- log2(differential_genes$FC)
differential_genes <- differential_genes[differential_genes$gene != target_gene, ]
differential_genes <- differential_genes[differential_genes$p.value < 0.05, ]
# 转换数值列类型
differential_genes[, 2:7] <- lapply(differential_genes[, 2:7], as.numeric)
write.xlsx(differential_genes, paste0(target_gene, "_result.xlsx"))
head(differential_genes)
cat("显著差异调控基因数:", nrow(differential_genes), "\n")
五、结果可视化
5.1 条形图：Top 10扰动最显著的基因
# 9. 条形图：Top 10扰动最显著的基因
top_n_genes <- 10
differential_genes_top <- differential_genes %>%
mutate(
abs_logfc = abs(log_fold_change),
# 添加绘图所需列
regulation = ifelse(log_fold_change > 0, "Up-regulated", "Down-regulated"),
label = sprintf("%.2f", log_fold_change),
label_hjust = ifelse(log_fold_change > 0, -0.2, 1.2)
) %>%
slice_max(abs_logfc, n = min(top_n_genes, nrow(.))) # 防止空数据
p1 <- ggplot(differential_genes_top,
aes(x = reorder(gene, log_fold_change),
y = log_fold_change,
fill = regulation)) +
geom_col(alpha = 0.8, width = 0.7) +
geom_text(aes(label = label, hjust = label_hjust),
size = 3.5, color = "black") +
scale_fill_manual(values = c("Up-regulated" = "#E74C3C",
"Down-regulated" = "#3498DB")) +
coord_flip() +
labs(title = paste("Top", top_n_genes, "Differentially Regulated Genes"),
subtitle = paste("KO gene:", target_gene),
x = "Gene",
y = "log2(Fold Change)",
fill = "Regulation") +
theme_classic() +
theme(
axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 12, face = "bold"),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray50"),
legend.position = "top",
legend.title = element_text(size = 10, face = "bold")
) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.15)))
ggsave("barplot_top_genes.pdf", p1, width = 10, height = 8)
解读：条形图展示受目标基因扰动影响最大的前10个基因。红色条表示敲除后上调的基因，蓝色条表示下调的基因。条越长，说明该基因受目标基因的调控影响越大。
5.2 火山图：扰动程度与显著性
# 10. 火山图
volcano_plot <- perturbation_results$diffRegulation %>%
mutate(
log2FC = log2(FC), # 预先计算，避免重复
significant = p.value < significant_threshold & gene != target_gene,
neg_log_pval = -log10(p.value),
perturbation_magnitude = case_when(
log2FC > 2 & significant ~ "High",
log2FC > 1 & significant ~ "Medium",
significant ~ "Low",
TRUE ~ "Not significant"
)
) %>%
mutate(
# 直接在管道中取显著基因中log2FC绝对值最大的10个
gene_label = if_else(
gene %in% (filter(., significant) %>% slice_max(abs(log2FC), n = 10) %>% pull(gene)),
gene,
NA_character_
)
)
p2 <- ggplot(volcano_plot, aes(x = log2FC, y = neg_log_pval)) +
# 散点层：按扰动程度着色
geom_point(aes(color = perturbation_magnitude),
alpha = 0.6, size = 2) +
# 颜色映射
scale_color_manual(
name = "Perturbation\nMagnitude",
values = c(
"High" = "#E74C3C", # 红色
"Medium" = "#F39C12", # 橙色
"Low" = "#3498DB", # 蓝色
"Not significant" = "gray80"
),
breaks = c("High", "Medium", "Low", "Not significant")
) +
# 显著性阈值线
geom_hline(yintercept = -log10(significant_threshold),
linetype = "dashed", color = "blue", alpha = 0.5) +
# log2FC阈值线
geom_vline(xintercept = c(-1, 1),
linetype = "dashed", color = "darkgreen", alpha = 0.4) +
# 基因标签
ggrepel::geom_text_repel(
data = subset(volcano_plot, !is.na(gene_label)),
aes(label = gene_label),
size = 3.5,
max.overlaps = 30,
box.padding = 0.5,
point.padding = 0.3,
segment.color = "gray50",
segment.alpha = 0.5,
fontface = "bold",
color = "#2C3E50"
) +
# 坐标轴和标题
labs(
title = paste("Virtual Gene Perturbation:", target_gene),
x = expression(log[2] ~ "(Fold Change)"),
y = expression(-log[10] ~ "(p-value)")
) +
# 主题
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 9),
axis.text = element_text(size = 11, color = "black"),
axis.title = element_text(size = 12, face = "bold"),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray50"),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", linewidth = 0.8)
)
ggsave("volcano.pdf", p2,width = 8, height = 6)
解读：火山图同时展示扰动程度（x轴） 和显著性（y轴）。右上角和左上角的点是最值得关注的——它们既有大的扰动效应，又有高的统计显著性。
5.3 MA图：扰动程度与表达量的关系
# 11. MA图
# 获取要标记的Top基因
genes_to_label <- differential_genes %>%
slice_max(log_fold_change, n = 5) %>%
pull(gene)
# 构建MA数据
MA_plot <- data.frame(
gene = rownames(expression_matrix),
mean_expression = log1p(rowMeans(expression_matrix))
) %>%
left_join(
perturbation_results$diffRegulation %>%
transmute(gene, perturbation_magnitude = log2(FC)),
by = "gene"
) %>%
mutate(
is_perturbed = gene %in% differential_genes$gene,
gene_label = if_else(gene %in% genes_to_label, gene, NA_character_),
perturbation_category = factor(
case_when(
is_perturbed & perturbation_magnitude > 2 ~ "High",
is_perturbed & perturbation_magnitude > 1 ~ "Medium",
is_perturbed ~ "Low",
TRUE ~ "Not perturbed"
),
levels = c("High", "Medium", "Low", "Not perturbed")
)
)
MA_plot <- MA_plot %>%
mutate(
significance_label = if_else(is_perturbed, "Significant", "Not significant")
)
p3 <- ggplot(MA_plot, aes(x = mean_expression, y = perturbation_magnitude)) +
# 散点层
geom_point(
aes(color = significance_label, size = mean_expression),
alpha = 0.5, shape = 16
) +
# 颜色映射
scale_color_manual(
name = "Significance",
values = c("Significant" = "#E74C3C", "Not significant" = "gray80")
) +
# 点大小映射
scale_size_continuous(
name = "Mean Expression",
range = c(0.5, 4),
guide = guide_legend(override.aes = list(alpha = 0.8))
) +
# LOESS趋势线
geom_smooth(
method = "loess",
color = "#2C3E50",
se = TRUE,
alpha = 0.2,
linewidth = 0.6,
fill = "#3498DB"
) +
# 基因标签
ggrepel::geom_text_repel(
data = subset(MA_plot, !is.na(gene_label)),
aes(label = gene_label),
size = 3.5,
max.overlaps = 30,
box.padding = 0.5,
point.padding = 0.3,
color = "#E74C3C",
fontface = "bold",
segment.color = "gray50",
segment.alpha = 0.4
) +
# 坐标轴和标题
labs(
title = paste("Perturbation Magnitude vs Mean Expression"),
subtitle = paste("KO gene:", target_gene),
x = expression(log ~ "(Mean Expression + 1)"),
y = expression(log[2] ~ "(Perturbation Magnitude)")
) +
# 主题
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 9),
axis.text = element_text(size = 11, color = "black"),
axis.title = element_text(size = 12, face = "bold"),
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray50"),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", linewidth = 0.8)
)
ggsave("MA_plot.pdf", p3, width = 8, height = 6)
解读：MA图展示基因的平均表达水平（x轴） 与扰动程度（y轴） 的关系。可以观察目标基因的扰动效应是否偏向于高表达或低表达的基因。
5.4 点图：Top扰动基因可视化
# 12. 点图
top_dotplot <- differential_genes %>%
slice_max(log_fold_change, n = 10)
p4 <- top_dotplot %>%
mutate(neg_log_pval = -log10(p.value)) %>%
ggplot(aes(reorder(gene, log_fold_change), log_fold_change,
size = neg_log_pval, color = log_fold_change)) +
geom_point(alpha = 0.8) +
scale_color_gradient2(low = "#3498DB", mid = "#F5F5F5", high = "#E74C3C") +
scale_size_continuous(range = c(3, 8)) +
coord_flip() +
labs(x = NULL, y = expression(log[2]~FC),
color = expression(log[2]~FC), size = expression(-log[10]~p)) +
theme_classic() +
theme(panel.grid.major.y = element_line(color = "gray90"))
ggsave("dotplot.pdf", p4, width = 8, height = 6)
解读：点图整合了扰动程度（x轴位置）、显著性（点的大小）和方向（颜色） 三个维度的信息。红色大点是最值得关注的候选基因。
六、常见问题与注意事项
6.1 scTenifoldKnk运行很慢怎么办？
-
减少
nc_nNet参数：从10减少到5（会降低稳定性） -
减少
nc_nCells参数：从500减少到300（每个子网络用更少的细胞） - 减少输入基因数：只使用高变基因（已做）
- 减少细胞数：抽样5000个细胞而非10000个
6.2 目标基因不在高变基因中怎么办？
在 union(variable_genes, target_gene) 中已经强制将目标基因加入分析列表，即使它不在高变基因中也会被保留。
6.3 结果中显著基因太少怎么办？
- 降低
significant_threshold（如从0.05放宽到0.1） - 减少
nc_nNet参数（网络数量减少，但单个网络的可靠性可能提高） - 检查目标基因在数据中是否有足够的表达量
6.4 虚拟敲除 vs 实验敲除
| 对比项 | 虚拟敲除 | 实验敲除 |
|---|---|---|
| 成本 | 极低（计算时间） | 高（实验成本） |
| 速度 | 快（数小时） | 慢（数周-数月） |
| 验证方式 | 计算预测 | 实验验证 |
| 适用范围 | 转录组层面 | 表型+转录组 |
| 局限性 | 仅预测调控网络变化 | 可观察完整表型 |
虚拟基因敲除是连接计算预测与实验验证的桥梁。在单细胞研究中，它帮助你从“这个基因在哪些细胞中表达”进一步追问“这个基因的功能是什么？敲掉它会有什么后果？”，从而推动研究从描述性分析向机制性探索迈进。
这里是两栖生物手册，中科院生物医学博士，持续记录生物医学实验和生物信息学笔记，干湿结合两不误～正努力成为最贴心负责的生信数据分析者和实验技术分享者，欢迎大家关注～
