scTenifoldKnk 虚拟基因敲除（In-Silico Knockout）
When to Use
- 想预测某个基因的功能，但没有敲除/干扰实验数据时——用野生型（WT）单细胞数据在计算上模拟敲除，看哪些基因被扰动。
- 已有一个野生型/对照条件的单细胞数据（不需要 KO 样本），想评估候选基因（如差异表达得到的 hub 基因、转录因子）的调控影响。
- 需要对多个候选基因做批量筛选（
transcriptomeWide模式，一次建网、逐个敲除）。 - 已有真实的 KO 实验数据，想用虚拟敲除做交叉验证或排序。
不适用：只有归一化/批次校正后的表达值（必须用原始 counts）；细胞数极少（< 数百，网络不稳）；想比较两个条件间的差异表达（那是 DE 分析的工作）。
Overview
scTenifoldKnk 建在 scTenifoldNet 之上，三步：
- 构建 WT scGRN：主成分回归为每个基因拟合其调控关系，重复
nc_nNet次（每次随机抽nc_nCells个细胞），用张量分解（CP/ALS，秩td_K）去噪，得到权重平均的调控网络。 - 虚拟敲除：把目标基因在邻接矩阵中的出度边全部置零（假设敲除消除了该基因对其靶基因的调控），得到 KO 网络。
- 流形对齐 + 差异调控：对 WT/KO 两个网络做非线性流形对齐，计算每个基因在两套坐标间的欧氏距离；距离经 Box-Cox 变换标准化后按卡方分布（1 自由度）算 p 值，BH 校正得到 p.adj。距离相对期望值放大最多的基因即虚拟敲除扰动基因。
注意：logFC = log2(FC) 中的 FC 是流形对齐距离相对期望值的倍数，不是表达量的差异倍数。
依赖
脚本运行时需要 R 环境中安装：scTenifoldKnk, scTenifoldNet, Seurat, Matrix, matrixStats, ggplot2, dplyr。
> scTenifoldKnk / scTenifoldNet 均不在 CRAN/Bioconductor，需从 GitHub 安装（cailab-tamu/scTenifoldKnk、cailab-tamu/scTenifoldNet）。见 abio-r-package-install skill。
Steps
可复用脚本：scripts/sctenifoldknk.R.txt（纯命令行参数，无硬编码路径）。
- 准备输入：确认是原始 counts（整数、未归一化），基因为行名、细胞为列名。
- 子集到高变基因（强烈建议，
--hvg TRUE默认）：基因数直接决定耗时与内存（张量是nGenes × nGenes × nNet）。可用--n-genes进一步限制到方差最大的前 N 个。 - 确认目标基因在矩阵内：
gKO必须在 rownames 中，否则报错。若目标基因不在 HVG 里，需--hvg FALSE或重新FindVariableFeatures把它纳入。 - 运行：
scTenifoldKnk(countMatrix, gKO, ...)（QC → 建网 → 张量分解 → 敲除 → 流形对齐 → 差异调控）。 - 整理结果：取
res$diffRegulation（列：gene, distance, Z, FC, p.value, p.adj），算logFC = log2(FC)，剔除被敲除基因自身，按 p 值排序。 - 筛选与可视化：p.value < 阈值（默认 0.05）为扰动基因；输出 top N 横条图 + 火山图；同时看 p.adj 作为保守参考。
- 保存：完整结果 RDS（含 WT/KO 网络与流形对齐坐标）+ CSV + markdown 报告。
Parameters to Customize
| 参数 | 说明 | 默认 |
|---|---|---|
--input |
输入 RDS（Seurat 对象，或 基因×细胞 counts 矩阵） | 必填 |
--output |
输出根目录（自动创建 Img/Stat/Rdata） | . |
--gene |
要敲除的基因符号（gKO） | 必填 |
--assay |
Seurat assay（对象里没有该 assay 时回落到 DefaultAssay） | RNA |
--hvg |
是否只用高变基因 | TRUE |
--n-genes |
再按方差取前 N 个基因（0 = 不限制） | 0 |
--min-lib |
QC：最小文库大小（qc_minLibSize） |
1000 |
--max-mt |
QC：线粒体比例上限（qc_maxMTratio） |
0.1 |
--min-pct |
QC：基因至少在多少比例的细胞中表达（qc_minPCT） |
0.05 |
--remove-outliers |
QC：剔除文库大小离群细胞（qc_removeOutlierCells） |
TRUE |
--n-net |
构建的网络数（nc_nNet） |
10 |
--n-cells |
每次建网抽样的细胞数（nc_nCells） |
500 |
--n-comp |
主成分回归的 PC 数（nc_nComp，须 >2 且 < 基因数） |
3 |
--k |
张量分解秩（td_K） |
3 |
--lambda |
网络稀疏化强度（nc_lambda，0–1；0 表示不稀疏化） |
0 |
--n-cores |
并行核数（硬性 ≤4） | 4 |
--p-cutoff |
扰动基因的 p 值阈值 | 0.05 |
--top |
报告中展示的 top 基因数 | 10 |
--empirical-null |
用 Efron 经验零分布代替卡方（需 locfdr 包） | FALSE |
--seed |
随机种子（建网抽样） | 42 |
Evidence Criteria / Expected Output
Stat/diffRegulation_all.csv：全部基因的 distance / Z / FC / logFC / p.value / p.adjStat/diffRegulation_p<cutoff>.csv：通过阈值的扰动基因Img/knk_top_genes_barplot.pdf|png：top 扰动基因的 logFC 横条图Img/knk_volcano.pdf|png：logFC vs -log10(p)Rdata/scTenifoldKnk_result_<gene>.rds：完整结果（tensorNetworksWT/KO 网络、manifoldAlignment、diffRegulation）Stat/knk_report.md：参数、规模、显著基因数与 Top 表- 判读：p.value < 0.05 的基因视为被扰动。p.value 基于卡方分布、偏激进，务必同时看
p.adj（BH）——若 p 显著但 p.adj 全不显著，说明扰动信号弱，应调大--n-net/--n-cells或谨慎下结论。
Pitfalls / Notes
- ⚠️ 参数名陷阱（原脚本常见错误）：正确写法是
qc_maxMTratio（不是qc_mtThreshold）和qc_minLibSize（不是qc_minLSize）。scTenifoldKnk()签名里没有...，写错名字会直接报unused argument而中断，不会静默忽略。 - ⚠️ 必须显式限制线程数，且旋钮不是
nCores：detectCores()与blas_get_num_procs()在容器里都返回宿主机核数（本机 64）。而 scTenifoldNet 的计算量主要走pcNetCoreRcpp编译后端（OpenMP），它完全不受nCores控制；nCores只在 Rcpp 后端失败时的furrr::future_map兜底路径生效。实测不限流单进程直接吃到 5766% CPU（≈57 核）。正确做法是设置OMP_NUM_THREADS/OPENBLAS_NUM_THREADS/MKL_NUM_THREADS（必须在加载相关包之前）+RhpcBLASctl::blas_set_num_threads()/omp_set_num_threads()，并把传给scTenifoldKnk()的nCores固定为 1（防止兜底路径再乘一层）。脚本已按此实现，--n-cores即 OpenMP 线程预算（硬上限 4）。实测：修复前单进程 5766%（≈57 核）；修复后单进程 CPU 峰值 403%（≈4 核，符合配额），进程树仅 1 个进程。注意 OpenMP/OpenBLAS 的线程池仍会按宿主核数建（/proc/<pid>/status里可见 64 线程），但活跃并发被限制在 4，空闲线程不耗 CPU。 - 输入必须是原始 counts：不要传归一化/log 化的表达值。scTenifoldKnk 内部自己做 QC 和 CPM 归一化。
- 规模决定可行性：耗时与内存约按
nGenes² × nNet增长。实测基准（631 基因 × 1253 细胞 QC 后，4 核，nc_nCells=500）：--n-net 3→ 3.6 分钟；--n-net 10→ 11.9 分钟（网络构建占 ~95%）。两版结果的 top 扰动基因完全一致、log10(p) 相关系数 r=1.000，说明该数据上结果对nNet稳健；但发表级分析仍建议 nNet≥10。5000 基因时张量约 2GB、耗时显著上升。脚本会打印预估张量内存并在 >20GB 时告警（上限 50GB）。 - QC 可能大幅削减基因数：
qc_minPCT按全部细胞计算表达比例，只在稀有细胞类型中表达的基因会被剔除。实测 1500 个 HVG + 7546 细胞，经 minLibSize=1000 / MT<0.1 / minPCT=0.05 后只剩 631 基因 × 1253 细胞。若目标基因或关键调控因子被 QC 掉，需放宽--min-pct/--min-lib，或先用--hvg FALSE检查。 nc_nCells要 ≤ 实际细胞数（脚本会自动截断并提示）。抽样细胞越多网络越稳，但单次回归更慢。nc_nComp必须 >2 且 < 基因数；td_K为 CP 张量分解秩，默认 3。- 可复现性：建网阶段有随机抽样（每次抽
nc_nCells个细胞），固定--seed才可复现。 - 线粒体基因识别：scQC 按基因名以
MT-开头（大小写不敏感）计算线粒体比例；基因名不是这种命名规范时该 QC 实际不生效。 - 结果是计算预测，不是实验证据：虚拟敲除只能给出候选基因排序和假设，需后续实验（qPCR / 敲低 / 敲除）验证。报告中不要把它写成"证实"。
- 单基因模式 vs 全转录组模式：脚本用单基因模式（
transcriptomeWide=FALSE，结果返回 WT/KO 两套网络）。若要批量敲多个基因，改为transcriptomeWide=TRUE会一次建网、逐个扰动，返回perturbationDistances矩阵（效率高得多，但不再返回 KO 网络）。 - 不要用结果反推"表达量变化"：
FC是流形对齐距离比，与表达 fold change 不同尺度，不能直接与其他工具的 logFC 比较。
> Packaging note. Scripts and tables in this skill carry a .txt suffix on their
> original extension (e.g. .R.txt, .tsv.txt) because the SkillHub upload allowlist
> does not accept .R or .tsv. This is cosmetic: Rscript is extension-agnostic, so
> Rscript scripts/<name>.R.txt ... behaves exactly like .R, and tab-separated data
> parses identically from .tsv.txt. Keep the suffix when repackaging.
