最新文章

20260912 111444 Seurat V5 数据整合分析教程

2026年09月12日

Seurat v5 数据整合分析教程 ¶

基于 Seurat v5.5.1 官方教程，系统学习多批次单细胞数据的整合分析。
使用 PBMC 3k 数据模拟批次效应，完整演示 CCA / RPCA 整合、保守 Marker、跨批次差异表达等核心流程。

目录 ¶

1. 背景与原理
2. 模拟批次效应
3. 未整合分析（展示批次效应）
4. CCA 整合
5. 整合前后对比
6. 细胞类型注释
7. 保守 Marker 基因（FindConservedMarkers）
8. 跨批次差异表达（FindMarkers）
9. RPCA 整合与 CCA 对比
10. 总结与核心知识点
附录: 完整脚本

1. 背景与原理 ¶

1.1 为什么需要数据整合？¶

在单细胞实验中，不同批次的样本由于技术差异（捕获效率、测序深度、试剂批次等）会产生批次效应。如果不进行整合，批次效应会导致：

相同细胞类型因批次不同而被分到不同聚类
批次特异的技术噪声干扰生物学信号
跨批次比较不可靠

1.2 Seurat v5 整合核心特性 ¶

Seurat v5 引入了 layers 机制，可以在一个 Seurat 对象中管理多个批次的数据，无需创建多个对象：

操作	函数	说明
按批次分割	`split(pbmc[["RNA"]], f = pbmc$batch)`	将 counts 分为多个 layer
整合	`IntegrateLayers(method = CCAIntegration)`	一行代码完成整合
合并	`JoinLayers(pbmc[["RNA"]])`	整合后合并回单 layer

1.3 两种整合方法对比 ¶

方法	函数	特点	适用场景
CCA	`CCAIntegration`	经典方法，整合力度大	批次差异大、需要彻底去除批次效应
RPCA	`RPCAIntegration`	更快更保守	批次差异小、希望保留更多生物学变异

2. 模拟批次效应 ¶

由于官方 ifnb 数据集（394MB）下载过慢，本教程使用已有的 PBMC 3k 数据，通过人为分割模拟两个「批次」。

2.1 数据准备 ¶

# 读取原始 10X 数据
pbmc_data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc_data, project = "pbmc3k",
                            min.cells = 3, min.features = 200)

# QC 过滤
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
# 过滤后细胞数: 2638

2.2 随机分批 + 添加技术偏差 ¶

# 随机分为两个批次
set.seed(42)
pbmc$batch <- sample(c("batch1", "batch2"), size = ncol(pbmc), replace = TRUE)
# batch1: 1348 个细胞 | batch2: 1290 个细胞

# 为 batch2 添加技术偏差：
# 1) 乘法噪声（0.7-1.3 的基因特异噪声因子，模拟不同捕获效率）
# 2) Dropout 500 个基因的 30% 细胞（模拟检测率差异）

2.3 分割 Layers（Seurat v5 核心）¶

# 按 batch 分割 RNA assay 的 counts layer
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)
# 分割后 pbmc 对象包含 counts.batch1 和 counts.batch2 两个 layer

3. 未整合分析（展示批次效应）¶

对每个 layer 独立进行标准流程（归一化、找高变基因、PCA），展示未整合时的批次效应。

pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)
pbmc <- FindNeighbors(pbmc, dims = 1:30, reduction = "pca")
pbmc <- FindClusters(pbmc, resolution = 0.8, cluster.name = "unintegrated_clusters")
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

3.1 未整合 UMAP - 按批次着色 ¶

可以清楚看到，两个批次在 UMAP 空间中明显分离，这就是批次效应。

未整合-按批次

3.3 未整合 - 按批次分割 ¶

用 split.by 参数分面展示，更直观地看到两个批次的分布差异。

未整合-按批次分割

4. CCA 整合 ¶

4.1 一行代码完成整合 ¶

Seurat v5 的 IntegrateLayers 函数封装了整个整合流程：

pbmc <- IntegrateLayers(
    object = pbmc,
    method = CCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.cca",
    verbose = FALSE
)

CCA 整合原理：
1. 在两个批次之间寻找典型相关分析（CCA）的「锚点」对
2. 基于锚点对计算批次校正向量
3. 将校正后的数据投影到新的降维空间（integrated.cca）

4.2 合并 Layers 并重新聚类 ¶

# 整合后必须合并 layers
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# 基于整合后的降维结果聚类
pbmc <- FindNeighbors(pbmc, reduction = "integrated.cca", dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5, cluster.name = "cca_clusters")
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "integrated.cca", reduction.name = "umap.cca")

4.3 整合后 UMAP - 按批次着色 ¶

整合后，两个批次在 UMAP 空间中充分混合，批次效应已消除。

CCA整合后-按批次

5. 整合前后对比 ¶

5.1 按批次着色对比 ¶

左：未整合（批次明显分离）| 右：CCA 整合后（批次充分混合）

整合前后-按批次

5.2 按聚类着色对比 ¶

左：未整合聚类（受批次影响）| 右：整合后聚类（反映真实细胞类型）

整合前后-按聚类

6. 细胞类型注释 ¶

6.1 寻找各聚类 Marker 基因 ¶

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top_markers <- pbmc.markers %>%
    group_by(cluster) %>%
    slice_max(n = 5, order_by = avg_log2FC)

各聚类 top 5 marker（示例）：

聚类	Top Markers	细胞类型
0	CD3D, CCR7, CD3E, LDHB, CD7	Naive CD4 T
1	IL32, CD3D, CD3E, IL7R, CD2	Memory CD4 T / T 细胞
2	LYZ, S100A9, S100A8	CD14+ Mono
3	MS4A1, CD79A, CD79B	B 细胞
4	CD8A, GZMK, CCL5	CD8 T
5	FCGR3A, MS4A7	FCGR3A+ Mono
6	GNLY, NKG7, KLRD1	NK
7	FCER1A, CST3, LST1	DC
8	PPBP, PF4	Platelet

6.3 Marker FeaturePlot ¶

经典 marker 基因在 UMAP 空间的表达分布：

FeaturePlot(pbmc, features = c("IL7R", "CD14", "MS4A1", "CD8A", "GNLY", "FCGR3A"),
            reduction = "umap.cca", ncol = 3)

IL7R: CD4 T 细胞
CD14: 单核细胞
MS4A1: B 细胞
CD8A: CD8 T 细胞
GNLY: NK 细胞
FCGR3A: FCGR3A+ 单核细胞

Marker FeaturePlot

7. 保守 Marker 基因（FindConservedMarkers）¶

7.1 原理 ¶

FindConservedMarkers 在两个批次中分别寻找某聚类的 marker 基因，然后用 meta-analysis 合并 p 值，找到不受批次效应影响的稳定 marker。

7.2 核心代码 ¶

# 注意：FindConservedMarkers 需要 layers 处于合并状态
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

Idents(pbmc) <- "cca_clusters"
# 选择最大的聚类进行分析
conserved_markers <- FindConservedMarkers(
    pbmc,
    ident.1 = target_cluster,  # 聚类1
    grouping.var = "batch",     # 按批次分组
    verbose = FALSE
)

7.3 结果解读 ¶

聚类1（T 细胞）的保守 marker 基因（前 15 个）：

基因	batch2 avg_log2FC	batch1 avg_log2FC	说明
IL32	1.457	1.153	T 细胞活化标志
CD3D	1.056	1.067	T 细胞受体组件
CD3E	1.052	1.052	T 细胞受体组件
IL7R	1.260	1.711	CD4 T 细胞标志
CD2	1.848	1.389	T 细胞粘附分子
LDHB	0.970	1.038	代谢酶（泛表达）
LTB	1.235	1.438	淋巴毒素
LCK	0.957	0.682	T 细胞信号转导
AQP3	1.826	2.623	水通道蛋白（T 细胞活化）
PRDX2	1.185	0.871	抗氧化蛋白

这些基因在两个批次中都显著高表达，说明它们是稳定的生物学 marker，不受批次效应影响。

基因	avg_log2FC	pct.1	pct.2	说明
B2M	-0.722	1.0	1.0	MHC-I 组件（技术偏差）
RPL15	-0.634	1.0	1.0	核糖体蛋白（技术偏差）
RPL13	-0.558	1.0	1.0	核糖体蛋白（技术偏差）
TMSB4X	-0.419	1.0	1.0	肌动蛋白结合（技术偏差）
RPL19	-0.487	1.0	1.0	核糖体蛋白（技术偏差）
RPL3	-0.468	1.0	1.0	核糖体蛋白（技术偏差）
RPS15A	-0.510	0.993	1.0	核糖体蛋白（技术偏差）
RPS19	-0.597	1.0	1.0	核糖体蛋白（技术偏差）

9. RPCA 整合与 CCA 对比 ¶

9.1 RPCA 整合 ¶

RPCA（Reciprocal PCA）是一种更快速、更保守的整合方法：

# RPCA 需要 layers 处于分割状态
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)

pbmc <- IntegrateLayers(
    object = pbmc,
    method = RPCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.rpca",
    verbose = FALSE
)

# 合并 layers
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

9.4 CCA vs RPCA 选择建议 ¶

指标	CCA	RPCA
整合力度	强（可能过度校正）	适中（保留更多变异）
计算速度	较慢	较快
适用场景	批次差异大	批次差异小
生物学信号保留	可能丢失部分	保留更好
推荐	常规首选	大数据集/小批次差异

10. 总结与核心知识点 ¶

10.1 核心流程图 ¶

原始数据 → QC 过滤 → 模拟/识别批次
                           ↓
                    split(pbmc[["RNA"]], f = pbmc$batch)
                           ↓
              ┌─────────────┴──────────────┐
              ↓                           ↓
        未整合分析                    CCA/RPCA 整合
        (展示批次效应)              IntegrateLayers()
              ↓                           ↓
        UMAP（批次分离）           JoinLayers() + 聚类
                                      ↓
                              UMAP（批次混合）
                                      ↓
                    ┌─────────────────┴──────────────────┐
                    ↓                                    ↓
           FindConservedMarkers                    FindMarkers
           (跨批次保守 marker)                   (跨批次差异基因)

10.2 关键注意事项 ¶

split vs JoinLayers 切换：
- IntegrateLayers 需要 layers 分割状态
- FindConservedMarkers 需要 layers 合并状态
- 脚本中需要反复切换
依赖包：
- FindConservedMarkers 需要 metap 包
- metap 依赖 multtest（Bioconductor 包）
- 安装：BiocManager::install("multtest") → install.packages("metap")
整合方法选择：
- 默认使用 CCA（整合力度大）
- 大数据集或小批次差异时优先使用 RPCA

10.3 输出文件清单 ¶

编号	文件名	内容
01	01_unintegrated_batch.png	未整合 UMAP - 按批次
02	02_unintegrated_cluster.png	未整合 UMAP - 按聚类
03	03_unintegrated_split.png	未整合 - 按批次分割
04	04_integrated_batch.png	CCA 整合后 - 按批次
05	05_integrated_cluster.png	CCA 整合后 - 按聚类
06	06_integrated_split.png	CCA 整合后 - 按批次分割
07	07_compare_batch.png	整合前后对比 - 按批次
08	08_compare_cluster.png	整合前后对比 - 按聚类
09	09_annotated_umap.png	注释后 UMAP
10	10_marker_featureplot.png	Marker FeaturePlot
11	11_conserved_markers_violin.png	保守 Marker 小提琴图
12	12_conserved_markers_dot.png	保守 Marker DotPlot
13	13_batch_de_violin.png	跨批次差异基因小提琴图
14	14_rpca_batch.png	RPCA 整合后 - 按批次
15	15_cca_vs_rpca.png	CCA vs RPCA 对比

附录: 完整脚本 ¶

integration_vignette.R ¶

###############################################################################
# Seurat v5 阶段3：数据整合分析（Integration）
# 基于官方教程:
#   - https://satijalab.org/seurat/articles/integration_introduction
#   - https://satijalab.org/seurat/articles/seurat5_integration
#
# 由于 ifnb 数据集（394MB）从 Seurat 官方服务器下载过慢，
# 本脚本使用已有的 PBMC 3k 数据，通过人为分割模拟两个"批次"，
# 完整演示 Seurat v5 整合流程的所有核心步骤：
#   1. 模拟批次效应
#   2. 未整合分析（展示批次效应）
#   3. CCA 整合
#   4. 整合前后对比
#   5. 保守 Marker 基因（FindConservedMarkers）
#   6. 跨条件差异表达（FindMarkers）
#   7. 整合方法对比（CCA vs RPCA）
###############################################################################

library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)

# === 路径配置 ===
base_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial"
data_dir <- file.path(base_dir, "data/filtered_gene_bc_matrices/hg19")
output_dir <- file.path(base_dir, "output")
integ_dir <- file.path(output_dir, "integration")
dir.create(integ_dir, showWarnings = FALSE, recursive = TRUE)

# ============================================================================
# 1. 加载数据 & 模拟批次效应
# ============================================================================
cat("\n========== 1. 加载数据 & 模拟批次 ==========\n")

# 读取原始 10X 数据
pbmc_data <- Read10X(data.dir = data_dir)

# 创建 Seurat 对象
pbmc <- CreateSeuratObject(counts = pbmc_data, project = "pbmc3k",
                            min.cells = 3, min.features = 200)

# 计算线粒体百分比
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# QC 过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
cat("过滤后细胞数:", ncol(pbmc), "\n")

# --- 模拟批次效应 ---
# 将细胞随机分为两个"批次"
set.seed(42)
cell_ids <- colnames(pbmc)
batch_labels <- sample(c("batch1", "batch2"), size = length(cell_ids), replace = TRUE)
names(batch_labels) <- cell_ids
pbmc$batch <- batch_labels

# 为 batch2 人为添加技术偏差：
# 1) 对 batch2 的部分基因表达量添加随机噪声（模拟批次效应）
# 2) 降低 batch2 中某些基因的检测率
cat("正在模拟批次效应...\n")

# 获取 counts 矩阵
counts_mat <- GetAssayData(pbmc, layer = "counts")

# 找到 batch2 的细胞
batch2_cells <- which(pbmc$batch == "batch2")
batch1_cells <- which(pbmc$batch == "batch1")
cat("batch1:", length(batch1_cells), "个细胞\n")
cat("batch2:", length(batch2_cells), "个细胞\n")

# 对 batch2 添加乘法噪声（模拟不同的捕获效率）
set.seed(123)
noise_factor <- 0.7 + 0.6 * runif(nrow(counts_mat))  # 0.7-1.3 的基因特异噪声
noise_matrix <- matrix(noise_factor, nrow = nrow(counts_mat), ncol = length(batch2_cells))
# 对 batch2 的 counts 乘以噪声因子（向下取整）
counts_batch2 <- ceiling(as.matrix(counts_mat[, batch2_cells]) * noise_matrix)
counts_mat[, batch2_cells] <- as(counts_batch2, "dgCMatrix")

# 为 batch2 随机 dropout 部分基因（模拟检测率差异）
set.seed(456)
dropout_genes <- sample(rownames(pbmc), size = 500)  # 随机选 500 个基因
for (gene in dropout_genes) {
    dropout_cells <- sample(batch2_cells, size = floor(length(batch2_cells) * 0.3))
    counts_mat[gene, dropout_cells] <- 0
}

# 写回修改后的 counts
pbmc[["RNA"]]$counts <- counts_mat
cat("批次效应模拟完成\n")

# ============================================================================
# 2. 按批次分割 Layers（Seurat v5 核心特性）
# ============================================================================
cat("\n========== 2. 分割 Layers ==========\n")

# Seurat v5 核心特性：在一个对象中通过 layers 管理多个批次
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)
cat("分割后 layers:\n")
print(pbmc)

# ============================================================================
# 3. 未整合分析（展示批次效应）
# ============================================================================
cat("\n========== 3. 未整合分析 ==========\n")

# 标准流程（对每个 layer 独立归一化和找高变基因）
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)

# 聚类
pbmc <- FindNeighbors(pbmc, dims = 1:30, reduction = "pca")
pbmc <- FindClusters(pbmc, resolution = 0.8, cluster.name = "unintegrated_clusters")

# UMAP（未整合）
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

# 3.1 未整合 UMAP - 按批次着色
p_unintegrated_batch <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "batch") +
    ggtitle("未整合 - 按批次着色")
ggsave(file.path(integ_dir, "01_unintegrated_batch.png"), plot = p_unintegrated_batch, width = 8, height = 7, dpi = 150)
cat("已保存: 01_unintegrated_batch.png\n")

# 3.2 未整合 UMAP - 按聚类着色
p_unintegrated_cluster <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "unintegrated_clusters") +
    ggtitle("未整合 - 按聚类着色")
ggsave(file.path(integ_dir, "02_unintegrated_cluster.png"), plot = p_unintegrated_cluster, width = 8, height = 7, dpi = 150)
cat("已保存: 02_unintegrated_cluster.png\n")

# 3.3 未整合 - split.by 批次
p_unintegrated_split <- DimPlot(pbmc, reduction = "umap.unintegrated", split.by = "batch") +
    ggtitle("未整合 - 按批次分割")
ggsave(file.path(integ_dir, "03_unintegrated_split.png"), plot = p_unintegrated_split, width = 14, height = 6, dpi = 150)
cat("已保存: 03_unintegrated_split.png\n")

# ============================================================================
# 4. CCA 整合
# ============================================================================
cat("\n========== 4. CCA 整合 ==========\n")
cat("正在运行 CCA 整合（可能需要几分钟）...\n")

# IntegrateLayers 一行代码完成整合
pbmc <- IntegrateLayers(
    object = pbmc,
    method = CCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.cca",
    verbose = FALSE
)
cat("CCA 整合完成\n")

# 重新合并 layers（整合后必须执行）
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])
cat("Layers 已重新合并\n")

# 基于整合后的降维结果聚类
pbmc <- FindNeighbors(pbmc, reduction = "integrated.cca", dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5, cluster.name = "cca_clusters")

# 整合后 UMAP
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "integrated.cca", reduction.name = "umap.cca")

# 4.1 整合后 UMAP - 按批次着色
p_integrated_batch <- DimPlot(pbmc, reduction = "umap.cca", group.by = "batch") +
    ggtitle("CCA 整合后 - 按批次着色")
ggsave(file.path(integ_dir, "04_integrated_batch.png"), plot = p_integrated_batch, width = 8, height = 7, dpi = 150)
cat("已保存: 04_integrated_batch.png\n")

# 4.2 整合后 UMAP - 按聚类着色
p_integrated_cluster <- DimPlot(pbmc, reduction = "umap.cca", group.by = "cca_clusters", label = TRUE) +
    ggtitle("CCA 整合后 - 按聚类着色")
ggsave(file.path(integ_dir, "05_integrated_cluster.png"), plot = p_integrated_cluster, width = 8, height = 7, dpi = 150)
cat("已保存: 05_integrated_cluster.png\n")

# 4.3 整合后 - split.by 批次
p_integrated_split <- DimPlot(pbmc, reduction = "umap.cca", split.by = "batch") +
    ggtitle("CCA 整合后 - 按批次分割")
ggsave(file.path(integ_dir, "06_integrated_split.png"), plot = p_integrated_split, width = 14, height = 6, dpi = 150)
cat("已保存: 06_integrated_split.png\n")

# ============================================================================
# 5. 整合前后对比图
# ============================================================================
cat("\n========== 5. 整合前后对比 ==========\n")

# 5.1 并排对比 - 按批次着色
p_compare_batch <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "batch") +
    ggtitle("未整合") +
DimPlot(pbmc, reduction = "umap.cca", group.by = "batch") +
    ggtitle("CCA 整合后")
ggsave(file.path(integ_dir, "07_compare_batch.png"), plot = p_compare_batch, width = 16, height = 6, dpi = 150)
cat("已保存: 07_compare_batch.png\n")

# 5.2 并排对比 - 按聚类着色
p_compare_cluster <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "unintegrated_clusters") +
    ggtitle("未整合聚类") +
DimPlot(pbmc, reduction = "umap.cca", group.by = "cca_clusters", label = TRUE) +
    ggtitle("整合后聚类")
ggsave(file.path(integ_dir, "08_compare_cluster.png"), plot = p_compare_cluster, width = 16, height = 6, dpi = 150)
cat("已保存: 08_compare_cluster.png\n")

# ============================================================================
# 6. 细胞类型注释
# ============================================================================
cat("\n========== 6. 细胞类型注释 ==========\n")

# 使用经典 marker 进行注释（与阶段1一致）
# 先寻找每个聚类的 marker 基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top_markers <- pbmc.markers %>%
    group_by(cluster) %>%
    slice_max(n = 5, order_by = avg_log2FC)
cat("各聚类 top 5 marker:\n")
print(top_markers)

# 保存 marker 表
write.csv(pbmc.markers, file.path(integ_dir, "integration_markers.csv"), row.names = FALSE)
cat("已保存: integration_markers.csv\n")

# 基于阶段1的知识进行注释
# 查看各聚类的 marker 来确定细胞类型
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B",
                      "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
# 根据实际聚类数调整
n_clusters <- length(unique(pbmc$cca_clusters))
cat("聚类数量:", n_clusters, "\n")

# 如果聚类数与注释数不匹配，使用通用注释
if (n_clusters == length(new.cluster.ids)) {
    names(new.cluster.ids) <- levels(pbmc)
    pbmc <- RenameIdents(pbmc, new.cluster.ids)
} else {
    cat("聚类数与预设注释数不一致，使用通用注释\n")
    # 使用数字标签
    new.cluster.ids <- paste0("Cluster", 0:(n_clusters - 1))
    names(new.cluster.ids) <- levels(pbmc)
    pbmc <- RenameIdents(pbmc, new.cluster.ids)
}
pbmc$cell_type <- Idents(pbmc)

# 6.1 注释后 UMAP
p_annotated <- DimPlot(pbmc, reduction = "umap.cca", label = TRUE, pt.size = 0.5) + NoLegend() +
    ggtitle("整合后 - 细胞类型注释")
ggsave(file.path(integ_dir, "09_annotated_umap.png"), plot = p_annotated, width = 8, height = 7, dpi = 150)
cat("已保存: 09_annotated_umap.png\n")

# 6.2 Marker FeaturePlot
p_markers_fp <- FeaturePlot(pbmc, features = c("IL7R", "CD14", "MS4A1", "CD8A", "GNLY", "FCGR3A"),
                             reduction = "umap.cca", ncol = 3)
ggsave(file.path(integ_dir, "10_marker_featureplot.png"), plot = p_markers_fp, width = 14, height = 10, dpi = 150)
cat("已保存: 10_marker_featureplot.png\n")

# ============================================================================
# 7. 保守 Marker 基因（FindConservedMarkers）
# ============================================================================
cat("\n========== 7. 保守 Marker 基因 ==========\n")

# FindConservedMarkers：在两个批次中分别寻找 marker，然后用 meta-analysis 合并 p 值
# 这可以找到不受批次效应影响的稳定 marker 基因

# 注意：FindConservedMarkers 需要 layers 处于合并状态
# （它在内部按 grouping.var 分组分析，不需要预先 split）
# 确保 layers 已合并
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# 以某个聚类为例（选择最大的聚类）
Idents(pbmc) <- "cca_clusters"
target_cluster <- as.character(which.max(table(pbmc$cca_clusters)))
cat("分析聚类", target_cluster, "的保守 marker\n")

conserved_markers <- FindConservedMarkers(
    pbmc,
    ident.1 = target_cluster,
    grouping.var = "batch",
    verbose = FALSE
)

cat("保守 marker 前 15 个:\n")
print(head(conserved_markers, 15))

# 保存保守 marker
write.csv(conserved_markers, file.path(integ_dir, "conserved_markers.csv"), row.names = TRUE)
cat("已保存: conserved_markers.csv\n")

# 7.1 保守 marker 小提琴图（按批次分割）
top_conserved <- rownames(head(conserved_markers, 6))
p_conserved_vln <- VlnPlot(pbmc, features = top_conserved, split.by = "batch", pt.size = 0, ncol = 3) +
    ggtitle("保守 Marker - 按批次分割")
ggsave(file.path(integ_dir, "11_conserved_markers_violin.png"), plot = p_conserved_vln, width = 14, height = 10, dpi = 150)
cat("已保存: 11_conserved_markers_violin.png\n")

# 7.2 保守 marker DotPlot（按批次分割）
# layers 已合并，直接使用
markers_to_plot <- rownames(head(conserved_markers, 12))
p_conserved_dot <- DotPlot(pbmc, features = markers_to_plot, split.by = "batch") + RotatedAxis() +
    ggtitle("保守 Marker DotPlot - 按批次对比")
ggsave(file.path(integ_dir, "12_conserved_markers_dot.png"), plot = p_conserved_dot, width = 12, height = 6, dpi = 150)
cat("已保存: 12_conserved_markers_dot.png\n")

# ============================================================================
# 8. 跨批次差异表达分析
# ============================================================================
cat("\n========== 8. 跨批次差异表达 ==========\n")

# 创建 celltype.batch 组合变量
pbmc$celltype.batch <- paste(pbmc$cell_type, pbmc$batch, sep = "_")
Idents(pbmc) <- "celltype.batch"

# 选择一个细胞类型，比较两个批次的差异
# 尝试找到 CD14+ Mono 或第一个细胞类型
unique_idents <- unique(pbmc$celltype.batch)
cat("可用的 celltype.batch 组合:\n")
print(table(pbmc$celltype.batch))

# 选择第一个细胞类型进行跨批次比较
# 找到同时存在于两个批次的细胞类型
batch1_types <- unique(pbmc$cell_type[pbmc$batch == "batch1"])
batch2_types <- unique(pbmc$cell_type[pbmc$batch == "batch2"])
common_types <- intersect(batch1_types, batch2_types)
cat("两个批次共有的细胞类型:\n")
print(common_types)

if (length(common_types) > 0) {
    target_type <- common_types[1]
    ident1 <- paste0(target_type, "_batch2")
    ident2 <- paste0(target_type, "_batch1")

    if (ident1 %in% unique_idents && ident2 %in% unique_idents) {
        de_genes <- FindMarkers(pbmc, ident.1 = ident1, ident.2 = ident2, verbose = FALSE)
        cat("\n", target_type, "的跨批次差异基因（前 15 个）:\n")
        print(head(de_genes, 15))

        write.csv(de_genes, file.path(integ_dir, "batch_de_genes.csv"), row.names = TRUE)
        cat("已保存: batch_de_genes.csv\n")

        # 8.1 差异基因小提琴图
        top_de <- rownames(head(de_genes, 6))
        if (length(top_de) > 0) {
            p_de_vln <- VlnPlot(pbmc, features = top_de, split.by = "batch", pt.size = 0, ncol = 3) +
                ggtitle(paste0("跨批次差异基因 - ", target_type))
            ggsave(file.path(integ_dir, "13_batch_de_violin.png"), plot = p_de_vln, width = 14, height = 10, dpi = 150)
            cat("已保存: 13_batch_de_violin.png\n")
        }
    }
}

# ============================================================================
# 9. RPCA 整合（更快速、更保守的整合方法）
# ============================================================================
cat("\n========== 9. RPCA 整合 ==========\n")

# RPCA 整合需要按批次分割 layers
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)

# RPCA 整合（比 CCA 更快、更保守）
cat("正在运行 RPCA 整合...\n")
pbmc <- IntegrateLayers(
    object = pbmc,
    method = RPCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.rpca",
    verbose = FALSE
)
cat("RPCA 整合完成\n")

# 合并 layers
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# RPCA 聚类和 UMAP
pbmc <- FindNeighbors(pbmc, reduction = "integrated.rpca", dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5, cluster.name = "rpca_clusters")
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "integrated.rpca", reduction.name = "umap.rpca")

# 9.1 RPCA 整合后 UMAP - 按批次着色
p_rpca_batch <- DimPlot(pbmc, reduction = "umap.rpca", group.by = "batch") +
    ggtitle("RPCA 整合后 - 按批次着色")
ggsave(file.path(integ_dir, "14_rpca_batch.png"), plot = p_rpca_batch, width = 8, height = 7, dpi = 150)
cat("已保存: 14_rpca_batch.png\n")

# 9.2 CCA vs RPCA 对比
p_cca_rpca <- DimPlot(pbmc, reduction = "umap.cca", group.by = "batch") +
    ggtitle("CCA 整合") +
DimPlot(pbmc, reduction = "umap.rpca", group.by = "batch") +
    ggtitle("RPCA 整合")
ggsave(file.path(integ_dir, "15_cca_vs_rpca.png"), plot = p_cca_rpca, width = 16, height = 6, dpi = 150)
cat("已保存: 15_cca_vs_rpca.png\n")

# ============================================================================
# 10. 保存结果
# ============================================================================
cat("\n========== 10. 保存结果 ==========\n")

saveRDS(pbmc, file.path(output_dir, "pbmc3k_integrated.rds"))
cat("已保存整合后对象: pbmc3k_integrated.rds\n")

# ============================================================================
# 总结
# ============================================================================
cat("\n")
cat("=======================================================\n")
cat("  阶段3 数据整合教程完成!\n")
cat("=======================================================\n")
cat("  输出目录:", integ_dir, "\n")
cat("  图片数量: 15 张\n")
cat("  覆盖内容:\n")
cat("    01-03: 未整合分析（展示批次效应）\n")
cat("    04-06: CCA 整合结果\n")
cat("    07-08: 整合前后对比\n")
cat("    09-10: 细胞类型注释\n")
cat("    11-12: 保守 Marker 基因\n")
cat("    13:    跨批次差异表达\n")
cat("    14-15: RPCA 整合与 CCA 对比\n")
cat("  核心知识点:\n")
cat("    - Seurat v5 layers 机制（split / Join）\n")
cat("    - IntegrateLayers 一行代码整合\n")
cat("    - CCA vs RPCA 两种整合方法\n")
cat("    - FindConservedMarkers 跨批次保守 marker\n")
cat("    - FindMarkers 跨批次差异表达\n")
cat("    - 整合前后对比可视化\n")
cat("=======================================================\n")

本教程基于 Seurat v5.5.1 官方文档整理
R 4.6.1 + Seurat 5.5.1 + Windows 环境
阶段3: 数据整合分析（CCA / RPCA / 保守 Marker / 跨批次差异表达）

阅读全文

20260912 093449 Seurat V5 数据可视化方法与 Sctransform 归一化教程

2026年09月12日

Seurat v5 数据可视化方法与 SCTransform 归一化教程 ¶

基于 Seurat v5.5.1 官方教程，系统学习单细胞数据的可视化技巧与 SCTransform 先进归一化方法。
阶段2-A：数据可视化方法（20张图，5大板块）
阶段2-B：SCTransform 归一化（9张图，7个步骤）

方法	函数	特点	适用场景
山脊图	`RidgePlot()`	展示表达密度分布	比较多个基因的表达模式
小提琴图	`VlnPlot()`	展示表达分布+密度	各聚类间基因表达比较
Feature Plot	`FeaturePlot()`	在低维空间展示表达	空间分布模式
气泡图	`DotPlot()`	大小=比例，颜色=表达量	大量基因×聚类的概览
热图	`DoHeatmap()`	颜色梯度展示表达	多基因多聚类综合展示

Part B: SCTransform 归一化 ¶

7. SCTransform 原理与流程 ¶

7.1 为什么需要 SCTransform？¶

传统标准流程使用 NormalizeData() + ScaleData() + FindVariableFeatures() 三步处理：
- LogNormalize：简单但对测序深度差异敏感
- 高变基因选择：可能受技术因素影响
- 标准化：无法完全去除批次效应

SCTransform（v2）使用负二项回归模型一步完成：
1. 对每个基因拟合 count 模型（考虑测序深度等协变量）
2. 计算 Pearson 残差（去除技术变异后的生物学信号）
3. 选择高变基因（默认 3000 个，比标准的 2000 个更多）
4. 生成「校正后」的 UMI counts（可理解为统一测序深度后的预期值）

7.2 核心代码 ¶

library(Seurat)
library(ggplot2)
library(sctransform)

# 从原始数据开始（独立于阶段1的标准流程）
pbmc_data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc_data)

# 计算线粒体百分比
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# QC 过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

# SCTransform 一步完成归一化 + 变量特征选择 + 标准化
# 同时回归线粒体百分比（去除技术变异）
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

7.3 SCT Assay 结构 ¶

运行 SCTransform() 后，创建了一个名为 SCT 的新 Assay 并设为默认：

层级	内容	用途
`scale.data`	Pearson 残差	PCA 输入
`counts`	校正后的 UMI counts	重新计算归一化值
`data`	校正 counts 的 log 归一化	可视化

> DefaultAssay(pbmc)
[1] "SCT"

> pbmc[["SCT"]]
SCTAssay data with 12519 features for 2638 cells, and 1 SCTModel(s)
Top 10 variable features:
 S100A9, GNLY, LYZ, S100A8, NKG7, FTL, GZMB, IGLL5, FTH1, CCL5

7.4 运行信息 ¶

细胞数：2,700 -> 2,638（QC 过滤后）
基因数：32,738 -> 12,519（过滤后参与建模）
耗时：约 5 分钟（未安装 glmGamPoi 加速包，使用原生实现）
变量特征：3,000 个（默认，比标准流程多 1,000 个）
模型版本：v2（Seurat v5 默认，固定斜率模型 + 排除 Poisson 基因）

亚群	标志基因	说明
Naive CD8 T	CD8A, CCR7	初始 T 细胞
Memory CD8 T	GZMK	记忆 T 细胞
Effector CD8 T	CCL5	效应 T 细胞

亚群	标志基因	说明
Naive CD4 T	CCR7, IL32	初始 T 细胞
Memory CD4 T	S100A4, ANXA1	记忆 T 细胞
IFN-activated	ISG15	干素激活

10. SCT vs 标准流程对比 ¶

10.3 对比总结 ¶

指标	标准流程	SCTransform
归一化方法	LogNormalize	负二项回归（v2）
变量特征数	2,000	3,000
推荐 PC 数	10	30
回归变量	可选	内置支持
测序深度影响	仍有残留	显著减弱
细胞亚群分辨率	基础	更精细（CD8 naive/memory/effector）
计算耗时	秒级	分钟级（约5分钟）
内存占用	较低	较高（非稀疏矩阵）
适用场景	快速探索	正式分析/发表

10.4 何时选择 SCTransform？¶

数据中测序深度差异较大时
需要检测精细细胞亚群时
需要同时回归多个协变量时
准备发表正式结果时
跨数据集整合分析时（SCT 是整合分析的前提）

10.5 一行代码完成全流程 ¶

SCTransform 与管道操作符配合，可以一行代码完成完整分析：

library(dplyr)
pbmc <- CreateSeuratObject(pbmc_data) %>%
    PercentageFeatureSet(pattern = "^MT-", col.name = "percent.mt") %>%
    SCTransform(vars.to.regress = "percent.mt") %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    RunUMAP(dims = 1:30) %>%
    FindClusters()

附录: 完整脚本 ¶

A. 可视化教程脚本（visualization_vignette.R）¶

###############################################################################
# Seurat v5 阶段2-A：数据可视化方法
# 基于官方教程: https://satijalab.org/seurat/articles/visualization_vignette
# 使用阶段1保存的 pbmc3k_final.rds 对象
###############################################################################

library(Seurat)
library(ggplot2)
library(patchwork)

# === 路径配置 ===
output_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output"
viz_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output/visualization"
dir.create(viz_dir, showWarnings = FALSE, recursive = TRUE)

# === 加载阶段1的最终对象 ===
pbmc <- readRDS(file.path(output_dir, "pbmc3k_final.rds"))
cat("已加载 pbmc3k_final.rds\n")
print(pbmc)

# 添加分组变量（用于演示 split.by 功能）
set.seed(42)
pbmc$groups <- sample(c("group1", "group2"), size = ncol(pbmc), replace = TRUE)

# 定义演示用基因列表
features <- c("LYZ", "CCL5", "IL32", "PTPRCAP", "FCGR3A", "PF4")


# ============================================================================
# 1. 五种可视化方法展示 Marker 基因表达
# ============================================================================

cat("\n========== 1. 五种可视化方法 ==========\n")

# 1.1 Ridge Plot（山脊图）- 展示各聚类的基因表达分布
p_ridge <- RidgePlot(pbmc, features = features, ncol = 2)
ggsave(file.path(viz_dir, "01_ridge_plot.png"), plot = p_ridge, width = 12, height = 10, dpi = 150)
cat("已保存: 01_ridge_plot.png\n")

# 1.2 Violin Plot（小提琴图）- 展示各聚类的基因表达分布
p_vln <- VlnPlot(pbmc, features = features)
ggsave(file.path(viz_dir, "02_violin_plot.png"), plot = p_vln, width = 14, height = 8, dpi = 150)
cat("已保存: 02_violin_plot.png\n")

# 1.3 Feature Plot - 在低维空间展示基因表达
p_feature <- FeaturePlot(pbmc, features = features)
ggsave(file.path(viz_dir, "03_feature_plot.png"), plot = p_feature, width = 14, height = 10, dpi = 150)
cat("已保存: 03_feature_plot.png\n")

# 1.4 Dot Plot（气泡图）- 点大小=表达比例，颜色=平均表达水平
p_dot <- DotPlot(pbmc, features = features) + RotatedAxis()
ggsave(file.path(viz_dir, "04_dot_plot.png"), plot = p_dot, width = 10, height = 6, dpi = 150)
cat("已保存: 04_dot_plot.png\n")

# 1.5 Heatmap（热图）- 降采样后展示基因表达
p_heat <- DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)
ggsave(file.path(viz_dir, "05_heatmap.png"), plot = p_heat, width = 10, height = 6, dpi = 150)
cat("已保存: 05_heatmap.png\n")


# ============================================================================
# 2. 自定义 FeaturePlot
# ============================================================================

cat("\n========== 2. 自定义 FeaturePlot ==========\n")

# 2.1 带图例的 FeaturePlot
p_fp1 <- FeaturePlot(pbmc, features = "MS4A1")
ggsave(file.path(viz_dir, "06_featureplot_legend.png"), plot = p_fp1, width = 7, height = 6, dpi = 150)
cat("已保存: 06_featureplot_legend.png\n")

# 2.2 调整对比度（设定 min/max cutoff）
p_fp2 <- FeaturePlot(pbmc, features = "MS4A1", min.cutoff = 1, max.cutoff = 3)
ggsave(file.path(viz_dir, "07_featureplot_contrast.png"), plot = p_fp2, width = 7, height = 6, dpi = 150)
cat("已保存: 07_featureplot_contrast.png\n")

# 2.3 分位数 cutoff（适合多基因同时绘图）
p_fp3 <- FeaturePlot(pbmc, features = c("MS4A1", "PTPRCAP"), min.cutoff = "q10", max.cutoff = "q90")
ggsave(file.path(viz_dir, "08_featureplot_quantile.png"), plot = p_fp3, width = 12, height = 5, dpi = 150)
cat("已保存: 08_featureplot_quantile.png\n")

# 2.4 混合展示两个基因的共表达（blend 模式）
p_fp4 <- FeaturePlot(pbmc, features = c("MS4A1", "CD79A"), blend = TRUE)
ggsave(file.path(viz_dir, "09_featureplot_blend.png"), plot = p_fp4, width = 14, height = 5, dpi = 150)
cat("已保存: 09_featureplot_blend.png\n")

# 2.5 按分组分割展示
p_fp5 <- FeaturePlot(pbmc, features = c("MS4A1", "CD79A"), split.by = "groups")
ggsave(file.path(viz_dir, "10_featureplot_split.png"), plot = p_fp5, width = 14, height = 10, dpi = 150)
cat("已保存: 10_featureplot_split.png\n")


# ============================================================================
# 3. 其他可视化功能
# ============================================================================

cat("\n========== 3. 其他可视化功能 ==========\n")

# 3.1 按分组分割的小提琴图
p_vln_split <- VlnPlot(pbmc, features = "percent.mt", split.by = "groups")
ggsave(file.path(viz_dir, "11_violin_split.png"), plot = p_vln_split, width = 8, height = 5, dpi = 150)
cat("已保存: 11_violin_split.png\n")

# 3.2 按分组分割的气泡图
p_dot_split <- DotPlot(pbmc, features = features, split.by = "groups") + RotatedAxis()
ggsave(file.path(viz_dir, "12_dotplot_split.png"), plot = p_dot_split, width = 12, height = 6, dpi = 150)
cat("已保存: 12_dotplot_split.png\n")

# 3.3 DimPlot 默认使用 UMAP
p_dim1 <- DimPlot(pbmc)
ggsave(file.path(viz_dir, "13_dimplot_umap.png"), plot = p_dim1, width = 8, height = 7, dpi = 150)
cat("已保存: 13_dimplot_umap.png\n")

# 3.4 DoHeatmap 带分组条
p_heat2 <- DoHeatmap(pbmc, features = VariableFeatures(pbmc)[1:100], cells = 1:500, size = 4, angle = 90) + NoLegend()
ggsave(file.path(viz_dir, "14_heatmap_grouped.png"), plot = p_heat2, width = 10, height = 8, dpi = 150)
cat("已保存: 14_heatmap_grouped.png\n")


# ============================================================================
# 4. 应用主题（Theme）
# ============================================================================

cat("\n========== 4. 应用主题 ==========\n")

baseplot <- DimPlot(pbmc, reduction = "umap")

# 4.1 添加标题
p_theme1 <- baseplot + labs(title = "Clustering of 2,700 PBMCs")
ggsave(file.path(viz_dir, "15_theme_title.png"), plot = p_theme1, width = 8, height = 7, dpi = 150)
cat("已保存: 15_theme_title.png\n")

# 4.2 暗色主题
p_theme2 <- baseplot + DarkTheme()
ggsave(file.path(viz_dir, "16_theme_dark.png"), plot = p_theme2, width = 8, height = 7, dpi = 150)
cat("已保存: 16_theme_dark.png\n")

# 4.3 链式主题
p_theme3 <- baseplot + DarkTheme() +
    theme(axis.title = element_text(size = 18), legend.text = element_text(size = 18)) +
    guides(colour = guide_legend(override.aes = list(size = 10)))
ggsave(file.path(viz_dir, "17_theme_chained.png"), plot = p_theme3, width = 10, height = 8, dpi = 150)
cat("已保存: 17_theme_chained.png\n")


# ============================================================================
# 5. 图组合并
# ============================================================================

cat("\n========== 5. 图组合并 ==========\n")

plot1 <- DimPlot(pbmc)
plot2 <- FeatureScatter(pbmc, feature1 = "LYZ", feature2 = "CCL5")

# 5.1 拼接两图
p_combine1 <- plot1 + plot2
ggsave(file.path(viz_dir, "18_combine_plots.png"), plot = p_combine1, width = 14, height = 6, dpi = 150)
cat("已保存: 18_combine_plots.png\n")

# 5.2 移除所有图例
p_combine2 <- (plot1 + plot2) & NoLegend()
ggsave(file.path(viz_dir, "19_combine_nolegend.png"), plot = p_combine2, width = 14, height = 6, dpi = 150)
cat("已保存: 19_combine_nolegend.png\n")


# ============================================================================
# 6. 聚类标签
# ============================================================================

cat("\n========== 6. 聚类标签 ==========\n")

# 6.1 标注聚类
p_label1 <- DimPlot(pbmc, reduction = "pca") + NoLegend()
p_label1 <- LabelClusters(plot = p_label1, id = "ident")
ggsave(file.path(viz_dir, "20_label_clusters.png"), plot = p_label1, width = 8, height = 7, dpi = 150)
cat("已保存: 20_label_clusters.png\n")


# ============================================================================
# 总结
# ============================================================================

cat("\n")
cat("=======================================================\n")
cat("  阶段2-A 可视化教程完成!\n")
cat("=======================================================\n")
cat("  输出目录:", viz_dir, "\n")
cat("  图片数量: 20 张\n")
cat("  覆盖内容:\n")
cat("    01-05: 五种核心可视化方法\n")
cat("    06-10: FeaturePlot 自定义技巧\n")
cat("    11-14: 分组分割与热图\n")
cat("    15-17: 主题定制\n")
cat("    18-19: 图组合并\n")
cat("    20:    聚类标签\n")
cat("=======================================================\n")

B. SCTransform 脚本（sctransform_vignette.R）¶

###############################################################################
# Seurat v5 阶段2-B：SCTransform 归一化分析
# 基于官方教程: https://satijalab.org/seurat/articles/sctransform_vignette
# SCTransform 替代 NormalizeData + ScaleData + FindVariableFeatures 三步
# 使用 PBMC 3k 数据从头创建新对象（独立于阶段1的标准流程）
###############################################################################

library(Seurat)
library(ggplot2)
library(sctransform)
library(patchwork)

# === 路径配置 ===
data_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/data/filtered_gene_bc_matrices/hg19"
output_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output"
sct_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output/sctransform"
dir.create(sct_dir, showWarnings = FALSE, recursive = TRUE)

# ============================================================================
# 1. 加载数据并创建 Seurat 对象
# ============================================================================
cat("\n========== 1. 加载数据 ==========\n")

pbmc_data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc_data)
cat("已创建 Seurat 对象\n")
print(pbmc)

# ============================================================================
# 2. 计算线粒体百分比并应用 SCTransform
# ============================================================================
cat("\n========== 2. SCTransform 归一化 ==========\n")

# 计算线粒体基因百分比
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
cat("已计算线粒体百分比\n")

# QC 过滤（与阶段1一致）
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
cat("QC 过滤后细胞数:\n")
cat(ncol(pbmc), "\n")

# 应用 SCTransform（同时回归线粒体百分比）
# SCTransform 一步完成：归一化 + 变量特征选择 + 标准化
cat("正在运行 SCTransform（可能需要1-2分钟）...\n")
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = TRUE)
cat("SCTransform 完成!\n")
cat("当前默认 assay:\n")
print(DefaultAssay(pbmc))

# ============================================================================
# 3. 降维分析：PCA + UMAP
# ============================================================================
cat("\n========== 3. 降维分析 ==========\n")

# PCA（使用 SCT residuals）
pbmc <- RunPCA(pbmc, verbose = FALSE)
cat("PCA 完成\n")

# UMAP（使用前30个PC，SCTransform 推荐使用更多PC）
pbmc <- RunUMAP(pbmc, dims = 1:30, verbose = FALSE)
cat("UMAP 完成\n")

# 3.1 PCA 可视化
p_pca <- DimPlot(pbmc, reduction = "pca")
ggsave(file.path(sct_dir, "01_pca.png"), plot = p_pca, width = 8, height = 7, dpi = 150)
cat("已保存: 01_pca.png\n")

# 3.2 UMAP 可视化
p_umap <- DimPlot(pbmc, reduction = "umap")
ggsave(file.path(sct_dir, "02_umap.png"), plot = p_umap, width = 8, height = 7, dpi = 150)
cat("已保存: 02_umap.png\n")

# 3.3 PCA Elbow Plot（选择PC数量参考）
p_elbow <- ElbowPlot(pbmc, ndims = 50)
ggsave(file.path(sct_dir, "03_elbow_plot.png"), plot = p_elbow, width = 8, height = 5, dpi = 150)
cat("已保存: 03_elbow_plot.png\n")

# ============================================================================
# 4. 聚类分析
# ============================================================================
cat("\n========== 4. 聚类分析 ==========\n")

pbmc <- FindNeighbors(pbmc, dims = 1:30, verbose = FALSE)
pbmc <- FindClusters(pbmc, verbose = FALSE)
cat("聚类完成\n")

# 4.1 带标签的 UMAP
p_umap_label <- DimPlot(pbmc, reduction = "umap", label = TRUE)
ggsave(file.path(sct_dir, "04_umap_labeled.png"), plot = p_umap_label, width = 8, height = 7, dpi = 150)
cat("已保存: 04_umap_labeled.png\n")

# ============================================================================
# 5. SCTransform 优势展示：更精细的细胞亚群
# ============================================================================
cat("\n========== 5. SCTransform 优势展示 ==========\n")

# 5.1 T 细胞亚群标志基因小提琴图
# SCTransform 能更清晰地区分 CD8 T 细胞亚群（naive/memory/effector）
p_vln_t <- VlnPlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7", "ISG15", "CD3D"),
                   pt.size = 0.2, ncol = 4)
ggsave(file.path(sct_dir, "05_violin_tcell_markers.png"), plot = p_vln_t, width = 16, height = 10, dpi = 150)
cat("已保存: 05_violin_tcell_markers.png\n")

# 5.2 T 细胞标志基因 FeaturePlot
p_fp_t <- FeaturePlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7"),
                       pt.size = 0.2, ncol = 3)
ggsave(file.path(sct_dir, "06_featureplot_tcell.png"), plot = p_fp_t, width = 14, height = 10, dpi = 150)
cat("已保存: 06_featureplot_tcell.png\n")

# 5.3 B 细胞与 NK 细胞标志基因 FeaturePlot
p_fp_bnk <- FeaturePlot(pbmc, features = c("CD3D", "ISG15", "TCL1A", "FCER2", "XCL1", "FCGR3A"),
                         pt.size = 0.2, ncol = 3)
ggsave(file.path(sct_dir, "07_featureplot_bnk.png"), plot = p_fp_bnk, width = 14, height = 10, dpi = 150)
cat("已保存: 07_featureplot_bnk.png\n")

# ============================================================================
# 6. SCTransform 与标准归一化对比
# ============================================================================
cat("\n========== 6. 对比分析 ==========\n")

# 6.1 SCTransform 变量特征数量（默认3000个，比标准流程的2000更多）
sct_features <- VariableFeatures(pbmc)
cat("SCTransform 变量特征数量:", length(sct_features), "\n")

# 6.2 展示前20个变量特征
p_var <- LabelPoints(plot = VariableFeaturePlot(pbmc), points = sct_features[1:20], repel = TRUE)
ggsave(file.path(sct_dir, "08_variable_features.png"), plot = p_var, width = 10, height = 8, dpi = 150)
cat("已保存: 08_variable_features.png\n")

# 6.3 SCT assay 结构说明
cat("\nSCT Assay 层级结构:\n")
print(pbmc[["SCT"]])

# 6.4 标准流程对比图：使用前10个PC（标准流程推荐）vs 前30个PC（SCT推荐）
p_umap_10 <- DimPlot(pbmc, reduction = "umap", group.by = "seurat_clusters") +
    ggtitle("SCT: 30 PCs (推荐)")
ggsave(file.path(sct_dir, "09_umap_sct_30pc.png"), plot = p_umap_10, width = 8, height = 7, dpi = 150)
cat("已保存: 09_umap_sct_30pc.png\n")

# ============================================================================
# 7. 保存 SCTransform 结果
# ============================================================================
cat("\n========== 7. 保存结果 ==========\n")

saveRDS(pbmc, file.path(output_dir, "pbmc3k_sct.rds"))
cat("已保存 SCTransform 对象: pbmc3k_sct.rds\n")

# ============================================================================
# 总结
# ============================================================================
cat("\n")
cat("=======================================================\n")
cat("  阶段2-B SCTransform 教程完成!\n")
cat("=======================================================\n")
cat("  输出目录:", sct_dir, "\n")
cat("  图片数量: 9 张\n")
cat("  覆盖内容:\n")
cat("    01-03: PCA + UMAP + Elbow Plot\n")
cat("    04:    聚类标注\n")
cat("    05-07: T/B/NK 细胞亚群标志基因可视化\n")
cat("    08:    SCTransform 变量特征\n")
cat("    09:    SCT UMAP（30 PCs）\n")
cat("  SCT vs 标准流程:\n")
cat("    - 变量特征: 3000（SCT） vs 2000（标准）\n")
cat("    - 推荐PC数: 30（SCT） vs 10（标准）\n")
cat("    - 归一化方法: 负二项回归 vs LogNormalize\n")
cat("    - 能揭示更精细的细胞亚群（CD8 naive/memory/effector）\n")
cat("=======================================================\n")

输出文件清单 ¶

Part A: 可视化教程（20张图）¶

编号	文件名	内容
01	01_ridge_plot.png	山脊图
02	02_violin_plot.png	小提琴图
03	03_feature_plot.png	Feature Plot
04	04_dot_plot.png	气泡图
05	05_heatmap.png	热图
06	06_featureplot_legend.png	带图例 FeaturePlot
07	07_featureplot_contrast.png	调整对比度
08	08_featureplot_quantile.png	分位数 cutoff
09	09_featureplot_blend.png	混合模式
10	10_featureplot_split.png	分组分割
11	11_violin_split.png	分组小提琴图
12	12_dotplot_split.png	分组气泡图
13	13_dimplot_umap.png	UMAP
14	14_heatmap_grouped.png	分组热图
15	15_theme_title.png	添加标题
16	16_theme_dark.png	暗色主题
17	17_theme_chained.png	链式主题
18	18_combine_plots.png	图组合并
19	19_combine_nolegend.png	移除图例
20	20_label_clusters.png	聚类标签

Part B: SCTransform（9张图）¶

编号	文件名	内容
01	01_pca.png	PCA
02	02_umap.png	UMAP
03	03_elbow_plot.png	Elbow Plot
04	04_umap_labeled.png	聚类标注 UMAP
05	05_violin_tcell_markers.png	T 细胞标志小提琴图
06	06_featureplot_tcell.png	T 细胞 FeaturePlot
07	07_featureplot_bnk.png	B/NK 细胞 FeaturePlot
08	08_variable_features.png	变量特征图
09	09_umap_sct_30pc.png	SCT UMAP（30 PCs）

本教程基于 Seurat v5.5.1 官方文档整理
R 4.6.1 + Seurat 5.5.1 + Windows 环境
阶段2-A: 数据可视化方法 | 阶段2-B: SCTransform 归一化

阅读全文

20260912 091245 Seurat V5 Pbmc 3K 单细胞分析完整教程

2026年09月12日

Seurat v5 PBMC 3k 单细胞分析完整教程 ¶

基于 Seurat v5.5.1 官方教程，从零开始跑通单细胞 RNA-seq 标准分析全流程。

1. 概述 ¶

1.1 数据集 ¶

项目	说明
数据来源	10X Genomics PBMC（外周血单核细胞）
细胞数	2,700
测序平台	Illumina NextSeq 500
参考基因组	hg19
数据格式	10X CellRanger 输出（barcodes.tsv, genes.tsv, matrix.mtx）

1.2 分析流程 ¶

读取数据 -> QC过滤 -> 归一化 -> 高变基因 -> 标准化 -> PCA
-> 确定维度 -> 聚类 -> UMAP -> 寻找Marker -> 细胞类型注释

1.3 环境信息 ¶

软件	版本
R	4.6.1 (2026-06-24, “Happy Hop”)
Seurat	5.5.1
SeuratObject	5.4.0
操作系统	Windows (D盘安装: D:\r)

PC	Positive（正加载）	Negative（负加载）
PC1	CST3, TYROBP, LST1, AIF1, FTL	MALAT1, LTB, IL32, IL7R, CD2
PC2	CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1	NKG7, PRF1, CST7, GZMB, GZMA
PC3	HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1	PPBP, PF4, SDPR, SPARC, GNG11
PC4	HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1	VIM, IL7R, S100A6, IL32, S100A8
PC5	GZMB, NKG7, S100A8, FGFBP2, GNLY	LTB, IL7R, CKB, VIM, MS4A7

聚类	Marker 数（log2FC > 1）
0	214
1	1,190
2	384
3	543
4	344
5	1,241
6	891
7	1,589
8	623

聚类	Marker 基因	细胞类型
0	IL7R, CCR7	Naive CD4+ T
1	CD14, LYZ	CD14+ Mono（单核细胞）
2	IL7R, S100A4	Memory CD4+ T
3	MS4A1	B（B细胞）
4	CD8A	CD8+ T
5	FCGR3A, MS4A7	FCGR3A+ Mono
6	GNLY, NKG7	NK（自然杀伤细胞）
7	FCER1A, CST3	DC（树突状细胞）
8	PPBP	Platelet（血小板）

4. 分析总结 ¶

指标	数值
输入细胞数	2,700
过滤后细胞数	2,638
基因数	13,714
高变基因数	2,000
PCA 维度	10
聚类数	9
细胞类型数	9
Marker 基因总数	7,019（log2FC > 1）

5. 输出文件清单 ¶

文件	说明
01_qc_violin.png	QC 小提琴图
02_qc_scatter.png	QC 散点图
03_variable_features.png	高变基因图
04_pca.png	PCA 散点图
05_pca_loadings.png	PCA 负载图
06_pca_heatmap.png	PCA 热图
07_elbow_plot.png	Elbow 拐点图
08_umap.png	UMAP 聚类图
09_marker_heatmap.png	Marker 热图
10_feature_plot.png	Feature 表达图
11_violin_markers.png	Marker 小提琴图
12_umap_labeled.png	注释后 UMAP 图
all_markers.csv	全部 Marker 基因表
pbmc_tutorial.rds	中间结果（Seurat 对象）
pbmc3k_final.rds	最终结果（Seurat 对象）

6. 关键知识点回顾 ¶

6.1 Seurat 对象结构（v5）¶

Seurat 对象
|-- assays
|   +-- RNA
|       |-- layers
|       |   |-- counts    # 原始 UMI 计数矩阵
|       |   |-- data       # 归一化后的矩阵
|       |   +-- scale.data # 标准化后的矩阵
|       +-- VariableFeatures  # 高变基因列表
|-- reductions
|   |-- pca               # PCA 降维结果
|   +-- umap              # UMAP 降维结果
|-- meta.data             # 细胞元数据（nCount, nFeature, percent.mt 等）
+-- active.ident          # 当前聚类标签

6.2 标准分析步骤要点 ¶

步骤	核心函数	关键参数	目的
读取数据	Read10X()	data.dir	读取 10X 输出
创建对象	CreateSeuratObject()	min.cells, min.features	初始过滤
QC	PercentageFeatureSet()	pattern=”^MT-“	计算线粒体比例
归一化	NormalizeData()	LogNormalize, scale.factor=10000	消除细胞间测序深度差异
高变基因	FindVariableFeatures()	vst, nfeatures=2000	识别信息量最大的基因
标准化	ScaleData()	features=all.genes	Z-score 标准化
降维	RunPCA()	features=VariableFeatures	线性降维
确定维度	ElbowPlot()	-	选择 PC 数量
聚类	FindNeighbors() + FindClusters()	dims, resolution	图聚类
可视化	RunUMAP() + DimPlot()	dims	2D 可视化
Marker	FindAllMarkers()	only.pos=TRUE	差异基因分析
注释	RenameIdents()	-	映射细胞类型

6.3 常用参数调优建议 ¶

resolution：0.4-1.2 适合约3K细胞；更大数据集可适当提高
PC 数量：建议 7-15 之间，宁多勿少（太少会丢失信号）
percent.mt 阈值：常用 5%-10%，不同组织类型可能需要调整
nFeature_RNA 过滤：下限 200-500，上限通常取 2500-6000

7. 完整脚本 ¶

以下是完整的 R 分析脚本，可直接运行：

###############################################################################
# Seurat v5 PBMC 3k 标准分析流程
# 基于官方教程: https://satijalab.org/seurat/articles/pbmc3k_tutorial
# 数据: 10X Genomics PBMC 2,700 个外周血单核细胞
###############################################################################
#
# 学习目标: 跑通单细胞 RNA-seq 标准分析全流程
# 流程链路: 读取数据 -> QC过滤 -> 归一化 -> 高变基因 -> 标准化 -> PCA
#           -> 确定维度 -> 聚类 -> UMAP -> 寻找marker -> 细胞类型注释
#
###############################################################################


# ============================================================================
# 0. 环境准备
# ============================================================================

library(dplyr)
library(Seurat)
library(patchwork)
library(ggplot2)

# 设置工作目录（请根据实际路径修改）
data_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/data/filtered_gene_bc_matrices/hg19"
output_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)


# ============================================================================
# 1. 读取数据 & 创建 Seurat 对象
# ============================================================================
# Read10X() 读取 10X CellRanger 输出目录（包含 barcodes.tsv, genes.tsv, matrix.mtx）
# 返回一个稀疏矩阵: 行=基因, 列=细胞, 值=UMI计数

cat("\n========== 1. 读取数据 ==========\n")
pbmc.data <- Read10X(data.dir = data_dir)

# 创建 Seurat 对象
# min.cells = 3:   至少在 3 个细胞中检测到的基因才保留
# min.features = 200: 每个细胞至少检测到 200 个基因才保留
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200)
print(pbmc)

# 查看计数矩阵的稀疏性（大多数值为 0，用稀疏矩阵存储节省内存）
cat("\n--- 稀疏矩阵示例（前30个细胞中3个基因）---\n")
print(pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30])


# ============================================================================
# 2. 质量控制（QC）& 细胞过滤
# ============================================================================
# 常用 QC 指标:
#   - nFeature_RNA: 每个细胞检测到的独特基因数（过低=空液滴/低质量，过高=doublet）
#   - nCount_RNA:   每个细胞检测到的分子总数（UMI总数）
#   - percent.mt:   线粒体基因占比（过高=低质量/凋亡细胞）

cat("\n========== 2. 质量控制 ==========\n")

# 计算线粒体基因表达比例（基因名以 MT- 开头）
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 查看 QC 统计
cat("\n--- QC 指标（前5个细胞）---\n")
print(head(pbmc@meta.data, 5))

# 可视化 QC 指标（小提琴图）
p_qc_vln <- VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ggsave(file.path(output_dir, "01_qc_violin.png"), plot = p_qc_vln, width = 12, height = 5, dpi = 150)
cat("已保存: 01_qc_violin.png\n")

# FeatureScatter: 查看指标间关系
p_scatter1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
p_scatter2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p_qc_scatter <- p_scatter1 + p_scatter2
ggsave(file.path(output_dir, "02_qc_scatter.png"), plot = p_qc_scatter, width = 12, height = 5, dpi = 150)
cat("已保存: 02_qc_scatter.png\n")

# 过滤细胞:
#   nFeature_RNA > 200  （已有 min.features=200 过滤，这里二次确认）
#   nFeature_RNA < 2500 （去除 doublet）
#   percent.mt < 5      （去除线粒体污染严重的细胞）
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
cat("过滤后细胞数:", ncol(pbmc), "\n")


# ============================================================================
# 3. 数据归一化（Normalization）
# ============================================================================
# LogNormalize: 每个细胞的表达量除以该细胞总UMI，乘以10000，再取log
# 结果存储在 pbmc[["RNA"]]$data

cat("\n========== 3. 归一化 ==========\n")
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
cat("归一化完成\n")


# ============================================================================
# 4. 识别高变基因（Feature Selection）
# ============================================================================
# FindVariableFeatures: 识别在不同细胞间表达变异最大的基因
# 默认选择 2000 个高变基因，用于下游 PCA 降维

cat("\n========== 4. 高变基因识别 ==========\n")
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 查看前10个高变基因
top10 <- head(VariableFeatures(pbmc), 10)
cat("Top 10 高变基因:", paste(top10, collapse = ", "), "\n")

# 可视化高变基因
p_var <- VariableFeaturePlot(pbmc)
p_var_label <- LabelPoints(plot = p_var, points = top10, repel = TRUE)
ggsave(file.path(output_dir, "03_variable_features.png"), plot = p_var + p_var_label, width = 10, height = 7, dpi = 150)
cat("已保存: 03_variable_features.png\n")


# ============================================================================
# 5. 数据标准化（Scaling）
# ============================================================================
# ScaleData: 对每个基因做 Z-score 标准化（均值=0, 方差=1）
# 使下游 PCA 分析中各基因权重一致
# 结果存储在 pbmc[["RNA"]]$scale.data

cat("\n========== 5. 标准化 ==========\n")
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
cat("标准化完成（", length(all.genes), "个基因）\n")


# ============================================================================
# 6. 线性降维（PCA）
# ============================================================================
# 对标准化后的高变基因做 PCA
# 每个 PC 代表一组相关基因的「元特征」

cat("\n========== 6. PCA 降维 ==========\n")
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# 查看前5个 PC 的 top 基因
cat("\n--- 前5个主成分的标记基因---\n")
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

# 可视化 PCA
p_pca <- DimPlot(pbmc, reduction = "pca") + NoLegend()
ggsave(file.path(output_dir, "04_pca.png"), plot = p_pca, width = 7, height = 6, dpi = 150)
cat("已保存: 04_pca.png\n")

# PCA 负载图（前2个 PC）
p_pca_load <- VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
ggsave(file.path(output_dir, "05_pca_loadings.png"), plot = p_pca_load, width = 10, height = 8, dpi = 150)
cat("已保存: 05_pca_loadings.png\n")

# PCA 热图（前15个 PC）
p_pca_heat <- DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
ggsave(file.path(output_dir, "06_pca_heatmap.png"), plot = p_pca_heat, width = 12, height = 15, dpi = 150)
cat("已保存: 06_pca_heatmap.png\n")


# ============================================================================
# 7. 确定数据维度
# ============================================================================
# 选择用于下游分析的主成分数量
# ElbowPlot: 按方差解释比例排序各 PC，拐点处即为建议维度

cat("\n========== 7. 确定维度 ==========\n")
p_elbow <- ElbowPlot(pbmc)
ggsave(file.path(output_dir, "07_elbow_plot.png"), plot = p_elbow, width = 7, height = 5, dpi = 150)
cat("已保存: 07_elbow_plot.png\n")
cat("拐点约在 PC 9-10，本教程选取 10 个 PC\n")


# ============================================================================
# 8. 细胞聚类（Clustering）
# ============================================================================
# FindNeighbors: 基于 PCA 空间构建 KNN 图（K近邻图）
# FindClusters:  使用 Louvain 算法做模块化优化聚类
#   resolution: 聚类粒度，0.4-1.2 适合约 3K 细胞的数据集

cat("\n========== 8. 聚类 ==========\n")
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# 查看聚类结果
cat("\n--- 前5个细胞的聚类 ID ---\n")
print(head(Idents(pbmc), 5))
cat("聚类数量:", length(levels(pbmc)), "\n")


# ============================================================================
# 9. 非线性降维（UMAP）
# ============================================================================
# UMAP: 将高维数据投射到 2D 空间用于可视化
# 注意: UMAP 用于可视化探索，不应仅凭图形下生物学结论

cat("\n========== 9. UMAP ==========\n")
pbmc <- RunUMAP(pbmc, dims = 1:10)

p_umap <- DimPlot(pbmc, reduction = "umap")
ggsave(file.path(output_dir, "08_umap.png"), plot = p_umap, width = 8, height = 7, dpi = 150)
cat("已保存: 08_umap.png\n")

# 保存中间结果（可重新加载，无需重跑前面的步骤）
saveRDS(pbmc, file.path(output_dir, "pbmc_tutorial.rds"))
cat("已保存中间结果: pbmc_tutorial.rds\n")


# ============================================================================
# 10. 寻找差异表达基因（Marker Genes）
# ============================================================================
# FindAllMarkers: 为每个聚类寻找特异性标记基因（与其他所有细胞比较）
# presto 包加速 DE 分析（Seurat v5 自动调用）

cat("\n========== 10. Marker 基因 ==========\n")
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)

# 筛选 log2FC > 1 的显著 marker
pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) -> top_markers

cat("\n--- 各聚类 marker 数量 ---\n")
print(table(top_markers$cluster))

# 保存 marker 表
write.csv(pbmc.markers, file.path(output_dir, "all_markers.csv"), row.names = FALSE)
cat("已保存: all_markers.csv\n")

# 每个聚类取 top 10 基因做热图
top10_all <- pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup()

p_heatmap <- DoHeatmap(pbmc, features = top10_all$gene) + NoLegend()
ggsave(file.path(output_dir, "09_marker_heatmap.png"), plot = p_heatmap, width = 12, height = 15, dpi = 150)
cat("已保存: 09_marker_heatmap.png\n")

# FeaturePlot: 在 UMAP 上展示经典 marker 基因表达
p_feature <- FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14",
                                              "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
ggsave(file.path(output_dir, "10_feature_plot.png"), plot = p_feature, width = 14, height = 12, dpi = 150)
cat("已保存: 10_feature_plot.png\n")

# VlnPlot: 小提琴图展示特定 marker
p_vln <- VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
ggsave(file.path(output_dir, "11_violin_markers.png"), plot = p_vln, width = 8, height = 5, dpi = 150)
cat("已保存: 11_violin_markers.png\n")


# ============================================================================
# 11. 细胞类型注释
# ============================================================================
# 根据已知 canonical markers 将聚类映射到细胞类型
#
# 聚类  Marker 基因          细胞类型
# 0    IL7R, CCR7           Naive CD4+ T
# 1    CD14, LYZ            CD14+ Mono（单核细胞）
# 2    IL7R, S100A4         Memory CD4+ T
# 3    MS4A1                 B（B细胞）
# 4    CD8A                  CD8+ T
# 5    FCGR3A, MS4A7         FCGR3A+ Mono
# 6    GNLY, NKG7            NK（自然杀伤细胞）
# 7    FCER1A, CST3          DC（树突状细胞）
# 8    PPBP                  Platelet（血小板）

cat("\n========== 11. 细胞类型注释 ==========\n")

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B",
                      "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)

# 带标签的 UMAP 图
p_umap_labeled <- DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
ggsave(file.path(output_dir, "12_umap_labeled.png"), plot = p_umap_labeled, width = 8, height = 7, dpi = 150)
cat("已保存: 12_umap_labeled.png\n")

# 保存最终结果
saveRDS(pbmc, file.path(output_dir, "pbmc3k_final.rds"))
cat("已保存最终结果: pbmc3k_final.rds\n")


# ============================================================================
# 12. 总结
# ============================================================================

cat("\n")
cat("=======================================================\n")
cat("  PBMC 3k 标准分析流程完成!\n")
cat("=======================================================\n")
cat("  输入细胞数: 2,700\n")
cat("  过滤后细胞数:", ncol(pbmc), "\n")
cat("  聚类数:", length(levels(pbmc)), "\n")
cat("  细胞类型: Naive CD4 T / CD14+ Mono / Memory CD4 T / B /\n")
cat("           CD8 T / FCGR3A+ Mono / NK / DC / Platelet\n")
cat("  输出目录:", output_dir, "\n")
cat("  输出文件:\n")
cat("    01_qc_violin.png       - QC小提琴图\n")
cat("    02_qc_scatter.png      - QC散点图\n")
cat("    03_variable_features.png - 高变基因图\n")
cat("    04_pca.png             - PCA散点图\n")
cat("    05_pca_loadings.png    - PCA负载图\n")
cat("    06_pca_heatmap.png     - PCA热图\n")
cat("    07_elbow_plot.png      - 拐点图\n")
cat("    08_umap.png            - UMAP聚类图\n")
cat("    09_marker_heatmap.png  - Marker热图\n")
cat("    10_feature_plot.png    - Feature图\n")
cat("    11_violin_markers.png  - Marker小提琴图\n")
cat("    12_umap_labeled.png    - 注释后UMAP图\n")
cat("    all_markers.csv        - 全部Marker表\n")
cat("    pbmc_tutorial.rds      - 中间结果\n")
cat("    pbmc3k_final.rds       - 最终结果\n")
cat("=======================================================\n")

8. 下一步学习 ¶

阶段	教程	技能目标
阶段2	Visualization Vignette	掌握专业级图表制作
阶段2	SCTransform	替代传统归一化的先进方法
阶段3	Integration / Mapping	跨批次整合 + 自动注释
阶段4	WNN / ATAC / Mixscape	多模态与扰动分析
阶段5	Sketch / BPCells	百万级细胞分析
阶段6	Visium HD / 空间组学	空间转录组全链路

本教程基于 Seurat v5.5.1 官方文档整理，数据来自 10X Genomics 公开数据集。
R 4.6.1 + Seurat 5.5.1 + Windows 环境

阅读全文

20260910 144229 蛋白质结构预测与虚拟筛选全流程实现详解

2026年09月10日

蛋白质结构预测与虚拟筛选全流程实现详解 ¶

AI生成 本文档由星辰超级智能体（TeleAgent）基于 bioSkills 技能库自动生成，内容来源于以下技能的 SKILL.md 文件，经组织编排为端到端工作流教程。

场景假设 ¶

你是一名药物化学研究员，收到一个新靶点蛋白——人源 BRD4 溴区结构域（BRD4-BD1, UniProt: O60885），需要：

获取该蛋白的三维结构（优先使用实验结构，无则用 AlphaFold 预测）
评估结构质量，判断哪些区域可信
制备受体结构（加氢、质子化、修补缺失原子）
检测结合口袋（如果是 apo 结构）
准备配体库（从 SMILES 生成 3D 构象）
进行层次化虚拟筛选（Vina → GNINA 重打分）
对 top hits 进行 PoseBusters 物理合理性验证
ADMET 性质预测与类药性筛选
对最终候选化合物进行自由能计算（FEP）精排序

全流程总览 ¶

靶点蛋白序列/PDB ID
        │
        ▼
┌──────────────────────────┐
│ Step 1: 结构获取           │  技能: bio-structural-biology-alphafold-predictions
│ · PDB 实验结构 → 直接使用  │
│ · 无实验结构 → AlphaFold DB│
│ · 需要复合物 → AF3/Boltz   │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 2: 结构质量评估       │  技能: bio-structural-biology-alphafold-predictions
│ · pLDDT 分段读取           │
│ · PAE 域间置信度           │
│ · 决策：AFDB vs 新预测     │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 3: 受体结构制备       │  技能: bio-structural-biology-structure-preparation
│ · PDBFixer 补缺失原子     │
│ · PROPKA 预测 pKa         │
│ · reduce 优化翻转/互变异构  │
│ · PDB2PQR 生成 PDBQT      │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 4: 结合口袋检测       │  技能: bio-structural-biology-binding-site-detection
│ · fpocket 几何枚举         │
│ · P2Rank ML 配体性评分     │
│ · 交叉验证，确定 docking 框 │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 5: 配体库制备         │  技能: bio-virtual-screening (ligand prep)
│ · RDKit SMILES → 3D       │
│ · ETKDGv3 构象嵌入         │
│ · MMFF94 能量最小化        │
│ · meeko → PDBQT           │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 6: 层次化虚拟筛选     │  技能: bio-virtual-screening
│ · Stage 1: Vina 筛全库     │
│ · Stage 2: GNINA CNN 重排  │
│ · Stage 3: top-1000 候选   │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 7: Pose 物理验证      │  技能: bio-pose-validation
│ · PoseBusters 全检         │
│ · 应变能量化               │
│ · 芳香环平面度             │
│ · 过滤 PB-invalid 构象     │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 8: ADMET 预测         │  技能: bio-admet-prediction
│ · ADMETlab 3.0 批量预测    │
│ · Lipinski/Veber 类药性    │
│ · hERG/CYP/AMES 毒性       │
│ · PAINS 结构警报            │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 9: 自由能精排序       │  技能: bio-free-energy-calculations
│ · RBFE 热力学循环          │
│ · MBAR 分析                │
│ · 闭环验证                 │
│ · 或 Boltz-2 快速替代      │
└──────────────────────────┘

Step 1: 结构获取 ¶

技能：`bio-structural-biology-alphafold-predictions`¶

目标 ¶

获取 BRD4-BD1 的三维结构坐标文件。优先策略：PDB 有实验结构则直接用；无则从 AlphaFold DB 获取预测模型。

决策表：AFDB 使用 vs 重新预测 ¶

情况	用 AFDB	重新预测
单条 UniProt 单体，折叠层面问题	是	否
需要复合物/界面	否（仅单体）	是 — AF3 / AF-Multimer
需要配体/离子/辅因子/PTM	否	是 — AF3 或对接到实验结构
设计突变体不在 UniProt 中	否	是
需要特定构象状态	否（单一状态）	是（子采样 MSA）

BRD4-BD1 决策：PDB 中有多个 BRD4-BD1 实验结构（如 3MXF, 4LYI），优先使用实验 holo 结构（含已知配体）。若需 AlphaFold 模型作为补充：

代码：从 AlphaFold DB 获取模型 ¶

import requests

def afdb_metadata(accession):
    """查询 AFDB 元数据接口，返回当前下载 URL"""
    url = f'https://alphafold.ebi.ac.uk/api/prediction/{accession}'
    r = requests.get(url)
    r.raise_for_status()
    return r.json()  # list；每个片段/异构体一个对象

def fetch_afdb(accession, out_dir='.'):
    """下载坐标文件和 PAE 矩阵"""
    entries = afdb_metadata(accession)
    if not entries:
        return None  # >2700aa 非人源蛋白或非 UniProt 序列常缺失
    entry = entries[0]
    cif_text = requests.get(entry['cifUrl']).text
    pae_json = requests.get(entry['paeDocUrl']).json()
    cif_path = f"{out_dir}/AF-{accession}.cif"
    with open(cif_path, 'w') as f:
        f.write(cif_text)
    return cif_path, pae_json

# BRD4 full-length UniProt: O60885
result = fetch_afdb('O60885')

关键注意点 ¶

陷阱	原因	正确做法
版本后缀 URL 404	AFDB 版本号会变（v4→v6）	从 `/api/prediction/{accession}` 动态获取 URL，不要手工拼
长蛋白只返回部分残基	人源 >2700aa 的蛋白会被分段（F1, F2…）	遍历所有 metadata 条目，不同片段间的相对位置不可信
AFDB 模型无配体/离子/辅因子	AFDB 是真空中单链预测	需要 holo 结构时用实验 PDB；需要复合物时用 AF3/Boltz
404/空列表	非人源大蛋白或非 UniProt 序列	重新预测或查 ESM Metagenomic Atlas

Step 2: 结构质量评估 ¶

技能：`bio-structural-biology-alphafold-predictions`¶

目标 ¶

读取 AlphaFold 模型的 pLDDT 和 PAE，判断哪些区域可信、哪些需裁剪。

pLDDT 置信度分级 ¶

pLDDT	等级	操作含义	可信用于
> 90	Very high	骨架和侧链方向均可信	侧链细节、催化几何假设、MR 核心
70-90	Confident	骨架基本正确	折叠、域拓扑、骨架级 MR；侧链需谨慎
50-70	Low	骨架不确定	仅粗略拓扑；不可信细节
< 50	Very low	占位符；通常是固有无序区	无序信号，不是构象

代码：提取 pLDDT 并分级 ¶

from Bio.PDB import MMCIFParser

def extract_plddt(cif_file):
    """从 B-factor 列读取 pLDDT（注意：与热运动 B-factor 极性相反）"""
    parser = MMCIFParser(QUIET=True)
    structure = parser.get_structure('af', cif_file)
    plddt = {}
    for residue in structure[0].get_residues():
        if residue.id[0] == ' ' and 'CA' in residue:
            # pLDDT 存在 B-factor 列中，但是置信度（0-100，越高越好）
            # 与热运动 B-factor 极性相反——切勿当作热运动读取
            plddt[residue.id[1]] = residue['CA'].get_bfactor()
    return plddt

def plddt_band(score):
    if score > 90:
        return 'very_high'
    if score >= 70:
        return 'confident'
    if score >= 50:
        return 'low'
    return 'very_low'   # <50: 通常是无序区（IDR），不是建模错误

plddt = extract_plddt('AF-O60885.cif')
mean_plddt = sum(plddt.values()) / len(plddt)
disordered = [res for res, s in plddt.items() if s < 50]  # 候选 IDR

代码：读取 PAE 判断域间可靠性 ¶

import numpy as np

def load_pae(pae_json):
    entry = pae_json[0] if isinstance(pae_json, list) else pae_json
    # 紧凑格式（2023+）: 2D num_res × num_res 数组
    # 遗留 1D 'distance'/'residue1'/'residue2' 字段已移除
    return np.array(entry['predicted_aligned_error'])

def interdomain_confidence(pae, domain_a, domain_b):
    """PAE 不对称，取两个 off-diagonal block 的平均"""
    block = np.concatenate([pae[np.ix_(domain_a, domain_b)].ravel(),
                            pae[np.ix_(domain_b, domain_a)].ravel()])
    return block.mean()  # 低（约 <5Å）= 域间相对位置可信

pae = load_pae(result[1])

核心原则 ¶

pLDDT 是局部逐残基置信度，高均值 pLDDT 不保证域间相对位置正确。 两个域各 95 pLDDT 仍可能以未知相对方向飘浮——读 PAE 才能判断域间几何。

低 pLDDT 长段通常是真实的固有无序区（IDR），不是建模失败。 裁剪掉它会丢失真实生物学信息。pLDDT 与专门的 IDR 预测器具有竞争力（Akdel 2022, Piovesan 2022）。

关键陷阱 ¶

症状	原因	修复
模型”按柔度着色”看起来反了	B-factor 列被当作热运动读	它是 pLDDT（置信度，越高越好）；按 pLDDT band 着色
删掉低 pLDDT 尾巴后已知 IDR 消失	把低 pLDDT 当”错误”	低 pLDDT 通常= 无序；保留并注释为 IDR
域间距离无意义但 pLDDT 都很高	信任了拼接卡通图	pLDDT 是局部的；读 PAE——高 off-diagonal PAE = 域间无约束
对接打分很好但验证差	对接到 apo AFDB 口袋	口袋缺配体/辅因子，rotamer 不可靠；用 holo 结构或柔性对接

Step 3: 受体结构制备 ¶

技能：`bio-structural-biology-structure-preparation`¶

目标 ¶

将原始 PDB/mmCIF 转化为加氢完全、质子化正确、缺失原子补齐的 docking-ready 受体。

核心原则 ¶

沉积或预测的结构不是模拟就绪的。 每个添加的原子、每个放置的质子、每段补建的 loop 都是假设而非测量。制备条件（pH、质子化模型、补建内容）必须随文件一起记录，否则下游 H-bond、盐桥、界面和对接结果不可复现。

工具选择决策表 ¶

任务	工具	最佳场景	失败/误导场景
补缺失重原子、短内部 loop	PDBFixer	截短侧链、1-3 残基缺口	长/末端缺口、域级缺失
加 pH 氢、替换非标准残基	PDBFixer `addMissingHydrogens(pH)`	快速一键制备	忽略微环境 pKa 偏移
优化 Asn/Gln/His 翻转和 His 互变异构	reduce `-build`	解决酰胺/环方向歧义	当作 pKa 预测器——它优化几何而非滴定
预测残基 pKa / pH 依赖质子化	PROPKA (`propka3`)	判断哪些可滴定残基偏离标准态	当作精确实验 pKa
指定状态 + 写 PQR	PDB2PQR `--with-ph --ff --titration-state-method propka`	APBS/泊松-玻尔兹曼静电	力场不匹配时用作通用加氢
力场一致氢 + 显式变体	OpenMM `Modeller.addHydrogens`	MD 系统构建	需要非标准态时变体留默认

质子化决策表：标准态 vs pKa 预测器 ¶

残基/情形	pH 7 标准态通常可用	需 pKa 预测器 + H-bond 优化
表面 Asp/Glu	去质子化 (-1)	仅当靠近金属或 H-bond 伙伴
埋藏或盐桥 Asp/Glu	—	pKa 可升数个单位 → 可能中性/质子化
表面 Lys/Arg	质子化 (+1)	埋藏 Lys 靠近酸性残基可中性
任意位置 His	默认歧义	几乎总是：按局部 H-bond 和金属配位选 HID/HIE/HIP
游离 Cys	中性硫醇	催化/金属配位 Cys 常为硫醇盐
二硫键 Cys	S 上无 H (CYX)	先检测 SS 键；不要质子化

一句话规则：标准态仅对体相溶剂环境中的残基可辩护；任何埋藏、电荷簇集、金属邻接或口袋内的可滴定残基都需要预测器（PROPKA/H++）和 H-bond 网络 pass（reduce）。

代码：PDBFixer 制备管线 ¶

from pdbfixer import PDBFixer
from openmm.app import PDBFile

# 注意：finder 调用顺序很重要——残基 → 非标准 → 原子
# hydrogens 必须最后加（重原子存在后才能加氢）
fixer = PDBFixer(filename='receptor.pdb')

fixer.findMissingResidues()           # 短内部缺口 + 末端
fixer.findNonstandardResidues()       # 如 MSE（硒代蛋氨酸）
fixer.replaceNonstandardResidues()    # 映射回标准母体
fixer.removeHeterogens(keepWater=False)  # 去缓冲离子/冷冻保护剂
fixer.findMissingAtoms()              # 截短侧链 + 上面找到的残基
fixer.addMissingAtoms()               # 补重原子；补建的 loop 是假设，需记录

# pH 7.0 是一个选择，不是安全默认——需声明并匹配实验/生物条件
# addMissingHydrogens 检测二硫键并将那些 Cys 留为 CYX
fixer.addMissingHydrogens(pH=7.0)

with open('receptor_prepared.pdb', 'w') as out:
    PDBFile.writeFile(fixer.topology, fixer.positions, out, keepIds=True)

# 记录 provenance
built = [(ci, pos, name) for (ci, pos), names in fixer.missingResidues.items() for name in names]
print(f'built {len(fixer.missingResidues)} missing-residue segment(s), {len(built)} residue(s); pH=7.0')

代码：reduce 优化翻转和互变异构 ¶

import subprocess

# -build 执行 -OH -ROTEXOH -HIS -FLIP：加氢并优化 OH/His 旋转 + Asn/Gln/His 翻转
# reduce 通过小探针全原子接触评分优化方向（Word 1999）
with open('receptor_reduced.pdb', 'w') as out:
    subprocess.run(['reduce', '-build', 'receptor_prepared.pdb'], stdout=out, check=True)

# Reduce2（CCTBX/Phenix）是维护中的继承者：mmtbx.reduce2 receptor_prepared.pdb

代码：PROPKA pKa 预测 + PDB2PQR 生成 PQR/PDBQT ¶

import subprocess

# propka3 写 receptor_prepared.pka；SUMMARY 列出预测 pKa vs 模型（标准）pKa
subprocess.run(['propka3', 'receptor_prepared.pdb'], check=True)

# PDB2PQR 通过 PROPKA 在 --with-ph 指定 pH 下指定质子化态，写出电荷/半径
# --ff 必须匹配下游力场；--with-ph 7.0 是声明的滴定条件
subprocess.run([
    'pdb2pqr30', '--ff=AMBER', '--with-ph', '7.0',
    '--titration-state-method', 'propka', '--keep-chain',
    'receptor_prepared.pdb', 'receptor.pqr',
], check=True)

# 再用 meeko mk_prepare_receptor.py 转 PDBQT
subprocess.run([
    'mk_prepare_receptor.py', '--read_pqr', 'receptor.pqr',
    '-o', 'receptor', '-p'
], check=True)

预测模型制备的特殊考虑 ¶

AlphaFold/ESMFold 模型制备时必须先裁剪低 pLDDT 区域，再加氢/质子化。口袋 rotamer 即使在骨架 pLDDT 高时也是最不可靠的原子。

# 1) 先裁剪低 pLDDT 残基（<50-70 = 不可靠）
#    phenix.process_predicted_model 做这个 + PAE 域分割
# 2) 然后对裁剪后模型运行 PDBFixer 管线
fixer = PDBFixer(filename='af_model_trimmed.pdb')
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(pH=7.0)  # 预测模型从不携带实验氢

关键陷阱 ¶

症状	原因	修复
制备文件没有氢	X-ray 模型极少分辨氢	显式调用 `addMissingHydrogens(pH=...)`
His H-bond/金属配位错	默认 HIE/HID 猜测未考虑局部网络	让 reduce 选中性互变异构（HID vs HIE）；PROPKA 仅决定电荷态（HIP vs 中性）
埋藏 Asp/Glu 去质子化但应为中性	假设标准 pKa 7	预测 pKa（PROPKA/H++）；质子化该残基（ASH/GLH）
`addMissingAtoms` 建出离谱 loop	长/末端缺口交给只做短 loop 的补全器	不在这里补长缺口；用 loop/同源建模器并标记
催化金属/辅因子被去掉	`removeHeterogens()` 清掉所有非水 heterogen	有目的地保留：显式过滤，不要批量移除
下游结果不可复现	制备文件未记录假设	记录 pH、质子化模型、每个补建原子/loop 作为 provenance

Step 4: 结合口袋检测 ¶

技能：`bio-structural-biology-binding-site-detection`¶

目标 ¶

在 apo 结构上检测候选口袋并评分，确定 docking 搜索框的中心坐标和大小。

核心原则 ¶

几何凹陷不自动等于功能或可成药位点。 fpocket 的 alpha 球或 CASTp 的表面拓扑会找到表面所有凹陷——包括结晶添加剂位点、非功能表面口袋和域间裂缝。口袋排名是假设而非结合位点。最可靠的位点是已被配体占据的位点（holo 结构）。

可成药性评分基于 holo 训练集，系统性低估 apo、浅口袋、别构位点和 PPI 口袋。 隐口袋在 apo 快照中完全不可见——仅在 MD 或构象系综中出现。

方法选择决策表 ¶

问题/情况	方法	最佳场景	失败/误导场景
快速枚举所有几何凹陷	fpocket 或 CASTp	apo 或 holo 全表面扫描	报告添加剂位点和非功能口袋
按配体性/可成药性排名	P2Rank, fpocket 药性分, DoGSiteScorer	球状结构域上多口袋优先排序	训练于 holo/可成药集 → 低估 apo/shallow/allosteric
检测隐口袋/瞬态口袋	mdpocket over MD trajectory	位点仅在结合时打开	单一静态结构——隐口袋不可见

代码：fpocket 检测和排名 ¶

import subprocess
from pathlib import Path

def run_fpocket(pdb_path):
    subprocess.run(['fpocket', '-f', pdb_path], check=True)
    stem = Path(pdb_path).stem
    return Path(pdb_path).with_name(f'{stem}_out')

def parse_fpocket_info(out_dir):
    info = next(Path(out_dir).glob('*_info.txt'))
    pockets, current = [], None
    for raw in info.read_text().splitlines():
        line = raw.strip()
        if line.startswith('Pocket'):
            current = {'pocket': int(line.split()[1])}
            pockets.append(current)
        elif ':' in line and current is not None:
            key, val = (part.strip() for part in line.split(':', 1))
            try:
                current[key] = float(val)
            except ValueError:
                pass
    return pockets

out_dir = run_fpocket('receptor_prepared.pdb')
pockets = sorted(parse_fpocket_info(out_dir),
                 key=lambda p: p.get('Druggability Score', 0), reverse=True)
# 0.5 药性分是分诊线，不是定律
druggable = [p for p in pockets if p.get('Druggability Score', 0) >= 0.5]
for p in pockets[:5]:
    print(f"pocket {p['pocket']}: drug={p.get('Druggability Score', 0):.2f}")

代码：P2Rank ML 配体性评分 ¶

import csv
import subprocess
from pathlib import Path

def run_p2rank(pdb_path, out_dir='p2rank_out'):
    subprocess.run(['prank', 'predict', '-f', pdb_path, '-o', out_dir], check=True)
    return Path(out_dir) / f'{Path(pdb_path).name}_predictions.csv'

def parse_p2rank(csv_path):
    with open(csv_path) as fh:
        reader = csv.DictReader(fh, skipinitialspace=True)  # P2Rank 列名有空格
        rows = [{k.strip(): v.strip() for k, v in row.items()} for row in reader]
    return [{'rank': int(r['rank']), 'score': float(r['score']),
             'probability': float(r['probability'])} for r in rows]

predictions = parse_p2rank(run_p2rank('receptor_prepared.pdb'))
top = predictions[0]
print(f"top pocket: probability={top['probability']:.2f} score={top['score']:.2f}")

正交验证信号 ¶

口袋检测排名需要以下信号交叉验证：

独立检测器一致性：fpocket 和 P2Rank 是否排名同一个口袋为 top
口袋内衬残基的进化保守性（ConSurf）：功能位点通常保守，添加剂/晶体裂缝不保守
热点/混合溶剂作图（FTMap, mixed-solvent MD）
同源 holo 结构
已知生物学信息

关键陷阱 ¶

症状	原因	修复
排名第一的口袋是表面添加剂位点	几何找到所有凹陷	交叉检查口袋原子 vs 输入中的 HETATM/添加剂
已知结合位点完全漏掉	位点是隐口袋——apo 快照中关闭	运行 mdpocket over MD/NMR 系综
真实靶点药性分全低	评分训练于 holo/可成药集	低分不等于”无位点”；用独立方法佐证
P2Rank parse KeyError	CSV 列名有前导空格	用 `csv.DictReader(skipinitialspace=True)` 并 strip keys
AlphaFold 口袋对接差	口袋内衬 rotamer 最不可靠	先验证/制备模型；优先用 holo 实验结构

Step 5: 配体库制备 ¶

技能：`bio-virtual-screening`（配体制备部分）¶

目标 ¶

从 SMILES 生成 3D docking-ready 配体文件，包含正确的质子化和构象。

代码：RDKit + meeko 配体制备 ¶

from rdkit import Chem
from rdkit.Chem import AllChem
from meeko import MoleculePreparation, PDBQTWriterLegacy

def prepare_ligand(smiles, pdbqt_out):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f'invalid SMILES: {smiles}')
    mol = Chem.AddHs(mol)

    # ETKDGv3 3D 构象嵌入
    embed_status = AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
    if embed_status != 0:
        raise RuntimeError('ETKDGv3 failed to generate a ligand conformer')

    # MMFF94 能量最小化
    if not AllChem.MMFFHasAllMoleculeParams(mol):
        raise ValueError('MMFF94 parameters unavailable for this ligand')
    optimization_status = AllChem.MMFFOptimizeMolecule(mol)
    if optimization_status != 0:
        raise RuntimeError('MMFF94 optimization did not converge')

    # meeko 0.5+ API：prepare() 返回 MoleculeSetup 列表
    # 用 PDBQTWriterLegacy.write_string() 生成 PDBQT
    mk_prep = MoleculePreparation()
    setups = mk_prep.prepare(mol)
    pdbqt_text, is_ok, err = PDBQTWriterLegacy.write_string(setups[0])
    if not is_ok:
        raise RuntimeError(f'meeko PDBQT export failed: {err}')
    with open(pdbqt_out, 'w') as f:
        f.write(pdbqt_text)
    return pdbqt_out

# 批量制备
library = [
    'CC(=O)c1ccccc1',      # acetophenone
    'O=C(O)c1ccccc1O',     # salicylic acid
    # ... 更多 SMILES
]
for i, smi in enumerate(library):
    prepare_ligand(smi, f'ligand_{i:04d}.pdbqt')

关键注意点 ¶

陷阱	原因	修复
`MolFromSmiles` 不预测 pH 7.4 质子化	它只解析提供的 SMILES 状态	需先用 pKa/protomer 工具生成正确质子化态
`Uncharger` 不预测质子化	它仅中和形式电荷	用专门工具预测 pH 7.4 下的质子化
meeko 0.5+ API 变化	`prep.write_pdbqt_file()` 已废弃	用 `PDBQTWriterLegacy.write_string(setups[0])`
Histidine 质子化不一致	受体用 pH 7.4，配体用 pH 7.0	统一 pH，用 PROPKA 检查催化残基

Step 6: 层次化虚拟筛选 ¶

Docking 工具分类 ¶

工具	评分函数	最佳场景	失败场景
AutoDock Vina 1.2	Vina（经验）	开放、表征良好的基线	交叉对接；隐口袋；金属中心
SMINA	Vina + 灵活 + 自定义	自定义评分；柔性侧链	Vina 评分同样的注意事项
GNINA 1.1	CNN 或 Vina 评分	CNN 辅助 pose 排序	新化学型 OOD
DiffDock-L	扩散生成	交叉对接 pose 采样	需 PoseBusters 验证
Boltz-2 + GNINA rescore	基础模型 + CNN	实验性多模型工作流	独立基准每个证据流

场景决策树 ¶

场景	推荐工作流
自对接（已知配体口袋）	GNINA `--cnn_scoring rescore`
交叉对接（apo 或相关靶点结构）	DiffDock-L pose + GNINA rescore + PoseBusters
超大库（10M+）	校准层次筛选：性质/警报分诊 → Vina → GNINA → 高成本后续
别构/未定义位点	P2Rank 口袋检测 → 对所有口袋做系综对接
金属配位配体	GOLD（商业）或手动参数化 Vina 金属评分
Hit-to-lead 精修	用共晶结构；MD 弛豫受体；FEP 亲和力

代码：Vina 单配体对接 ¶

from vina import Vina

def dock_single(receptor_pdbqt, ligand_pdbqt, center, box_size,
                exhaustiveness=8, n_poses=10):
    """
    center: (x, y, z) 口袋中心坐标
    box_size: (x, y, z) 搜索框大小
    """
    v = Vina(sf_name='vina')
    v.set_receptor(receptor_pdbqt)
    v.set_ligand_from_file(ligand_pdbqt)
    v.compute_vina_maps(center=center, box_size=box_size)
    v.dock(exhaustiveness=exhaustiveness, n_poses=n_poses)
    return v.energies(), v.poses()

代码：GNINA CNN 重打分 ¶

gnina -r receptor.pdb -l ligand.sdf \
      --autobox_ligand reference_ligand.sdf \
      --cnn_scoring rescore \
      -o poses.sdf.gz \
      --num_modes 9 --exhaustiveness 8

--cnn_scoring 模式说明：

模式	行为	计算开销
`none`	无 CNN；全程经验评分	最低
`rescore`（默认）	经验搜索 + CNN 重排最终 pose	CNN 选项中最低
`refinement`	CNN 精修 pose 并排名	~10x rescore
`metrorescore`	CNN 在 Metropolis 搜索中评分	中等
`metrorefine`	CNN 在搜索和精修中均评分	高
`all`	全程 CNN	极高，不推荐

代码：层次化虚拟筛选管线 ¶

import pandas as pd
from concurrent.futures import ProcessPoolExecutor
from functools import partial

def vs_pipeline(library_smi, receptor_pdbqt, center, box, output_dir, n_workers=16):
    # Stage 0: 读取化合物库
    df = pd.read_csv(library_smi)

    # Stage 1: 类药性过滤（Lipinski + Veber + PAINS）
    df_stage1 = drug_like_filter(df)  # 委托给 admet-prediction 技能

    # Stage 2: Vina 对接全过滤后库
    worker = partial(vina_dock, receptor_pdbqt=receptor_pdbqt,
                     center=center, box=box)
    with ProcessPoolExecutor(max_workers=n_workers) as ex:
        affinities = list(ex.map(worker, df_stage1['smiles']))
    df_stage1['vina_affinity'] = affinities

    # 取 Vina top 1%
    df_stage2 = df_stage1.nsmallest(int(len(df_stage1) * 0.01), 'vina_affinity')

    # Stage 3: GNINA CNN 重打分 top 1%
    df_stage2['gnina_affinity'] = df_stage2['smiles'].apply(
        lambda smi: gnina_rescore(smi, receptor_pdbqt, center, box))

    # 取 GNINA top 1000
    df_stage3 = df_stage2.nsmallest(1000, 'gnina_affinity')

    # Stage 4: Pose 物理验证
    return pose_validate(df_stage3)  # 委托给 pose-validation 技能

Vina vs GNINA 分歧调解表 ¶

Vina top pose	GNINA top pose	行动
同一 pose，分数相似	同一 pose，分数相似	视为支持证据；仍需物理验证
不同 pose	同口袋，不同方向	保留两者，对照靶点相关对照比较
Vina 极好，GNINA 一般	不同 pose，分数差异大	检查两个 pose；不要仅从分数分歧推断谁对
两者都差	大量配体分数相似偏低	可能口袋/蛋白构象错；重新考虑受体

关键陷阱 ¶

症状	原因	修复
Vina segfault	PDBQT 损坏（原子名）	用 meeko 重新制备
GNINA 挂起	GPU OOM	减少并发；`--num_modes 5`
所有色预测结合能很差（-3 到 -5）	质子化错误；配体对框太大	重检 pKa；扩大框
配体间亲和力完全相同	受体 grid 未计算	dock 前调 `v.compute_vina_maps()`
GPU 模式慢	Vina 仅 CPU；仅 GNINA 用 GPU	GNINA 用 GPU；Vina 用多核 CPU

Step 7: Pose 物理验证 ¶

目标 ¶

用 PoseBusters 对 docked poses 进行物理合理性验证，过滤 PB-invalid 构象。

核心原则 ¶

Pose 准确性和化学合理性是不同的维度。 在 PoseBusters Benchmark 上，DiffDock 的 top-pose RMSD ≤ 2Å 成功率为 72%，但 RMSD+PB-valid 联合成功率仅 47%——25 个百分点的差距。不要将此百分比泛化到其他数据集或更新模型版本。

PoseBusters 检查组 ¶

检查组	测试内容	2024 基准标准
Sanity	配体化学合理性	RDKit sanitization 通过
Bond lengths	键长在参考范围内	0.75–1.25 × RDKit 距离几何下界
Bond angles	1-3 距离在参考范围内	0.75–1.25 × RDKit 下界
Internal steric	无配体内 clash	对距离 > 0.70 × RDKit 下界
Aromatic planarity	芳香环平面	最大偏离拟合平面 ≤ 0.25 Å
Double-bond stereo	Z/E 保持	匹配输入 SMILES
Internal energy	能量 vs 生成构象	UFF 能量比 ≤ 100 × 50 个生成弛豫构象均值
Volume overlap	与蛋白 vdW 重叠	< 7.5% 配体 vdW 体积
Minimum distance	无严重蛋白-配体 clash	距离 ≥ 0.75 × vdW 半径和
Chirality	R/S 保持	匹配输入 SMILES

通过所有检查 = “PB-valid”。PB-valid + RMSD ≤ 2Å = 现代标准。

代码：PoseBusters 验证 ¶

from posebusters import PoseBusters
import pandas as pd

bust = PoseBusters(config='dock')  # 'dock' = 全检除 RMSD 参考

results = bust.bust(
    mol_pred='/path/to/docked_poses.sdf',
    mol_cond='/path/to/receptor.pdb',
)

# 聚合所有 bool 列为 pb_valid 标志
check_cols = [
    col for col in results.select_dtypes(include='bool').columns
    if not col.lower().startswith('rmsd')
]
results['pb_valid'] = results[check_cols].all(axis=1)
valid = results[results['pb_valid']]
print(f'{len(valid)} / {len(results)} poses are PB-valid')

配置选择 ¶

配置	包含检查	使用场景
`redock`	全检 + RMSD vs 参考 + 蛋白 vdW 重叠	自对接基准、回溯验证
`dock`	全检除 RMSD 参考	盲对接、前瞻性虚拟筛选
`mol`	仅配体内（sanity, bonds, angles, rings, stereo, energy）	构象 QC；无蛋白上下文

代码：应变能量化 ¶

from rdkit import Chem
from rdkit.Chem import AllChem

def ligand_strain(docked_sdf, n_ref=20):
    """量化 docked pose 偏离最低能量自由构象的程度"""
    suppl = Chem.SDMolSupplier(docked_sdf, removeHs=False)
    strains = []
    for docked in suppl:
        if docked is None:
            continue
        smi = Chem.MolToSmiles(docked)
        ref = Chem.MolFromSmiles(smi)
        if ref is None:
            strains.append({'strain': None, 'note': 'reference_parse_failed'})
            continue
        ref = Chem.AddHs(ref)
        props_ref = AllChem.MMFFGetMoleculeProperties(ref)
        if props_ref is None:
            strains.append({'strain': None, 'note': 'no_mmff_params'})
            continue
        conf_ids = list(AllChem.EmbedMultipleConfs(
            ref, numConfs=n_ref, params=AllChem.ETKDGv3()))
        if not conf_ids:
            strains.append({'strain': None, 'note': 'embedding_failed'})
            continue
        AllChem.MMFFOptimizeMoleculeConfs(ref)

        ref_energies = []
        for c in conf_ids:
            ff = AllChem.MMFFGetMoleculeForceField(ref, props_ref, confId=c)
            if ff is not None:
                ref_energies.append(ff.CalcEnergy())
        min_ref = min(ref_energies) if ref_energies else None

        # 加 H 坐标，弛豫 H 同时固定重原子
        docked_h = Chem.AddHs(Chem.Mol(docked), addCoords=True)
        if docked_h.GetNumAtoms() != ref.GetNumAtoms():
            strains.append({'strain': None, 'note': 'atom_mismatch'})
            continue
        props_docked = AllChem.MMFFGetMoleculeProperties(docked_h)
        docked_ff = AllChem.MMFFGetMoleculeForceField(docked_h, props_docked) \
            if props_docked is not None else None
        if docked_ff is not None:
            for atom in docked_h.GetAtoms():
                if atom.GetAtomicNum() != 1:
                    docked_ff.AddFixedPoint(atom.GetIdx())
            docked_ff.Minimize(maxIts=200)
            docked_e = docked_ff.CalcEnergy()
            strains.append({
                'min_ref_energy': min_ref,
                'docked_energy': docked_e,
                'strain': docked_e - min_ref,
                'note': 'ok'
            })
        else:
            strains.append({'strain': None, 'note': 'docked_ff_failed'})
    return strains

RMSD vs PB-valid 调解表 ¶

RMSD ≤ 2Å	PB-valid	行动
是	是	物理合理且接近参考；仍需验证下游任务适用性
是	否	接近参考但某项合理性检查未过；检查失败项
否	是	物理合理但与参考不同；调查对齐、蛋白状态和替代结合模式
否	否	与参考不同且不合理；检查两方面再决定是否拒绝

关键陷阱 ¶

症状	原因	修复
期望的检查缺失	输入加载失败或配置省略	检查返回 DataFrame、加载状态列和输入格式
全检通过但 pose 无效	受体文件格式错	用带氢的 PDB；PDBQT 可能不工作
vdW 重叠假阳性（共价）	共价键被计为 clash	用共价对接专用验证
应变计算慢	参考构象太多	`n_ref` 减到 5-10

Step 8: ADMET 预测 ¶

目标 ¶

对 top hits 进行 ADMET 性质预测和类药性筛选，排除有毒性或药代性质差的化合物。

ADMET 模型分类 ¶

工具	端点数	架构	不确定性	访问方式	失败场景
ADMETlab 3.0	119 特征（77 模型 + 34 计算性质 + 8 规则）	多任务 DMPNN + 描述符	证据不确定性	Web 服务/ API	OOD；金属；大环
ADMET-AI	TDC 衍生	Chemprop D-MPNN	版本相关	Python 包	v2 与 v1/服务器不同
chemprop（自建）	用户定义	D-MPNN ± 描述符	集成	Python 包	需充足训练数据
ProTox-3.0	61 毒性端点	RF/DNN + 相似性	置信分	Web 服务	仅毒性

场景决策树 ¶

场景	工作流	理由
库分诊，无内部数据	ADMETlab 3.0 API 批量	广覆盖 + 模型不确定性
单端点，有内部数据	chemprop D-MPNN 基准	按部署相关验证选择
需校准概率	chemprop 集成 + Platt	DL 原生概率很少校准
心脏毒性倾向	ADMETlab hERG + ProTox-3.0 + 文献	多模型三角验证
药物-药物相互作用（CYP）	CYP1A2/2C9/2C19/2D6/3A4	标准 5 个 CYP

代码：Lipinski/Veber 类药性筛选 ¶

from rdkit.Chem import Descriptors, Lipinski, QED
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams

def druglike_score(mol):
    """计算类药性指标"""
    mw = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    tpsa = Descriptors.TPSA(mol)
    rotbonds = Lipinski.NumRotatableBonds(mol)
    qed = QED.qed(mol)

    lipinski_violations = sum([mw > 500, logp > 5, hbd > 5, hba > 10])
    veber_pass = rotbonds <= 10 and tpsa <= 140
    bbb_simple_screen = tpsa <= 90 and mw <= 500 and hbd <= 3  # BBB 启发式

    return {'MW': mw, 'LogP': logp, 'HBD': hbd, 'HBA': hba,
            'TPSA': tpsa, 'RotBonds': rotbonds, 'QED': qed,
            'Lipinski_violations': lipinski_violations,
            'Veber_pass': veber_pass,
            'BBB_simple_screen': bbb_simple_screen}

def alerts(mol, catalogs=('PAINS_A', 'BRENK', 'ZINC')):
    """PAINS/BRENK/REOS 结构警报"""
    params = FilterCatalogParams()
    for cat in catalogs:
        params.AddCatalog(getattr(FilterCatalogParams.FilterCatalogs, cat))
    catalog = FilterCatalog(params)
    hits = catalog.GetMatches(mol)
    return [h.GetDescription() for h in hits]

ADMETlab 3.0 端点样例 ¶

类别	端点示例
吸收	Caco-2 通透性、HIA (%)、Pgp 抑制剂/底物、MDCK
分布	BBB+、PPB (%)、VDss (L/kg)、Fu（游离分数）
代谢	CYP1A2/2C9/2C19/2D6/3A4 抑制剂/底物
排泄	CL (mL/min/kg)、T1/2 (h)
毒性	hERG、AMES、肝毒性(DILI)、致癌性、免疫毒性、致突变性
类药性	Lipinski, Veber, Ghose, Egan, Muegge, QED, SAscore

hERG 心脏毒性端点 ¶

模型	架构	训练数据	AUC
Shan et al. D-MPNN + MOE	D-MPNN + 206 MOE	7,889 化合物	0.956
CardioTox-net	五 DL + 神经元集成	BindingDB + ChEMBL	0.930
ADMETlab 3.0 hERG	DMPNN 多任务	内部	0.92
ProTox-3.0 心脏毒性	RF	5,252 ChEMBL	0.86 CV / 0.95 外部

单模型概率 > 0.5 不是 kill signal。 需多模型三角验证 + 文献检索 + 考虑暴露量 vs hERG 效力。

PAINS 作为 kill filter 的陷阱 ¶

PAINS 校准的是 HTS 分析干扰，不预测药物开发失败。将 PAINS 匹配作为排除标准会从库中清除有效先导物（如姜黄素类似物、多酚天然产物）。正确做法是标记 PAINS 命中以便正交分析确认，而非预先排除。

模型间调解 ¶

当 ADMETlab、ProTox-3.0 和 chemprop 在 hERG 上不一致：

情况	行动
全预测 hERG+	更高关注；计划体外膜片钳
结果分歧	检查适用域、活性阈值和测定定义再决定
全预测 hERG-	较低关注，但临床候选和新化学型仍考虑体外筛选
不要将相关模型当作独立证据	仅仅因为托管在不同服务上

OECD QSAR 五支柱（监管级）¶

明确定义端点 — 具体生物测定、单位、条件
明确算法 — 可复现模型 + 代码
定义适用域（AD） — 模型有效的范围
适当统计验证 — 外部测试集、交叉验证
机制解释 — 生物/化学原理

关键陷阱 ¶

症状	原因	修复
ADMETlab API 超时	服务负载/配额	遵循批量限制；退避重试
chemprop 过拟合	随机拆分	用骨架拆分 `--split scaffold_balanced`
hERG 预测 50/50	OOD	检查适用域（kNN 距离、集成方差）
QED 计算失败	分子缺失/未 sanitize	拒绝解析失败；sanitize 输入
预测跨运行不一致	chemprop 集成随机种子	`--seed 42` 并复用模型
校准不匹配	DL 原生概率未校准	在验证集上 Platt scaling

Step 9: 自由能精排序 ¶

技能：`bio-free-energy-calculations`¶

目标 ¶

对 top 候选化合物进行相对结合自由能（RBFE）计算，以物理方法精确排序类似物。

FEP 方法分类 ¶

方法	成本/对	精度	适用场景	失败场景
FEP+（Schrödinger）	系统/协议相关 GPU	商业 RBFE	商业先导优化	许可和复现约束
OpenFE RBFE	系统/协议相关 GPU	开源 RBFE	开源项目	映射/设置/采样需审查
GROMACS/AMBER RBFE	实现相关	自定义	专家控制	手动验证负担
MM/PBSA / MM/GBSA	低成本端点	近似端点分	系列内探索排序	熵、采样、模型依赖
Boltz-2 affinity	秒级 GPU	基准相关	ML 替代；≥1000x 更快	新化学型

场景决策树 ¶

场景	推荐工作流
排序近似物（R-group SAR）	RBFE via OpenFE（循环：lig1↔lig2↔lig3）
跨骨架排序	ABFE 每个配体；或协调 RBFE 星形网络
同源先导优化集	RBFE + 连通冗余感知扰动图
单配体亲和力	ABFE（无需参考化合物）
低成本探索排序	项目验证端点或 ML 方法 + 正交确认
选择性（靶点 vs 脱靶）	两个蛋白上 RBFE；报告 delta-delta-G
别构 vs 正构	ABFE 可比；用 MD 检查 pose 稳定性

RBFE 热力学循环 ¶

ΔΔG_binding = ΔG_lig1→lig2(bound) - ΔG_lig1→lig2(solvent)

代码：OpenFE RBFE 设置 ¶

from openfe import SmallMoleculeComponent, ProteinComponent, SolventComponent
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol

protein = ProteinComponent.from_pdb_file('receptor_prepared.pdb')
ligA = SmallMoleculeComponent.from_sdf_file('ligand_A.sdf')
ligB = SmallMoleculeComponent.from_sdf_file('ligand_B.sdf')
solvent = SolventComponent()

protocol = RelativeHybridTopologyProtocol(
    RelativeHybridTopologyProtocol.default_settings()
)
# 需创建/检查原子映射（Kartograf 是 OpenFE 1.7 CLI 默认映射器）
# 构建 bound 和 solvent Transformation 对象
# 通过 openfe quickrun 或 Python 执行接口运行

Lambda 窗口调度 ¶

阶段	Lambda 值	目的
解耦（vdW）	0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.0	关闭配体 vdW
充电（Coulomb）	0.0, 0.25, 0.5, 0.75, 1.0	关闭配体部分电荷
约束（仅 ABFE）	0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0	Boresch 式约束

代码：MBAR 分析 ¶

from alchemlyb import concat
from alchemlyb.parsing import gmx
from alchemlyb.estimators import MBAR
from alchemlyb.postprocessors.units import to_kcalmol

u_nks = []
for window in range(12):
    df = gmx.extract_u_nk(f'window_{window}.xvg', T=300)
    u_nks.append(df)

u_nk = concat(u_nks)
mbar = MBAR().fit(u_nk)
delta_g = to_kcalmol(mbar.delta_f_).iloc[0, -1]
d_delta_g = to_kcalmol(mbar.d_delta_f_).iloc[0, -1]
print(f'ΔG: {delta_g:.2f} ± {d_delta_g:.2f} kcal/mol')

代码：闭环验证 ¶

def cycle_closure_residual(cycle):
    """cycle = [(lig_i, lig_j, delta_g, sd), ...]"""
    total = sum(d_g for _, _, d_g, _ in cycle)
    total_var = sum(sd**2 for _, _, _, sd in cycle)
    return total, total_var ** 0.5

# 残差应与传播不确定性一致（在 ~2σ 内）

力场选择 ¶

力场	用于	备注
OPLS4（Schrödinger）	FEP+ 默认	商业；测试充分
OpenFF 2.1.1 (Sage)	OpenFE 1.7 默认	检查序列化设置
GAFF2	AMBER FEP	仅配体；蛋白 FF14SB
CGenFF	CHARMM 系 FEP	CHARMM 力场族

Boltz-2 快速替代 ¶

Boltz-2 报告在评估基准上亲和力精度接近 FEP，至少 1000 倍速度优势。但结果是基于特定基准的，不建立对任意 ChEMBL 数据的通用 RMSE 或相关性。新化学型或模态超出演示的训练/基准域时不应单独依赖。

何时用 Boltz-2：hit 发现分诊用 affinity_probability_binary；hit-to-lead 用 affinity_pred_value
何时不依赖 Boltz-2：新化学型/模态超出版本域，或无靶点相关验证的生产决策

关键陷阱 ¶

症状	原因	修复
Lambda 窗口模拟发散	初始 pose 差	先用 MM 最小化弛豫 pose
闭环残差与不确定性不一致	采样/映射/力场/相关边问题	检查有符号残差、重叠、映射和独立重复
MBAR 返回 NaN	窗口间重叠不足	添加中间 lambda 窗口
ABFE 系统偏移	约束/标准态/采样/力场问题	检查协议文档化的约束/自由能项和符号
MM/GBSA 排序与对接不符	不同轨迹帧	在 MD 弛豫 pose 上计算 MM/GBSA
Boltz-2 亲和力全 0	输入格式错	检查 SMILES 有效性；先标准化

全流程总结图 ¶

靶点蛋白 BRD4-BD1 (O60885)
    │
    ├─ PDB 实验结构 (h3MXF) ──────────────────────────────┐
    │                                                      │
    └─ AlphaFold DB (O60885) ─┐                           │
         │                     │                           │
         ▼                     │                           │
    pLDDT/PAE 评估              │                           │
    裁剪低置信区域               │                           │
         │                     │                           │
         ▼                     ▼                           ▼
    ┌────────────────────────────────────────────────────────┐
    │           受体结构制备 (PDBFixer + PROPKA + reduce)     │
    │           补缺失原子 → 加氢 → 质子化 → PDBQT            │
    └────────────────────────┬───────────────────────────────┘
                             │
    ┌────────────────────────▼───────────────────────────────┐
    │  口袋检测 (fpocket + P2Rank)  ← 如果是 apo             │
    │  或从 holo 共晶配体中心定义搜索框                       │
    └────────────────────────┬───────────────────────────────┘
                             │
                  ┌──────────┴──────────┐
                  │                     │
            配体库制备               受体 PDBQT
          (RDKit + meeko)              │
                  │                     │
                  └──────────┬──────────┘
                             │
    ┌────────────────────────▼───────────────────────────────┐
    │  Stage 1: 类药性过滤 (Lipinski + Veber + PAINS)        │
    │  Stage 2: Vina 对接 (exhaustiveness=8, 并行)           │
    │  Stage 3: GNINA CNN rescore top 1%                    │
    │  Stage 4: PoseBusters 物理验证 top 1000               │
    │  Stage 5: ADMET 预测 (ADMETlab 3.0 + hERG + CYP)      │
    └────────────────────────┬───────────────────────────────┘
                             │
                    top ~50 候选
                             │
    ┌────────────────────────▼───────────────────────────────┐
    │  自由能精排序                                          │
    │  · 同系列: RBFE (OpenFE, MBAR, 闭环验证)              │
    │  · 跨骨架: ABFE 或 Boltz-2 快速分诊                    │
    └────────────────────────┬───────────────────────────────┘
                             │
                    top 5-10 候选
                             │
                    ▼
              实验验证 (SPR / ITC)

涉及技能一览表 ¶

步骤	技能名称	主要工具	核心功能
Step 1	bio-structural-biology-alphafold-predictions	requests, Bio.PDB	AFDB 模型获取与置信度解读
Step 2	bio-structural-biology-alphafold-predictions	numpy	pLDDT 分级、PAE 域间分析
Step 3	bio-structural-biology-structure-preparation	PDBFixer, PROPKA, reduce, PDB2PQR	受体制备：加氢、质子化、补原子
Step 4	bio-structural-biology-binding-site-detection	fpocket, P2Rank	结合口袋检测与可成药性评分
Step 5	bio-virtual-screening (ligand prep)	RDKit, meeko	SMILES → 3D PDBQT
Step 6	bio-virtual-screening	AutoDock Vina, GNINA	层次化虚拟筛选
Step 7	bio-pose-validation	PoseBusters, RDKit	Pose 物理合理性验证 + 应变能
Step 8	bio-admet-prediction	ADMETlab 3.0, RDKit, chemprop	ADMET 预测与类药性筛选
Step 9	bio-free-energy-calculations	OpenFE, alchemlyb, MBAR	RBFE/ABFE 自由能精排序
补充	bio-ml-docking-rescoring	DiffDock-L, Boltz-2, Chai-1	ML 对接与重打分（替代/补充）

量化阈值速查表 ¶

指标	阈值	含义	来源
pLDDT > 90	Very high	骨架+侧链方向可信	Jumper 2021
pLDDT 70-90	Confident	骨架基本正确	AFDB 定义
pLDDT 50-70	Low	仅粗略拓扑	AFDB 定义
pLDDT < 50	Very low	通常是无序区（IDR）	Akdel 2022
PAE < 5Å（域间）	可信	域间相对位置有约束	AFDB
Druggability Score ≥ 0.5	可成药分诊线	药物分子可结合	Schmidtke 2010
Lipinski violations ≤ 1	类药	MW≤500, LogP≤5, HBD≤5, HBA≤10	Lipinski
Veber pass	可药	RotBonds≤10, TPSA≤140	Veber
PoseBusters PB-valid	物理合理	通过全部 ~20 项检查	Buttenschoen 2024
RMSD ≤ 2Å + PB-valid	现代标准	准确且合理	Buttenschoen 2024
芳香环偏离 ≤ 0.25Å	平面	芳香环非平面 = 不合理	PoseBusters
vdW 体积重叠 < 7.5%	无严重 clash	与蛋白的体积重叠	PoseBusters
内部应变能 ≤ 3 kcal/mol	低应变	~70% 蛋白结合配体满足	Boström 1998
Vina affinity < -7 kcal/mol	可能活性	经验阈值（非保证）	Vina 手册
闭环残差 < 2σ	一致	热力学循环闭合	最佳实践
hERG prob < 0.5	较低关注	需多模型三角验证	ADMETlab 3.0
BBB simple screen	可过 BBB	TPSA≤90, MW≤500, HBD≤3	启发式
Exhaustiveness = 8	Vina 默认	可增到 16/32/64	Vina 手册

关键决策点汇总 ¶

1. 实验结构 vs 预测模型 ¶

PDB 有 holo 实验结构 → 直接用实验结构
PDB 只有 apo → 用 apo + 口袋检测，或考虑 MD 开口袋袋
无实验结构 → AlphaFold DB → 检查 pLDDT/PAE → 裁剪后制备
需要复合物 → AF3/Boltz-1/Chai-1

2. 对接策略选择 ¶

有 holo 受体 + 已知口袋 → GNINA --cnn_scoring rescore
Apo 或预测受体 → DiffDock-L + GNINA rescore + PoseBusters
超大库 → Vina → top 1% GNINA → top 0.1% Boltz-2/FEP
金属/辅因子 → 需要专门力场（ZAFF, MCPB.py）

3. 自由能方法选择 ¶

同系列类似物（R-group SAR）→ RBFE
跨骨架 → ABFE
快速分诊 → Boltz-2 affinity
-低成本探索 → MM/GBSA（仅系列内排序）

4. ML vs 经典对接 ¶

Vina PB-valid + RMSD ≤ 2Å = 58%（PoseBusters Benchmark）
DiffDock PB-valid + RMSD ≤ 2Å = 12%（原始版本）
ML 模型 pose 必须 PoseBusters 验证（不是可选）
保留采样器置信度、独立评分和物理验证作为分离的证据列

参考文献 ¶

Jumper J, et al. (2021) Highly accurate protein structure prediction with AlphaFold. Nature 596:583-589.
Varadi M, et al. (2024) AlphaFold Protein Structure Database in 2024. Nucleic Acids Res 52:D368-D375.
Akdel M, et al. (2022) A structural biology community assessment of AlphaFold2 applications. Nat Struct Mol Biol 29:1056-1067.
Word JM, et al. (1999) Asparagine and glutamine: using hydrogen atom contacts in the choice of side-chain amide orientation. J Mol Biol 285:1735-1747.
Olsson MHM, et al. (2011) PROPKA3. J Chem Theory Comput 7:525-537.
Dolinsky TJ, et al. (2004) PDB2PQR. Nucleic Acids Res 32:W665-W667.
Le Guilloux V, et al. (2009) Fpocket. BMC Bioinformatics 10:168.
Krivak R, Hoksza D (2018) P2Rank. J Cheminform 10:39.
Trott O, Olson AJ (2010) AutoDock Vina. J Comput Chem 31:455-461.
Eberhardt J, et al. (2021) Vina 1.2 features. J Chem Inf Model 61:3891-3898.
McNutt AT, et al. (2021) GNINA 1.0 CNN docking. J Cheminformatics 13:43.
Buttenschoen M, et al. (2024) PoseBusters. Chem Sci 15:3130-3139.
Corso G, et al. (2024) DiffDock-L. ICLR.
Passaro S, et al. (2025) Boltz-2. bioRxiv.
Abramson J, et al. (2024) AlphaFold 3. Nature 630:493-500.
Fu et al. (2024) ADMETlab 3.0. Nucleic Acids Res 52:W422-W431.
Capuzzi et al. (2017) PAINS reality check. J Chem Inf Model 57:417.
Wang L, et al. (2015) FEP+ method. J Am Chem Soc 137:2695-2703.
Mey ASJS, et al. (2020) Alchemical free-energy best practices. Living J Comput Mol Sci 2:18378.
Shirts MR, Chodera JD (2008) MBAR. J Chem Phys 129:124105.
Boström J, et al. (1998) Conformational energy penalties of protein-bound ligands. J Comput-Aided Mol Des 12:383-396.
Schmidtke P, Barril X (2010) Understanding and predicting druggability. J Med Chem 53:5858-5867.
Cimermancic P, et al. (2016) CryptoSite: cryptic binding sites. J Mol Biol 428:709-719.
Lyu et al. (2019) Ultralarge virtual-screening proof of concept. Nature 566:224-229.

AI生成 本文档由星辰超级智能体（TeleAgent）基于 bioSkills 技能库中以下 SKILL.md 文件编写：bio-structural-biology-alphafold-predictions、bio-structural-biology-structure-preparation、bio-structural-biology-binding-site-detection、bio-virtual-screening、bio-pose-validation、bio-admet-prediction、bio-free-energy-calculations、bio-ml-docking-rescoring。所有命令、参数和陷阱警告均来源于技能文件的原始内容。

AI生成

阅读全文

20260910 143546 Crispr 筛选 Hits 鉴定全流程实现详解

2026年09月10日

CRISPR 筛选 Hits 鉴定全流程实现详解 ¶

基于 bio-research-suite 技能库中 15 个 CRISPR 筛选相关技能的实际内容编写
涵盖从 sgRNA 库设计 → FASTQ 计数 → QC → 拷贝数校正 → Hit Calling → 通路富集的端到端流程

场景假设 ¶

你在研究一种新型靶向抗癌药物 compound-X，作用于某癌症细胞系 HCT116。需要设计并执行一个全基因组 CRISPR-Cas9 敲除筛选，找出：
- 致敏基因（sensitizer）：敲除后使细胞对药物更敏感（合成致死）
- 抗性基因（resistance）：敲除后使细胞耐受药物

实验设计：药物处理组 vs. 载体对照组，各 2 个生物学重复， brunello 文库（4 条 sgRNA/基因，约 76,000 条 sgRNA），Day 14 收集细胞。

全流程概览 ¶

 sgRNA库设计 → FASTQ → 计数 → QC → 拷贝数校正 → Hit Calling → 通路富集 → 可视化
     |           |        |       |         |            |           |          |
  Step 1      Step 2   Step 3  Step 4    Step 5       Step 6      Step 7    Step 8

Step 1: sgRNA 文库设计 ¶

涉及技能 ¶

技能	功能
`bio-crispr-screens-library-design`	文库设计全流程（化学体系选择、on-target 打分、control 设计）
`bio-genome-engineering-grna-design`	单条 gRNA 设计

文库 QC 要点 ¶

Gini 指数 < 0.1（质粒库均匀度）
零计数 sgRNA < 0.5%
每条 sgRNA 覆盖 ≥ 500x 细胞（MOI ≈ 0.3）
非靶向对照（NTC）≥ 500 条

# 检查文库均匀度
import pandas as pd
import numpy as np

def library_qc(counts_df):
    """文库质量控制指标"""
    out = pd.DataFrame(index=counts_df.columns)
    out['n_detected'] = (counts_df > 0).sum()
    out['pct_zero'] = (counts_df == 0).sum() / len(counts_df) * 100
    out['pct_lowcount'] = (counts_df < 30).sum() / len(counts_df) * 100
    out['gini'] = counts_df.apply(gini_coefficient)
    out['median_count'] = counts_df.median()
    out['p90_p10_ratio'] = counts_df.quantile(0.9) / counts_df.quantile(0.1).replace(0, np.nan)
    return out

def gini_coefficient(x):
    """计算 Gini 指数"""
    x = np.sort(x[x > 0])
    n = len(x)
    cumsum = np.cumsum(x)
    return (2 * np.sum((np.arange(1, n+1)) * x) - (n + 1) * cumsum[-1]) / (n * cumsum[-1])

关键陷阱 ¶

CRISPRi/a 必须用 FANTOM5 CAGE TSS：靶向没有 CAGE 证据的 TSS 会丢失大部分敲低效率
Cas9 KO 对旁系同源基因不敏感：42% 组成型表达基因因旁系同源缓冲而无法检出（Dede 2020），需用 Cas12a 多重文库
统计预测的 guide 效率不等于实验效率：Azimuth/DeepSpCas9 打分仅作初筛，仍需实验验证

Step 2: FASTQ 预处理 ¶

涉及技能 ¶

技能	功能
`bio-read-qc-fastp-workflow`	fastp 一站式预处理
`bio-read-qc-quality-reports`	FastQC + MultiQC

# fastp 预处理
fastp \
    --input_plasmid_R1.fq.gz --input_plasmid_R2.fq.gz \
    --output_Day0_R1.fq.gz --output_Day0_R2.fq.gz \
    --detect_adapter_for_pe \
    --correction \
    --cut_front --cut_tail \
    --thread 16

# 生成质量报告
fastqc *.fq.gz -t 16 -o qc_reports/
multiqc qc_reports/ -o multiqc_report/

关键陷阱 ¶

fastp 自动接头检测不通用：若使用自定义文库构建 kit，需显式指定 --adapter_sequence
CRISPR 文库的 sgRNA 区段通常固定长度：不依赖 fastp 自动 trim，直接在 MAGeCK count 步骤用 --trim-5 AUTO

Step 3: sgRNA 计数 ¶

涉及技能 ¶

技能	功能
`bio-crispr-screens-mageck-analysis`	MAGeCK count：FASTQ → sgRNA 计数矩阵

MAGeCK count 命令 ¶

mageck count \
    --list-seq brunello_library.csv \
    --sample-label Plasmid,Day0,Veh_r1,Veh_r2,Drug_r1,Drug_r2 \
    --fastq Plasmid.fq.gz Day0.fq.gz Veh_r1.fq.gz Veh_r2.fq.gz Drug_r1.fq.gz Drug_r2.fq.gz \
    --norm-method median \
    --output-prefix screen \
    --trim-5 AUTO

文库文件格式（CSV，列顺序固定）：

sgRNA,Sequence,Gene
BRCA1_1,ATGGATTTATCTGCTCTTCG,BRCA1
BRCA1_2,CAGCAGATACTTGATGCATC,BRCA1
NTC_0001,GACGCATCGAATCAATAGCC,NonTargeting_0001

输出文件 ¶

文件	内容
`screen.count.txt`	原始计数矩阵
`screen.count_normalized.txt`	中位数归一化计数
`screen.countsummary.txt`	每样本 QC（Gini、reads、mapping rate）

归一化方法选择 ¶

`--norm-method`	适用场景	机制	何时失效
`median`（默认）	标准筛选，<40% sgRNA 变化	缩放到所有 sgRNA 中位数	重度筛选（>40% 变化）时中位数偏移
`control`	重度筛选，有 ≥500 条 NTC	用 NTC 中位数缩放	必须提供 `--control-sgrna ntcs.txt`
`total`	等测序深度假设	缩放到总 reads	对 PCR jackpot 敏感

诊断流程：先用 median 跑一遍，如果 QC 显示 CEGv2 PR-AUC < 0.5 但 Gini 正常，换用 control 重新归一化。

Step 4: 质量控制 ¶

涉及技能 ¶

技能	功能
`bio-crispr-screens-screen-qc`	全方位 QC（Gini、replicate correlation、CEGv2 PR-AUC）

六阶段瓶颈检测 ¶

阶段	关键指标	合格阈值	不及格后果
质粒库	Gini, skew, %zero	Gini < 0.1, skew < 2, zero < 0.5%	sgRNA 丢失无法区分
Day-0 感染	文库覆盖, MOI	≥99% guide 检出, MOI 0.3	奠基者效应
筛选 (puro/blast)	存活率, Gini 漂移	30-40% 存活, ΔGini < 0.05	筛选假象
终点	Replicate 相关性	Pearson ≥ 0.8, Spearman > 0.7	噪声主导, FDR 膨胀
生物学信号	CEGv2 PR-AUC	> 0.7 at FDR 5%	筛选无信号
拷贝数假象	扩增区富集	sgRNA 切割数与 depletion 无相关	扩增区假阳性

# Replicate 一致性检查
import pandas as pd
import numpy as np
from scipy.stats import pearsonr, spearmanr

def replicate_concordance(counts_df, groups):
    """检查生物学重复的一致性"""
    results = []
    for group_name, samples in groups.items():
        if len(samples) >= 2:
            log_counts = np.log2(counts_df[samples] + 1)
            for i in range(len(samples)):
                for j in range(i+1, len(samples)):
                    r_pearson, p_p = pearsonr(
                        log_counts[samples[i]], log_counts[samples[j]])
                    r_spearman, p_s = spearmanr(
                        log_counts[samples[i]], log_counts[samples[j]])
                    results.append({
                        'group': group_name,
                        'pair': f"{samples[i]} vs {samples[j]}",
                        'pearson_r': r_pearson,
                        'spearman_r': r_spearman
                    })
    return pd.DataFrame(results)

# 使用
groups = {
    'Vehicle': ['Veh_r1', 'Veh_r2'],
    'Drug': ['Drug_r1', 'Drug_r2']
}
concordance = replicate_concordance(counts, groups)
print(concordance)
# 要求: Pearson >= 0.8, Spearman > 0.7

关键陷阱 ¶

Plasmid Gini > 0.1 → 文库合成或扩增偏差，筛选不可用
Endpoint Gini > 0.30 → 需区分是强生物学选择（药物筛选可接受）还是瓶颈效应
CEGv2 PR-AUC < 0.7 → 筛选质量不达标，hit 列表不可信

Step 5: 拷贝数校正（癌症细胞系必备）¶

涉及技能 ¶

技能	功能
`bio-crispr-screens-copy-number-correction`	CRISPRcleanR / CERES / Chronos

为什么必须做 ¶

机制（Aguirre 2016, Munoz 2016）：
1. 局部扩增 → 某区域有 4-50+ 个拷贝
2. 每条 sgRNA 同时切割所有拷贝
3. 多重切割 → DNA 损伤反应 → G2 期阻滞
4. 细胞停滞 → sgRNA 表现为”必需”（假阳性）

影响范围：所有 Cas9-KO 癌症细胞系筛选，普遍性的，不可通过增加测序深度或重复数修复。

典型案例：ERBB2 在 HER2 扩增的 SK-BR-3 中表现为必需（假阳性）

方法选择 ¶

可用数据	推荐方法	理由
无匹配 CN 数据	CRISPRcleanR	无监督，仅用基因组位置
单细胞系 + 匹配 CN	CRISPRcleanR 或 Chronos	Chronos 更严谨
DepMap 规模（1000+ 系）	Chronos	群体动力学 + CN + 筛选质量

CRISPRcleanR 预校正（R）¶

library(CRISPRcleanR)

# 加载文库注释
data(KY_Library_v1.0)

# 归一化 + 计算折变化
norm_counts <- ccr.NormfoldChanges(
    filename = 'counts.txt',
    min_reads = 30,
    EXPname = 'hct116_screen',
    libraryAnnotation = KY_Library_v1.0
)

# 基因组位置校正
cleaned <- ccr.GWclean(
    norm_counts,
    libraryAnnotation = KY_Library_v1.0,
    saveToPath = 'corrected/'
)

# 校正后的计数可用于 MAGeCK 下游分析

关键陷阱 ¶

CRISPRi/CRISPRa 无需 CN 校正：不切割 DNA，无 DNA 损伤假象
CRISPRcleanR 输出的是校正后的计数，兼容任何下游工具（MAGeCK/BAGEL2/drugZ）
MYC 扩增的结直肠细胞系（10+ 拷贝）、FGFR1 扩增的头颈鳞癌细胞系均为典型假阳性

Step 6: Hit Calling ¶

涉及技能 ¶

技能	功能
`bio-crispr-screens-mageck-analysis`	MAGeCK RRA（两条件）或 MLE（多条件）
`bio-crispr-screens-drugz-chemogenomic`	drugZ（药物筛选专用）
`bio-crispr-screens-bagel-essentiality`	BAGEL2（必需性分类）
`bio-crispr-screens-hit-calling`	跨方法决策树 + 共识策略

方法选择决策树 ¶

实验设计	首选方法	原因
两条件（药物 vs 载体），单细胞系	drugZ	药物筛选专用，双向 Z，对小效应更敏感
多时间点（Day 0→7→14→21）	MAGeCK MLE	RRA 无法建模多时间点
多细胞系面板	Chronos	联合建模 CN + 筛选质量
必需性分类（Day 0 → 终点）	BAGEL2	贝叶斯分类，CEGv2/NEGv1 校准
多筛联合分析	JACKS	共享 sgRNA 效率

本场景：drugZ + MAGeCK RRA 共识策略 ¶

为什么用 drugZ 而非 MAGeCK 作为药物筛选首选：
- drugZ 是双向 Z-score，同时检测致敏（synth）和抗性（supp）
- 对小效应（药物筛选典型 10-30% sgRNA 变化）比 RRA 更敏感
- Colic et al. 2019 基准测试：drugZ hits 对预期 DDR 通路富集远超 MAGeCK/STARS/edgeR

drugZ 分析 ¶

python drugz.py \
    -i screen.count.txt \
    -o drugz_results.txt \
    -c Veh_r1,Veh_r2 \
    -x Drug_r1,Drug_r2

drugZ 输出列：

列	含义
`gene`	基因名
`sumZ`	sgRNA 级 Z-score 之和
`normZ`	归一化后的 per-gene Z
`synth_pval` / `synth_fdr`	致敏方向 p 值 / FDR
`supp_pval` / `supp_fdr`	抗性方向 p 值 / FDR

MAGeCK RRA 分析（对照验证）¶

mageck test \
    --count-table screen.count.txt \
    --treatment-id Drug_r1,Drug_r2 \
    --control-id Veh_r1,Veh_r2 \
    --norm-method median \
    --output-prefix drug_vs_veh \
    --gene-lfc-method median \
    --sort-criteria pos

MAGeCK 输出解读：

列	含义	判定标准
`neg\\|fdr`	负选择（depletion）FDR	< 0.05 → 敏感基因
`neg\\|lfc`	负选择 log2FC	< -1 → 2 倍以上消耗
`pos\\|fdr`	正选择（enrichment）FDR	< 0.05 → 抗性基因
`pos\\|lfc`	正选择 log2FC	> 1 → 2 倍以上富集

BAGEL2 必需性分析（补充）¶

# 1. 计算折叠变化
BAGEL.py fc \
    -i screen.count_normalized.txt \
    -o bagel_fc.txt \
    -c Day0 \
    -e Veh_r1,Veh_r2,Drug_r1,Drug_r2

# 2. 计算 Bayes Factor
BAGEL.py bf \
    -i bagel_fc.txt \
    -o bagel_bf.txt \
    -e CEGv2.txt \
    -n NEGv1.txt

# 3. 精确率-召回率曲线
BAGEL.py pr \
    -i bagel_bf.txt \
    -o bagel_pr \
    -e CEGv2.txt \
    -n NEGv1.txt

BAGEL2 校准：BF > 6 ≈ 90% 后验概率，~5% FDR（Hart 2017 惯例）

共识 Hit 调用 ¶

import pandas as pd

def consensus_hits(drugz_path, mageck_path, fdr_threshold=0.05):
    """跨方法共识 hit 调用"""
    # drugZ hits
    drugz = pd.read_csv(drugz_path, sep='\t')
    drugz_synth = drugz[drugz['synth_fdr'] < fdr_threshold]['gene'].tolist()
    drugz_supp = drugz[drugz['supp_fdr'] < fdr_threshold]['gene'].tolist()

    # MAGeCK hits
    mageck = pd.read_csv(mageck_path, sep='\t')
    mageck_neg = mageck[mageck['neg|fdr'] < fdr_threshold]['id'].tolist()
    mageck_pos = mageck[mageck['pos|fdr'] < fdr_threshold]['id'].tolist()

    # 共识
    consensus_sensitizer = set(drugz_synth) & set(mageck_neg)  # 敏感基因
    consensus_resistance = set(drugz_supp) & set(mageck_pos)   # 抗性基因

    # 单方法 hits（需正交验证）
    drugz_only_synth = set(drugz_synth) - set(mageck_neg)
    mageck_only_neg = set(mageck_neg) - set(drugz_synth)

    return {
        'consensus_sensitizers': consensus_sensitizer,
        'consensus_resistance': consensus_resistance,
        'drugz_only': drugz_only_synth,
        'mageck_only': mageck_only_neg
    }

跨方法分歧时的处理 ¶

分歧模式	可能原因	处理方式
MAGeCK FDR 显著, BAGEL2 不显著	参考集不同；BAGEL2 训练于 CEGv2/NEGv1	必需性判断信 BAGEL2；一般 LFC 信 MAGeCK
MAGeCK 显著, drugZ 不显著	RRA 对小效应不敏感	药物筛选中信 drugZ
MAGeCK 和 Chronos top 100 一致, 100-300 分歧	Chronos 校正了 CN，MAGeCK 没有	癌症细胞系信 Chronos
RRA 和 MLE 在同一数据上分歧	RRA 对离群值更鲁棒；MLE 更敏感	两方法都排名靠前 = 高可信

Step 7: 通路富集 ¶

涉及技能 ¶

技能	功能
`bio-pathway-go-enrichment`	GO 富集分析
`bio-pathway-kegg-pathways`	KEGG 通路分析
`bio-pathway-gsea`	GSEA（fgsea）

分析方式 ¶

library(clusterProfiler)

# 对 hit 基因做 GO 富集
go_result <- enrichGO(
    gene = hit_genes,
    OrgDb = org.Hs.eg.db,
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05,
    qvalueCutoff = 0.2
)

# KEGG 通路富集
kegg_result <- enrichKEGG(
    gene = hit_genes,
    organism = "hsa",
    pAdjustMethod = "BH"
)

# GSEA：用 MAGeCK 排名（neg|score 或 neg|lfc）
library(fgsea)
ranks <- mageck_results$`neg|lfc`
names(ranks) <- mageck_results$id
fgsea_result <- fgsea(
    pathways = kegg_pathways,
    stats = ranks,
    minSize = 15,
    maxSize = 500
)

关键陷阱 ¶

GSEA 用排序列表，不用 p 值过滤列表：按 log2FC 或 signed -log10(p) 排序所有基因
对 hit 基因做超几何检验时使用背景基因 = 文库中所有基因，不是全基因组

Step 8: 可视化 ¶

涉及技能 ¶

技能	功能
`bio-data-visualization-volcano-and-ma-plots`	火山图
`bio-data-visualization-heatmaps-clustering`	热图
`bio-data-visualization-rank-plots`	排名图

import matplotlib.pyplot as plt
import numpy as np

def crispr_volcano(gene_summary_path, direction='neg', fdr_threshold=0.05, lfc_threshold=1.0):
    """CRISPR 筛选火山图"""
    df = pd.read_csv(gene_summary_path, sep='\t')
    lfc_col = f'{direction}|lfc'
    fdr_col = f'{direction}|fdr'

    fig, ax = plt.subplots(figsize=(10, 8))
    sig = (df[fdr_col] < fdr_threshold) & (np.abs(df[lfc_col]) > lfc_threshold)

    ax.scatter(df.loc[~sig, lfc_col],
               -np.log10(df.loc[~sig, fdr_col].clip(lower=1e-10)),
               c='lightgray', alpha=0.3, s=8, label='n.s.')
    ax.scatter(df.loc[sig, lfc_col],
               -np.log10(df.loc[sig, fdr_col].clip(lower=1e-10)),
               c='red' if direction == 'neg' else 'blue',
               alpha=0.7, s=20, label=f'FDR<{fdr_threshold}, |LFC|>{lfc_threshold}')

    top = df.loc[sig].nsmallest(20, fdr_col)
    for _, r in top.iterrows():
        ax.annotate(r['id'], (r[lfc_col], -np.log10(max(r[fdr_col], 1e-10))),
                    fontsize=7, ha='center', va='bottom')

    ax.axhline(-np.log10(fdr_threshold), ls='--', c='black', lw=0.5)
    ax.axvline(0, c='gray', lw=0.5)
    ax.axvline(-lfc_threshold, ls=':', c='gray', lw=0.5)
    ax.axvline(lfc_threshold, ls=':', c='gray', lw=0.5)
    ax.set_xlabel(f'{direction} LFC (log2 fold change)')
    ax.set_ylabel('-log10(FDR)')
    ax.legend()
    return fig

全流程总结图 ¶

 ┌────────────┐   ┌────────────┐   ┌─────────────────┐
 │  Step 1    │   │  Step 2    │   │   Step 3        │
 │ 文库设计    │──►│ FASTQ 预处理│──►│  MAGeCK count   │
 │ Brunello   │   │ fastp      │   │  FASTQ→计数矩阵  │
 └────────────┘   └────────────┘   └────────┬────────┘
                                            │
                    ┌───────────────────────▼──────────────┐
                    │         Step 4: QC                   │
                    │  Gini / replicate / CEGv2 PR-AUC     │
                    └───────────────┬──────────────────────┘
                                    │
                    ┌───────────────▼──────────────────────┐
                    │    Step 5: 拷贝数校正 (癌症细胞系)    │
                    │    CRISPRcleanR / Chronos             │
                    └───────────────┬──────────────────────┘
                                    │
              ┌─────────────────────▼──────────────────────┐
              │            Step 6: Hit Calling              │
              │  ┌─────────┐  ┌──────────┐  ┌───────────┐  │
              │  │ drugZ   │  │ MAGeCK   │  │ BAGEL2    │  │
              │  │ (药物)  │  │ RRA/MLE  │  │ (必需性)  │  │
              │  └────┬────┘  └────┬─────┘  └─────┬─────┘  │
              │       └─────┬──────┘              │        │
              │         共识 hits                │        │
              └─────────────────┬─────────────────────────┘
                                │
              ┌─────────────────▼──────────────┐
              │     Step 7: 通路富集            │
              │  GO / KEGG / GSEA               │
              └─────────────────┬──────────────┘
                                │
              ┌─────────────────▼──────────────┐
              │     Step 8: 可视化              │
              │  Volcano / Heatmap / Rank plot  │
              └────────────────────────────────┘

涉及技能一览表 ¶

步骤	技能	库	主要功能
Step 1	`bio-crispr-screens-library-design`	bio	sgRNA 文库设计
Step 1	`bio-genome-engineering-grna-design`	bio	gRNA 设计
Step 2	`bio-read-qc-fastp-workflow`	bio	fastp 预处理
Step 2	`bio-read-qc-quality-reports`	bio	FastQC/MultiQC
Step 3	`bio-crispr-screens-mageck-analysis`	bio	MAGeCK count
Step 4	`bio-crispr-screens-screen-qc`	bio	全方位 QC
Step 5	`bio-crispr-screens-copy-number-correction`	bio	CN 校正
Step 6	`bio-crispr-screens-mageck-analysis`	bio	MAGeCK RRA/MLE
Step 6	`bio-crispr-screens-drugz-chemogenomic`	bio	drugZ 药物筛选
Step 6	`bio-crispr-screens-bagel-essentiality`	bio	BAGEL2 必需性
Step 6	`bio-crispr-screens-hit-calling`	bio	跨方法决策树
Step 7	`bio-pathway-go-enrichment`	bio	GO 富集
Step 7	`bio-pathway-kegg-pathways`	bio	KEGG 通路
Step 7	`bio-pathway-gsea`	bio	GSEA
Step 8	`bio-data-visualization-volcano-and-ma-plots`	bio	火山图
Step 8	`bio-data-visualization-heatmaps-clustering`	bio	热图
Pipeline	`bio-workflows-crispr-screen-pipeline`	bio	端到端流程

量化阈值速查 ¶

阈值	值	来源
Hit FDR（基因级）	< 0.05	MAGeCK (Li 2014)
Hit LFC 绝对值	> 1（2 倍）	生物学效应大小
中位数归一化适用条件	< 40% sgRNA 变化	中位数假设
文库覆盖	≥ 500x/sgRNA	Joung 2017
MOI	~0.3	低 MOI 确保单插入
Plasmid Gini	< 0.1	MAGeCK-VISPR
Endpoint Gini	< 0.30	MAGeCK-VISPR
Replicate Pearson	≥ 0.8	MAGeCK-VISPR
CEGv2 PR-AUC	> 0.7	Hart 2017
BAGEL2 BF	> 6 ≈ 90% posterior	Hart 2017
NTC 数量	≥ 500	归一化可靠性
每基因 sgRNA 数（MLE）	≥ 3 非零	MAGeCK 收敛要求

本文由 AI 生成，基于 bio-research-suite 技能库中 15 个 CRISPR 筛选相关技能的真实内容（SKILL.md）编写。所有命令、参数、阈值均来源于技能文档原文。

AI生成

阅读全文

20260910 113051 单细胞转录组细胞类型鉴定与轨迹分析全流程实现详解

2026年09月10日

20260910 113033 宏基因组微生物分类鉴定全流程实现详解

2026年09月10日

宏基因组微生物分类鉴定全流程实现详解 ¶

基于 bioSkills 技能库（561 个生物信息学技能）的案例 3 实现
涉及技能：5 个 bio-metagenomics- / bio-microbiome- 技能 | 场景：肠道宏基因组 shotgun 测序数据

场景假设 ¶

研究者拿到一批人类肠道宏基因组（shotgun 全基因组鸟枪法）测序数据（双端 150bp），目标是：

回答”样本里有哪些微生物”（分类鉴定，物种水平）
估计各物种的相对丰度（细胞分数 vs 读段分数）
解析群落的功能潜能（基因家族与代谢通路）
比较不同分组（健康 vs 疾病）的群落多样性

数据特征：paired-end reads，含宿主污染（人类 DNA）、接头与低质量碱基，需要先质控再去宿主。

第一步：质控与宿主去除（前置必备）¶

技能：bio-fastq-quality-control + bio-contamination-screening（或 bio-contamination-controls）
工具：fastp + BBDuk / KneadData

1.1 质控 ¶

fastp -i R1.fastq.gz -I R2.fastq.gz \
    -o clean_R1.fastq.gz -O clean_R2.fastq.gz \
    --qualified_quality_phred 20 --length_required 50 \
    --detect_adapter_for_pe --thread 8

关键点：
- --length_required 50 保留至少 50 bp 的读段（Kraken2 对短读分类效果差）
- 去接头是分类的前置条件，接头会制造”宿主 DNA”假信号

1.2 宿主去除 ¶

bowtie2 -x human_genome_index -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \
    --un-conc-gz host_removed_R%.fastq.gz --very-sensitive

关键点：宿主 reads 占低生物量样本的绝大部分，不去除会淹没微生物信号；且 GenBank 中 >200 万条记录带有人类/载体污染序列（Steinegger & Salzberg 2020），必须去宿主。

二、Kraken2 分类 + Bracken 丰度重估 ¶

技能：bio-metagenomics-kraken（primary_tool: Kraken2）+ bio-abundance-estimation（Bracken）
工具：Kraken2 + Bracken + KrakenTools

2.1 基础分类命令 ¶

kraken2 --db "$KRAKEN_DB" \
    --paired --gzip-compressed --threads 8 \
    --confidence 0.1 \            # 从默认0提高到0.1-0.4，抑制单k-mer假阳性
    --minimum-hit-groups 2 \      # 要求>=2个独立hit区域（临床建议3）
    --report out.kreport \
    --output out.kraken \
    R1.fq.gz R2.fq.gz

核心洞察——Kraken 报告是数据库条件的相似性账本，不是样本清单：
1. 分类 ≠ 存在：--confidence 0 时一个共享 k-mer 就把读分类到数据库最近亲属，即使真实生物不存在
2. 读数 ≠ 丰度：counts 随基因组长度和拷贝数变化，报告百分比是片段分数，不是细胞分数。Bracken 修复的是错误分类层级问题，不修复这个
3. 报告尾部的分类单元是假设，不是发现：单区域 hit、哈希碰撞、污染参考序列会填充长尾

2.2 假阳性控制：独特 minimizer ¶

kraken2 --db "$KRAKEN_DB" --paired --confidence 0.1 \
    --report-minimizer-data \     # 插入2列：总 minimizer + 独特 minimizer
    --report out.krep --output /dev/null \
    R1.fq.gz R2.fq.gz
# 高reads但低独特minimizer的物种 = 假阳性（单一保守区域反复命中）

不要硬编码 “>=1024 unique k-mers”——那是临床数据/特定数据库的经验值，不是常量。用阴性对照标定自己的阈值。

2.3 自建数据库 ¶

kraken2-build --download-taxonomy --db custom_db
for lib in bacteria archaea viral human UniVec_Core; do
    kraken2-build --download-library "$lib" --db custom_db
done
kraken2-build --build --db custom_db --threads 16
kraken2-build --clean --db custom_db   # 删除 library/ + taxonomy/ 减小体积
bracken-build -d custom_db -t 16 -k 35 -l 150   # -k 必须等于 Kraken2 的 k(35)；-l = 读长

大坑：
- --max-db-size 会随机下采样 k-mers 来凑容量，这就是预建 8GB/16GB 数据库的做法，也是为什么高 confidence 会崩溃——读太少无法通过高阈值
- 数据库版本是最大的实验变量：数据库里没有的分类单元，无论怎么调参都检测不到

2.4 交接给 Bracken ¶

bracken -d "$KRAKEN_DB" -i out.kreport -o out.bracken -w out.bracken.kreport \
    -r 150 \   # 必须匹配 built databaseRLENmers.kmer_distrib 且等于真实读长
    -l S -t 10 # species level；-t 是重分配下限，不是置信度

Bracken 修复的只有”错误层级”问题——它把滞留在 genus 的 reads 按基因组先验重新分配到物种。它不修复基因组大小偏差，也不修复假阳性（甚至可能把不存在物种的 reads 重新分配给最近的近亲来放大它）。所以 先做假阳性控制再跑 Bracken。

2.5 常见错误对照表 ¶

错误/症状	原因	解决
成百上千个不可能物种	`--confidence 0`，无 hit-group 下限	提高 confidence；`--minimum-hit-groups >=2`
Bracken “kmer_distrib file not found”	`-r` 没有匹配的内置分布	`bracken-build -l <readlen>`
报告列错位	`--report-minimizer-data` 插入了2列	用 8 列布局解析
小型数据库近乎零分类	截断数据库 + 高 confidence	用大数据库，或降低 confidence

三、MetaPhlAn 4：标记基因法（互补视角）¶

技能：bio-metagenomics-metaphlan（primary_tool: MetaPhlAn）
工具：MetaPhlAn 4.1+ / Bowtie2 / minimap2

3.1 基本分类 ¶

# paired-end 作为 ONE 逗号分隔参数（MetaPhlAn 按单端处理，不用配对信息）
metaphlan reads_R1.fastq.gz,reads_R2.fastq.gz \
    --input_type fastq \
    --index mpa_vJun23_CHOCOPhlAnSGB_202403 \   # 钉住 index，数据库版本是批次变量
    --nproc 8 \
    --mapout sample.map.bz2 \                    # 缓存 read->marker 映射（4.2 前是 --bowtie2out）
    --output_file profile.txt

核心洞察——MetaPhlAn 百分比是细胞分数，不是读段分数：
- 标记基因测 WHO 在（细胞）；k-mer 测 HOW MUCH DNA 在（读段）
- MetaPhlAn 的百分比是基因组大小归一化的分类丰度；Kraken/Bracken 的百分比是 reads 的序列丰度
- 两者不可合并：它们衡量不同量纲的东西（Sun 2021 Nat Methods 18:618），合并成一张表互相比对是常见错误

3.2 SGB 单位 ¶

MetaPhlAn 4 的原子分类单位是 SGB（物种级基因组分箱，~95% ANI 簇）：
- kSGB 有培养参考基因组、有拉丁名
- uSGB 只由 MAG 定义（需要 >=5 个 MAG），报告为占位 ID 没有名字——不要丢弃 uSGB 行，这是 MetaPhlAn 4 的核心创新，能解释 ~20% 更多肠道读段

3.3 从映射缓存重剖分 ¶

metaphlan sample.map.bz2 --input_type mapout \   # 前4.2是 bowtie2out
    --tax_lev s \            # k,p,c,o,f,g,s,t
    --stat_q 0.2 \          # 截断均值：去掉上下20%，平均中间60%
    --output_file profile_species.txt

保存一次 --mapout，之后不同层级/估计器设置无需重新比对——重比对是昂贵步骤，下游全部免费。

3.4 未知分数的缩放效应 ¶

未知估计关闭（4.2 前默认）：已知 taxa 吸收 100%，数据库缺失群落不可见，高估每个已知分类单元
未知估计开启（4.2 默认）：出现 UNCLASSIFIED 行，每个已知丰度按比例缩小
在土壤/海洋/瘤胃中，未知分数可能是样本中最大的”taxon”
4.2 两个 flag 重命名：--bowtie2out → --mapout；--bowtie2db → --db_dir

3.5 MetaPhlAn 常见错误 ¶

错误	原因	解决
与 Kraken 百分比合并	分子分数 vs 读段分数，量纲不同	分开放；如需统一用基因组长度换算
低映射率 = QC 失败	只打标记基因，低映射正常	低映射正常；大未知 = 数据库缺失
合并后行不齐	不同 index 构建的 profiles	钉住一个 index 全部重建

四、HUMAnN 3：功能谱分析 ¶

技能：bio-metagenomics-functional-profiling（primary_tool: HUMAnN）
工具：HUMAnN 3.6+ / MetaPhlAn / DIAMOND / Bowtie2

4.1 分层搜索算法 ¶

HUMAnN 不是对每条读段暴力翻译。三个层级逐层过滤：

分类学预筛（MetaPhlAn）：检测物种，然后只对高于 --prescreen-threshold（默认 0.01%）的物种构建样本特异 ChocoPhlAn 泛基因组
核苷酸泛基因组搜索（Bowtie2）：比对到泛基因组的读段得到 UniRef90 家族 + 高置信物种标签
翻译蛋白搜索（DIAMOND）：未通过的读段 6 框翻译后比对全 UniRef90/50；只有家族，物种是推断的 |unclassified

4.2 运行命令 ¶

# 预QC：接头/质量trim + 去宿主（例如 KneadData）。宿主读段膨胀 UNMAPPED 且浪费 DIAMOND 时间
cat sample_R1.fq.gz sample_R2.fq.gz > sample.fq.gz   # 配对读段没有原生配对，拼接成单文件
humann --input sample.fq.gz --output out/ \
    --taxonomic-profile sample_metaphlan.tsv \   # 复用 MetaPhlAn profile；不要 --remove-temp-output
    --threads 8   # 默认：prescreen 0.01, translated-id 80, gap-fill on, minpath on

# 归一化：跨样本统计前先 PER SAMPLE 归一化（RPK 是深度依赖的）
humann_renorm_table -i out/sample_genefamilies.tsv -o out/sample_cpm.tsv -u cpm
humann_join_tables -i out -o merged_pathabundance.tsv --file_name pathabundance
humann_regroup_table -i merged_pathabundance.tsv -g uniref90_ko -o merged_ko.tsv   # 会加 UNGROUPED 行
humann_split_stratified_table -i merged_pathabundance.tsv -o .   # 统计跑在 UNSTRATIFIED 文件上

4.3 不要丢弃 UNMAPPED/UNINTEGRATED ¶

import pandas as pd

df = pd.read_csv('merged_pathabundance_unstratified.tsv', sep='\t', index_col=0)
meta = pd.read_csv('metadata.tsv', sep='\t', index_col=0)
unmapped = df.loc[['UNMAPPED', 'UNINTEGRATED']]     # 分母，检查但不要丢弃
g1 = meta.index[meta['condition'] == 'healthy']
g2 = meta.index[meta['condition'] == 'disease']
unmapped_shift = unmapped[g2].mean(axis=1) - unmapped[g1].mean(axis=1)
print(unmapped_shift.round(1))
# 然后交给 MaAsLin2/ANCOM-BC（对组分与零膨胀建模），而不是裸 t-test/Mann-Whitney

为什么不能丢：丢弃后重新归一化会把已分配的特征重新缩放到和为 1，凭空制造丰度并抹掉数据库覆盖信号。两个样本 UNMAPPED 20% vs 60% 会看起来一模一样。

4.4 功能谱 = 潜力，不是活性 ¶

一个基因家族或”完整通路”是该群落编码的能力——从来不是速率、通量或表达证据。RNA（宏转录组）才说表达。叙述”疾病微生物组上调 X”是绝对错误——它只是更丰富地携带了该基因。

五、多样性分析（α 和 β 多样性）¶

技能：bio-microbiome-diversity-analysis（primary_tool: phyloseq）
工具：QIIME2 / phyloseq / vegan / scikit-bio

5.1 三个”把手”决定答案 ¶

一个多样性值不是测量，而是三个选择的输出：
1. 抽平深度（rarefaction depth）
2. 系统发育树（tree）
3. 度量（metric）

--p-sampling-depth 是”披着归一化外衣的样本删除旋钮”：core-metrics 把每个样本稀平到该深度，并静默丢弃总计数低于它的样本——不警告，只是 PCoA 上少几个点。丢弃的通常是最低产出的样本（最低生物量的拭子、病重患者、失败提取），丢失几乎从不随机。

5.2 选择抽样深度 ¶

qiime feature-table summarize --i-table table.qza --o-visualization table.qzv  # 深度信息在这
qiime diversity alpha-rarefaction \
    --i-table table.qza --i-phylogeny rooted-tree.qza \
    --p-max-depth 20000 \
    --m-metadata-file metadata.tsv --o-visualization alpha-rarefaction.qza
qiime diversity core-metrics-phylogenetic \
    --i-phylogeny rooted-tree.qza --i-table table.qza \
    --p-sampling-depth 10000 \    # 在 observed-features 平台上；会静默丢弃低深度样本
    --m-metadata-file metadata.tsv --output-dir core-metrics-results

规则：深度要在 observed-features 平台的饱和区；不要用 min(sample_sums)——一个小文库会把所有样本拖到欠饱和噪声；报告丢弃样本列表；在邻近深度确认结论。

5.3 建树（建模选择）¶

优先 SEPP 片段插入到全长参考树（或 Greengenes2），而不是 de novo MAFFT+FastTree。~250bp 读长的 de novo 树分辨率差、中点根任意，会让 unweighted UniFrac 不稳定（Janssen 2018）。

qiime fragment-insertion sepp \
    --i-representative-sequences rep-seqs.qda \
    --i-reference-database sepp-refs-gg-13-8.qda \
    --p-threads 4 \
    --o-tree insertion-tree.qda --o-placements insertion-placements.qda
qiime fragment-insertion filter-features \
    --i-table table.qda --i-tree insertion-tree.qda \
    --o-filtered-table table-sepp.qda --o-removed-table removed-table.qda

5.3 α 多样性（R）¶

library(phyloseq); library(vegan)
ps_rare <- rarefy_even_depth(ps, sample.size = chosen_depth, rngseed = 42, replace = FALSE)
alpha <- estimate_richness(ps_rare, measures = c('Observed', 'Shannon', 'InvSimpson'))
alpha$Group <- sample_data(ps_rare)$Group
alpha$Shannon_eff <- exp(alpha$Shannon)   # 有效种数（Hill q=1），跨工具可比
kruskal.test(Shannon ~ Group, data = alpha)

QIIME2 的 Shannon 是 log2（bits），R 是自然对数（nats）——数值差一个常数因子，直接比较看起来像真实效应。报告 exp(H') 可绕过基数问题。

5.4 β 多样性（必须加权 + 非加权都要报告）¶

wu  <- UniFrac(ps_rare, weighted = TRUE)    # 优势谱系视角
uwu <- UniFrac(ps_rare, weighted = FALSE)   # 稀有谱系 + 拓扑视角
gu  <- as.dist(GUniFrac::GUniFrac(t(as(otu_table(ps_rare), 'matrix')), phy_tree(ps_rare), alpha = 0.5)$unifracs[, , 'd_0.5'])

meta <- data.frame(sample_data(ps_rare))
adonis2(wu ~ Group, data = meta, permutations = 999)   # 显著 = 质心位置差异
permutest(betadisper(wu, meta$Group))                  # 强制：是离散度差异还是位置差异？

如果不做 betadisper，adonis2 的显著可能只是”一组变异更大”而不是”群落组成不同”——尤其不平衡设计（Anderson & Walsh 2013）。unweighted 与 weighted 结论相反时，说明变化在稀有谱系 vs 优势谱系，必须都报告并说明各自暗示的谱系轴。

六、可视化与群落统计 ¶

技能：bio-metagenomics-visualization（primary_tool: phyloseq）
工具：phyloseq / vegan / microViz / pandas + scikit-bio / matplotlib

6.1 诚实的组成图 ¶

Krona 是”里面有什么”的最诚实答案——完整可钻取层级，而不是预先折叠到 top-10 的柱状图：

krepdeport2krona.py -r kraken.kreport -o krona.txt && ktImportText krona.txt -o krona.html

柱状图规则：把稀有 taxa 折叠成”Other”时必须标注隐藏了多少 taxa 和多少百分比；每 n 展示（绝不叠加一个代表样本）；≤12 种色盲安全色；如果有总生物量数据就同时展示绝对负载（相对丰度会隐藏群落 bloom/collapse 的真实幅度）。

6.2 α 多样性：Hill 数 + 丰富度诚实性 ¶

alpha <- estimate_richness(ps_counts, measures = c('Shannon', 'InvSimpson'))
hill_q1 <- exp(alpha$Shannon)   # 有效种数 q=1
hill_q2 <- alpha$InvSimpson     # q=2 优势度加权

shotgun richness 是数据库 + 置信度阈值的读出，不是生物学：25-70% 的 shotgun 物种可能是假阳性。Observed/Chao1 只能用在整数 counts（Bracken 可以，MetaPhlAn 百分比不行）。

6.3 β 多样性：组成学正确选择 ¶

距离	组成学正确	需要树	何时
Bray-Curtis	否	否	领域默认、直观；标注不协调
Aitchison（CLR 上的欧氏）	是	否	组成学正确选择，需处理零
RPCA (DEICODE)	是	否	稀疏无伪计数
UniFrac（加权/非加权）	部分	是	仅当有树（MetaPhl SGB 树）

library(vegan)
dist_bc <- vegdist(otu_matrix, method = 'bray')
pm <- adonis2(dist_bc ~ group, permutations = 999, by = 'terms')
bd <- betadisper(dist_bc, group); permutest(bd, permutations = 999)  # 必须配对

6.4 Python 组成学排序 ¶

import pandas as pd
from skbio.stats.composition import clr, multi_replace
from sklearn.decomposition import PCA

ab = pd.read_csv('merged_abundance.txt', sep='\t', index_col=0)
ab = ab[ab.index.str.contains(r'\|s__') & ~ab.index.str.contains(r'\|t__')]  # 物种行
proportions = (ab.T.values / ab.T.values.sum(axis=1, keepdims=True))
clr_mat = clr(multi_replace(proportions))      # 零替换 + CLR（不是 StandardScaler）
pca = PCA(n_components=2).fit(clr_mat)
coords = pca.transform(clr_mat)

6.5 差异丰度：共识，不是单工具 ¶

DA 方法跨数据集中差异巨大（Nearing 2022）：报告的显著 taxa 更依赖工具而非生物学。至少跑两个组成学感知的方法，交集 = 高置信，并集 = 探索性，全部报告。ALDEx2 与 ANCOM-II 最保守/一致；LinDA 与 ANCOM-BC FDR 控制好。禁止：未校正的 Wilcoxon/t-test 直接作用于原始相对丰度——一行代码同时错在组成学和多重检验上。

工具	模型	引用
ALDEx2	Dirichlet 蒙特卡洛 → CLR	Fernandes 2014 Microbiome
ANCOM-BC	log-linear + 抽样分数偏倚校正	Lin & Peddada 2020 Nat Commun
MaAsLin2	一般线性模型 + 协变量	Mallick 2021 PLoS Comput Biol
LinDA	CLR 上线性模型 + 偏倚校正	Zhou 2022 Genome Biol

七、全流程总结 ¶

原始 shotgun reads（R1+R2）
  │
  ├─ Step 1: fastp 质控（去接头、去低质量）→ Bowtie2 宿主去除
  │
  ├─ Step 2: Kraken2 分类（--confidence 0.1 --minimum-hit-groups 2）
  │          → Bracken 丰度重估计（-r 读长匹配）
  │          → 假阳性控制（--report-minimizer-data / KrakenUniq）
  │
  ├─ Step 3: MetaPhlAn 4 标记基因剖面（细胞分数，uSGB 量化未知 taxa）
  │          （注意：与 Kraken 分数不可合并）
  │
  ├─ Step 4: HUMAnN 3 功能谱（MetaPhl 预筛 → Bowtie2 泛基因组 → DIAMOND）
  │          → 保留 UNMAPPED/UNINTEGRATED → 归一化到 CPM
  │
  ├─ Step 5: 多样性（α: Hill q=0/1/2；β: adonis2 + betadisper 配对）
  │
  └─ Step 6: 可视化（Krona / 堆叠柱 / CLR-PCA / 差异丰度共识）

涉及技能一览 ¶

步骤	技能名	工具	核心作用
1	`bio-metagenomics-kraken`	Kraken2 + Bracken	k-mer 分类 + 丰度重估计
1-2	`bio-metagenomics-metaphlan`	MetaPhlAn 4	标记基因分类（细胞分数）
3	`bio-metagenomics-functional-profiling`	HUMAnN 3	功能基因 + 通路丰度
4	`bio-microbiome-diversity-analysis`	QIIME2/phyloseq	α/β 多样性
5	`bio-metagenomics-visualization`	phyloseq/vegan/microViz	可视化 + 群落统计

三步走总结 ¶

“谁在里面”：Kraken2 + Bracken（读段分数）+ MetaPhlAn（细胞分数）双视角
“能做什么”：HUMAnN 功能谱（潜力，不是活性）
“怎么比较”：多样性 + 排序 + 差异丰度共识

阅读全文

20260910 110123 Rna Seq 差异表达分析全流程实现详解

2026年09月10日

基于 bioSkills 技能库（561 个生物信息学技能）的案例 1 实现
涉及技能：9 个 bio-* 技能 | 场景：肿瘤 vs 正常组织 bulk RNA-seq

场景假设 ¶

研究者在肿瘤组织（treat）和正常组织（ctrl）各取 3 个生物学重复样本，进行 bulk RNA-seq 测序（Illumina NovaSeq，PE150），目标是找出两组间的差异表达基因并定位受影响的生物学通路。

第一步：原始数据质控 ¶

技能：bio-read-qc-quality-reports + bio-fastq-quality
工具：FastQC + MultiQC + Biopython

1.1 逐文件 QC（FastQC）¶

# 对6个样本的R1/R2共12个FASTQ文件运行FastQC
fastqc -t 8 -o qc/raw/ raw_data/*.fastq.gz

FastQC 生成每个文件的 HTML 报告，包含 10 个模块（per-base quality、GC content、duplication、adapter content 等）。

关键注意点：

FastQC 的 pass/warn/fail 红绿灯是为 WGS DNA 校准的，RNA-seq 数据必然在 per-base content（随机六聚体偏倚）和 duplication（高表达转录本）模块上报红，这是正常的，不是数据问题。
NovaSeq/NextSeq 是 2-color chemistry，质量值是分箱的（只有 4 个值），箱线图看起来”方块化”是正常的。

1.2 汇总报告（MultiQC）¶

multiqc qc/raw/ -o qc/multiqc/ -f

MultiQC 把所有样本的报告聚合到一张总览表，方便横向比较。先读 General Statistics 表，看每个样本的 reads 数、Q30 比例、GC 含量是否有离群值。

1.3 编程级质量分析（Biopython）¶

from Bio import SeqIO
from collections import defaultdict

# 确认编码格式：2011年后的Illumina数据均为Phred+33（'fastq'）
# 千万不要猜编码——从仪器元数据确认

# 计算每个read位置的平均质量，识别3'端衰减
position_quals = defaultdict(list)
for record in SeqIO.parse('sample1_R1.fastq', 'fastq'):
    for i, q in enumerate(record.letter_annotations['phred_quality']):
        position_quals[i].append(q)

for pos in sorted(position_quals)[:20]:
    quals = position_quals[pos]
    print(f'Position {pos}: mean={sum(quals)/len(quals):.1f}')

关键判断标准：

Phred 分数	错误率	准确率
Q20	1/100	99%
Q30	1/1000	99.9%
Q40	1/10000	99.99%

Q30 比例是常规质量目标。

第二步：数据预处理 ¶

技能：bio-read-qc-fastp-workflow
工具：fastp（一步完成接头修剪 + 质量过滤 + poly-G 去除 + QC 报告）

# NovaSeq RNA-seq 预处理：轻修剪（aligner 会 soft-clip，不要 --dedup）
fastp -i raw_data/sample1_R1.fq.gz -I raw_data/sample1_R2.fq.gz \
      -o clean/sample1_R1.fq.gz -O clean/sample1_R2.fq.gz \
      --detect_adapter_for_pe \
      --trim_poly_g \
      -q 20 -l 50 \
      -w 8 \
      -h clean/sample1.html -j clean/sample1.json

关键参数解释：

--detect_adapter_for_pe：fastp 通过 R1/R2 overlap 分析自动检测接头，不需要手动指定接头序列，连只有 1 个碱基的 read-through 也能切掉。
--trim_poly_g：NovaSeq 是 2-color chemistry，G 是无信号碱基，poly-G 尾巴质量过滤器看不到但必须去除。
-q 20 -l 50：质量低于 Q20 的碱基过滤，过滤后短于 50bp 的 read 丢弃。
绝对不加 --dedup：fastp 的 --dedup 是序列级去重，会删掉高表达转录本的真实 reads，RNA-seq 中绝不能用。

预处理后再 QC ¶

fastqc -t 8 -o qc/trimmed/ clean/*.fastq.gz
multiqc qc/ -o qc/compare/ -f   # 同时对比 raw vs trimmed

第三步：基因组比对 ¶

技能：bioinformatics-align-star
工具：STAR（splice-aware aligner）

3.1 构建 STAR 索引（一次性）¶

# sjdbOverhang = readlength - 1 = 150 - 1 = 149
STAR --runMode genomeGenerate \
     --genomeDir star_index/ \
     --genomeFastaFiles GRCh38.fa \
     --sjdbGTFfile gencode.v44.annotation.gtf \
     --sjdbOverhang 149 \
     --runThreadN 16

关键点：--sjdbOverhang 应设为 read 长度减 1，决定剪接位点侧翼序列长度，直接影响 junction 检测灵敏度。

3.2 比对 ¶

for sample in sample1 sample2 sample3 sample4 sample5 sample6; do
  STAR --genomeDir star_index/ \
       --readFilesIn clean/${sample}_R1.fq.gz clean/${sample}_R2.fq.gz \
       --readFilesCommand zcat \
       --outFileNamePrefix aligned/${sample}_ \
       --outSAMtype BAM SortedByCoordinate \
       --quantMode GeneCounts \
       --twopassMode Basic \
       --runThreadN 8
done

关键参数：

--outSAMtype BAM SortedByCoordinate：直接输出按坐标排序的 BAM。
--quantMode GeneCounts：STAR 顺便生成 per-gene 计数，第 2 列是无链计数、第 3 列正向链、第 4 列反向链——可据此判断 library strandedness。
--twopassMode Basic：第一遍比对发现新剪接位点，第二遍用扩展后的 junction 注释重新比对，提高灵敏度。

STAR 特性：唯一比对的 read 获得 MAPQ=255，这在 SAM 规范中意为”MAPQ 不可用”，但对 STAR 而言就是”唯一比对”。对后续 featureCounts 计数无影响，但做 RNA variant calling 时 GATK 可能对 255 有不同处理。

第四步：基因水平定量 ¶

技能：featurecounts-counting
工具：featureCounts（Subread）

featureCounts -a gencode.v44.annotation.gtf \
              -o counts.txt \
              -s 2 \
              -p --countReadPairs \
              -T 8 \
              aligned/sample1_Aligned.sortedByCoord.out.bam \
              aligned/sample2_Aligned.sortedByCoord.out.bam \
              aligned/sample3_Aligned.sortedByCoord.out.bam \
              aligned/sample4_Aligned.sortedByCoord.out.bam \
              aligned/sample5_Aligned.sortedByCoord.out.bam \
              aligned/sample6_Aligned.sortedByCoord.out.bam

关键参数：

-s 2：反向链特异性。大部分 dUTP/TruSeq 文库是反向链特异性建库。如果链方向选错，计数会减半。不确定时用 RSeQC infer_experiment.py 先验证，或看 STAR GeneCounts 的第 3/4 列哪个数字大。
-p --countReadPairs：paired-end 数据按 fragment（read pair）计数，而不是按单条 read 计数（不加 -p 会导致计数翻倍）。

验证输出：

cat counts.txt.summary

正常 poly-A 文库的 Assigned 比例应在 70-90%。如果只有 ~50%，很可能链方向选错了。

第五步：计数矩阵质控 ¶

技能：count-matrix-qc
工具：DESeq2 (R) + Python

5.1 加载并检查 ¶

library(DESeq2)
library(edgeR)

counts <- read.csv('counts.txt', sep='\t', row.names=1, skip=1)
counts <- counts[, 6:11]  # 去掉前5列元数据列，只保留计数列
coldata <- data.frame(condition = factor(c('ctrl','ctrl','ctrl','treat','treat','treat')),
                      row.names = colnames(counts))

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~condition)

colSums(counts(dds))        # 每样本总reads
colSums(counts(dds) > 0)    # 每样本检出基因数

5.2 过滤低表达基因 ¶

# 用 edgeR 的 filterByExpr 做有原则的过滤，优于硬编码 CPM>1
keep <- filterByExpr(counts(dds), group = dds$condition)
dds <- dds[keep, ]
cat("保留基因数:", sum(keep), "，过滤基因数:", sum(!keep), "\n")

5.3 归一化和变换 ¶

# VST（方差稳定变换）用于 QC 可视化，不是用于 DE 检验
# n>30 用 vst()，n<30 用 rlog()——本例 n=6 偏小，但 vst 通常也能用
vsd <- vst(dds, blind = TRUE)   # blind=TRUE 用于无监督QC

关键区分：CPM/TPM 不校正组成偏倚，只按样本内总量缩放，只能用于可视化，不能用于 DE 归一化。DESeq2 用 median-of-ratios 估计 size factor 做归一化，这在 DESeq() 内部自动完成。

5.4 PCA 和样本关系 ¶

plotPCA(vsd, intgroup = 'condition')

# 样本距离热图
pheatmap::pheatmap(as.matrix(dist(t(assay(vsd)))), main = 'Sample distance (VST)')

判断标准：

PCA 模式	含义	行动
按 condition 分开	信号强	继续
按批次分开	批次主导	在 design 中加入 batch
某个样本离群	异常或样本错配	查QC、查size factor

5.5 离群值检测 ¶

dds <- DESeq(dds)
cooks <- assays(dds)[['cooks']]   # 注意：不是 results(dds)$cooksd
boxplot(log10(cooks), las = 2, main = "Cook's distance per sample")

5.6 批次效应处理 ¶

# 已知批次：放入 design 公式，让引擎在原始计数层面校正
design(dds) <- ~ batch + condition   # condition 在最后 = 感兴趣的对比

# 绝不要对 DE 输入数据使用 removeBatchEffect() 或 ComBat 后再喂给 DESeq2
# 那样会双重校正

致命情况：如果 batch 与 condition 完全混杂（比如所有 treat 在 batch A，所有 ctrl 在 batch B），统计上无法补救。

第六步：差异表达检验 ¶

技能：deseq2-basics
工具：DESeq2 (R)

library(DESeq2)

# 重新构建（如果上面改了 design）
dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~condition)
dds <- DESeq(dds)

# 提取结果：显式命名对比，alpha=0.05
res <- results(dds, contrast = c('condition', 'treat', 'ctrl'), alpha = 0.05)
summary(res)

# LFC 收缩（apeglm）：降噪但不改变 p 值
res_shrunk <- lfcShrink(dds, coef = 'condition_treat_vs_ctrl', type = 'apeglm')

两个关键陷阱：

默认 results() 静默使用最后一个系数——一定要用 contrast= 显式命名想要的对比，否则可能拿到错误的系数。
收缩后的 LFC 使用的是未收缩 Wald 检验的 p 值——p 值不变，不要在收缩后的值上做基于阈值的 FDR 声明。LFC 收缩是为了 ranking 和可视化，不是为了重新检验。

Wald vs LRT 选择：

检验	适用场景	特点
Wald	两组对比（默认）	检验单个系数是否为 0
LRT	多水平因子/综合检验	需要 `reduced` 模型，检验”任何差异”

七步：结果可视化 ¶

library(DESeq2)
library(EnhancedVolcano)
library(ggplot2)

7.1 离散度诊断 ¶

plotDispEsts(dds)

看基因级离散度（黑点）是否跟随拟合趋势线（红线），收缩后（蓝点）是否被拉向趋势——健康的表现。

7.2 p 值直方图 ¶

res_df <- as.data.frame(res)
ggplot(res_df, aes(x = pvalue)) + geom_histogram(bins = 50) + theme_bw()

零假设下应为均匀分布 + 0 附近一个尖峰。如果是 U 形，说明有隐藏批次效应需要加入协变量。

7.3 MA 图 ¶

# 未收缩 vs 收缩对比
par(mfrow = c(1, 2))
plotMA(res, ylim = c(-5, 5), main = 'Unshrunken')
plotMA(res_shrunk, ylim = c(-5, 5), main = 'apeglm shrunk')

收缩后低计数基因的噪声被拉向 0，云更紧凑。

7.4 火山图 ¶

library(EnhancedVolcano)
EnhancedVolcano(res_shrunk,
    lab = rownames(res_shrunk),
    x = 'log2FoldChange',
    y = 'pvalue',          # 注意：用未收缩的 p 值
    pCutoff = 0.05,
    FCcutoff = 1,
    max.overlaps = Inf)

关键标注：x 轴标注”shrunken log2FC (apeglm)”，y 轴标注”unshrunken Wald p-value”——收缩不重新计算 p 值。

7.5 差异基因热图 ¶

# 取 top 50 差异基因
top_genes <- head(order(res$padj), 50)
mat <- assay(vsd)[top_genes, ]
pheatmap::pheatmap(mat, scale = 'row',     # 行z-score：看pattern
                   annotation_col = coldata,
                   show_rownames = FALSE)

7.6 单基因表达图 ¶

# 查看某个关键基因的表达模式
plotCounts(dds, gene = 'TP53', intgroup = 'condition')

第八步：GSEA 通路富集 ¶

技能：pathway-gsea
工具：clusterProfiler (R)

8.1 构建排序列表 ¶

library(clusterProfiler)
library(org.Hs.eg.db)

# 用 DESeq2 的 stat（Wald 统计量）作为排名指标
# 关键：用带符号的、方差校准的统计量，不用原始 p 值
gl <- res$stat
names(gl) <- rownames(res)
gl <- sort(gl, decreasing = TRUE)

# 去掉 NA
gl <- gl[!is.na(gl)]

关键原则：排名指标必须是有符号（signed）且方差校准的统计量（如 DESeq2 的 stat、limma 的 t），绝不能用原始 p 值——p 值没有方向性。

8.2 运行 GSEA ¶

# GO Biological Process
set.seed(42)   # 必须设种子，保证可重复
gse_go <- gseGO(geneList = gl,
                OrgDb = org.Hs.eg.db,
                ont = 'BP',
                minGSSize = 10,
                maxGSSize = 500,
                seed = TRUE)
gse_go <- setReadable(gse_go, OrgDb = org.Hs.eg.db)

# KEGG 通路
gse_kegg <- gseKEGG(geneList = gl,
                    organism = 'hsa',
                    minGSSize = 10,
                    maxGSSize = 500,
                    seed = TRUE)

8.3 查看结果 ¶

head(gse_go@result, 10)

列名	含义
`enrichmentScore` (ES)	富集得分
`NES`	归一化富集得分（跨基因集可比）
`pvalue`	原始 p 值
`p.adjust`	FDR 校正后 p 值
`leading_edge`	核心驱动基因
`core_enrichment`	leading edge 中的基因列表

8.4 可视化 ¶

# 气泡图
dotplot(gse_go, showCategory = 20) + ggtitle('GSEA GO BP')

# GSEA 运行得分曲线
library(enrichplot)
gseaplot(gse_go, by = 'count', geneSetID = 1)

全流程总结 ¶

原始 FASTQ
  │
  ├─ Step 1: FastQC + MultiQC 质控报告
  │          （RNA-seq 在 per-base content / duplication 报红是正常的）
  │
  ├─ Step 2: fastp 一步预处理
  │          （自动接头修剪 + poly-G去除 + 质量过滤，绝不 --dedup）
  │
  ├─ Step 3: STAR 比对
  │          （sjdbOverhang=149, twopassMode, --quantMode GeneCounts）
  │
  ├─ Step 4: featureCounts 定量
  │          （-s 2 反向链, -p 按fragment计数, 验证Assigned>70%）
  │
  ├─ Step 5: 计数矩阵 QC
  │          （filterByExpr过滤, VST变换, PCA聚类, Cook's离群值, 批次检查）
  │
  ├─ Step 6: DESeq2 差异检验
  │          （Wald检验, alpha=0.05, apeglm LFC收缩, 显式contrast）
  │
  ├─ Step 7: 可视化
  │          （离散度图, p值直方图, MA图, 火山图, 热图, plotCounts）
  │
  └─ Step 8: GSEA 通路富集
             （用 stat 排名, clusterProfiler gseGO/gseKEGG, NES+FDR判定）

涉及技能一览 ¶

步骤	技能名	工具	核心作用
1	`bio-fastq-quality`	Biopython	编程级质量分析
1	`bio-read-qc-quality-reports`	FastQC/MultiQC	全局质控报告
2	`bio-read-qc-fastp-workflow`	fastp	一步预处理
3	`bio-read-alignment-star-alignment`	STAR	剪接感知比对
4	`bio-rna-quantification-featurecounts-counting`	featureCounts	基因水平定量
5	`bio-rna-quantification-count-matrix-qc`	DESeq2	矩阵质控+PCA
6	`bio-differential-expression-deseq2-basics`	DESeq2	差异检验+LFC收缩
7	`bio-differential-expression-de-visualization`	DESeq2+ggplot2	诊断与结果图
8	`bio-pathway-gsea`	clusterProfiler	GSEA通路富集

阅读全文

20260901 104042 抖音作品采集工具提示词

2026年09月01日

请开发一个基于 Flask 的抖音作品采集工具，满足以下全部需求：

## 1. 核心功能
- 用户输入抖音**作者主页链接**或**单个作品链接**，系统自动采集该链接下的所有作品（视频/图集）并下载保存到本地文件夹。
- 对每个作品调用 **DeepSeek API** 进行内容分析，生成 120 字以内的内容总结，并将总结结果展示在网页上。
- 提供简单的 Web 操作界面，支持任务提交、状态实时查看、结果展示。

## 2. 技术栈
- Python 3.8+
- Flask（Web 框架）
- requests（HTTP 请求）
- python-dotenv（环境变量管理）
- 不使用数据库，任务状态存储于内存（简单字典即可）
- 抖音解析：使用**第三方解析 API**（需要在代码中预留接口，用户可自行替换为实际可用的服务，例如聚合数据、快解析等，本提示词中默认使用一个虚构的 API `https://api.example.com/douyin/parse`，开发者需自行替换）
- AI 总结：调用 **DeepSeek Chat API**（`https://api.deepseek.com/v1/chat/completions`），模型使用 `deepseek-chat`，API Key 从环境变量 `DEEPSEEK_API_KEY` 读取。

## 3. 项目结构

douyin_collector/
├── app.py # Flask 主程序（路由、任务管理）
├── config.py # 配置文件（下载目录、API Key 默认设置等）
├── requirements.txt # 依赖列表
├── utils/
│ ├── init.py
│ ├── douyin_parser.py # 抖音解析模块
│ ├── downloader.py # 下载器模块
│ └── ai_summarizer.py # DeepSeek 总结模块
├── templates/
│ ├── index.html # 首页：输入链接的表单
│ └── result.html # 结果页：展示任务状态和结果（含 AJAX 轮询）
├── static/
│ └── style.css # 可选，简单样式（可省略）
└── downloads/ # 下载文件保存目录（运行时自动创建）

## 4. 模块详细要求

### 4.1 抖音解析模块 (`utils/douyin_parser.py`)
- 函数 `parse_douyin_url(url: str) -> List[Dict]`，接收用户输入的链接，返回**作品信息列表**（即使只有一个作品也返回列表）。
- 该函数内部调用第三方解析 API（POST 或 GET 均可，示例使用 GET 请求，参数 `url`）。
- 解析 API 返回 JSON 格式，假设包含字段：
  - `is_author`: 布尔值，表示是否为作者主页（若为作者主页，则 `works` 字段是一个包含多个作品信息的列表）
  - 若 `is_author` 为 `True`，则从 `data["works"]` 中提取作品列表；
  - 若为 `False`，则直接返回一个作品对象（包装成列表）。
- 每个作品对象包含以下字段（根据实际 API 可能有所不同，但需在代码中做兼容处理）：
  ```python
  {
      "type": "video" 或 "image",
      "author": "作者昵称",
      "title": "作品标题或描述",
      "video_url": "无水印视频直链（type=video时必填）",
      "image_urls": ["图片链接1", "图片链接2", ...] （type=image时必填）,
      "cover": "封面图链接（可选）",
      "duration": 15,  # 视频时长（秒，可选）
      "music": "背景音乐名称（可选）",
      "likes": 100,    # 点赞数（可选）
      "comments": 20,  # 评论数（可选）
      "shares": 5,     # 分享数（可选）
      "create_time": "2023-01-01 12:00:00"  # 发布时间（可选）
  }
  ```
- 注意：第三方 API 可能没有所有这些字段，需要给出默认值（例如 `duration` 默认 0，`music` 默认 "无"，等等）。
- 如果解析失败（网络错误、API 返回错误等），抛出异常，由上层捕获。

### 4.2 下载器模块 (`utils/downloader.py`)
- 函数 `download_work(work: Dict, save_dir: str) -> str`，接收单个作品信息和保存目录，返回下载后的本地文件路径（视频文件路径或图集文件夹路径）。
- 根据 `work["type"]` 分两种情况：
  - **视频**：从 `work["video_url"]` 下载视频文件，文件名格式建议为 `{作者}_{标题前20字符}.mp4`（需处理非法字符），使用流式下载（requests.get(stream=True)）。
  - **图集**：创建子文件夹，文件夹名同样为 `{作者}_{标题前20字符}`，依次下载所有图片，保存为 `image_1.jpg`, `image_2.jpg`, ...。
- 确保 `save_dir` 存在，若不存在则创建。
- 如果下载失败，抛出异常。

### 4.3 AI 总结模块 (`utils/ai_summarizer.py`)
- 函数 `summarize_work(work: Dict) -> str`，接收作品信息，返回内容总结字符串（不超过 120 字）。
- 使用 `requests.post` 调用 DeepSeek API：
  - URL: `https://api.deepseek.com/v1/chat/completions`
  - 请求头：`Content-Type: application/json`，`Authorization: Bearer {DEEPSEEK_API_KEY}`
  - 请求体 JSON：
    ```json
    {
      "model": "deepseek-chat",
      "messages": [
        {"role": "system", "content": "你是一个严谨的内容总结助手。"},
        {"role": "user", "content": "这里是完整的提示词模板（见下）"}
      ],
      "temperature": 0.7,
      "max_tokens": 300,
      "stream": false
    }
    ```
  - API Key 从环境变量 `DEEPSEEK_API_KEY` 读取，若未设置则返回错误信息。
- **完整提示词模板**（用户消息内容，使用 f-string 填充作品信息）：
  ```
  你是一个专业的内容分析师，擅长从抖音作品的元数据中提炼核心信息。请根据以下作品信息，生成一段简洁、准确、有吸引力的内容总结（不超过120字），帮助用户快速了解作品主题和亮点。

  【作品信息】
  作者：{author}
  标题/描述：{title}
  类型：{'视频' if type == 'video' else '图集'}
  时长：{duration}秒
  背景音乐：{music}
  点赞数：{likes}
  评论数：{comments}
  分享数：{shares}
  发布时间：{create_time}

  【总结要求】
  1. 概括作品的主要内容或主题；
  2. 突出作品的特色或可能的爆点（如情感、搞笑、知识、颜值等）；
  3. 语言流畅，避免冗余，直接输出总结文本，不要任何前缀或解释。
  ```
- 如果 API 调用失败（如网络问题、密钥错误等），捕获异常并返回错误信息字符串，格式如 `"DeepSeek API 调用失败：{错误详情}"`。

### 4.4 Flask 主程序 (`app.py`)
- 定义以下路由：
  - `GET /`：渲染 `index.html`，显示一个表单，包含一个输入框（name="url"）和提交按钮。
  - `POST /submit`：接收表单中的 `url`，创建任务 ID（使用 `os.urandom(8).hex()`），将任务信息存入全局字典 `tasks`（结构见下），启动后台线程执行 `process_task(task_id, url)`，然后重定向到 `/task/<task_id>`。
  - `GET /task/<task_id>`：根据任务 ID 从 `tasks` 字典获取任务，渲染 `result.html`，传入任务对象和任务 ID。
  - `GET /api/task/<task_id>`：返回该任务的 JSON 数据（用于前端 AJAX 轮询），若任务不存在返回 404。
- **任务字典结构**：
  ```python
  tasks[task_id] = {
      'status': 'pending',   # pending / parsing / downloading / completed / error
      'message': '任务已提交',
      'results': []           # 最终结果为列表，每个元素为字典，包含 work, local_path, summary
  }
  ```
- `process_task(task_id, url)` 函数逻辑：
  1. 更新状态为 `parsing`，调用 `parse_douyin_url(url)` 获取作品列表。
  2. 如果列表为空或解析失败，设置状态为 `error`，消息为对应错误，返回。
  3. 更新状态为 `downloading`，遍历作品列表，对每个作品：
     - 调用 `download_work(work, config.DOWNLOAD_DIR)` 下载，得到本地路径。
     - 调用 `summarize_work(work)` 生成 AI 总结。
     - 将 `{'work': work, 'local_path': local_path, 'summary': summary}` 追加到 `tasks[task_id]['results']`。
     - 更新消息为 `下载中 {idx+1}/{total}`。
  4. 全部完成后，设置状态为 `completed`，消息为 `处理完成，共 {total} 个作品`。
  5. 若过程中出现异常，设置状态为 `error`，消息包含错误详情。

### 4.5 配置文件 (`config.py`)
- 定义下载目录 `DOWNLOAD_DIR = os.path.join(os.path.dirname(__file__), 'downloads')`
- 尝试从 `.env` 文件加载环境变量（使用 `python-dotenv` 的 `load_dotenv()`）
- 不直接写死 API Key，全部从环境变量读取。

### 4.6 前端模板
- `index.html`：简单的 HTML 表单，包含一个输入框和提交按钮，表单 action 指向 `/submit`，method 为 POST。
- `result.html`：显示任务状态（status 和 message），并使用 JavaScript 每 2 秒通过 `fetch('/api/task/{{ task_id }}')` 轮询任务状态，当状态变为 `completed` 或 `error` 时停止轮询，并渲染结果列表（如果有 results）。结果列表中每一项显示：作品标题、作者、本地路径、AI 总结。页面底部提供返回首页的链接。

## 5. 依赖文件 `requirements.txt`

flask==2.3.2
requests==2.31.0
python-dotenv==1.0.0

## 6. 运行说明
- 用户需要在项目根目录创建 `.env` 文件，内容为 `DEEPSEEK_API_KEY=sk-xxxx`（实际密钥）。
- 用户需要将 `utils/douyin_parser.py` 中的第三方解析 API URL 替换为实际可用的服务，并确保返回格式符合上述约定（或自行调整解析逻辑）。
- 运行 `python app.py` 启动服务，默认端口 5000。

## 7. 注意事项
- 代码中要包含适当的注释，特别是关键步骤（如解析、下载、API 调用）。
- 请确保所有模块都能正确导入，路径处理正确（使用 `os.path`）。
- 需要处理文件名中的非法字符（如 `/ \ : * ? " < > |`），替换为下划线。
- 异常处理要完善，避免进程崩溃，所有错误都应反馈到任务状态中。
- 生成的代码应该是完整可运行的，不需要用户再手动补充核心逻辑。

请严格按照以上要求生成完整的项目代码，包括所有目录下的文件内容。

阅读全文

20260831 152705 Ai时代落地攻坚者Fde前沿部署工程师全景解析

2026年08月31日

AI时代落地攻坚者：FDE前沿部署工程师全景解析 ¶

全球企业在生成式AI领域投入超300亿美元，但95%的项目无法产生正向利润——技术早已不是瓶颈，技术与真实业务之间的最后一公里彻底断裂。FDE（前沿部署工程师）作为AI时代唯一连接标准化技术与非标场景的复合型岗位，正以年薪百万的姿态成为科技招聘市场的核心焦点。本文从产品经理视角，拆解FDE的底层逻辑与完整闭环。

很多AI小白都会有一个困惑：如今大模型、AI算法层出不穷，开源模型随处可得，企业AI团队不断扩招、预算持续加码，为什么绝大多数企业的AI项目，最终都沦为“演示版产物”，无法真正落地产生业务价值？

麻省理工学院NANDA团队2025年发布的《生成式AI的鸿沟：2025年商业AI现状》给出了扎心答案：全球企业在生成式AI领域累计投入超300-400亿美元，但95%的AI项目无法对企业利润表产生可衡量的正向影响。

行业核心瓶颈早已不是“AI技术不够强”，而是技术与真实业务之间的最后一公里彻底断裂——实验室里完美运行的模型，落地到企业复杂、非标、老旧的真实业务场景中，立刻失效、水土不服。

在传统互联网、软件行业，产品落地的闭环相对顺畅：标准化产品+远程迭代+线上更新即可完成价值交付。但AI行业完全不同，AI价值的核心是“场景化赋能”，没有通用万能的AI产品，只有适配特定业务、特定流程、特定数据的AI解决方案。

这种行业变革下，传统研发、产品、交付岗位的分工壁垒彻底失效，一个全新的核心岗位应运而生——FDE（Forward Deployed Engineer，前沿部署工程师）。

尤其在近几年科技招聘市场上，全球顶尖AI企业疯狂扩招FDE，岗位薪资远超普通研发。字节跳动开出月薪3.5万到7万元、全年15薪的待遇，年薪最高预计达105万元；阿里云智能给出2万到5万元的月薪、全年16薪的方案；在硅谷，前沿实验室资深FDE的年总薪酬中位数达到48.5万美元。

与此同时，OpenAI联手19家PE投入超40亿美元成立部署公司；Anthropic与黑石、高盛等合作成立估值15亿美元的企业服务公司；亚马逊云科技宣布投入10亿美元组建数千人规模的FDE团队。FDE相关岗位从2025年4月的643个飙升到2026年4月的5330个，同比增长近729%。

FDE成为AI行业最稀缺、最核心的岗位之一。但绝大多数AI小白、职场从业者、甚至企业管理者，都无法清晰定义FDE，分不清它和普通开发、实施、售前的区别，更看不懂企业高薪招聘FDE的底层逻辑。

这究竟是AI时代催生的全新黄金职业，还是“新瓶装旧酒”的概念炒作？

作为一名产品经理，我想用一套完整的产品思维框架，帮AI小白真正理解FDE——它是什么、做什么、为什么火、值不值得。

第一章什么是FDE？——从产品经理的视角重新定义 ¶

1.1 一句话说清楚FDE ¶

如果用一句话向AI小白解释FDE，我会这样说：

FDE是“带着工具箱上前线的工程师”——他驻扎在客户现场，既懂业务又写代码，把一个模糊的AI想法拆解成可运行、可评估、可持续维护的生产系统。

拆开来看，这个定义包含了四个关键词：

“上前线” ——不在总部写通用代码，而是到客户那里去
“带工具箱” ——拥有完整的工程技术能力
“既懂业务又写代码” ——横跨业务理解与技术实现
“拆解成生产系统” ——不只做演示（Demo），而是交付真正能用的东西

从字面拆解，Forward Deployed（前沿部署）是军事术语，核心含义是脱离后方总部，深入一线战场，就地决策、就地执行、就地迭代。这也是FDE和所有传统技术岗位最核心的区别：传统工程师是“后方办公、远程交付”，FDE是“前线驻扎、现场攻坚、闭环落地”。

结合产品经理的价值闭环理论，我给FDE下一个小白能彻底看懂的精准定义：

FDE是AI时代专属的复合型价值落地岗位，是连接企业标准化AI技术能力与客户非标真实业务场景的唯一核心桥梁。它打破了传统“产品定需求、研发做功能、交付做部署”的割裂分工，一人承担「需求挖掘、现场落地、问题攻坚、产品反哺、价值验证」全链路工作，核心目标不是“交付产品”，而是“落地价值”。

简单来说：传统工程师负责“把产品做出来”，FDE负责“把产品用起来、用出效果、迭代更好用”。

1.2 FDE的起源：从Palantir的“不得已”到AI时代的“必然”¶

FDE并非AI时代的凭空发明。要理解它，必须回到它的源头——Palantir。

2003年，Palantir成立，最初的服务对象是美国情报与国防机构。这些客户面临一个根本性难题：他们的工作内容高度机密，无法向外部工程师详细描述需求。Palantir的创始人给潜在客户展示Demo后，得到的反馈是：“这个产品太糟糕了，和我们做的事情完全没关系。”

传统的软件开发模式在这里彻底失效——不可能先写需求文档、再开发、再交付。Palantir的解法是：让工程师直接进驻客户现场。工程师不只是技术支持人员，而是真正嵌入客户组织、理解业务痛点、直接动手写程序解决问题的人。

这个模式逐渐成为Palantir的核心竞争优势，并催生了“Forward Deployed Engineer”这个职称。在Palantir内部，FDE被称为“Delta”，与负责产品开发的“Dev”（Software Engineer）形成互补。直到2016年左右，Palantir的FDE人数甚至超过了传统软件工程师。

Palantir官方将这种模式描述为“工程版本的反向传播”（the human equivalent of backpropagation）——工程师尽可能贴近问题本身，与核心工程团队协作，不断将现场反馈合成并转化为新功能。

AI浪潮兴起后，OpenAI、Anthropic等公司开始广泛采用这个模式。为什么？因为生成式AI的落地面临同样的困境——Demo很漂亮，但进了真实企业环境就撞上各种现实问题。FDE这个“老办法”，在AI时代找到了全新的用武之地。

1.3 FDE的核心底层逻辑（产品二元性理论）¶

市面上绝大多数FDE科普都只罗列岗位职责，却忽略了其最核心的二元性本质，这也是90%AI小白看不懂FDE的根本原因。

从产品经理岗位迭代、价值沉淀视角来看，FDE天生具备交付落地+产品迭代双向闭环能力，两股力量相互制衡、相互赋能，缺一不可：

左手·交付端（即时价值）：直面客户真实业务场景，解决非标落地问题，适配客户老旧系统、复杂流程、特殊规则，快速完成AI功能部署、调试、试运行，对单客户业务价值落地全权负责，解决“能用、好用”的问题。

右手·产品端（长期价值）：在现场落地过程中，持续挖掘通用痛点、共性需求、落地坑点，将单个客户的非标解决方案，抽象、沉淀为可复用的产品能力、通用模板、标准化流程，反哺后方产品与研发团队，对产品整体迭代升级负责，解决“可复制、可规模化”的问题。

这组二元性逻辑，是FDE区别于售前、实施、外包、普通研发的核心护城河，也是其高薪、稀缺的底层原因。普通技术岗位只承担单一维度工作，而FDE承担“短期交付价值+长期产品价值”的双重闭环，是企业AI商业化规模化的核心枢纽。

1.4 一张图看懂FDE的定位 ¶

这张图的核心逻辑：FDE是连接“客户的模糊需求”与“可运行的生产系统”之间的桥梁，同时将一线洞察反馈给产品团队，形成正向循环。

Palantir官方曾给出一个精准且经典的定位：FDE的工作模式与创业公司CTO高度一致。在小型前线团队中，FDE对高风险、高难度的AI落地项目端到端全权负责，没有繁琐的分工流转，无需层层审批，可现场快速决策、快速落地、快速迭代。

不同于传统岗位“各司其职、分段交付”，FDE是全链路负责人：需求听不懂自己拆解、系统跑不通自己调试、流程不匹配自己改造、产品有缺陷自己反馈、效果不达标自己优化。这种“一站式闭环能力”，正是AI复杂落地场景最需要的核心能力。

第二章 FDE的职责范围——从进场到离场的完整闭环 ¶

2.1 职责全景图 ¶

一个FDE的完整工作周期，可以拆解为五个阶段：

2.2 各阶段详解 ¶

阶段一：进场与发现

FDE进驻客户现场后的第一件事，不是写代码，而是理解。他需要绘制客户现有的系统地图——哪些数据源、哪些遗留系统、什么权限结构、什么合规要求。更重要的是，他要理解业务人员真正的工作方式——那些写在流程文档之外的、实际发生的操作。

AI落地最大的坑，是“伪需求落地”。客户往往无法清晰表达自身痛点，只会提出模糊的功能诉求，传统售前、产品容易被表面需求误导，最终导致“做出来的功能没人用，用户真正的痛点没解决”。

FDE的首要职责，就是深入业务一线，完成需求去伪存真、深度拆解：

沉浸式驻场调研：全程扎根客户业务场景，跟随业务人员完成全流程工作，看懂真实业务逻辑、流程卡点、效率痛点，而非依靠客户口头描述、书面文档判断需求；
模糊需求翻译：将客户“想要智能提效、想要智能分析”等模糊诉求，翻译为可落地、可量化、可验证的AI技术需求；
需求优先级判定：结合业务价值、落地难度、投入成本，筛选高价值刚需需求，剔除无效伪需求，避免资源浪费；
stakeholder对齐：对接客户业务、IT、管理层多角色，统一落地预期，明确AI落地的核心指标、验收标准。

一个典型的场景是：客户说“我想要一个AI客服系统”。FDE进场后发现，真正的痛点是多个业务系统的客户数据完全不通，客服效率低下只是表层症状。如果直接做AI客服，核心问题根本没解决。FDE要先穿透表层，找到真正值得解决的问题。

阶段二：方案设计

找到真实问题后，FDE要将模糊的业务诉求转化为具体的技术方案。这需要将“我想用AI降本增效”这类空泛想法，拆解成可落地的技术路径。

标准化AI产品只能适配通用场景，而企业真实业务场景90%都是非标场景：老旧系统不兼容、数据格式混乱、业务流程特殊、安全权限严苛、定制化规则繁多。这也是普通AI产品落地失败的核心原因。

FDE的核心工程职责，就是解决非标场景的最后一公里落地难题：

现场系统适配：对接客户本地部署系统、老旧ERP、CRM、工控系统，解决接口不通、数据隔离、格式不兼容等问题；
快速原型搭建与迭代：基于公司核心AI模型、标准化平台，快速开发定制化适配模块，搭建适配客户场景的落地原型，无需后方研发远程排期等待；
生产环境部署调试：独立完成AI系统的现场部署、压力测试、漏洞调试、安全适配，保障系统稳定接入客户生产环境，真正实现业务可用；
问题即时攻坚：现场解决落地突发问题，比如数据异常、模型精度偏差、流程适配冲突、权限报错等，杜绝项目卡顿、烂尾。

阶段三：快速构建

FDE的核心交付物不是PPT或方案文档，而是可运行的生产级代码。他们通常需要在几天内从模糊的客户问题推进到可工作的原型（Proof of Concept），在几周内推进到生产级部署。

不同于销售以“签单”为核心目标，FDE以“价值深化”为核心目标。通过沉浸式驻场落地，FDE最懂客户业务痛点与潜在需求，是企业挖掘存量客户价值、提升客户粘性的核心角色：

解答客户技术疑问，解决长期使用问题，建立深度信任；
挖掘客户潜在业务场景，推动AI能力跨场景复用，实现客户项目扩容、续费增购；
收集客户长期战略需求，为企业产品长期规划、行业布局提供一线依据。

FDE是企业前线与后方的核心枢纽，需要打通销售、产品、研发、算法、售后全团队：向前对接客户业务，向后同步落地进度、问题、需求，协调后方资源解决前线无法攻克的核心技术难题，实现“前线落地+后方迭代”的高效联动闭环。

阶段四：部署上线

这是最容易被低估的环节。AI模型在Demo环境里完美运行，但进了真实企业环境就撞上：遗留系统没有标准接口、数据质量参差不齐、权限体系不兼容、合规审计不通过。FDE要逐一解决这些问题，让系统真正跑起来。

产品经理核心准则：没有价值验证的落地，都是无效交付。很多AI项目完成部署后就宣告结束，但从未验证是否真正提升业务效率、降低成本，这也是AI项目“看似落地、实则无效”的关键。

FDE区别于普通工程师的核心职责之一，就是对业务结果全权负责：

业务人员培训：指导客户员工使用AI系统，梳理标准化使用流程，解决用户操作卡点；
落地效果监控：实时跟踪AI系统运行数据、业务使用数据，统计效率提升、成本降低、错误率下降等核心指标；
价值闭环验证：对照前期约定的验收标准，验证AI落地的真实业务价值，形成可量化的落地报告；
体验优化调整：根据用户使用反馈、业务运行数据，现场微调系统功能、流程逻辑，提升产品适配度。

阶段五：沉淀与反馈

FDE的职责不止于交付。他还要培训客户团队，让他们在FDE离场后能自主运维；同时将现场的经验教训反馈给产品团队，推动产品迭代。

这是FDE最核心、最稀缺、最容易被忽略的职责，也是区分“真FDE”和“假换皮FDE”的核心标准。普通交付岗位只服务单客户，做完即结束，而FDE承担企业产品规模化迭代的核心数据源职责。

FDE需要将单客户的非标落地经验，转化为企业可复用的通用产品能力：

共性痛点提炼：汇总多个客户的落地坑点、高频需求、适配难点，筛选行业通用问题；
非标方案抽象：将单个客户的定制化适配方案、接口改造、流程优化逻辑，抽象为通用模板、可配置组件、标准化适配工具；
产品需求反哺：向后方产品、研发团队输出精准的迭代需求，推动核心产品优化模型精度、新增适配功能、完善行业解决方案；
行业方法论沉淀：总结不同行业的AI落地流程、适配标准、避坑指南，搭建企业行业落地知识库，降低后续客户的落地成本。

2.3 FDE工作的时间分配 ¶

根据行业实践，FDE的工作时间大致可以这样分配：

25% ：纯写代码
50% ：做集成和调试（对接客户系统、处理数据、解决兼容性问题）
25% ：沟通和项目管理（与客户对齐需求、协调资源、推进决策）

第三章 FDE与其他工程师的区别——一张表看懂 ¶

很多AI小白、职场从业者最容易混淆的问题：FDE到底和后端开发、算法工程师、实施工程师、售前架构师、解决方案工程师有什么区别？

很多企业为了招人，会把传统实施、售前岗位包装成FDE，导致市场上“假FDE”泛滥。本节通过产品分工理论、价值导向、工作场景、能力模型、薪资逻辑五大维度，做精准区分，同时附上真假FDE鉴别标准，彻底杜绝认知混淆。

3.1 核心岗位差异化对照表 ¶

3.2 关键岗位深度差异化解读 ¶

3.2.1 FDE vs 算法工程师：模型“能用”和“好用”的分工

AI小白最容易混淆的两个岗位：算法工程师负责造出好模型，FDE负责用好模型。

算法工程师的核心工作是在实验室优化模型参数、提升通用精度、训练数据集，追求“模型理论性能最大化”，完全脱离具体业务场景。这就导致很多模型实验室效果满分，但落地到真实场景中，因为数据混乱、流程不适配、场景特殊化，完全无法使用。

而FDE不做模型训练、不调核心算法，核心工作是把算法工程师做好的通用模型，适配到具体业务场景中，解决模型落地的适配问题、流程问题、数据问题、体验问题。

简单总结：算法工程师解决“技术上限”，FDE解决“落地下限”。没有FDE，再好的AI模型都是空中楼阁。

3.2.2 FDE vs 普通开发工程师：后方造产品 vs 前线做价值

Palantir内部有一个非常精辟的概括：

Dev（传统软件工程师）的焦点是 “一项能力，服务多个客户” （one capability， many customers）
Delta（FDE）的焦点是 “一个客户，多项能力” （one customer， many capabilities）

普通后端/前端开发工程师，是典型的“后方执行者”，所有工作都基于产品经理输出的需求文档，只负责功能开发，无需关注需求真伪、业务适配、落地效果。开发完成、测试上线，工作即结束，完全不接触客户、不了解业务。

FDE是“前线决策者+执行者”，没有固定需求文档，需要自己挖掘需求、判断需求、定制方案、落地实施、迭代优化。普通开发是“为海量客户做一个标准化功能”，FDE是“为一个客户打通全链路业务能力”，并将能力复用给海量客户。

3.2.3 FDE vs 实施工程师：交付闭环 vs 执行交付

这是最容易被混淆、也是市面上“假FDE”最多的岗位。很多企业把实施工程师改名FDE，本质是岗位包装。

实施工程师的核心是标准化执行：产品定型后，按照固定流程完成部署、安装、培训，只能解决标准化问题，无法做定制化开发、无法迭代产品、无法处理复杂非标场景。一旦遇到系统不兼容、流程冲突、需求变更，只能反馈后方团队，无法现场解决。

而FDE的核心是非标攻坚+价值闭环：具备独立开发、现场改造、即时迭代的能力，能解决所有非标落地难题，且核心目标不是“交付完成”，而是“业务产生价值”，同时承担产品迭代沉淀的核心职责。

3.2.4 FDE vs 售前架构师：落地结果负责 vs 签单结果负责

售前架构师的所有工作，核心目标都是辅助销售签单，工作集中在项目前期，签单后立刻退出项目，不参与落地、不负责效果。薪资结构包含销售提成、业绩分红，和签单金额强绑定。

FDE的工作集中在项目签单后，对落地效果、业务价值、长期迭代负责，不参与签单、无销售提成，薪资结构为“底薪+绩效+股票”，和落地价值、产品迭代效果强绑定。

3.3 一个更宏观的视角：FDE是“三栖人才”¶

综合来看，FDE是一个高度复合的角色：

20%销售：理解客户需求、建立信任、推动合作
30%产品：需求分析、方案设计、产品思维
50%工程：写代码、做集成、部署上线

它处于产品管理、软件工程和战略咨询三者的交汇点。

基于行业头部企业JD与落地模式，总结出5条压力测试标准，小白可直接用来辨别FDE岗位真伪：

1.薪资结构（最硬核标准）：真FDE无销售提成，薪资为底薪+绩效+股权；有业绩提成、签单分红的均为售前/实施换皮岗。

2.工作内容：真FDE70%以上时间写生产代码、做现场开发调试；仅写方案、做PPT、演示原型的是架构师/售前。

3.产品回流通道：真FDE有完善的机制，将现场经验沉淀为产品通用能力；无产品反哺、做完即走的是外包/实施。

4.进场时机：真FDE签单后进场，对落地效果兜底；签单前进场辅助谈单的是售前。

5.需求决策权：真FDE可自主判断需求真伪、优化方案、拒绝无效需求；只能被动执行客户要求的是基础交付岗。

第四章企业为什么需要FDE——痛点、必要性、价值 ¶

4.1 产品战略视角：解决AI行业核心痛点 ¶

从产品经理痛点-解决方案-价值匹配理论来看，FDE的诞生不是企业跟风，而是AI行业发展到现阶段的必然战略选择。传统岗位分工模式，已经无法解决AI商业化落地的三大核心痛点，而FDE是唯一适配的解决方案。

4.1.1 痛点一：AI技术标准化与业务场景非标化的天然矛盾

传统互联网产品（社交、电商、办公软件）的用户需求、使用场景高度标准化，一套产品可以服务千万级用户，依靠“后方研发迭代+线上更新”即可完成价值交付。

但AI产品完全不同，尤其是企业级AI、行业AI，没有标准化场景，只有非标化需求。不同行业、不同企业、不同部门的业务流程、数据规范、系统架构、核心痛点完全不同。通用AI模型必须经过场景适配、定制改造、流程融合，才能产生价值。

传统分工模式下，后方研发不了解一线业务，前线交付不懂技术迭代，导致“研发做的功能用不上，客户需要的功能没人做”，形成落地死循环。而FDE作为前线技术枢纽，完美消解这一矛盾，实现标准化技术与非标场景的精准适配。

4.1.2 痛点二：传统分工链路过长，AI落地效率极低、成本极高

传统AI项目落地链路：客户提需求→销售记录→产品拆解→研发排期→开发测试→交付部署→售后运维。整条链路涉及5-6个岗位，层层传递、层层损耗，需求偏差、沟通滞后、迭代缓慢是常态。

AI落地的核心要求是快速试错、快速迭代、快速验证，模糊需求、动态场景是常态，过长的分工链路会直接导致项目延期、烂尾、投入白费。

FDE的全链路闭环模式，彻底砍掉中间冗余环节，一人完成需求、开发、落地、迭代、反馈，大幅提升落地效率，降低企业沟通成本、人力成本、试错成本。

4.1.3 痛点三：落地经验无法沉淀，企业无法形成规模化壁垒

绝大多数AI企业的致命问题：每个客户落地都是“从零开始”，做完一个客户，经验全部流失，下一个客户重新踩坑、重新开发、重新适配，永远无法实现规模化复制，企业只能靠人力堆砌做项目，无法形成产品壁垒。

FDE的核心战略价值，就是将一次性项目交付，转化为持续性产品资产。通过前线落地沉淀通用能力、行业模板、适配方案，让企业产品越做越强、落地成本越来越低、复制速度越来越快，构建行业核心壁垒。

4.1.4 痛点四：AI项目重交付、轻价值，客户留存率极低

传统交付模式只关注“是否部署完成”，不关注“是否产生业务价值”。很多AI项目交付后，因为适配度低、使用复杂、无法解决核心痛点，被客户闲置、弃用，最终导致客户不续费、不增购，企业营收无法持续增长。

FDE以“业务价值落地”为核心目标，全程跟进使用效果、优化体验、深挖场景，确保AI能力真正赋能业务，大幅提升客户满意度、留存率、复购率，实现企业商业价值的持续闭环。

4.2 企业什么时候必须启用FDE ¶

FDE不是通用岗位，并非所有软件、互联网项目都需要。基于产品场景适配理论，我总结出6类必须启用FDE的核心场景，也是当前AI企业FDE扩招的核心赛道，精准区分“需要FDE”和“不需要FDE”的业务场景。

4.2.1 场景一：企业级AI、行业AI非标落地项目（核心刚需场景）

针对工业制造、金融、政务、医疗、能源、交通等传统行业的AI改造项目，行业壁垒高、业务流程复杂、系统老旧、数据非标、安全合规要求严苛，标准化AI产品完全无法适配，必须派驻FDE驻场攻坚，完成场景适配、定制迭代、价值落地。这是目前FDE最核心、最主流的应用场景。

4.2.2 场景二：前沿AI模型商业化落地项目

大模型、多模态AI、AI Agent、智能决策等前沿AI技术，目前行业无成熟落地标准、无固定使用场景、需求高度不确定。这类项目无法依靠后方标准化研发落地，需要FDE深入一线，探索落地场景、验证商业模式、沉淀落地标准，为产品规模化铺路。OpenAI、Anthropic大规模扩招FDE，核心就是服务这类前沿模型商业化落地。

4.2.3 场景三：高风险、高涉密、高安全要求的政企项目

政务、国防、军工、金融核心系统等涉密、高安全场景，数据无法外传、系统无法对外开放、远程调试受限，后方研发团队无法远程落地，必须由FDE现场驻场，在客户内网、涉密环境中完成部署、开发、调试、迭代，同时保障数据安全与合规。这也是Palantir首创FDE模式的核心场景。

4.2.4 场景四：新产品、新赛道的0-1探索落地项目

企业布局全新AI赛道、孵化全新AI产品时，没有成熟需求、没有落地案例、没有标准化方案，需要FDE深入一线，完成0-1场景探索、需求验证、模式打磨、案例沉淀，帮企业跑通商业模式和产品落地路径。

4.2.5 场景五：需要深度客户共创的战略级项目

企业核心大客户、战略合作伙伴的AI定制项目，需要深度绑定客户业务、共创解决方案、长期迭代优化，需要FDE长期驻场，全程对接、持续优化、深挖价值，打造标杆客户和标杆案例，带动行业规模化拓展。

4.2.6 场景六：AI项目落地失败率高、迭代卡顿的存量项目

企业原有AI项目频繁烂尾、落地效果差、客户投诉多、无法规模化复制，本质是缺失前线落地迭代能力，引入FDE可打通落地闭环，解决卡点问题，沉淀标准化能力，盘活存量项目。

4.2.7 无需FDE的场景（反向对标）

标准化C端互联网产品、通用SaaS软件、轻量化工具类产品，场景统一、需求标准、远程可迭代、落地无门槛，无需派驻FDE，传统研发+交付模式即可满足需求。

4.3 FDE的核心价值 ¶

结合产品价值分层理论，将FDE的价值分为三个维度，从表层落地价值到深层战略价值，全方位解释其高薪、稀缺、不可替代的核心原因。

4.3.1 表层价值：解决AI落地最后一公里，实现业务价值闭环

这是FDE最直观的价值：彻底解决AI项目“落地难、效果差、烂尾多”的行业痛点。传统模式下，AI技术只能停留在演示、试点阶段，无法融入真实业务流程；FDE通过现场驻场、非标适配、即时迭代，让AI系统真正跑在客户生产环境中，真正替代人工、提升效率、降低成本，完成从“技术demo”到“业务价值”的最终转化。

传统软件销售是“我给你软件，你自己用”。FDE模式是“我帮你把问题解决掉”。

对客户来说，FDE的价值在于：

更快的投产时间：FDE处理的集成工作，如果交给客户自己的工程团队，可能在内部 backlog 里排好几个月
更高的采用率：根据客户实际工作方式定制的部署才会被真正使用，通用配置往往被束之高阁
更低的部署风险：FDE在部署前就发现并解决环境特定的问题——认证缺口、数据治理缺口、合规障碍

4.3.2 中层价值：降低企业落地成本，提升商业化效率

1.降本：FDE一人承担多岗位复合工作，替代传统产品+交付+售前+运维的碎片化分工，大幅降低企业人力成本、沟通成本、试错成本；现场即时迭代，避免跨团队沟通损耗，减少项目延期、返工带来的资源浪费。

2.提效：砍掉冗余链路，实现需求到落地的快速闭环，AI项目落地周期缩短50%以上；同时通过经验沉淀，让后续同类项目落地无需重复造轮子，规模化落地效率持续提升。

3.增收：通过价值落地提升客户满意度，实现高留存、高复购、高扩容；同时挖掘客户潜在场景，拓展新的商业增长点，提升单客户ARPU值。

对提供FDE服务的公司来说，FDE的价值在于：

把一线洞察沉淀回平台产品，降低未来项目的定制化成本与时间。

Palantir的模式就是典型：FDE在空客解决了供应链数据整合的问题，这个解决方案后来被沉淀为Foundry平台的能力，再去服务航空公司、制造商时就不需要从头开始了。

这形成了一个正向循环：

4.3.3 深层价值：沉淀产品核心资产，构建企业长期壁垒

这是FDE最核心、最稀缺的战略价值，也是普通技术岗位不具备的价值。FDE是企业产品迭代的一线触角、行业数据的核心来源、落地方法论的沉淀载体。

每一次客户落地，都是一次产品迭代的机会。FDE将零散的非标需求、落地坑点、适配方案，抽象为通用功能、行业模板、标准化流程，让企业AI产品从“零散定制项目”升级为“可规模化复制的标准化产品”，持续积累产品资产、行业经验、客户口碑，构建同行无法复制的核心壁垒。

FDE不仅仅在“交付”，还在“发现”。他们工作直接嵌入客户工作流，发现真实的运营约束、组织动态和架构层面的瓶颈。

一个优秀的FDE会做五件事：

绘制客户的工作流，包括每一步背后的隐藏逻辑
发现隐藏的约束——技术的、组织的、流程的
教育和协调利益相关方，把模糊需求转化为清晰的问题陈述
构建能立即产生可衡量结果的“砂石路”方案
向客户展示产品能解锁的新可能性
将模式反馈给总部，由产品团队固化为可复用的功能

Think of an FDE as a product discovery loop embodied as a person——FDE就是“产品发现循环”的人格化存在。

4.3.4 行业价值：重构AI商业化落地模式，破解行业悖论

FDE模式的普及，彻底重构了AI行业的落地逻辑，让AI行业从“重技术研发、轻落地应用”的畸形模式，转向“技术研发+场景落地+价值迭代”的健康闭环模式，破解了“技术过剩、落地贫瘠”的行业悖论，推动AI技术真正赋能千行百业，实现产业数字化、智能化升级。

FDE这个角色对个人职业生涯也有独特价值。

一位Palantir的FDE这样描述自己的工作：“我选择成为FDE，是因为从创造方案到看到效果之间的反馈周期极短。”

这种“发现问题→设计方案→写代码→部署上线→看到效果”的完整闭环，让FDE积累了远超普通工程师的实战经验。难怪如此多的前FDE后来成为创始人：Kalshi的Tarek Mansour、Sourcegraph的Quinn Slack、Anduril的Matt Grimm等。

FDE的工作本质上是在经营一家“一人公司” ——你独立面对客户、定义问题、交付方案、对结果负责。

第五章真实案例深度拆解 ¶

5.1 案例一：金融银行业AI智能风控落地案例 ¶

5.1.1 项目背景：金融AI落地的极致痛点

银行业是AI商业化落地要求最高、壁垒最强、最能体现FDE核心价值的核心场景。相较于工业场景，银行AI落地具备强监管合规、数据高度敏感、业务逻辑极度复杂、系统架构老旧、容错率趋近于零五大核心特征，也是90%通用AI产品落地失败的重灾区。

本次案例以国内某股份制商业银行零售信贷智能风控AI升级项目为原型，完整还原行业真实落地现状。该银行原有传统信贷风控体系依赖人工审核+静态规则引擎，存在严重业务痛点：小额信贷申请量大、人工审核效率极低、夜间无人工值守导致用户体验差、传统规则无法识别新型欺诈行为，每年坏账损失超千万元。

该银行采购国内头部AI企业通用智能风控大模型系统，该模型在测试数据集、标准化金融样本中准确率达98.7%，理论上可完全替代传统人工风控审核、拦截欺诈订单、降低坏账率。

但项目初期进场后，传统研发、交付团队完全无法推进，项目直接停滞，濒临作废，核心落地痛点全部为非标场景适配问题，也是所有金融AI项目的共性难题：

1.数据非标与合规壁垒：银行核心信贷数据全部隔离在内网涉密环境，无法外传、无法对接云端模型，通用AI模型的云端调用架构完全无法适配银行本地化、私有化、内外网隔离的部署要求；同时金融数据受银保监会监管，数据脱敏、字段调用、日志留存有严格合规标准，通用产品无适配能力。

2.老旧系统兼容壁垒：银行核心信贷系统、用户征信系统、交易流水系统均为十年以上老旧架构，无标准化API接口，数据格式混乱、字段不统一，通用AI风控系统无法完成数据对接与业务联动。

3.业务规则高度非标：该银行针对小额信贷、青年客群、下沉客群有专属差异化风控规则，区别于行业通用标准，通用AI模型的风控判定逻辑过于标准化，无法适配银行差异化授信、拦截规则，直接落地会出现大量误判、错拒、错放问题，引发金融风险。

4.迭代链路过长，无法适配金融实时风控：传统后方研发迭代模式，一次规则优化、模型微调需要1-2个月，而金融欺诈手段实时迭代、新型诈骗套路层出不穷，缓慢迭代完全无法满足风控实时拦截需求。

5.业务价值无法落地：单纯部署通用模型只能实现简单数据展示，无法真正融入信贷审批流程、无法自动完成风控判定、无法降低坏账率，仅能作为演示工具，无实际业务价值。

5.1.2 FDE进场后的全流程攻坚

为盘活项目、规避金融风险、落地真实业务价值，AI厂商紧急派驻资深金融FDE专项团队长期驻场银行总行科技部门与零售信贷部，全程沉浸式攻坚6个月，完整落地FDE「需求校验-非标适配-工程落地-价值验证-产品沉淀」全链路能力，完美诠释FDE二元性核心逻辑，具体工作拆解如下：

第一步：沉浸式驻场调研，拆解金融真伪需求（产品核心能力）

区别于传统售前远程对接、看文档判需求的模式，FDE团队全程驻场，跟随银行信贷审核员、风控专员、合规专员轮岗办公，吃透银行完整信贷业务流程、监管规则、风控痛点，完成需求去伪存真：

剔除客户表层伪需求：银行初期要求“提升模型整体准确率”，FDE调研后判定该需求无业务价值，金融风控核心不是通用准确率，而是降低优质用户错拒率、提升恶意欺诈拦截率，双指标平衡才是核心；
挖掘隐性刚需需求：除基础风控审核外，银行存在合规日志自动留存、风控报表自动生成、异常订单溯源、新型欺诈规则快速迭代、监管报备数据自动输出五大隐性刚需；
对齐多方预期：同步银行业务部门、科技部门、合规部门三方诉求，明确落地核心指标：坏账率下降15%以上、人工审核效率提升60%、合规零违规、无批量误判事件。

第二步：金融非标场景定制适配，破解落地壁垒（工程核心能力）

针对银行私有化部署、老旧系统、合规管控、差异化风控四大非标痛点，FDE团队现场自主开发、快速迭代，无需后方研发排期，当场解决所有落地卡点：

私有化合规改造：彻底重构模型调用架构，摒弃云端部署模式，完成AI风控模型全本地化内网部署，搭建银行专属数据脱敏模块、操作日志全留存模块、权限分级管控模块，100%满足银保监会金融数据合规要求；
老旧系统打通适配：自主编写适配代码，开发银行老旧信贷系统专属兼容接口，完成多源异构数据清洗、字段统一、格式适配，实现AI风控系统与银行核心业务系统的无缝联动，做到信贷申请数据实时同步、风控结果实时回传；
差异化风控规则定制：基于银行专属客群规则，现场微调模型权重、新增自定义风控策略引擎，支持银行灵活配置差异化授信、拦截、复核规则，解决通用模型“一刀切”的误判问题，适配银行个性化业务逻辑；
实时迭代攻坚：针对市场新型金融诈骗套路，FDE现场快速更新风控特征、迭代判定规则，实现24小时内完成策略上线，彻底解决传统研发迭代滞后、风控失效的行业痛点。

第三步：业务赋能与可量化价值验证（运营落地能力）

FDE不只是完成系统部署，更对银行最终业务结果全权负责，落地后持续做业务赋能与效果优化：

分层赋能培训：针对风控专员开展AI系统操作、异常订单处理、规则配置培训，针对合规部门开展合规校验、日志核查培训，搭建标准化使用流程；
全维度数据监控：实时监控信贷审核效率、订单通过率、欺诈拦截率、坏账率、误判率等核心业务指标，每日输出AI风控运行报表；
精细化效果优化：针对初期少量误判订单，现场复盘优化模型策略与规则，持续打磨适配效果，保障系统稳定运行。

最终落地量化成果：银行零售小额信贷人工审核效率提升68%，夜间无人值守场景实现全自动审核，用户授信时长从平均5分钟缩短至10秒；恶意欺诈订单拦截率提升22%，年度信贷坏账损失下降18.3%，直接为银行节约千万级成本，且全程合规零违规。

第四步：金融行业能力沉淀与产品反哺（FDE核心战略价值）

项目稳定落地、价值验证完成后，FDE团队并未离场，而是完成单次项目经验向企业通用产品资产的转化，这也是FDE区别于普通交付工程师的核心壁垒：

沉淀金融合规通用模块：将银行私有化部署、数据脱敏、日志留存、权限管控、监管适配的非标方案，抽象为银行业AI合规标准化组件，可一键复用至所有银行、金融机构项目；
沉淀老旧金融系统适配模板：整理国有行、股份行、城商行老旧核心系统的对接方案、数据清洗规则，形成标准化适配工具包；
沉淀金融风控策略引擎：将银行差异化客群风控规则、动态迭代逻辑，打磨为可配置、可自定义的通用风控策略引擎，纳入企业核心AI风控平台；
输出金融AI落地方法论：形成《银行业AI风控落地合规手册》《金融非标系统适配避坑指南》，搭建金融行业落地知识库。

通过FDE的沉淀赋能，该AI厂商后续服务全国数十家银行的AI风控项目，落地周期从原来的4-6个月缩短至1个月，落地成本降低75%，同时产品合规性、适配性行业领先，快速构建起金融AI落地的核心竞争壁垒。

5.1.3 金融案例深度复盘：全方位看懂FDE的不可替代性

通过本次银行AI风控落地案例，可精准区分FDE与传统研发、交付、实施岗位的核心差异，彻底印证FDE在高端商业AI场景的刚需价值：

1.传统交付团队只能“搬砖落地”：仅能完成标准化部署，无法处理金融合规非标、老旧系统兼容、差异化业务规则，最终只能交付一个无法使用的演示系统，无法产生业务价值；

2.后方算法/研发团队只能“做标准化产品”：不了解金融监管规则、银行非标业务逻辑，打磨的通用模型无法适配真实金融场景，远程迭代完全跟不上金融风控的实时需求；

3.FDE实现“落地+价值+迭代+沉淀”四重闭环：既现场解决了银行所有落地卡点、落地真实降本增效价值，又将单次定制化落地经验转化为企业可复用的产品能力，帮企业实现从“单一项目交付”到“行业规模化复制”的升级。

同时，该案例完美印证前文的产品理论：AI行业的核心壁垒从来不是技术模型，而是场景适配与落地能力。同等模型精度下，有无FDE团队，直接决定AI项目是“摆设演示版”还是“盈利工具”，这也是头部金融AI企业高薪扩招FDE的核心原因。

本次项目最终价值总结：

客户（银行）层面：在满足强监管合规的前提下，实现信贷风控智能化升级，大幅提升审核效率、降低坏账损失、优化用户体验，落地实打实的金融业务价值；
企业（AI厂商）层面：打造金融行业标杆案例，沉淀专属金融AI标准化产品能力与落地方法论，大幅降低后续行业落地成本，构建行业壁垒；
产品层面：完善金融行业AI产品矩阵，补齐合规适配、老旧系统兼容、差异化风控迭代能力，实现产品从“通用化”到“行业专业化”的升级。

5.2 案例二：工业制造业AI质检落地案例 ¶

5.2.1 项目背景

某大型汽车零部件工厂，传统零部件质检完全依赖人工，存在效率低、漏检率高、人工成本高三大痛点，常年困扰生产。某AI企业为其提供AI视觉智能质检系统，核心能力是通过图像识别，自动检测零部件划痕、变形、瑕疵等缺陷。

项目初期，企业后方算法团队打磨的模型，实验室检测准确率可达99.2%，理论上完全可以替代人工。但项目进场后，立刻出现大面积落地问题，传统交付团队无法解决，项目濒临烂尾：

工厂生产环境光线复杂、设备震动、拍摄角度不固定，导致现场采集的图片和实验室数据集差异极大，模型现场准确率暴跌至85%，无法满足生产标准；
工厂老旧生产系统无标准化接口，AI系统无法对接生产流水线，无法实现自动化采集、自动化质检、数据同步；
工厂质检流程有特殊校验规则、分级标准，标准化AI系统无法适配，和员工工作流程冲突；
传统交付团队只能完成基础部署，无法现场修改代码、适配场景、优化模型调用逻辑，反馈后方研发迭代需要1-2个月，项目严重卡顿。

5.2.2 FDE进场后的全流程攻坚（核心工作还原）

企业紧急派驻10人FDE团队驻场9个月，全程扎根生产车间，完成全链路落地攻坚，完整还原FDE二元性工作模式：

第一步：沉浸式需求挖掘与问题定位（产品能力）

FDE团队全程跟随车间工人轮岗，吃透生产流程、质检标准、环境痛点，精准定位落地失败的核心原因：不是模型精度不足，而是场景适配缺失、数据环境不符、流程不匹配。同时梳理出客户未明确表达的隐性需求：AI质检数据需同步至生产管理系统、支持缺陷分类统计、生成质检报表、适配多型号零部件检测。

第二步：现场非标适配与快速迭代（工程能力）

现场优化图像采集逻辑：根据车间光线、设备环境，调整拍摄参数、图像预处理算法，适配现场复杂场景；
开发专属对接接口：自主编写代码，打通AI系统与工厂老旧生产设备、管理系统的数据壁垒，实现自动化采集、质检、同步；
定制化流程改造：根据工厂质检分级规则，优化AI缺陷判定逻辑、新增自定义校验模块，适配车间真实工作流程；
现场数据迭代：采集车间真实生产数据，筛选有效样本，同步后方算法团队微调模型，将现场检测准确率从85%提升至99.1%，达到生产标准。

第三步：业务赋能与价值验证（运营能力）

FDE团队对车间员工开展分层培训，梳理标准化使用流程、异常处理流程；实时监控系统运行数据，统计核心价值指标：质检效率提升70%、漏检率从3.2%降至0.3%、每年为工厂节省人工成本超200万，完成业务价值闭环验证。

第四步：经验沉淀与产品反哺（战略能力）

这是FDE区别于传统交付的核心动作。项目落地稳定后，FDE团队没有直接离场，而是完成核心经验沉淀：

提炼工业车间复杂光线、震动环境下的图像采集通用方案；
开发老旧工业系统通用适配接口模板；
沉淀零部件质检分级自定义配置模块；
输出工业AI质检落地避坑指南、标准化流程手册。

所有沉淀的能力全部同步至后方产品、研发团队，迭代进企业核心AI质检平台，成为标准化通用能力。后续企业服务同类工厂客户时，无需重新开发、重新适配，落地周期从3个月缩短至2周，落地成本降低80%，实现规模化复制。

5.2.3 项目最终价值总结

客户层面：成功落地AI智能质检系统，实现降本增效，完成智能化升级；
企业层面：交付成功标杆案例，提升客户粘性，同时沉淀通用产品能力，打通行业规模化落地路径；
产品层面：完善工业AI质检产品矩阵，优化场景适配能力，构建行业竞争壁垒。

5.3 案例三：Palantir国防涉密AI落地案例 ¶

5.3.1 项目背景

Palantir早期核心业务是为美国国防、情报机构提供AI数据分析、智能决策系统，服务反恐、反扩散、后勤保障等涉密 mission 场景。这类场景具备三大核心特点：需求高度非标、数据高度涉密、任务动态变化、无通用落地标准。

传统远程研发、交付模式完全无法适配：涉密数据无法外传、业务需求无法书面化、任务场景动态调整、系统适配要求极高。为此Palantir首创Echo-Delta双团队FDE模式，成为全球FDE模式的标杆。

5.3.2 Echo-Delta核心协同模式

1.Echo团队（业务洞察团队）：由行业 domain 专家组成，深耕国防业务场景，熟悉现有工作流程、核心痛点、任务短板，负责挖掘真实业务需求、对齐高层预期、筛选高价值落地场景，解决“做什么”的问题；

2.Delta团队（FDE前线工程团队）：即核心FDE团队，负责快速搭建落地原型、现场开发适配、迭代优化、系统部署，在真实涉密场景中验证解决方案，快速跑出“可用、有效”的落地雏形（行业称为“砂石路方案”），解决“怎么做”的问题；

3.后方核心产品团队：观察前线FDE落地的砂石路方案，提炼通用逻辑、共性能力，将一次性定制化方案，打磨为标准化、可复用的产品能力（行业称为“柏油路能力”），实现从单点落地到规模化产品的升级。

5.3.3 FDE落地核心价值

解决涉密场景落地难题：FDE驻场涉密环境，现场完成系统开发、部署、迭代，规避数据外传风险，满足安全合规要求；
适配动态非标需求：针对反恐、后勤等不同mission的差异化需求，现场快速迭代方案，适配动态变化的业务场景；
构建产品迭代闭环：将国防场景的非标落地经验，持续沉淀为Palantir核心平台Foundry的通用能力，让产品适配更多政务、国防细分场景，成为全球政企AI落地的标杆产品。

5.4 案例四：某科技公司在某国有银行海外分行实施案例 ¶

5.4.1 背景：一家出海银行的三重困境

某国有银行在东南亚开设了海外分行。业务做得不错——跨境中小企业开户、资产托管、投资理财、债券交易——但IT系统已经快撑不住了。

第一个困境：技术债堆积如山

分行存量的IT系统里，有大量基于VB（Visual Basic）开发的老旧模块，搭载着超长的SQL存储过程。业务文档早就丢失了，代码没有统一规范。整个分行有70多套相互孤立的系统，像70多个互不通话的孤岛。

第二个困境：业务全链路低效

债券交易、跨境开户、理财审批这些核心业务流程，大量依赖邮件和电话线下流转。业务部门和技术部门之间信息不通，需求传递靠口口相传，交付周期拖得很长。

第三个困境：双重合规压力

这家分行要同时满足两套监管要求——总行的集团管控和海外属地的金融监管。两套标准对操作留痕、权限管控、变更追溯都有严格要求。传统的“按人天计费”外包模式，交付权责模糊，根本无法同时满足两套标准。

传统IT服务商的解决方案是什么？派几个人过去，按天收钱，写完代码就走。但这家银行需要的不是一个“写代码的人”——它需要的是能同时搞定技术、业务、合规的复合型团队。

5.4.2 FDE进场：不是“派人”，是“共建”

某科技公司拿下了这个项目，但用的不是传统外包模式。他们采用了一套组合方案：

FDE深度驻场：组建专属FDE团队长期驻扎客户现场
AI全链路赋能：用AI工具覆盖需求梳理、代码开发、自动化测试全流程
双监管底层架构：一套数据底座同时输出两套监管报表

项目从2024年3月启动，历时两年。

5.4.3 FDE团队做什么？——不是“写代码”，是“做业务”

这家科技公司的FDE团队不是传统意义上的“驻场程序员”。他们兼具三重能力：

业务顾问的能力：理解银行的业务逻辑——跨境开户的流程是什么、债券交易的审批链条有多长、理财产品的合规要求有哪些。FDE团队要能跟业务人员用同一个语言对话。

全栈研发的能力：能写代码、能重构架构、能对接70多套遗留系统、能做云迁移。

产品沉淀的能力：不只是一个项目做一个方案，而是把每次解决的问题沉淀为可复用的标准化组件。

最关键的差异在这里：传统外包是“客户说要什么，就做什么”。FDE团队是“主动挖掘业务痛点、输出定制化解决方案”。

5.4.4 结果：两年交出什么成绩单？

系统层面：

2024年3-10月，FDE团队切入7大业务领域，承接12个关键系统开发
2024年12月，办公区扩容，增设120个专属工位及独立安全运维中心
2024年底，团队规模达到74人
2025年全年，累计开发项目超过50个，累计代码量超过1500万行，上云重构项目16个

架构层面：

将传统整体式开发模式升级为松耦合的微服务架构
搭建了海外系统与总行新一代组件的融合桥梁
以一套数据底座同时输出总行、海外属地两套监管报表

9套核心系统（企业应用、电子渠道、投资理财等）顺利上线。强积金系统、网上银行、债券交易系统等核心应用不仅提升了业务处理效率，还强化了系统安全性与可扩展性。

第六章总结：FDE的本质与未来 ¶

6.1 用一个类比理解FDE ¶

如果把AI公司比作一家武器制造商：

AI研究员：研究更先进的火药配方
软件工程师：在工厂里造标准化的武器
销售：向客户展示武器有多厉害
FDE：跟着部队上前线，根据实际战场情况调整武器、教士兵使用、现场修理故障

AI能力再强，如果没人把它送到战场上、让它真正发挥作用，它就只是一堆躺在实验室里的技术。FDE就是那个把技术送上战场的人。

6.2 FDE的核心本质 ¶

回顾全文，FDE的本质可以归结为三点：

第一，FDE是“结果导向”的工程师。他不对自己的代码量负责，不对交付的功能数量负责，只对客户的业务结果负责。

第二，FDE是“闭环思维”的工程师。他不只是把东西做出来就结束，而是要把一线的经验教训反馈回产品，让产品变得更好。

第三，FDE是“两栖作战”的工程师。他既要有扎实的工程能力（能写生产级代码），又要有业务理解力（能跟非技术主管沟通），还要有在模糊环境中独立推进的能力。

企业愿意为FDE支付高昂溢价，正是因为这种“两栖作战”能力极为稀缺。

6.3 FDE的未来 ¶

FDE正在从一个独立岗位变成ToB团队的标配能力。杭州已有大厂推行“新人先做三个月FDE”的制度，甚至直接派驻产品经理驻场做需求调研与Demo开发。

在AI技术快速迭代的背景下，客户对“既懂技术、又懂业务”的现场工程支持需求正在急剧增长。FDE这个角色，不会消失，只会越来越重要。

但需要警惕的是：FDE不是一个可以随意复制的模式。它需要高价值客户、成熟的产品平台、以及真正具备“两栖作战”能力的人才。盲目跟风招FDE，只会变成一场昂贵的错误。

6.4 给AI小白的建议 ¶

如果你是一名AI小白，想理解FDE，记住这三句话就够了：

FDE是“带着工具箱上前线的工程师” ——不在总部写通用代码，而是到客户现场解决问题
FDE对结果负责，不对代码量负责 ——目标是让客户真正用起来、产生价值
FDE是“产品发现循环的人格化” ——在一线发现问题、解决问题、把经验带回产品

如果你正在考虑是否要从事FDE这个职业，问自己三个问题：

我能不能接受长期在客户现场工作（而非坐在总部舒适的办公室里）？
我能不能在信息不充分、目标模糊的情况下推进工作？
我是不是既有工程能力，又愿意深入理解业务？

如果三个答案都是“是”，FDE可能是最适合你的职业之一。

6.5 结语 ¶

FDE的本质，是在最复杂的环境中、用最短的反馈周期、把最模糊的问题变成可运行的系统、再把现场的经验变成可复用的产品能力。

AI行业已经彻底告别“拼技术、拼模型参数”的野蛮生长时代，进入拼落地、拼场景、拼价值、拼规模化的精细化竞争时代。

过去，企业的核心竞争力是“能不能做出先进的AI技术”；未来，企业的核心竞争力是“能不能把先进的AI技术落地、用好、做出价值、规模化复制”。而FDE，正是掌握这一核心竞争力的关键岗位。

对于AI小白而言，看懂FDE的本质，就看懂了AI行业的底层商业逻辑：技术是基础，落地是核心，价值是终极目标。

FDE的崛起，不是岗位的迭代，而是整个AI行业从“技术噱头”走向“产业实用”的标志性变革。掌握FDE的核心认知，才能真正看懂AI商业化的未来趋势，规避行业认知误区，精准把握AI时代的职业与商业机遇。

在AI时代，FDE这个角色只会越来越重要。大家加油！

阅读全文

文库	化学体系	每基因 sgRNA	特点
Brunello	SpCas9 KO	4	Doench Rule Set 2 打分，on-target 效率高
Avana	SpCas9 KO	6	DepMap 使用，更密覆盖
TKOv3	SpCas9 KO	4	加入了增强子靶向 sgRNA
Dolcetto	CRISPRi	6	不切割 DNA，适合剂量敏感基因

指标	参考值	说明
`min_genes`	200	低于此多为空液滴/碎片
`log1p_total_counts`	5 MAD	在对数尺度过滤（深度右偏）
`pct_counts_in_top_20_genes`	5 MAD	高值标记低复杂度/垂死细胞
`pct_counts_mt`	3 MAD + 硬 >8%	组织依赖：心肌/肝细胞/肌肉天然 20-40% 且健康；细胞核 ~0-2%
`min_cells`（基因）	3	去掉稀有基因

方法	模型/假设	何时用	何时失效
shifted-log (CP10k)	常数总 mRNA	通用默认	组成差异大的细胞类型（浆细胞、增殖细胞）
scran deconvolution	pooled size factors	低深度、高 dropout	R-only，需预聚类
sctransform v1/v2	NB 正则回归	Seurat 深度去除	慢，v1 过拟合
Pearson residuals	解析残差	UMI 计数	实验性

细胞类型	Markers
CD4 T	CD3D, CD4, IL7R
CD8 T	CD3D, CD8A, CD8B
B	MS4A1, CD79A, CD19
CD14 Mono	CD14, LYZ, S100A8
DC	FCER1A, CST3

方法	模型/假设	使用场景	失效时
PAGA	cluster 图连通性	判断连续谱是否存在（强制第一步）	只给连通性不给定方向
Slingshot	簇中心 MST + 主曲线	已知树/分支拓扑	依赖输入聚类
Monocle3	反向图嵌入主图	树/环形/断连拓扑	在 UMAP 空间学习，继承畸变
DPT	可逆扩散随机游走	线性/简单分支	根敏感
Palantir	有向马尔可夫链	分支命运	根敏感
CellRank 2	不可逆马尔可夫链 + GPCCA	有向命运映射	kernel 质量依赖

混杂因素	如何制造虚假交互	缓解
Ambient RNA	高表达分泌基因泄漏到每个 cluster	先 SoupX 去污染
细胞丰度	大 cluster 给更紧凑 null → 成为网络枢纽	下采样，检查 counts vs 大小
测序深度	深度差异混淆技术与生物	整合计数，先归一化深度
解离应激	FOS/JUN/HSP 是真实配体	打分应激模块

步骤	技能名	工具	核心作用
1	`bio-single-cell-data-io`	Scanpy	数据读取
1-2	`bio-single-cell-preprocessing`	Scanpy/Seurat	质控、归一化、HVG
3	`bio-single-cell-batch-integration`	Harmony/scVI	批次整合
4	`bio-single-cell-clustering`	Scanpy/Seurat	降维+聚类
5	`bio-single-cell-markers-annotation`	Scanpy/Seurat	Marker+注释
6	`bio-single-cell-trajectory-inference`	PAGA/Slingshot/Monocle3/CellRank	轨迹+伪时间
7	`bio-single-cell-cell-intercommunication`	LIANA/CellChat/NicheNet	细胞通讯

文章归档

最新文章