企业官网建设流程全解析

写单细胞数据整合代码，这件事最容易被低估的其实是开头那两步：你以为要解决的是算法选型，实际上先卡住你的是数据格式和批次效应的判断。我在跑一个多来源的小鼠免疫细胞数据时，第一批次和第二批次无论是细胞类型组成还是测序深度都差得很远，直接丢进整合管线出来的结果非常难看。后来把“用什么思路整合”“在哪个层面整合”这两个问题先想清楚，代码写起来就顺多了。

这篇内容适合正在处理多批量单细胞数据、准备把 Seurat 或 Scanpy 流程真正跑起来的人。不管你是想校正批次效应，还是想在多来源数据里做细胞类型注释，甚至只是想把整合结果输出成一张能放进论文的 UMAP，下面这些代码和思路都能直接用。我尽量把每一步为什么这么做讲明白，也把我踩过的坑原样写出来。

1. 动手写整合代码前，先把“要整合什么”这件事想清楚

1.1 批次效应、技术差异与真实生物学信号的边界

单细胞数据整合，本质上是区分“技术带来的差异”和“生物学本身的差异”。同一个细胞类型，在不同文库制备批次里，检出基因数、线粒体比例、测序饱和度都可能不同，这些技术因素会掩盖真实的细胞状态差异。整合代码要做的不是把所有差异都抹平，而是只抹掉和批次标签相关的技术成分，尽量保留细胞类型、发育轨迹、疾病状态这些生物学信号。

我自己的判断标准很简单：如果某个差异可以用批次标签完全解释，而和细胞类型没有明确对应关系，那它大概率是技术噪声。比如某次整合两组数据，一组的 CD4 T 细胞整体高表达某些核糖体基因，另一组总体表达量偏低，但两组里都存在 naive 和 memory 亚群，这种整体偏移就是典型的批次效应。相反，如果只在某一组里出现了明确的罕见细胞亚群，这个差异就必须保留。

所以在写任何整合代码之前，我会先画三张图：按批次着色的 UMAP、按主要细胞类型着色的 UMAP、每个批次的质控指标分布。如果批次之间的聚类完全分离，且细胞类型混杂在不同批次里，说明整合是必要的。如果每个批次本身就有独特且稳定的细胞类型结构，这时候强行整合反而会破坏真实信号。这一步判断花的十分钟，通常能省掉后面调试整合脚本的半天时间。

1.2 数据输入形态：从表达矩阵到 AnnData / Seurat 对象

整合代码的第一步是把原始数据处理成算法能吃的标准对象。不同工具链的入口不一样：Scanpy 体系通常用anndata，Seurat 体系用Seurat对象。我从一开始就会把所有批次统一成同一种输入格式，避免后面反复转换。

import scanpy as sc import anndata as ad import numpy as np # 假设每个批次的表达矩阵是 csv 或者 h5ad # 我习惯先把每个样本读进来，加上 batch 标签，再 concat adata_list = [] for batch_name, path in [ ("batch1", "data/batch1.h5ad"), ("batch2", "data/batch2.h5ad"), ("batch3", "data/batch3.h5ad"), ]: adata = sc.read_h5ad(path) adata.obs["batch"] = batch_name adata_list.append(adata) adata = ad.concat(adata_list, join="outer", label="sample", index_unique="-") adata.obs["batch"] = adata.obs["batch"].astype("category") print(adata.shape)

这里有一个关键细节：join="outer"是取所有基因的并集，缺失的基因位填 0；join="inner"则只保留共有基因。大部分整合算法要求基因特征一致，如果用 inner join，可能丢掉一些低表达但重要的标记基因。我通常先用 outer，再做一次基础的过滤，而不是一开始就 inner。

R 的 Seurat 流程则简单一些，直接merge多个 Seurat 对象，通过add.cell.ids区分来源：

library(Seurat) obj_list <- list() obj_list[["batch1"]] <- Read10X("data/batch1/") |> CreateSeuratObject(min.cells = 3, min.features = 200) obj_list[["batch2"]] <- Read10X("data/batch2/") |> CreateSeuratObject(min.cells = 3, min.features = 200) merged_obj <- merge(obj_list[[1]], y = obj_list[2:length(obj_list)], add.cell.ids = names(obj_list)) merged_obj$batch <- gsub("_.*", "", colnames(merged_obj))

很多新手会在这里直接进入归一化和高变基因选择。但我建议先完成一个基本质控：过滤低质量细胞、剔除线粒体基因比例过高的细胞、去掉在极少细胞中表达的基因。这一步做不好，后面所有整合代码都会把“死细胞”或“空液滴”当成一类新的细胞群来校正。

2. 基于锚点的整合：Seurat CCA 流程的代码落地与细节

2.1 为什么要先做 SCTransform 再做整合

Seurat 的经典整合路线是：先对每个样本分别做 SCTransform 归一化，再执行FindIntegrationAnchors，最后用IntegrateData生成一个批次校正后的表达矩阵。这里我特别强调 SCTransform 要先在每个样本内部单独跑，而不是在合并后的对象上跑。因为 SCTransform 本质上是在拟合每个细胞的技术噪声模型，如果跨样本一起跑，模型会把批次间的表达量差异误当成基因表达特征。

# 对每个批次单独做 SCTransform for (b in unique(merged_obj$batch)) { obj_sub <- subset(merged_obj, subset = batch == b) obj_sub <- SCTransform(obj_sub, vars.to.regress = c("percent.mt"), verbose = FALSE) obj_list[[b]] <- obj_sub }

如果样本数量多，这一步会非常慢，可以考虑开启BPPARAM并行参数。但要注意，不是所有电脑都能顺利并行跑 SCTransform，小数据集反而串行更快。我的经验是，单样本细胞数在 5 万以下时，并行提升有限，还容易因为内存不够崩掉。

2.2 锚点集成代码的完整结构与参数语义

做完 SCTransform 后就可以寻找锚点了。FindIntegrationAnchors的核心逻辑是：先在降维空间里找到不同批次之间互相匹配的细胞对（锚点），再根据这些锚点计算校正向量。锚点数量由k.anchor和k.filter控制。k.anchor越大，每对批次之间的候选锚点越多，整合越平滑，但也越容易把不相似的细胞强行拉在一起；k.filter用于过滤掉那些可信度低的锚点。

anchors <- FindIntegrationAnchors( object.list = obj_list, anchor.features = 3000, reduction = "cca", k.anchor = 5, k.filter = 200 ) integrated <- IntegrateData( anchorset = anchors, new.assay.name = "integrated", k.weight = 50 )

跑完IntegrateData后，默认的整合表达矩阵存储在integrated[["integrated"]]这个 assay 里。很多人在这一步之后直接用原始 RNA assay 做后续分析，导致看到的 UMAP 根本没有批次校正效果。正确做法是后续 PCA、聚类、UMAP 都指定DefaultAssay(integrated) <- "integrated"，只在找标记基因或者做表达可视化时切回RNAassay。

2.3 锚点整合里最容易翻车的三个点

第一，锚点太少。如果两个批次之间细胞类型组成差异很大，或者技术差异过于剧烈，FindIntegrationAnchors会提示锚点数量不足，整合结果几乎等于没整合。此时可以适度调大k.anchor到 10 或 20，但更推荐的做法是先用 Harmony 这类线性校正方法做一个快速整合，再在这个基础上做锚点整合，而不是死磕参数。

第二，IntegrateData内存暴涨。锚点整合会构造一个稀疏校正矩阵，当细胞数超过 20 万时，内存占用非常可观。我遇到过一次 30 万细胞的数据，整合阶段直接吃掉 128 GB 内存。后来改成先分面整合：先把类似批次合并成子集，整合完再合并，才把内存降下来。

第三，不要把整合后的表达矩阵直接当作真实表达量去做差异表达分析。integratedassay 里的值是校正后的残差，不是原始 UMIs，用它做FindMarkers会得到一些看似显著但实际不可解释的结果。正确的做法是差异表达永远回到RNAassay，这算是单细胞整合代码里最常见的一个误用。

3. 快速线性整合：Harmony 在整合代码里的定位

3.1 Harmony 为什么快，以及它适用的场景

Harmony 的核心思路是在 PCA 降维空间里做迭代聚类和线性校正。它首先把细胞嵌入到低维 PCA 空间，然后通过一个“软分配”过程把每个细胞分配到多个批次相关的聚类中，再计算每个聚类的批次校正向量，反复迭代直到不同批次的细胞分布趋于一致。因为全程是在低维空间里做矩阵运算，不需要构建细胞两两之间的锚点关系，所以速度比基于锚点的方法快得多。

Harmony 特别适合两个场景：一是作为整合代码的“快速预跑”，先看看批次效应能不能被基本校正；二是批量处理的样本数量很多、且你并不需要精确保留罕见亚群的数据集。它对这个罕见亚群的容忍度没有锚点整合好，但大多数做细胞图谱的项目，Harmony 的表现已经足够。

在 Scanpy 流程里，Harmony 的接入非常顺滑：

import scanpy as sc sc.pp.pca(adata, n_comps=50) sc.external.pp.harmony_integrate(adata, key="batch", basis="X_pca", adjusted_basis="X_pca_harmony") sc.pp.neighbors(adata, use_rep="X_pca_harmony") sc.tl.umap(adata)

关键是把use_rep指向X_pca_harmony，而不是原来的X_pca。很多初学者在这里容易漏掉，结果 UMAP 画出来还是分群严重的模样，以为是 Harmony 没用，实际上只是没有在后续最近邻计算中使用校正后的嵌入了。

3.2 把 Harmony 写进 Scanpy 流程的代码片段

如果只是默认参数跑，Harmony 一般不会出错，但要写出稳定可复现的整合代码，有几个参数值得调。核心是theta和lambda。theta控制每个批次簇的多样性惩罚，值越大，批次之间的校正力度越大；lambda控制岭回归惩罚，值越大，校正越保守。我通常把theta设为 0 到 2 之间，lambda保持默认 1。当批次效应很强、且你知道某些批次里的细胞类型确实缺失时，过高的theta会把缺失的细胞类型也强行拉平，这反而不好。

import harmonypy # 在 scanpy 外部用 harmonypy 直接把 PCA embedding 校正 pca_emb = adata.obsm["X_pca"] harmony_out = harmonypy.run_harmony( pca_emb, adata.obs, "batch", theta=1.5, lambda_reg=1.0, max_iter_harmony=20, ) adata.obsm["X_pca_harmony"] = harmony_out.Z_corr.T

harmonypy 返回的矩阵方向容易让人困惑，转置后才是细胞 × 主成分的形状。我在第一次用的时候在这里卡了半小时，所以特意记下来：Z_corr是 主成分 × 细胞 的矩阵，要.T之后才能放进obsm。

3.3 Harmony 与 Seurat 锚点整合的结果能不能混用

我的建议是不要混用，至少不要在同一份交付结果里说“我用 Harmony 做了批次校正，然后用 Seurat 的锚点整合结果来做聚类”。这两种方法对数据分布假设不同，得到的嵌入空间也不可互换。如果你的下游分析已经明确要使用某个整合矩阵，那就从一而终地基于它做 PCA、聚类和标记基因分析。

不过有一套流程是可以接受的：先用 Harmony 快速评估整合难度，观察整合后的 UMAP 里稀有细胞群是否还在、主要细胞类型是否清晰，再决定是否切换到更重的锚点整合。这样做的好处是避免在数据还没有被理解清楚时，花大量计算资源跑一次FindIntegrationAnchors。我最近一个项目就是先用 Harmony 跑通完整分析管线，确认真实生物学信号都存在，才回头用锚点整合做正式结果。

4. 深度生成模型整合：scVI 的代码实践与调参心得

4.1 scVI 在整合代码中的角色

scVI 是把整合问题转化为一个条件变分自编码器的建模问题。它的基本假设是：每个细胞的观测表达量由一个隐变量（代表真实细胞状态）和一个条件变量（代表批次标签）共同生成。训练时，模型学习去除掉批次条件的影响，把隐变量作为整合后的表示。

和 Harmony、Seurat 锚点这类方法相比，scVI 最大的优势是能够利用全部基因的表达信息，而不是局限于高变基因或 PCA 空间。代价是训练时间长、需要 GPU，而且对数据规模和批次配置更敏感。在 5 万以下的细胞数据上，scVI 的优势不明显；但在几十万细胞、且存在技术差异很大的多模态数据时，scVI 能保留更细腻的亚群结构。

用 scVI 写整合代码并不复杂，核心在训练前后对AnnData的设定：

import scvi import scanpy as sc scvi.model.SCVI.setup_anndata( adata, batch_key="batch", layer=None, categorical_covariate_keys=None, ) model = scvi.model.SCVI(adata, n_latent=30, n_layers=2) model.train(max_epochs=200, early_stopping=True, early_stopping_patience=20) adata.obsm["X_scVI"] = model.get_latent_representation()

注意setup_anndata必须在建模之前调用，它负责把批次标签编码成模型需要的格式。如果你之后修改了adata.obs["batch"]的值，必须重新 setup，否则模型会用到旧的索引，结果非常隐蔽地出错。

4.2 核心模型代码与训练参数

n_latent通常设为 20 到 50 之间。这个值决定隐变量的维度，维度太高容易过拟合噪声，太低则可能丢失细胞类型细节。我一般先用 30，然后观察整合后的聚类效果，再决定要不要增减。

max_epochs不是越大越好。我从实际操作里发现，当每个批量内细胞数超过 10 万时，200 个 epoch 通常是足够的；如果细胞数少，100 个 epoch 反而可能已经开始过拟合。这时候要依赖early_stopping，监控的是训练集上的损失，而不是验证集，所以看到训练损失一直降并不代表模型越来越好。

训练完模型后，还可以把它当做一个生成模型来使用，而不只是拿隐变量。比如model.get_normalized_expression()可以用模型把原始表达量转换成校正后的表达值。但是我不建议用这个值做差异表达分析，原因和 Seuratintegratedassay 类似：生成模型给出的期望表达量是模型推断出来的，不是真实测量值，后续统计检验的解释性会受影响。

4.3 推断批次标签如何影响整合质量

scVI 假设不同批次共享同一个隐空间，只是通过条件变量校正批次偏移。如果某个细胞类型只出现在一个批次中，模型仍然能把这个细胞放回对应的隐空间位置，但它对“这个位置的细胞在所有批次都应该存在”的约束是天然的。换句话说，scVI 不会像 Harmony 那样强制所有批次在所有区域都有细胞，因此在处理批次间细胞组成明显不平衡的数据时，scVI 的表现更稳。

但这里有个隐含的风险：当批次标签和某个生物学变量高度相关时，比如第一批次全是患病组、第二批次全是健康组，模型可能把疾病相关的生物差异也当成批次效应抹掉。为了避免这种情况，我在设置categorical_covariate_keys时会加入一个“疾病状态”变量，让模型知道这个差异不是批次噪声。这是从项目一开始就容易忽略的代码细节，建议在做任何深度整合前，先检查批次标签是否与分组变量存在强相关。

5. 整合效果怎么用代码评估：不要只看 UMAP 图

5.1 三种常用指标的计算代码：混合熵、批次熵、生物学保守性

很多教程把整合效果等同于 UMAP 上不同颜色是否混在一起，这在正式分析里不够。我需要一套可以量化的指标。用得比较多的有三种：批次平均熵、细胞类型平均纯度、以及整合前后基因表达差异的保留程度。

批次平均熵衡量的是在每个局部邻域中，不同批次细胞的混合程度。我用scanpy的sc.metrics相关实现或者自己写一个简单版本。简单思路是：对每个细胞取它在整合嵌入空间中的 k 近邻，计算这 k 个近邻里批次标签的香农熵，然后求所有细胞的平均熵。熵越高说明混合得越好。

from sklearn.neighbors import NearestNeighbors import numpy as np def batch_entropy(embedding, batch_labels, k=30): nbrs = NearestNeighbors(n_neighbors=k, metric="cosine").fit(embedding) _, indices = nbrs.kneighbors(embedding) entropies = [] for idx_list in indices: counts = np.bincount(batch_labels[idx_list], minlength=len(np.unique(batch_labels))) probs = counts / counts.sum() ent = -np.sum(probs * np.log(probs + 1e-12)) entropies.append(ent) return np.mean(entropies)

细胞类型纯度则反过来看，它主要检查已知细胞类型的标签在邻域内的纯度。如果整合后同一细胞类型的细胞在嵌入空间里聚在一起，纯度就高；如果被过度校正打散，纯度就下降。好的整合应该是在两个指标之间取平衡：批次混合度较高，同时细胞类型纯度不能掉得太厉害。

5.2 用代码对比多个整合结果的表格

我会把不同整合方法的结果放在同一套指标下做横向对比。为此，需要一个统一的函数，接收嵌入矩阵、批次标签、细胞类型标签，输出所有指标：

def evaluate_integration(embedding, batch_labels, celltype_labels, k=30): return { "batch_entropy": batch_entropy(embedding, batch_labels, k), "celltype_purity": celltype_purity(embedding, celltype_labels, k), "n_clusters": len(set(celltype_labels)), }

跑完 Harmony、Seurat 锚点、scVI 三套流程后，我会把结果整理成一个小表格。这个表格在汇报给合作者时非常有用：它能直观看出一种方法虽然混合度最高，但细胞类型纯度掉了很多，说明它牺牲了真实生物学信号。而从实际操作来看，scVI 通常在几十万细胞的数据集上能同时保住两个指标，但计算时间也最长。这个决策不是单纯选“最好”的算法，而是选“在可接受时间和资源内能达到最高质量”的算法。

5.3 UMAP 可视化的一些代码习惯

UMAP 不是评估指标，但它是最直观的交付物。我的习惯是同时输出四张图：按批次着色、按细胞类型着色、按某个关键标记基因表达着色、按细胞周期评分着色。前两张用来判断整合混合度，后两张用来判断生物学信号是否保留。如果批次图混得很好，但生物学标记基因的梯度也消失了，那整合就是过校正。

sc.pl.umap(adata, color=["batch", "celltype"], show=False, save="_overview.png") sc.pl.umap(adata, color=["Cd4", "Cd8a", "Foxp3"], show=False, save="_markers.png")

UMAP 本身也有随机性，我在sc.tl.umap之前会固定随机种子，并且在代码里写清楚版本号。比起肉眼判断，我更信任指标表格，但 UMAP 图仍然是给合作者看的第一印象，所以宁可多花几十秒钟把图调好。

6. 实战中踩过的整合报错与性能瓶颈

6.1 内存溢出与稀疏矩阵处理

单细胞数据是典型的稀疏矩阵，但很多整合代码在中间环节会把矩阵转成稠密格式，导致内存暴涨。我遇到过最典型的一次：FindIntegrationAnchors在构造锚点矩阵时默认用了稠密矩阵，10 万细胞 × 2 万基因的规模，直接让一台 128 GB 内存的服务器触发 OOM。

解决思路有两个方向。一是尽可能对特征做裁剪，比如先只保留高变基因再进行锚点搜索，减少矩阵维度。二是显式使用稀疏格式，在 Scanpy 流程里确保adata.X是csr_matrix而不是 numpy 数组，在 Seurat 流程里避免不必要的as.matrix调用。很多函数会隐式转稠密，所以我会在关键步骤前打印type(adata.X)来确认。

另一个内存陷阱在于scvi模型训练时，如果不做setup_anndata的层选择，它默认会用原始计数矩阵，这是没有问题的。但如果你把layer设为一个稠密去噪矩阵，那么数据在拷贝进模型时会被加载为稠密张量，大数据的显存瞬间被吃满。所以我一般保持layer=None，让模型直接从稀疏表达矩阵中读取。

6.2 批次标签不平衡导致的整合崩溃

批次不均衡是另一个高频问题。假设第一个批次有 12 万个细胞，第二个批次只有 5 千个细胞，Harmony 和锚点整合都会倾向把小批次的数据当作噪声处理。从代码层面看，批次标签类别不均衡，会导致模型学习的批次校正向量更多偏向大批次的特征。

对这种问题，我会先做一步直接的数据观察：按批次打印细胞数和主要细胞类型比例。如果发现某个批次的细胞数过少，我会先尝试对这个批次单独做二次抽样或者加权处理，而不是盲目地用整合代码去硬拉。具体做法是在 Harmony 的harmonypy里传入weights，在小批次样本上给更高的权重；在 Seurat 锚点整合里则可以通过samples参数调整不同样本的权重。

还有一点值得注意：批次和高变基因选择之间的耦合。如果你先合并所有批次再统一选高变基因，高变基因列表会被细胞数最多的批次主导，导致低细胞数批次里的特异性基因被过滤掉。这会让本就脆弱的整合雪上加霜。正确做法是每个批次单独选高变基因，再取并集，或者直接用 SCTransform 的anchor.features参数自动处理。

6.3 整合代码的可复现性：随机种子与依赖版本

单细胞整合代码里潜伏着一个很容易被忽略的问题：可复现性。UMAP、Harmony、scVI 都包含随机初始化或迭代过程，不同运行环境下结果可能完全不同。要保证自己第二天跑同一段代码能拿到同样结果，我一般在脚本开头固定三样东西：Python 随机种子、NumPy 随机种子、PyTorch 随机种子。

import random import numpy as np import torch random.seed(42) np.random.seed(42) torch.manual_seed(42)

但这不代表完全可复现。Seurat 的 CCA 整合底层有 C++ 实现，有些步骤的随机性不受 R 层面种子控制，所以即使固定了种子，不同机器上的结果仍可能有小幅波动。我的应对方式是：整个分析流程用版本锁定的环境来跑。我习惯把整合代码放进一个环境文件里，记录scanpy、scvi、harmonypy、Seurat的具体版本。别人拿到这个环境文件，才能复现出基本一致的整合结果。

在可复现性之外，我还建议把整合流程写成一个独立脚本，不要和分析可视化代码混在一起。这样一旦某个环节报错，或者合作者要求换一种整合方法，可以只替换接入的那一层。我之前就是先把整合写成一个可供其他分析复用的模块，现在每次开新项目只需要把批次标签和文件路径传进去，整套流程能很快跑通。

单细胞数据整合代码看上去只是几个函数调用，但每个函数背后的参数、数据假设和应用边界，都能决定结果是否可信。每次开启一份新数据，我最先看的不是算法能跑到什么程度，而是先确认批次效应和生物学信号是否真的有可能被区分。做整合这件事，最终交付的不仅是几张图和几个指标，而是一整套能解释数据来源、参数选择和性能取舍的判断过程。