ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

单细胞数据格式转换:rds/mtx转h5ad实操与避坑指南

单细胞数据格式转换:rds/mtx转h5ad实操与避坑指南 1. 转换前先想清楚rds、mtx、h5ad 分别是给谁用的做单细胞分析的人估计都遇到过这种场景合作方发来一个.rds文件或者从 GEO 下载到一坨.mtx文件而你手里主力分析工具是 Python 的 scanpy或者你需要把数据交给只会用 Python 的同事去跑流程。这时候你面临的第一道坎就是数据格式转不转得动、转得准不准。先说清楚这三样东西到底什么来头。.rds是 R 语言里用saveRDS()存出来的单个 R 对象文件通常里面是一个Seurat对象。Seurat 对象是一个极度复杂的嵌套结构除了表达矩阵还背着 meta.data细胞注释信息、降维结果tSNE、UMAP、PCA、marker 基因列表、以及各种自定义的 assay 信息。如果你只是想把表达矩阵弄出来用Read10X()或者Seurat::GetAssayData()就能拆出来但你要是想连注释信息一起转过去就得摸清楚 Seurat 对象内部组织方式。.mtx是 Matrix Market 格式本身是一种通用的稀疏矩阵存储格式单细胞领域最常见的形态是 10X Genomics 的 Cell Ranger 输出目录里面包含matrix.mtx、barcodes.tsv、features.tsv或genes.tsv三个文件。其中matrix.mtx只存稀疏矩阵里的非零元素三个文件合在一起才能真正还原一个完整的、带行名基因和列名细胞的表达式矩阵。.h5ad是 Python 生态里 scanpy 和 anndata 的核心格式基于 HDF5。这个格式的好处非常直观它不仅能存表达矩阵还能把obs细胞信息、var基因信息、obsm、varm、uns等各种元数据全部塞进一个文件里并且支持随机读取子集不用一次性全部加载进内存。转录组数据动辄几十万细胞h5ad 在 IO 效率和存储灵活性方面目前没有对手。一句话概括**rds 是 R 世界的全量快递箱mtx 是只装矩阵的真空袋h5ad 是 Python 世界的全量快递箱。**转换这件事本质上就是把数据从一个完整容器里安全地搬运到另一个完整容器里。顺便说一句搜索“rds 转换”的时候常会看到“阿里云RDS”之类的词那是关系型数据库服务跟单细胞分析里的.rds文件完全是两码事别混淆了。下面我讲的都是生信里面的 R 数据对象文件。2. 转换前先自查三个关键问题能帮你省半天时间2.1 你的数据源结构到底是什么接到一个.rds文件第一步千万别急着写转换代码先在 R 里看看这个对象到底长什么样。我见过不少人把test.rds读进来直接用as.matrix()去取数据结果报错半小时才发现里面不是 Seurat 对象而是自己当年随手存的一个 list。打开 R 或者 RStudio跑这几行library(Seurat) obj - readRDS(your_data.rds) class(obj)输出可能是Seurat也可能是SeuratObject、list、dgCMatrix。如果是 Seurat 对象继续查# 查看 Seurat 对象的内部结构 str(obj, max.level 2) # 查看默认 assay 名称和内容 DefaultAssay(obj) Assays(obj) # 查看 meta.data 列名 colnames(objmeta.data) # 查看降维结果名称 reductions(obj) # 查看不同 assay 里的表达矩阵 GetAssayData(obj, assay RNA, slot counts)[1:5, 1:5]这里有个特别容易翻车的点Seurat 对象的表达数据分为两个 slotcounts和data。counts是原始整数计数data是归一化后的数据。如果你下游要用 scanpy 的pp.normalize_total()和log1p()那推荐导出counts如果你只想做可视化或者直接用已经归一化好的数据那就导出data。也别天真地以为GetAssayData()默认拿到的就是counts——它默认拿的是dataslot。2.2 mtx 目录下文件是否齐全且对得上号10X 输出的 mtx 目录通常长这样filtered_feature_bc_matrix/ ├── barcodes.tsv.gz ├── features.tsv.gz └── matrix.mtx.gz旧版是genes.tsv代替features.tsv。读取之前先确认三个文件的行数是否匹配matrix.mtx注释里有一行定义了维度比如10000 5000 400000表示 10000 个基因、5000 个细胞、40 万个非零值。features.tsv有多少列新版至少有三列ensembl ID、symbol、type第三列通常标注是Gene Expression、Antibody Capture还是CRISPR Guide Capture。如果你只想保留基因表达数据读取后要做一次过滤。barcodes.tsv里的细胞条码可能带后缀-1也可能不带这本身不是问题但转换到 h5ad 后要保持一致不然之后 merge 不同样本会很痛苦。2.3 下游目标流程决定你怎么转这一点我觉得是最容易被忽视的。h5ad 只是一个数据容器你在里面装什么完全取决于下游要怎么用。如果你下游要做 cell type annotationobs里必须带上细胞类型注释、样本分组、批次信息。如果你要从uns里读 marker 基因列表那就要把 Seurat 里的markers结果以结构化方式导进去。如果你只是拿这个 h5ad 做 sanity check那只要保证表达矩阵无误即可别的不重要。先想清楚这个问题能帮你省很多后面修 bug 的时间。3. 实操rds 转 h5ad 的三种路子3.1 最稳妥的方式R 中转 AnnData再用 scanpy 读在高版本 R4.0里可以直接用anndata这个 R 包library(Seurat) library(anndata) obj - readRDS(your_data.rds) # 导出原始计数或归一化数据 counts - GetAssayData(obj, assay RNA, slot counts) meta - objmeta.data # 用 anndata 包创建 h5ad adata - AnnData( X Matrix::t(counts), # 注意scanpy 是 cell x geneR 是 gene x cell obs meta, var data.frame(gene_symbol rownames(counts), row.names rownames(counts)) ) ad$write_h5ad(your_data.h5ad)这里有一个极其重要的维度问题。Seurat 的对象是基因在行、细胞在列而 scanpy 的 AnnData 是细胞在obs、基因在var所以表达矩阵必须转置即Matrix::t(counts)。写反会导致后续所有分析全乱套而且因为 h5ad 本身是稀疏矩阵存储写反了不太会报错只会得到错误结果。这个坑我见过至少 5 个朋友踩过。一个更大的坑是anndataR 包的安装。CRAN 上虽然有但要求 Python 环境里有 anndata并且要通过reticulate连接。建议在 R 里先确认library(reticulate) reticulate::conda_list()如果你本机根本没有配 Python 环境推荐先创建一个独立的 conda 环境把scanpy装上然后让 R 的 reticulate 指向这个环境再装 R 版 anndata。具体操作conda create -n scvi -y python3.10 conda activate scvi pip install scanpy anndata然后回到 Rlibrary(reticulate) use_condaenv(scvi, required TRUE) library(anndata)这么一套走通之后转换速度会快很多因为底层调用的是 Python 的 anndata 写入逻辑性能没问题。但问题是很多人的 Seurat 对象里除了主 assay 之外还有SCT、integrated等 assay如果只转RNA的 counts那些信息就丢了。如果你需要保留多个 assay建议逐个GetAssayData()后放进layers。3.2 更通用的方式先导出 mtx再用 scanpy 读如果你不想折腾 R 和 Python 的环境互通或者目标对象不是标准 Seurat 对象可以先从 R 里把矩阵和元数据导出成常规文件再到 Python 里拼装。这个方法虽然多了一步但是每一步都比较可控。R 侧操作obj - readRDS(your_data.rds) counts - GetAssayData(obj, assay RNA, slot counts) # 写成 mtx 三件套 Matrix::writeMM(counts, file matrix.mtx) # 细胞信息 write.table(colnames(counts), file barcodes.tsv, row.names FALSE, col.names FALSE, quote FALSE) # 基因信息 gene_df - data.frame( gene_id rownames(counts), gene_symbol rownames(counts), type Gene Expression ) write.table(gene_df, file features.tsv, row.names FALSE, col.names FALSE, quote FALSE, sep \t) # 元数据单独存 write.csv(objmeta.data, file metadata.csv, row.names TRUE)Python 侧操作import scanpy as sc import pandas as pd adata sc.read_mtx(matrix.mtx).T # 又见转置 adata.var pd.read_csv(features.tsv, sep\t, headerNone, names[gene_id, gene_symbol, type]) adata.var_names adata.var[gene_symbol].values adata.obs_names pd.read_csv(barcodes.tsv, headerNone)[0].values meta pd.read_csv(metadata.csv, index_col0) adata.obs meta adata.write_h5ad(your_data.h5ad)这一步里最值得注意的细节是sc.read_mtx()读进来的维度是基因 x 细胞所以必须.T转置成细胞 x 基因。另外features.tsv里的基因名如果存在重复scanpy 在设置var_names时会自动加上后缀不会报错但后续做 marker 分析时需要留意最好提前去重。3.3 直接从 mtx 三件套转 h5ad 的完整流程如果你手里的数据本身就是 10X 标准输出压根不需要经过 R直接用 scanpy 就能搞定import scanpy as sc adata sc.read_10x_mtx(filtered_feature_bc_matrix/, var_namesgene_symbols, make_uniqueTrue) adata.var_names_make_unique() adata.write_h5ad(your_data.h5ad)这里var_namesgene_symbols表示使用features.tsv里的第二列symbol作为基因名如果你想要 ensembl ID就改用var_namesgene_ids。读进来之后我会先做一个快速 QCprint(adata) print(adata.var.head()) print(adata.obs.head())检查一下矩阵维度、非零元素数量、细胞数、基因数是否符合预期。如果barcodes.tsv里的细胞条码带-1后缀obs_names会是类似AAACCTGGTAT-1的样子。一些下游工具对此敏感建议统一去掉后缀adata.obs_names [x.split(-)[0] if x.endswith(-1) else x for x in adata.obs_names]但注意如果有多个样本合并后导出去掉后缀很大概率造成条码撞车。所以这个操作要放在合并之后、导出之前做或者干脆保留后缀。4. mtx 转 h5ad 时最容易翻车的四个细节点4.1 稀疏矩阵的存储类型10X 的 mtx 文件默认是dgCMatrix格式读进来的scanpy 读取后默认是scipy.sparse.csr_matrix。这两种格式本质上都能存稀疏矩阵但在做矩阵切片、行归一化这些操作时性能差异很大。如果你要频繁按基因取行用csr_matrix更合适如果按细胞取列csc_matrix更快。写入 h5ad 时anndata 会自动处理稀疏矩阵的编码不需要我们手动转换但如果你把adata.X设成了numpy.ndarray稠密矩阵导出文件体积会瞬间膨胀十倍甚至几十倍。40 万细胞的矩阵稠密存储分分钟几十 GB稀疏存储可能只有几百 MB。所以在构造 AnnData 时务必确保X保持稀疏格式。有个判断小技巧import scipy.sparse as sp if not sp.issparse(adata.X): print(警告X 不是稀疏矩阵)假如读取时不小心用了pd.DataFrame之类的结构可以用scipy.sparse.csr_matrix(adata.X)强行转回来。4.2 基因名的唯一化与大小写10X 的 features 文件里经常出现重复基因名比如MIR600HG出现两次或者 ensembl ID 和 symbol 之间的对应关系有错位。scanpy 的make_unique()会把重复名改成MIR600HG-1、MIR600HG-2这种形式看起来很丑但总比直接报错强。另一个容易忽略的问题是大小写。人的基因 symbol 通常按官方名首字母大写、其余小写如TP53但 Mouse 的基因 symbol 通常是全大写如Trp53是错的正确是Trp53也就是首字母大写后面小写实际上 mouse 基因 symbol 官方是首字母大写、其余小写如Trp53。如果你从某个老旧注释文件里拿到的基因名是全大写或全小写下游做富集分析时会有大量匹配不上。建议在转换后统一做一次 gene symbol 标准化# 以人基因为例把符号统一成官方格式太复杂这里只做简单去空格 adata.var[gene_symbol] adata.var[gene_symbol].str.strip().str.upper() # 谨慎使用这里提个醒**转换成全大写只适用于特定物种比如斑马鱼不要无脑用。**正确的做法是弄一个物种对应的基因注释文件做一次 symbol 映射。4.3 obs/var 索引的类型一致性AnnData 对索引类型极其敏感。如果你obs_names里面混入了整数和字符串或者var_names里有空值、NA后面做子集操作时会出一堆诡异错误比如“KeyError: None”或者维度对不上。写入 h5ad 之前务必检查两件事# 检查是否有空值 print(adata.obs_names.isnull().sum()) print(adata.var_names.isnull().sum()) # 检查是否有重复 print(adata.obs_names.duplicated().sum()) print(adata.var_names.duplicated().sum()) # 统一转为字符串 adata.obs_names adata.obs_names.astype(str) adata.var_names adata.var_names.astype(str)这个操作我在转一个 GEO 下载的 mtx 数据时踩过坑当时features.tsv最后多了一个空行导致读进来 33693 个基因名里有 1 个是NaN后续 scanpy 全流程一跑就炸。后来写进转换前检查清单再没出过问题。4.4 元数据列的 dtypes 问题从 R 导出的metadata.csv列类型经常是乱的。比如orig.ident是字符串nCount_RNA是整数percent.mt是小数但read_csv()默认推断类型偶尔会出错尤其是那些全列为0的列会被推断成整数导致后续运行sc.pp.calculate_qc_metrics时发生类型冲突。建议读取后手动指定meta pd.read_csv(metadata.csv, index_col0) # 数值型列统一转 float num_cols [nCount_RNA, nFeature_RNA, percent.mt] for col in num_cols: if col in meta.columns: meta[col] pd.to_numeric(meta[col], errorscoerce)另外 R 里 factor 型列导出后到 Python 就是普通字符串没问题。但 R 里rownames如果是一些数字开头的条码CSV 读进来后index_col0可能会被自动转成整数一定要dtypestr兜住meta pd.read_csv(metadata.csv, index_col0, dtypestr)5. 大矩阵转换的内存管理与加速技巧单细胞数据越来越卷动不动就是百万级细胞。我最近转一个 78 万细胞的样本原始 mtx 文件 3.6 GB内存 64 GB 的机器差点跑死。分享几个实操经验。5.1 用chunked方式写入和读取如果你用 R 的 anndata 包或者 Python 的 anndata 直接写超大文件默认一次性加载全部稀疏矩阵并编码非常吃内存。anndata 支持分块写入adata.write_h5ad(large_data.h5ad)这是常规方式适合 30 万细胞以内的数据。超过这个量级可以考虑用adata.write_h5ad(large_data.h5ad, compressiongzip)加compressiongzip能让文件体积小很多但写入速度会慢一些读取时也会增加解压开销。对于需要长期归档的数据压缩是值得的对要频繁读写的中转文件不压缩体验更好。5.2 不要用 pandas 做超大矩阵的中间载体很多人习惯把矩阵读成pd.DataFrame再拼装我强烈不建议这么做。DataFrame在列名、行名复杂时的内存开销远高于纯稀疏矩阵。正确姿势是matrix sc.read_mtx(matrix.mtx).X # 保持稀疏格式或者from scipy.io import mmread matrix mmread(matrix.mtx).T.tocsr()记得调用.tocsr()转成稀疏 CSR 格式后续操作效率高很多。5.3 利用 scanpy 的 backed 模式做超大数据集scanpy 提供backed模式不必一次性把全部数据加载进内存adata sc.read_h5ad(large_data.h5ad, backedr) print(adata[:10000, :].X.shape) # 按需切片读取写入大 h5ad 时也可以adata.write_h5ad(large_data.h5ad, chunksTrue)chunksTrue会按照数组分块写入减少峰值内存。实测下来70 万细胞、3 万基因的数据默认方式写入大约需要 20 GB 内存chunksTrue能降到 10 GB 左右。5.4 能瘦身的先瘦身再转换很多 rds 文件里存了多个冗余 assay比如原始数据里既有RNA又有integrated还有一个没删干净的SCT。如果下游只需要RNA的 counts建议先在 R 里删掉多余的 assay 再导出DefaultAssay(obj) - RNA obj[[integrated]] - NULL obj[[SCT]] - NULL另外 Seurat 对象里的reductionstSNE、UMAP、PCA如果不打算保留也一并清掉能进一步缩小内存占用。只保留表达矩阵和 meta 数据转换效率会有质的提升。6. 转换后的验证清单别急着删原始文件转换完成不代表万事大吉。我习惯写一个验证脚本花两分钟检查五个项目确认无误后才考虑删掉原始文件。import scanpy as sc import numpy as np adata sc.read_h5ad(transformed_data.h5ad) # 1. 维度检查细胞数、基因数是否和源数据一致 print(fShape: {adata.shape}) # 2. 表达矩阵非零元素数量是否合理 print(fNon-zero: {adata.X.nnz if hasattr(adata.X, nnz) else (adata.X 0).sum()}) # 3. obs 和 var 索引是否有空值 print(fNaN in obs_names: {adata.obs_names.isnull().sum()}) print(fNaN in var_names: {adata.var_names.isnull().sum()}) # 4. 抽查几个已知基因的表达值是否和 R 侧对得上 gene TP53 if gene in adata.var_names: cell_idx 0 print(f{gene} expression at first cell: {adata.X[cell_idx, adata.var_names.get_loc(gene)]}) # 5. 保存后是否能正常读取 reloaded sc.read_h5ad(transformed_data.h5ad) print(reloaded)挑一个在 R 侧已经查过表达量的基因对比它在 h5ad 里对应细胞的值能快速发现行列是否写反、矩阵是否转置错误这类致命问题。比如 R 里某个基因在第一个细胞的 counts 是 3那在 h5ad 里也应该同一个细胞读到 3如果差得离谱大概率是转置搞错了。另外一个容易被忽略的验证点是稀疏矩阵的对角线。当你打印adata.X[:5, :5]时如果发现对角线几乎全是非零值暂时不用慌因为稀疏矩阵的存储是按行索引的打印的结果看起来像对角线其实是显示效果。关键还是做数值比对。7. 避坑经验与工作流建议做格式转换这两年我个人遇到最多、最隐蔽的问题一是转置二是元数据丢失三是基因名不统一。转置问题前面已经强调过了这里再补充一个检查技巧如果你从 R 导出的matrix.mtx第一行是%%MatrixMarket matrix coordinate integer general用sc.read_mtx()读进来直接.T后var_names对应的是 R 里的行名obs_names对应 R 里的列名这个对应关系可以帮你核对是否写反。元数据丢失是最防不胜防的。比如 Seurat 对象里很多注释是存在objmeta.data里的但如果你用的导出路径是GetAssayData(obj, slot data)而不是完整 Seurat 对象那些注释就丢了。所以我在团队内部定了一条规矩转换流程必须写清楚“保留哪些元数据”并且转换后抽查两三个样本的注释信息是否完整。基因名不统一的问题建议不要在转换阶段解决而是在拿到最终数据后用一个标准注释文件统一映射。转换阶段只需保证var_names和原始符号一致不要顺手做大小写转换或者添加前缀后缀否则查原始文件时会疯掉。关于工作流设计我现在比较推荐的方式是所有原始数据先统一转成 h5ad后续无论是 R 还是 Python 分析都从 h5ad 出发。R 侧读 h5ad 可以用zellkonverter包转成SingleCellExperiment对象再用Seurat::as.Seurat()转回 Seurat。这样团队协作时R 和 Python 之间不再被格式卡脖子谁的活儿谁去读同一个中间文件数据版本也更整齐。最后分享一个我踩过的真实教训有一个 30 万细胞的样本从 rds 转 h5ad 后发现obs_names全部变成了1-1、2-1这种编号排查了半天才发现是 R 侧导出colnames(counts)时矩阵被隐式转换成了普通 matrix自带的行名全部丢了。从那之后我每次导出前都会强制检查stopifnot(!is.null(colnames(counts)), !is.null(rownames(counts)))任何时候都不要轻信原始文件格式的稳定性转换前检查、转换后验证这两步做到位单细胞数据跨语言流转其实是很顺手的事。
返回列表