rds与h5ad互转全攻略：Seurat到Scanpy的数据迁移指南

rds和h5ad互转这个话题我上一篇写过基础流程今天这篇算是姊妹篇但我要把视角拉得更深一点。在单细胞转录组这个领域混久了几乎每个项目都会遇到这种局面合作方用R的Seurat跑完了整套QC和聚类交付给你一个.rds文件而你自己的主分析管线在Python里后面想做velocity、做跨平台整合或者维护一个大型Atlas都离不开h5ad。两边数据格式不通最直接的办法就是做格式互转。这个需求对谁最有用两类人。一类是偏R的拿到Python生态出来的h5ad不知道如何塞回Seurat继续跑另一类是偏Python的被塞了一个.rds却要融入Scanpy流程。不管你是哪一类核心诉求其实都一样——把对象完整、无损地搬到对方那边尤其是metadata、降维坐标、邻接图这些“分析结论”不能丢。这篇我重点讲清楚底层的结构对应关系、手动转换的完整代码、几个容易翻车的隐蔽问题以及怎么设计一个可靠的转换验证环节。1. 转换思路先搞清楚两边数据结构在“装”什么1.1 rds到底装了什么rds是R的序列化格式本质上是把一个对象用saveRDS打包。一个常见的Seurat对象并不是一个简单的二维表而是一个S4类型的复合结构。我拆解它的时候重点会看这么几个区域assaysRNA、蛋白、ATAC等不同的测定单元每个assay内部又有多个表达矩阵典型的包括counts、data通常是log1p后的归一化数据、scale.data用于热图的标准化矩阵。meta.data细胞注释一行对应一个barcode。聚类编号、样本分组、批次信息都在这里。reductions降维结果比如PCA、UMAP、tSNE。每个reduction内部有cell.embeddings细胞在主成分或低维空间的坐标和feature.loadingsPCA的基因载荷UMAP一般没有这个字段。graphs细胞间的邻接图比如RNA_snn理论上可以理解成一个细胞×细胞的稀疏矩阵。commandSeurat运行过的命令记录有版本和时间信息属于辅助内容。misc其他杂项可能塞了不少自定义数据。我拆这些槽位的时候最常犯的错是只盯着表达量矩阵忽略了reductions和graphs。尤其是别人交付的rds聚类注释基本都是画在UMAP上的结果你只转了表达矩阵换个环境重新算一遍UMAP坐标对不上注释全部错位。互转的核心不是“把表达矩阵搬过去”而是“把分析结论完整搬过去”。1.2 h5ad的结构与命名习惯h5ad是anndata的HDF5存储格式它的结构设计得比较规整有明确的硬盘读取接口。它通常包含X细胞×基因的主矩阵可以是稠密、稀疏或分块存储。obs细胞注释DataFrame行索引与X的行一一对应。var基因注释DataFrame列索引与X的列一一对应。obsm / varm细胞和基因的低维嵌入比如obsm里有X_pca、X_umap。obsp / varp细胞对细胞、基因对基因的pairwise结构connectivities邻接图通常放在obsp里。layers多层表达矩阵可以同时存放counts、log1p、scale.data或者去批次后的校正矩阵。uns非结构化metadata聚类颜色、marker基因列表、差异分析结果都可以丢这里。两边的对应关系在逻辑上其实相当直接Seurat的meta.data对应adata.obsassay内部的counts和data对应adata.layers里的counts和lognorm或者干脆用data作为Xreductions对应obsm。只有当你用这种视角去拆一个Seurat对象时代码写起来才顺手——不是想着“用一个函数一键转”而是思考“每一块信息究竟要搬去哪个目标结构”。1.3 索引从1到0方向从基因×细胞到细胞×基因这里额外强调两个最容易被忽略的基础细节。一个是索引规则。R是1-based索引Python是0-based索引所以转换之后任何基于位置的矩阵切片边界都要重新核对。最直观的表现是某些在R里看着正常的代码搬到Python里发现最后一行或最后一列永远对不上而且这类问题极难定位它在向量长度不一致时才会爆出来如果恰好遇到长度相同的情况错误值会悄悄混进下游分析。另一个更关键的是维度方向。Seurat的assay矩阵默认是基因×细胞而anndata要求X是细胞×基因。简单说表达量矩阵必须做转置转置后原矩阵的列名barcode进入adata.obs_names原矩阵的行名基因进入adata.var_names。我第一次用SeuratDisk时以为它内部已经处理好了结果拿到h5ad后一检查有个assay的行列方向是反的。手动组装最大的好处就在这里每一块数据的方向都是由你自己控制的。2. 工具选型解析为什么我“嫌弃”SeuratDisk却还要提它大概所有人第一反应都是用SeuratDisk。它提供Convert和SaveH5Seurat系列函数理论上可以把Seurat对象保存为.h5seurat再转成.h5ad。我早期用这条路径做过不少数据转换但项目做多了慢慢总结出几条绕不开的痛点。第一个痛点是版本兼容混乱。SeuratDisk在社区维护上明显没有跟上Seurat自身演进的节奏Seurat版本迭代一轮之后默认assay结构已经变化了好几次老版本的SeuratDisk写出来的h5seurat文件在新环境下打不开或者新版本Seurat写出的文件在老SeuratDisk上读不进去。遇到这种情况第一反应是升级或降级某一个包但升级之后其他依赖又要适配问题就像滚雪球一样膨胀。对于长期项目来说这种不确定性是很耗精力的。第二个痛点也是我最介意的是信息丢失。它有固定的默认映射规则并不会自动把你没显式告诉它的字段放进去。实际测试中我遇到过一次转换后的h5ad里logs矩阵不在预期位置某个降维坐标直接缺失还有一个图的邻接信息放进了uns而不是obsp。这种“看上去转了但数据不完整”的结果最坑人——后续分析不是立即报错而是在某个隐性环节才出问题排查起来比一开始就手动组装要痛苦得多。第三个痛点是控制力不足。不同来源的Seurat对象里meta.data字段风格千差万别有的是factor型有的列名带空格有的列名是中文。SeuratDisk并不会自动清洗这些信息也不支持你中途介入处理。所以就算走SeuratDisk路线清洗工作终归要自己做那还不如直接走手动管线至少每一步都知道发生了什么。但我也承认它在简单场景下的价值。比如快速试读一下数据、粗略确认基因名和细胞数量用SeuratDisk转一次是很快的。只是如果你追求可靠性和完整性我强烈建议直接手动组装。下面这套方案就是我目前在实际项目中稳定使用的转换流程。3. 实操流程从Seurat对象到h5ad的完整转换方案3.1 R侧准备把Seurat拆成“标准零件”我这里以最常见的单细胞RNA Seurat对象为例代码逻辑可以直接推广到其他assay。第一步确认RNA assay存在导出三个关键矩阵library(Seurat) library(Matrix) seu - readRDS(seurat_object.rds) # 表达量矩阵 counts - GetAssayData(seu, assay RNA, slot counts) lognorm - GetAssayData(seu, assay RNA, slot data) # 检查维度R侧通常是基因 × 细胞 dim(counts)第二步处理meta.data。有两个问题在这里必须提前处理一是factor列二是列名规范。factor列如果直接导出内部存储是整数加levels属性序列化过程如果没把levels信息带出去对方拿到的就是一堆整数非常容易误读。列名里如果带空格、括号、减号scanpy里有些操作会直接因为语法问题报错。meta - seumeta.data # 把factor列统一转成字符串避免整数编码 meta[sapply(meta, is.factor)] - lapply(meta[sapply(meta, is.factor)], as.character) # 列名规范化 names(meta) - make.names(names(meta))第三步导出降维结果和邻接图pca_emb - seureductions$pcacell.embeddings umap_emb - seureductions$umapcell.embeddings pca_loadings - seureductions$pcafeature.loadings snn - seugraphs$RNA_snn第四步做矩阵转置并保存中间结果。这一步很关键Seurat矩阵是基因×细胞anndata需要的是细胞×基因所以转置是必须的。counts_t - t(counts) lognorm_t - t(lognorm) # 确保是稀疏矩阵 counts_t - as(counts_t, CsparseMatrix) lognorm_t - as(lognorm_t, CsparseMatrix) save(counts_t, lognorm_t, meta, pca_emb, umap_emb, pca_loadings, snn, file seurat_exported.rda)这里值得多说一句为什么强调CsparseMatrix。R中dgCMatrix就是CSC格式转置之后再显式转一次CsparseMatrix是为了确保下游Python端拿到的是稀疏列压缩格式。anndata底层用CSC和CSR都支持但scanpy读入后很多算法对CSC做切片和矩阵乘法更顺手而且h5ad在写入时对CSC的支持也比较成熟。当然如果你在R侧本来就习惯了dgCMatrixt()之后本身就是CSC不需要额外转换但千万不要用as.matrix()把它变密。几万细胞乘几万基因的矩阵一旦变稠密内存瞬间爆掉。3.2 Python侧组装构造AnnData并写入h5ad接下来是Python端的接收脚本import anndata as ad import scipy.sparse as sp import numpy as np import pandas as pd # 读取R侧保存的中间结果这里用rpy2或者转换成csv/h5格式都可以 # 这里我用load_rda做示意实际项目中可以直接用h5文件交换 adata ad.AnnData( Xsp.csc_matrix(lognorm_t), obsmeta_df, varpd.DataFrame(indexgene_names), layers{counts: sp.csc_matrix(counts_t)} ) # 索引统一转成字符串 adata.obs_names [str(x) for x in adata.obs_names] adata.var_names [str(x) for x in adata.var_names] adata.var_names_make_unique() # 降维结果写入obsm/varm adata.obsm[X_pca] pca_emb adata.obsm[X_umap] umap_emb if pca_loadings is not None: adata.varm[loadings] pca_loadings # 邻接图写入obsp adata.obsp[connectivities] snn # 可选的uns信息把聚类注释也放进去 adata.uns[orig_seurat_clusters] adata.obs[seurat_clusters].tolist() # 写文件 adata.write_h5ad(converted.h5ad)组装时还有几个细节要单独讲一下。第一obs和var的索引类型。R的barcode和基因名到了Python里通常是字符串但偶尔会出现全数字barcode。scanpy内部很多函数是按字符串索引去匹配obs/var的如果全数字barcode被读成int后面取子集会非常难受。所以我习惯在组装前统一用str()包一下。第二var_names_make_unique不是万能药。它对于少量重复基因名确实有效但如果原始数据里有一堆空的、NA的基因名它的拼接逻辑会产生大量带“-1”“-2”后缀的名字后面做marker基因筛选时会让你头疼不已。所以最佳策略是在源头清洗基因名不要让空值进入组装环节。第三scale.data要不要放进h5ad。我的建议是如果后续要复用Seurat的scale.data做差异分析再考虑放layers[scale]否则不建议放因为它是dense矩阵体积大而scanpy的大多数流程不会主动调用它。如果你确实要放记得转float32不然内存开销会上涨得很快。3.3 大对象的backed模式别一次性load进内存再来专门说说大规模数据。如果R端导出的矩阵是5万细胞×3万基因这种级别直接全量读进内存再转Python中途很容易崩掉。我的处理习惯是两步走R侧不要用saveRDS一次性打包分块把稀疏矩阵写入一个独立的h5文件每个assay对应一组数据集细胞注释单独写一个CSV。这样Python端读取时可以先读细胞注释再按需抽取矩阵块。Python侧用anndata.read_h5ad(..., backedr)先以只读模式打开文件不急着把全部矩阵load进内存等真正需要计算的时候再对某个层做to_memory()或只抽取需要的子集。backed模式有两个隐藏的坑必须说明。第一backed模式是只读的不能直接在读取的backed对象上修改obs、obsm等属性。如果你想往里加一列注释先得把它读进内存或者用write_h5ad生成新文件后再继续操作。第二某些scanpy函数是不支持backed对象的比如运行harmony整合时你发现它正在悄悄把整个矩阵load进来内存又爆了。所以backed模式适合“轻量探索”不适合“直接跑全套分析”这个边界要心里有数。4. 反向转换从h5ad回rds的实操要点4.1 读取h5ad的两种路线从h5ad回rds我也试过两条路线。路线一是用R的anndata包底层靠reticulate调Python环境里的anndata库直接读h5ad然后拆给Seurat。这个思路最直接前提是你的R环境里有可用的Python环境而且anndata的Python版本与R包的reticulate兼容。我用下来感觉它适合临时性读取不适合做正式交付——因为你的R脚本多了一个Python运行依赖换一台机器可能就起不来了。路线二是用zellkonverter先把h5ad读成SingleCellExperiment再用Seurat::as.Seurat(sce)转成Seurat对象。zellkonverter的底层依赖管理更干净basilisk会自动搭建好独立Python环境不用你手动去指定路径。而且SingleCellExperiment作为中间格式可以在后期做很多质控检查转换过程也更透明。我推荐路线二下面代码以zellkonverter为例。4.2 从h5ad到SingleCellExperiment再到Seuratlibrary(zellkonverter) library(SingleCellExperiment) library(Seurat) sce - readH5AD(converted.h5ad) sce seu - as.Seurat(sce)但这条路线有几个隐含细节需要处理。第一个是X和layers的映射。readH5AD默认把adata.X读成assay(sce, X)layers会变成其他assays。如果你希望Seurat对象的counts是原始counts而data是log1p建议在as.Seurat之前手动指定seu - as.Seurat(sce, counts counts, data X)前提是sce的assays里确实存在counts。如果源h5ad里只有X没有layers那你得根据X的类型决定是否能用expm1等函数逆推counts。这里我要特别提醒不同来源的h5ad对X的定义不一致有的X是counts有的是log1p还有的甚至做了z-score标准化。转回R之前必须先确认X到底是什么。否则配出来的Seurat对象看着没问题跑FindMarkers时结果却完全不符合预期因为数据口径根本不是同一种。第二个是obsm到reductions的映射。as.Seurat会把reducedDims映射成Seurat的reductions但名字可能是X_pca、X_umap而不是pca、umap。如果后续代码里按名字取reduction建议转完后统一重命名。另外某些Scanpy里常见的嵌入如X_diffmap放进Seurat里没有对应的原生类型Seurat会保留为未识别的reduction不影响对象完整性但你可能需要在分析前决定是否保留。第三个是obs中factor和字符串的处理。h5ad里的obs列是字符串或数值到了SingleCellExperiment变成DataFrame再转Seurat列名中带特殊字符或全角符号的问题依然存在建议转完后用make.names统一清洗。另外注意一个高频现象很多h5ad在导出时把缺失值写成了字符串“NA”在R里会被当作普通字符串而不是真正的NA。如果你后面要用is.na()做过滤这种字符串会漏掉需要先判断是否要替换回NA。4.3 逆变换和稀疏矩阵转置反向转换同样要处理维度方向。h5ad的X是细胞×基因R的assay要求是基因×细胞所以要用Matrix::t()做转置而不是base::t()。base::t()会把稀疏矩阵转成稠密矩阵几万细胞的数据直接爆内存这个坑我踩过不止一次。counts_gene_cell - Matrix::t(counts_cell_gene)如果是counts层也要做一样的转置。转完之后如果数据变成了dgTMatrix或其他格式建议统一用as(counts_gene_cell, dgCMatrix)转回dgCMatrix因为Seurat的很多函数默认对dgCMatrix做处理。这一步在代码里只多一行但对后续的稳定性影响非常大。5. 常见问题与排查实录5.1 基因名变成日期格式这个坑有极强的迷惑性。某次我在R侧把基因名写成CSV再读回来发现“1-Mar”“2-May”这类基因名全部被读成了日期格式。其实readRDS不会出现这问题但一旦你在某个环节里走了CSV或者某个同事习惯性把metadata写成文本文件再导入就会碰到。解决方法是read.csv时加check.namesFALSE读入后强制row.names as.character(row.names)。这种问题排查起来比较隐蔽因为只有部分基因名受影响而且在R里可能不会立即报错等到下游比对基因名时才露馅。5.2 factor导出后变成整数编码这个在前面已经提到。再展开一个场景你从某个公共数据集的Seurat对象里取meta.data发现上面有seurat_clusters、cell_type、sample_id等列其中很多是factor。如果直接构造AnnDataobs里的值会变成0、1、2这样的整数等你用scanpy画图时会发现图例是数字而不是标签再做分组时也会出错。建议在R侧加一行保险代码meta[sapply(meta, is.factor)] - lapply(meta[sapply(meta, is.factor)], as.character)Python侧也可以事后补救for col in adata.obs.columns: if str(adata.obs[col].dtype) category: adata.obs[col] adata.obs[col].astype(str)5.3 h5ad读取后var_names对不上有一回我拿到别人转好的h5ad用scanpy读取后var_names末尾多了一堆“-1”“-2”明显是之前做var_names_make_unique时把原本唯一的基因名也改了。排查发现原始数据里有少量空基因名和重复的“NA”字符串。这类问题在做marker基因筛选时会让结果非常难以对应。我的建议是如果源数据来自公共数据库基因注释一般没问题如果是自己拼的矩阵一定在源头清洗基因名不要让空值进入组装环节。如果已经进入h5ad了再清洗的成本就高很多。5.4 降维坐标丢失或放错位置这个归因于SeuratDisk的默认映射问题。实际案例中某同事用SeuratDisk把一个带UMAP坐标和聚类注释的Seurat对象转成h5ad结果读取后scanpy画图坐标疯狂报错。最终定位到原因是UMAP嵌入没有被正确写入obsm而是被塞进了uns。手动组装方案不会出现这个问题因为每一块都是你显式放进去的。转换后检查几个关键字段是否存在应该是流程的一部分而不是等到画图时报错了才发现。5.5 邻接图方向搞反Seurat里的graph例如RNA_snn是细胞×细胞anndata里的obsp也是细胞×细胞方向一致按理说不需要转置。但如果某些自定义图结构是“基因×细胞”或者从其他软件导出的图矩阵方向就可能颠倒。一旦方向错了下游用connectivities做聚类或可视化邻接关系时会出现非常离谱的结果。建议转换后在两边各抽一个小子集手动核对某几个细胞之间的关系值确认以后再继续下游工作。这一步虽然繁琐但在我的经验里值得做因为方向搞反带来的问题往往不是报错而是静默的坏结果。5.6 backed模式导致的只读误解如果使用了backedr读取h5ad要特别注意它只能读、不能原地改。很多新手会直接对backed对象赋值obs列运行时报错“object is read-only”然后一头雾水。正确做法是先copy到内存或者用write_h5ad生成新文件后再继续操作。这一点对大规模数据尤其重要因为很多人当初就是为了省内存才用backed模式结果赋值obs时又不知不觉把整个对象load进来内存直接炸了。记住一个模式backed模式适合“只读探索”不适合“修改属性”。6. 验证环节转换之后如何确保数据没丢转换本身只是开始真正决定交付质量的是验证。我每次交付转换脚本都会在脚本里加入一个“验证区块”检查三件事第一表达量矩阵是否一致。在R端和Python端各抽同一个细胞子集用同样的基因子集做sum和mean对比。数值不要求完全一样但浮点误差应该控制在极小范围内。如果log1p矩阵来自Seurat而Python端用的是scanpy自己的标准化函数重算的结果那数值可能有系统性偏差这不是转换的锅而是算法口径不一致。第二降维坐标与注释是否对齐。抽几个已知的cluster查看UMAP坐标和聚类标签是否在同一批细胞上。如果一个cluster中某个细胞的坐标在R端和Python端差异很大往往是样本顺序错位而不是数值精度问题。这种错位最容易发生在obs或var的排序被某些中间步骤打乱的情况下。第三邻接图结构是否完整。检查obsp里connectivities矩阵的维度、非零元素数量以及稀疏度是否和R端保持一致。非零元素数量差异过大说明边缘或者细胞过滤信息丢了。我每次都会在转换脚本最后打印一份简单的“转换报告”记录细胞数、基因数、矩阵稀疏度、降维坐标维度、邻接图非零数。这份报告随h5ad一起交付协作者拿到数据后只要不是直接报错都可以用它做初步核验。虽然看起来只是几行print但在实际协作中能省下不少沟通时间。再补充一点反向转换回Seurat时验证方向是一样的只是对调。我会在R端重新检查adata.obs里某几个关键字段是否完整UMAP坐标是否是源h5ad里的值。整个过程不需要写特别复杂的逻辑难的是养成“转完必查”的习惯。数据格式互转这件事数值精度可以接受微小误差但结构完整性必须百分之百对齐。7. 结尾里的一点实操体会做的次数多了我越来越觉得rds和h5ad互转本质上不是格式翻译而是理解两个生态对“分析结果”这个概念的封装方式不同。Seurat侧重把对象做成一个自包含的容器所有中间产物都堆在里面anndata则更偏向把不同层次的数据物理分开放置靠命名规范把它们关联起来。理解了这层差异手动组装反而比黑盒工具更可靠。最后再分享一个小技巧如果转换是团队内部的高频需求不要每次都在临时脚本里复制粘贴代码可以把R侧导出和Python侧组装各自封装成两三个函数接口只留输入rds路径和输出h5ad路径。这样内部处理逻辑统一以后遇到新的metdata字段类型或者新的assay结构只需要在一个地方做扩展。我在实际项目中就是这么做下来的一年多没再因为格式转换问题被人凌晨喊醒过。数据搬运这件事虽然不够炫酷但稳定性做好了能让上下游合作顺畅很多。