前阵子帮某实验室的一位同学收拾数据他手上只剩一份cellranger aggr产出的 merged 矩阵、一个aggr.csv各样本原始的 counts 目录早被清理干净。偏偏他要做的是分样本比较而 merged 那份矩阵里几万个 barcode 混在一块儿谁来自哪个样本完全看不出来。我花了大半天把这份合并文件按gem_group拆回原样过程中踩了几个不显眼、但会让你对不上数的坑这里完整记一遍。这篇东西讲的就是多个单细胞测序 aggr 合并文件的分解一份 merged 的filtered_feature_bc_matrix或者h5ad怎么还原成若干个独立的单样本矩阵还原的时候依据是什么、有哪些参数会埋雷、拆完怎么验证是对的。它解决的是数据被合并后失去样本标签这个很具体的痛点——上游跑完 aggr 把几个样本拼成一张大表下游突然要做分样本的差异分析、分样本的质控、或者只想拿其中一个样本画图这时候没有原始的 per-sample 输出就只能从合并文件里往回拆。适合谁看手上正拿着 aggr 输出、但又需要样本级信息的分析人员做单细胞下游流程、需要兼容别人递来的只剩 merged数据的人以及想搞清楚 barcode 后缀到底怎么生成、为什么拆出来数量对不上的同学。文中代码以 Python 为主R 端给等价实现两条路都走通。1. aggr 到底把什么东西合在了一起拆分这件事本质上是把合并时加进去的信息再剥出来。所以在动手之前得先弄明白合并动作到底做了什么、留下了什么痕迹。很多人拆不干净不是代码写错是没搞清楚合并时的信息流向随手按顺序切一刀结果样本边界切偏了。1.1 输出目录里哪些文件带样本信息aggr跑完之后outs/目录大致是这么个结构不同版本略有出入但主干一致outs/ ├── filtered_feature_bc_matrix.h5 ├── filtered_feature_bc_matrix/ │ ├── barcodes.tsv.gz │ ├── features.tsv.gz │ └── matrix.mtx.gz ├── raw_feature_bc_matrix.h5 ├── raw_feature_bc_matrix/ ├── molecule_info.h5 # 部分流程会保留 ├── analysis/ │ ├── clustering/graphclust/clusters.csv │ ├── pca/10_components/projection.csv │ ├── tsne/2_components/projection.csv │ └── differential_expression/*.csv ├── metrics_summary.csv └── cloupe.cloupe这里头真正带样本身份的只有两个地方一是barcodes.tsv.gz里每个 barcode 尾巴上的-N后缀二是一些分析产物比如clusters.csv、projection.csv里同样带着后缀的 barcode 列。features.tsv.gz是全体样本共用的所有样本用的是同一套基因注释这一列天然一致不需要拆。matrix.mtx.gz是基因乘细胞的坐标格式样本信息完全没写在里面只能靠列号去和 barcodes 对齐。换句话说样本标签只存在于 barcodes 的字符串里。你拿到 merged 文件之后如果 barcodes 那一列被某个工具自动去掉了后缀那就等于把唯一的分组依据给删了再想拆就只能靠细胞数比例去猜非常难受。所以第一件要提醒的事是任何读取步骤都要确认后缀还在。metrics_summary.csv也是好东西它按 library 列了每个样本的细胞数、reads 数、中位基因数。拆分完成后拿这张表对账比自己拿原始输出到处找要方便得多。1.2 barcode 后缀是唯一可靠的拆分依据aggr给 barcode 加后缀的规则并不复杂按aggr.csv里的gem_group没写这一列时就按行序自动编号从 1 开始给每个细胞挂一个-gem_group。比如第一个 library 的所有 barcode 都变成AAACCTGAGAAACCAT-1第二个是...-2依此类推。这里的顺序很关键后缀数字对照的是gem_group不是library_id的字母顺序也不是文件名顺序。我在一次复盘里就见过有人按library_id排序去切结果样本顺序和gem_group不一致整份数据的样本标签全体错位还照样跑出了显著差异非常危险。反过来如果你拿不到aggr.csv也可以从后缀反推样本数把 barcodes 全部提出来用正则抓最后一个-后面的数字看有多少个不同的值就是有多少个独立的分组。数量对得上metrics_summary.csv的行数基本就能确认没拆错。注意取后缀一定要用右锚定正则比如-(\d)$。某些平台的自定义 barcode 或者加上样本前缀之后字符串中间也可能出现-用split(-)取最后一段虽然结果一样但用正则更明确、更不容易被后续拼接破坏。1.3 aggr.csv 与 gem_group 的对应关系aggr.csv是个小文件长这样library_id,molecule_h5,gem_group sample_A,/path/A/molecule_info.h5,1 sample_B,/path/B/molecule_info.h5,2 sample_C,/path/C/molecule_info.h5,3gem_group这一列是可选的。如果三个样本各自跑了独立的芯片通道那就三个不同的组一一对应。但如果 sample_A 和 sample_B 是同一张芯片上两个不同的 lane规范写法是两行都写gem_group1这样它们会被当成同一个 GEM well 处理计数在相同 barcode 上合并。这就引出一个非常重要的边界情况同一个gem_group下的多个 library拆不出来。因为它们的 barcode 加了同一个后缀信息在合并阶段就已经被抹平了。这种情况只能回到molecule_info.h5或者各样本的原始 barcode 白名单去区分。所以在动手拆之前第一件事是翻开aggr.csv看一眼gem_group列有没有重复值有重复就提前死心别浪费时间。2. 三条拆分路线怎么选拆分在技术上不难难在选对入口。根据你手上的文件形态大致有三条路从 mtx 目录切、从 h5/h5ad 切、从 R 侧的 Seurat 对象切。三条路的结果应该完全一致但内存占用、依赖环境和出错概率差别不小。2.1 路线 A直接在 mtx 坐标文件上切matrix.mtx.gz是纯文本的坐标格式第一行是注释第二行是行数 列数 非零元素数之后每行是基因行号 细胞列号 计数值。理论上完全可以按列号做流式过滤内存占用极低几个 G 的文件也能扛。但真要在文本层面切需要先算出每个细胞列属于哪个组而这个映射来自 barcodes 的顺序所以你还是得把 barcodes 读进内存这一列很小几十万行也就几 MB。然后做一次列号到组号的映射表再用脚本按行扫描、分流写出。麻烦在于 mtx 要求坐标按列升序排列分流写入之后还得各自排序手写容易出错。我的实践是文件不大非零元素在五亿以内就用scipy.io.mmread一路读进来切简单可靠文件特别大再考虑文本流式方案。别为了省那点内存把流程搞复杂中途出错排查成本更高。2.2 路线 B读 h5 或 h5ad 再切filtered_feature_bc_matrix.h5是 10x 的 HDF5 格式scanpy.read_10x_h5一行就能读成 AnnDatabarcode 落在obs_names上后缀原样保留。这条路最舒服因为切完之后直接write_h5ad就是标准下游输入还能顺便做基因名去重、类型转换。要注意read_10x_h5默认会执行var_names_make_unique遇到重名基因会加-1、-2这种后缀容易和 barcode 后缀混淆。我的习惯是读进来之后立刻把原始gene_ids存到var里后续需要严格匹配再用 id别用被改过的名字。如果手上已经是h5ad比如别人转好的逻辑一样只是 barcode 可能被处理过。先跑一行print(adata.obs_names[:5])看看尾巴还在不在这是所有拆分动作的第一步。2.3 路线 CR 侧 Seurat 读入后 SplitObjectR 用户更习惯Read10X加SplitObject。这条路可行但有两个版本差异要注意Read10X有个strip.suffix参数默认是FALSE也就是保留-1后缀某些旧版本或者某些封装函数会默认剥掉剥掉之后就再也拆不出来了。所以我一般显式写Read10X(data.dir ..., strip.suffix FALSE)把意图钉死。SplitObject需要一个 metadata 列做分组依据所以流程是先建对象、算出gem列、再拆。拆完记得把 barcode 上的后缀去掉并检查是否有跨样本重名否则后面merge的时候会出问题。2.4 三条路线的取舍对比维度路线 Amtx 直切路线 Bh5/h5ad路线 CR/Seurat内存峰值低可流式中高全量载入高对象开销大上手难度中要自己管坐标排序低几行代码低但要熟悉 Seurat输出格式标准 10x 三件套h5adrds适合场景超大矩阵、只要矩阵下游走 Python 生态下游走 R 生态主要风险排序错、列号偏移基因名被改写后缀被自动剥掉实际选择上我基本按下游用什么就用什么入口来定下游是 Python 就选 B是 R 就选 C只有矩阵大到内存扛不住才退到 A。三条路交叉验证一次比如 B 和 C 各跑一遍比细胞数和总 counts是最省心的正确性保障。3. 从零开始的完整拆分实操下面按 Python 主线走一遍完整流程中间穿插 R 的等价写法。所有路径都假设你在outs/的上一级目录执行。3.1 第一步环境与前置校验依赖就三个numpy、pandas、scipy走 h5 再加scanpy和h5py。先做三件校验确认aggr.csv里gem_group有没有重复值有重复就没法按后缀拆。确认barcodes.tsv.gz里的 barcode 是否带后缀。确认metrics_summary.csv里每个 library 的细胞数作为后面拆完的对账基线。import gzip, shutil, re from pathlib import Path import numpy as np import pandas as pd import scipy.io import scipy.sparse as sp OUT Path(outs) MAT OUT / filtered_feature_bc_matrix # 1) 看 aggr.csv 的分组情况 aggr pd.read_csv(aggr.csv) print(aggr) if gem_group in aggr.columns: dup aggr[gem_group].duplicated(keepFalse) if dup.any(): print([警告] 存在共用 gem_group 的 library这些样本无法按后缀拆分) print(aggr[dup]) # 2) 看 barcode 后缀 bcs pd.read_csv(MAT / barcodes.tsv.gz, headerNone)[0].astype(str) gem bcs.str.extract(r-(\d)$)[0] print(后缀分布) print(gem.value_counts(dropnaFalse).sort_index()) print(无后缀细胞数, gem.isna().sum())gem.isna().sum()这个数如果大于零说明有一部分 barcode 没有后缀。这种情况一般出现在混合来源上——比如有人手工拼过矩阵或者用了非标准流程。处理策略后面第 5 节细说。3.2 第二步按 gem_group 切矩阵并回写 10x 目录读取 mtx 时要注意方向matrix.mtx是基因行 × 细胞列所以切样本是切列。按列切片用 CSC 格式效率最高别用 CSR。def read_10x_mtx(d: Path): with gzip.open(d / matrix.mtx.gz, rb) as f: m scipy.io.mmread(f).tocsc() # genes x cells bcs pd.read_csv(d / barcodes.tsv.gz, headerNone)[0].astype(str).tolist() feats pd.read_csv(d / features.tsv.gz, headerNone, sep\t) assert m.shape (feats.shape[0], len(bcs)), 维度对不上先查文件是否配套 return m, bcs, feats维度断言这一步很有必要。我遇到过有人把新版本的features.tsv.gz和旧版本的matrix.mtx.gz混着拷肉眼看不出来切片之后基因全错位做出来的热图怎么看怎么怪。回写函数要保证三件事mtx 的 header 行格式符合 10x 要求、features.tsv.gz三列不丢、barcodes不带换行残留。def write_10x_mtx(m, barcodes, feats, outdir: Path): outdir.mkdir(parentsTrue, exist_okTrue) tmp outdir / matrix.mtx scipy.io.mmwrite(str(tmp), m, fieldinteger) with open(tmp, rb) as fin, gzip.open(outdir / matrix.mtx.gz, wb) as fout: shutil.copyfileobj(fin, fout) tmp.unlink() with gzip.open(outdir / barcodes.tsv.gz, wt, encodingutf-8) as f: f.write(\n.join(barcodes) \n) with gzip.open(outdir / features.tsv.gz, wt, encodingutf-8) as f: for row in feats.itertuples(indexFalse): f.write(\t.join(str(x) for x in row) \n)主流程m, bcs, feats read_10x_mtx(MAT) bc_s pd.Series(bcs, namebc) gem bc_s.str.extract(r-(\d)$)[0] assert gem.notna().all(), 存在无后缀 barcode先处理再拆 total_before int(m.sum()) split_dir Path(split) report [] for g in sorted(gem.unique(), keyint): cols np.where(gem.values g)[0] sub m[:, cols].tocsc() sub_bcs bc_s.iloc[cols].str.replace(r-\d$, , regexTrue).tolist() write_10x_mtx(sub, sub_bcs, feats, split_dir / fgem{g}) report.append({gem_group: g, n_cells: len(sub_bcs), n_counts: int(sub.sum())}) report pd.DataFrame(report) print(report) print(总 counts 校验, total_before, -, report[n_counts].sum())这里有个刻意设计拆出来的 barcode 去掉了后缀。理由是 10x 目录作为标准输入交给下游时带-1后缀容易被某些工具当异常处理但去掉之后跨样本必然出现 barcode 重名所以每个 gem 单独一个目录别急着合并。如果后面要重新合并用add.cell.ids加样本前缀别裸拼。打印出来的total counts校验是必须做的拆分前后总 counts 必须严格相等前提是原始合并没做下采样这一点第 5 节会讲。数值不相等就说明切片有遗漏最常见的是某个gem组没被遍历到。3.3 第三步走 h5 或 h5ad 的更短路径如果矩阵不大用 AnnData 更省事还能顺手把基因名整理干净import scanpy as sc adata sc.read_10x_h5(outs/filtered_feature_bc_matrix.h5) print(adata.shape) print(adata.obs_names[:3].tolist()) # 先确认后缀在 gem pd.Index(adata.obs_names).str.extract(r-(\d)$)[0].values for g in sorted(set(gem), keyint): mask gem g sub adata[mask].copy() sub.obs_names [b[: -(len(str(g)) 1)] for b in sub.obs_names] sub.obs[gem_group] str(g) sub.write_h5ad(fsplit/gem{g}.h5ad) print(g, sub.shape)注意sub adata[mask].copy()里的.copy()少了它后面改名会触发视图警告甚至在旧版本里改不动。3.4 第四步R 端的等价实现R 侧写法有版本差异下面这段在 Seurat v4 和 v5 都能跑library(Seurat) library(Matrix) mat - Read10X(data.dir outs/filtered_feature_bc_matrix, strip.suffix FALSE) obj - CreateSeuratObject(counts mat, project merged, min.cells 0, min.features 0) bcs - colnames(obj) obj$gem_group - ifelse(grepl(-\\d$, bcs), sub(^.*-(\\d)$, \\1, bcs), NA_character_) print(table(obj$gem_group, useNA ifany)) lst - SplitObject(obj, split.by gem_group) dir.create(split, showWarnings FALSE) for (nm in names(lst)) { x - lst[[nm]] colnames(x) - sub(-\\d$, , colnames(x)) saveRDS(x, file file.path(split, paste0(gem, nm, .rds))) cat(nm, ncol(x), sum(x$nCount_RNA), \n) }min.cells 0, min.features 0这两项一定要显式写。默认值是min.cells 0, min.features 0看着没问题但有些人为了顺手过滤一下随手改成 3 或 200拆完的细胞数就和原始对不上还以为是拆分出错其实是建对象时被过滤了。这类顺手造成的误差最难查。4. 拆分结果的校验与对齐拆完不校验等于没拆。这一步看着枯燥但它决定了你后面所有分析结论能不能站住。我给自己定的最低标准是三项对账细胞数、总 counts、特征行数。4.1 细胞数对账的三个层级第一层是组内自洽每个gem组的细胞数之和必须等于 merged 的总细胞数。这一层最容易过因为就是切片除非漏组。第二层是外部分组对账和metrics_summary.csv里每个 library 的细胞数比。这一层是真正的问题探测器。数值对不上通常有三种可能——gem_group存在共用多个 library 拼进同一组该组细胞数偏大、原始输出本身被过滤过、aggr时启用了下采样导致部分细胞被丢弃。第三层是原始文件对账如果各样本的outs/filtered_feature_bc_matrix/barcodes.tsv.gz还在直接数行数。这一层最权威但往往也是最缺的——数据都递到你手上了原始目录大概率已经清了。metrics pd.read_csv(outs/metrics_summary.csv) print(metrics.columns.tolist()) # 常见列名Library ID / Estimated Number of Cells chk metrics[[Library ID, Estimated Number of Cells]].copy() chk.columns [library_id, n_cells_expected] chk[n_cells_expected] chk[n_cells_expected].astype(int) # 按 aggr.csv 把 library_id 映射到 gem_group m aggr[[library_id, gem_group]].merge( report, left_ongem_group, right_ongem_group, howleft) print(m.merge(chk, onlibrary_id))三个数一一对上基本可以放心。对不上但差值很小个位数级别通常是边界细胞在合并阶段被去重或者被 QC 门槛卡掉可以接受差值大到百分之几就要回去查gem_group和归一化参数了。4.2 总 counts 与稀疏结构校验细胞数对上不代表计数对得上。我习惯再加两项import scipy.sparse as sp sub scipy.io.mmread(split/gem1/matrix.mtx.gz).tocsc() print(非零元素, sub.nnz) print(总 counts, int(sub.sum())) print(每细胞中位 counts, float(np.median(np.asarray(sub.sum(axis0)).ravel()))) print(每细胞中位基因数, float(np.median(np.asarray((sub 0).sum(axis0)).ravel())))这两个中位数拿去和metrics_summary.csv里的Median Genes per Cell、Median UMI Counts per Cell比是非常灵敏的探针。如果细胞数对得上但中位数明显偏小八成是合并时做过深度归一化计数被缩放了。特征维度也要查每个拆分目录的features.tsv.gz行数必须完全一致且和 merged 的一致。不一致说明写回的时候丢了行或者被某个工具做了基因过滤。4.3 分析产物也要跟着拆outs/analysis/下的东西同样是 merged 尺度的不是拆完矩阵就完事了。clusters.csv、projection.csv、differential_expression/这些文件的 barcode 同样带后缀可以按同样的规则拆tsne pd.read_csv(outs/analysis/tsne/2_components/projection.csv) tsne[gem] tsne[Barcode].str.extract(r-(\d)$)[0] for g, sub in tsne.groupby(gem, sortTrue): sub.drop(columnsgem).to_csv(fsplit/gem{g}_tsne.csv, indexFalse)但要说清楚这些是 merged 尺度下的降维和聚类结果拆到单样本只对看看每个样本在合并空间里落在哪有用不能当单样本自己的 UMAP 用。想拿单样本的降维图得在拆分后的矩阵上重跑。同理differential_expression/是拿library_id做分组算出来的组间差异拆开后没有意义直接忽略。5. 踩过的坑与排查速查表下面这几条都是我在实际项目里真踩过的按坑的隐蔽程度排序越靠前越容易被忽略。5.1 后缀被悄悄剥掉最典型的一次某同学拿Read10X读完、建好对象、存在中间 rds 里准备第二天接着做第二天用了另一个封装函数重新读 rds结果那个函数内部调了类似make.names的处理把-1变成了...1。再想拆的时候正则全失效gem全变成NA。这类问题的排查方法很简单在任何中间产物落盘前后各打印一次head(barcodes)。发现后缀变形优先怀疑读取函数的默认参数而不是自己的正则。Read10X记得写strip.suffix FALSEPython 侧read_10x_h5不动obs_names但要小心别自己跑obs_names_make_unique把重复 barcode 改成新的后缀。5.2 归一化参数埋的雷aggr有个--normalize参数可选none默认、mapped、hm。走默认不重采样拆出来的 counts 和原始一致一旦选了mapped每个 library 的 reads 会被下采样到同一水平合并矩阵里的 counts 是缩过的整数这时候拆出来的不是原始 counts做绝对定量、做饱和度评估全都不准。怎么判断看metrics_summary.csv里各 library 的Mean Reads per Cell是否被人为拉平或者直接看拆分后每个样本的总 counts 是否和原始文件如果还在差一截。差得整齐就是归一化过。这种情况不是不能拆而是要在结论里写清counts 经过深度归一化仅用于相对比较。我一般会在项目 README 里加一行备注免得半年后自己都忘了。5.3 共用 gem_group拆不开的死结回到第 1.3 节那个问题。如果aggr.csv里两个 library 共用同一个gem_group它们的 barcode 在合并时已经混成一锅按后缀只能拆出组拆不出样本。这时候唯一可行的路子是如果还留着各样本的molecule_info.h5用分子级别的信息重建分样本矩阵如果还留着各样本的barcodes.tsv.gz白名单哪怕只有 barcode 列表可以用它做交集反查如果两者都没有那就接受这组样本无法区分的事实别硬拆。顺便提一个反例有人为了拆得干净手工改aggr.csv给每个 library 都写独立gem_group然后重跑合并这是对的但如果只是把 merged 的 barcode 后缀按比例重新编号那就是纯粹的数据造假做出来的样本间差异全是人为构造的非常危险。5.4 常见问题速查表现象最可能原因排查动作处理后缀解析全为 NA读取时后缀被剥离/变形打印 barcodes 头几行换读取方式禁用剥后缀参数拆分细胞数之和小于总数遍历漏组或正则不匹配检查gem.value_counts补齐分组用右锚定正则某组细胞数明显偏大多个 library 共用 gem_group查aggr.csv重复值回到分子文件重建或接受counts 总量对不上启用了 depth normalization比对 metrics 的中位 reads标注为归一化数据仅做相对比较基因数不一致features 版本混杂比对行数与基因 id统一参考版本后重做R 中 barcode 变成...1中间函数做了名字规范化head(colnames(obj))关掉规范化或从源重读拆完后再 merge 报重复去后缀后跨样本重名检查交集大小merge 时加样本前缀h5ad 读入后 obs 少一列写入时被覆盖对比shape与源文件重新拆分并保留分组列5.5 一个容易被忽略的细节空组与单细胞组有些gem_group在 merged 里可能只有寥寥几个细胞甚至是零。零细胞的情况意味着那个 library 在合并阶段就被全部过滤了。拆分脚本遍历sorted(gem.unique())时不会遍历到空组所以split/下不会生成对应目录这是正常的但要在报告里显式标出来否则业务方会以为你漏了一个样本。细胞数极少的组比如只有十几个也要警惕这类样本单独做 QC 和聚类基本没有统计意义建议在报告里标注细胞数不足建议合并分析或不单独解读。6. 什么时候不该拆以及拆完怎么收场写到这里拆分的技术路径已经完整了。但我更想说的是判断层面的东西有些场景下硬拆反而是错的。这一节没有什么代码全是经验。6.1 三种更适合重做而不是回拆的情况第一种原始molecule_info.h5还在。那最优解不是拆合并文件而是直接用分子级别信息重建各样本矩阵或者干脆按各自的molecule_info.h5重跑一遍定量。拆合并文件本质上是有损恢复能不做就不做。第二种合并时做过深度归一化。拆出来的 counts 是缩过的任何基于绝对值的指标都不能信。这种情况下我更倾向于找到原始数据重做实在找不到就在所有下游结论前加一句基于归一化 counts。第三种共用gem_group且没有其他线索。前面说过这是死结。与其编一套看似合理的样本标签不如老实告诉需求方这几个样本合在一起了分不开。还有一类特殊情况你只是想看某个样本的 QC 指标并不需要矩阵。这时候完全不用拆metrics_summary.csv里已经有每个 library 的细胞数、reads、中位基因数直接拿来看就行。别为了一个数字去切一个几 G 的矩阵。6.2 拆完之后重新整合的正确姿势拆完再合是很多流程的必经之路比如先分样本 QC、再统一整合。这里有个高频翻车点去掉后缀之后不同样本的 barcode 会重名。10x 的 barcode 是 16bp 随机序列几万个细胞里跨样本撞车几乎是必然的。Python 侧用scanpy.concat时带上label和keys它会自动加前缀import anndata as ad parts [sc.read_h5ad(fsplit/gem{i}.h5ad) for i in (1, 2, 3)] merged ad.concat(parts, labelsample, keys[fgem{i} for i in (1, 2, 3)], index_unique-, joinouter) print(merged.obs[sample].value_counts())index_unique-是关键参数它会把样本标签拼到 barcode 后面避免重名。R 侧对应的是merge(x, y, add.cell.ids c(gem1,gem2))一样的作用。整合完成之后建议再跑一次总 counts 对账确认拼接过程没有丢细胞。这一步很多人省掉结果在整合后的对象上发现细胞数比拆分前少了几千个回头一查是join参数写成了inner把基因并集之外的东西全丢了。我个人在实际操作中的体会是拆分这件事的技术难度远低于它的判断难度。代码就是几十行写一次能复用很久真正花时间的是想清楚这份合并文件里到底还残留多少样本信息、哪些结论在拆分后依然成立。每次动手前先把aggr.csv和metrics_summary.csv翻一遍把分组关系和细胞数基线记在纸上后面所有的对账和排查都会顺很多。