如何Vibe Coding生信分析：以scRNAseq为例

这次用自然语言驱动完成了一轮单细胞转录组分析。整个过程不是把命令丢给模型生成完就算结束，而是把需求、规则、检查点和产物都固定下来：数据怎么导入、哪些细胞留下、参数为什么这样选、每个 cluster 怎么注释、最后输出哪些对象和图。回头看，它更像在搭一条可以解释、可以复跑的分析流水线。

1. 先立规矩，再跑分析

Prompt：请检查我的单细胞数据结构、样本信息和现有文件；为 10x 数据、TCR/BCR 文件与 canonical sample IDs 建立显式映射；所有脚本、日志、中间对象和图表都保存到独立输出目录；不要覆盖原始数据；在执行分析前先说明数据检查结果和目录结构。

单细胞分析最容易失控的地方不是某个函数报错，而是每一步参数都“看起来合理”，最后却很难解释为什么这样选。所以开工前先约定几件事：

所有日志、脚本、中间对象和图表都写进独立输出目录，不覆盖原始数据；
每个样本的 10x 目录、TCR/BCR 文件和 canonical sample ID 先建立显式映射；
QC、PCA、resolution 都要看分布和结果再决定，不套固定阈值；
注释必须回到 marker gene 和生物学背景，不只依赖聚类编号。

导入环节的核心结构其实很短：

sample_map <- c(
  sample_1 = "path/to/sample_1",
  sample_2 = "path/to/sample_2"
)

objects <- lapply(names(sample_map), function(sample_id) {
  counts <- Read10X(sample_map[[sample_id]])
  obj <- CreateSeuratObject(counts, project = sample_id)
  obj$sample_id <- sample_id
  obj
})

TCR/BCR 这类辅助数据没有立即用于聚类，但同样保留 barcode 和 sample 的映射。这样后续想回溯克隆型信息时，不会因为中途转换对象或合并数据而失去追踪。

2. QC 不是设阈值，而是看分布

Prompt：请为每个样本计算并可视化 nFeature、nCount 和 mitochondrial percentage；根据实际数据分布推荐 QC thresholds；输出过滤前后统计、QC plots 和过滤依据；不要使用固定阈值。

QC 我没有直接写死 nFeature、nCount 和 mitochondrial percentage 的阈值，而是先对每个样本做分布检查，再结合整体数据形态设定过滤标准。这样可以避免把某个组织或批次本身的生物学差异误删掉。

obj <- PercentageFeatureSet(obj, pattern = "^mt-", col.name = "percent.mt")

VlnPlot(
  obj,
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
  group.by = "sample_id",
  ncol = 1
)

最终保留了约 2.9 万个细胞。这个数字本身不是目标，而是 QC 后的自然结果。重要的是过滤前后都留下统计表和图，能看清楚每个样本损失了多少细胞、高线粒体比例细胞集中在哪些批次。

QC violin plots

3. 整合的关键是去批次，不是抹平差异

Prompt：请对过滤后的样本执行 normalization、HVG 和 PCA；使用 Seurat v5 进行多样本整合；系统比较不同 PCA dimensions 和 clustering resolutions；输出 ElbowPlot、不同 resolution 的比较和最终 UMAP；根据聚类稳定性、UMAP 结构和生物学合理性选择参数，避免 over-correction。

多样本整合时，我更关心的是能否保留真实生物学信号。流程上先分样本做 normalization、HVG 和 PCA，再用 Seurat v5 的整合层做 batch correction：

obj_list <- SplitObject(obj, split.by = "sample_id")

obj_list <- lapply(obj_list, function(x) {
  x <- NormalizeData(x)
  x <- FindVariableFeatures(x)
  x <- ScaleData(x)
  RunPCA(x)
})

obj <- merge(obj_list[[1]], obj_list[-1])
obj <- IntegrateLayers(
  obj,
  method = CCAIntegration,
  new.reduction = "integrated.cca"
)

obj <- RunUMAP(obj, reduction = "integrated.cca", dims = 1:30)
obj <- FindNeighbors(obj, reduction = "integrated.cca", dims = 1:30)

选 dimensions 时看 ElbowPlot 和 PC 解释率，也检查不同 PC 数对 UMAP 和聚类结构的影响。选 resolution 时更直接：跑一组 resolution，比较 cluster 数量、UMAP 形态和小 cluster 是否有稳定 marker 支持。这里最终使用 30 PCs 和 resolution 0.4，得到 26 个 clusters。

ElbowPlot

4. 注释前先让 marker gene 说话

Prompt：请计算每个 cluster 的 marker genes，输出 Top 10 markers 和完整 marker 表；结合 canonical markers 与必要生物学背景进行 cell-type annotation；先将相近 clusters 合并成更大的细胞分类；输出 annotation 表、最终 UMAP 和更新后的 Seurat object。

聚类只是给出结构，注释还是要回到表达证据。每个 cluster 都输出 Top 10 markers，再结合 canonical marker 和组织背景判断身份。为了避免条目过碎，我再把相近细胞群合并成大类，例如把 macrophage 和 DC 归入 Macrophage/DC，把 T cell 和 NK 归入 T/NK。

markers <- FindAllMarkers(
  obj,
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25
)

top_markers <- markers %>%
  group_by(cluster) %>%
  slice_max(avg_log2FC, n = 10)

celltype_map <- c(
  "0" = "Macrophage/DC",
  "1" = "Fibroblast",
  "2" = "T/NK"
)

obj$celltype_coarse <- celltype_map[as.character(obj$seurat_clusters)]

最终把 26 个 clusters 归并成 10 个大细胞类型。这一层 annotation 就是后续 CellChat 分析的基础，也让结果图比直接展示 cluster 编号更容易读。

Final annotated UMAP

5. CellChat：把细胞类型放回通信网络

Prompt：请使用注释后的合并细胞大类运行 CellChat；保留 cell、sample、barcode 与 TCR/BCR 的映射关系；输出 interaction strength、interaction number、主要 signaling pathway 热图和通信数据表；必要时比较不同分组间的通信模式。

注释完成后，用大类做 CellChat 分析，而不是直接用 26 个小 cluster。这样能降低零碎群体带来的噪声，也更容易观察主要细胞群之间的通信模式。

cellchat <- CreateCellChat(obj, group.by = "celltype_coarse")
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
cellchat <- computeCommunProb(cellchat)
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- aggregateNet(cellchat)

这次数据里共识别出 50 多条 signaling pathway 和 1200 多条 communication。分析输出两类总览图：interaction strength 和 interaction number。前者看通信权重，后者看通信事件数量，两者结合可以避免只凭一个指标过度解读。

CellChat interaction strength

CellChat interaction number

通路层面的热图进一步展示了不同细胞类型之间哪些信号更强、哪些通路集中在一组 sender-receiver 关系里。

netVisual_heatmap(cellchat, signaling = pathway_name)

6. 这套 vibe coding 流程真正省掉的是什么

Prompt：请总结本次单细胞分析：数据导入情况、QC 标准与结果、整合方法、最终 dimensions/resolution 及选择依据、最终 clusters 数量、每个 cluster 的 annotation 与主要 markers、CellChat 主要结论、输出文件位置和可复现实验环境。

省掉的不是判断，而是重复劳动。自然语言先把意图和约束说清楚，模型负责展开脚本；真正需要人工介入的地方变成数据检查、参数选择和生物学解释。整个过程中最有价值的不是某一行代码，而是下面这些产物：

可复跑的 R scripts；
QC 前后的统计表和图；
PCA、UMAP、dimension 与 resolution 的比较结果；
cluster marker genes 和 annotation 表；
最终 Seurat object 与 CellChat object；
每一步的分析日志。

换句话说，vibe coding 的关键仍然是把分析决策写清楚。模型可以很快生成代码，但 QC 为什么这么过滤、整合为什么不能过度校正、cluster 为什么这样命名，这些还是必须由数据和生物学判断来回答。