之前试过用第三方工具做RDS和H5AD的格式互转，结果踩了不少坑：各种莫名其妙的bug层出不穷，调试起来费时费力，最后还没搞定。这次我们换个思路——直接从底层数据结构入手，用几行简洁的代码自己实现，彻底摆脱对外部工具的依赖。
R需要的包：
library(Seurat)
library(Matrix)
python需要的模块
import scanpy as sc
import anndata as ad
import scipy.sparse as sp
from scipy.sparse import coo_matrix, csr_matrix, csc_matrix
import pandas as pd
import sys
from pathlib import Path
任务1：实现rds向h5ad的转换
情景1：提取rds的重要信息
首先读取rds文件，查看相关信息
library(Seurat)
library(Matrix)
seurat_obj <- readRDS("pbmc.rds")
seurat_obj
An object of class Seurat
13714 features across 2638 samples within 1 assay
Active assay: RNA (13714 features, 2000 variable features)
3 layers present: counts, data, scale.data
2 dimensional reductions calculated: pca, umap
接下来提取一些信息（counts,datas,metadata等），写个函数优化提取过程（避免某些信息不存在报错）
#安全提取数据
safe_get <- function(expr, name) {
tryCatch(expr, error = function(e) {
cat(sprintf("Warning: Could not extract %s (%s)\\n", name, e$message))
NULL
})
}
提取counts数据
counts <- safe_get(
GetAssayData(seurat_obj, assay = "RNA", layer = "counts"),
"counts"
)
提取data数据
normalized <- safe_get(
GetAssayData(seurat_obj, assay = "RNA", layer = "data"),
"normalized data"
)
提取元数据
obs <- seurat_obj@meta.data
# 确保行名是细胞ID
if (is.null(rownames(obs))) {
rownames(obs) <- colnames(seurat_obj)
}
提取基因信息
var <- data.frame(
gene = rownames(seurat_obj),
row.names = rownames(seurat_obj),
stringsAsFactors = FALSE
)
counts 元数据 基因信息是必须提取的，接下来保存一下
# 保存矩阵（转置为cell x gene以兼容Python），outdir是输出路径哈
writeMM(t(counts), file.path(outdir, "counts.mtx"))
write.csv(obs, file.path(outdir, "obs.csv"), row.names = TRUE)
write.csv(var, file.path(outdir, "var.csv"), row.names = TRUE)
额外的，如果存在已经处理ok的降维数据，比如pca,umap,tsne，也可以提取出来。下面编写一个函数处理一下：
# 提取并保存降维结果
extract_dimred <- function(seurat_obj, reduction, filename) {
embed <- safe_get(
Embeddings(seurat_obj, reduction = reduction),
reduction
)
if (!is.null(embed)) {
write.csv(embed, file.path(outdir, filename), row.names = TRUE)
cat(sprintf("Saved %s (%d dimensions)\\n", reduction, ncol(embed)))
return(TRUE)
}
FALSE
}
extract_dimred(seurat_obj, "pca", "pca.csv")
extract_dimred(seurat_obj, "umap", "umap.csv")
extract_dimred(seurat_obj, "tsne", "tsne.csv")
在提取的过程中，可以发现tsne是不存在的（也就是rds中没有做tsne非线性降维分析哈），不过不影响程序的执行。
rds中关键数据提取完成之后，接下来就是转换成h5ad的步骤操作了。
情景2：将rds中提取的数据构建h5ad文件
##导入模块
import scanpy as sc
import anndata as ad
import scipy as sp
from scipy.sparse import coo_matrix,csr_matrix,csc_matrix
import pandas as pd
from pathlib import Path
上面的基于rds提取的重要数据存放到一个特定文件夹indir目录，判断核心文件存在否。
required = ['counts.mtx', 'obs.csv', 'var.csv']
indir = Path(indir)
missing = [f for f in required if not (indir / f).exists()]
if missing:
raise FileNotFoundError(f"Missing required files: {missing}")
读取表达矩阵
counts = sp.mmread(indir / "counts.mtx").T.tocsr()
读取元数据
obs = pd.read_csv(indir / "obs.csv", index_col=0)
var = pd.read_csv(indir / "var.csv", index_col=0)
创建AnnData对象
# 验证维度匹配
if counts.shape[1] != len(obs):
raise ValueError(f"Counts rows ({counts.shape[0]}) != obs rows ({len(obs)})")
if counts.shape[0] != len(var):
raise ValueError(f"Counts cols ({counts.shape[1]}) != var rows ({len(var)})")
# 创建AnnData对象
adata = ad.AnnData(X=counts, obs=obs, var=var)
adata
AnnData object with n_obs × n_vars = 2638 × 13714
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.5', 'seurat_clusters'
var: 'gene'
添加降维结果
dimred_files = {
"X_pca": "pca.csv",
"X_umap": "umap.csv",
"X_tsne": "tsne.csv"}
for obsm_key, filename in dimred_files.items():
filepath = indir / filename
if filepath.exists():
print(f"Adding {obsm_key}…")
df = pd.read_csv(filepath, index_col=0)
# 确保索引匹配
if df.index.equals(adata.obs_names):
adata.obsm[obsm_key] = df.values
else:
# 尝试重新索引
df = df.reindex(adata.obs_names)
if df.isna().any().any():
print(f"Warning: {filename} has mismatched indices, skipped")
else:
adata.obsm[obsm_key] = df.values
adata
AnnData object with n_obs × n_vars = 2638 × 13714
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.5', 'seurat_clusters'
var: 'gene'
obsm: 'X_pca', 'X_umap'
结果保存
adata.write("adata.h5ad")
任务2：实现h5ad向rds的转换
情景1：提取h5ad的重要信息
import scanpy as sc
import scipy as sp
from scipy.sparse import coo_matrix, csr_matrix, csc_matrix
import pandas as pd
import numpy as np
import os
读取h5ad文件
adata = sc.read_h5ad("adata.h5ad")
adata
AnnData object with n_obs × n_vars = 2638 × 13714
obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent.mt', 'RNA_snn_res.0.5', 'seurat_clusters'
var: 'gene'
obsm: 'X_pca', 'X_umap'
layers: 'log_normalized'
提取表达矩阵
if 'counts' in adata.layers:
counts = adata.layers['counts']
else:
counts = adata.X
##转置
if sp.sparse.issparse(counts):
counts = counts.T.tocsr() # 转置：R是gene x cell，Python是cell x gene
else:
counts = csr_matrix(counts).T
提取标准化数据
if 'log_normalized' in adata.layers:
normalized = adata.layers['log_normalized'].T.tocsr()
elif 'normalized' in adata.layers:
normalized = adata.layers['normalized'].T.tocsr()
else:
normalized = None
提取元数据和基因信息
# 提取obs和var（Seurat的meta.data和gene信息）
obs = adata.obs
var = adata.var
提取降维结果
pca = pd.DataFrame(adata.obsm['X_pca'], index=adata.obs_names) if 'X_pca' in adata.obsm else None
umap = pd.DataFrame(adata.obsm['X_umap'], index=adata.obs_names) if 'X_umap' in adata.obsm else None
tsne = pd.DataFrame(adata.obsm['X_tsne'], index=adata.obs_names) if 'X_tsne' in adata.obsm else None
结果保存（outdir是输出路径）
sp.io.mmwrite(f"{outdir}/counts.mtx", counts)
if normalized is not None:
sp.io.mmwrite(f"{outdir}/normalized_counts.mtx", normalized)
obs.to_csv(f"{outdir}/obs.csv")
var.to_csv(f"{outdir}/var.csv")
if pca is not None:
pca.to_csv(f"{outdir}/pca.csv")
if umap is not None:
umap.to_csv(f"{outdir}/umap.csv")
if tsne is not None:
tsne.to_csv(f"{outdir}/tsne.csv")
情景2：将h5ad中提取的数据构建rds文件
library(Seurat)
library(Matrix)
首先读取上述情景1中h5ad输出的重要数据，indir是重要数据存放目录
## 读取表达矩阵
counts <- readMM(paste0(indir, "/counts.mtx"))
counts <- as(counts, "CsparseMatrix")
## 读取元数据
obs <- read.csv(paste0(indir, "/obs.csv"), row.names = 1, stringsAsFactors = FALSE)
var <- read.csv(paste0(indir, "/var.csv"), row.names = 1, stringsAsFactors = FALSE)
## 设置行列名（关键！Matrix Market格式不保存名称）
rownames(counts) <- var$gene
colnames(counts) <- rownames(obs)
## 创建Seurat对象
seurat_obj <- CreateSeuratObject(counts = counts, meta.data = obs)
seurat_obj
An object of class Seurat
13714 features across 2638 samples within 1 assay
Active assay: RNA (13714 features, 0 variable features)
2 layers present: counts, data
还可以添加标准化数据呦
norm_file <- paste0(indir, "/normalized_counts.mtx")
if(file.exists(norm_file)){
normalized <- readMM(norm_file)
normalized <- as(normalized, "CsparseMatrix")
rownames(normalized) <- var$gene
colnames(normalized) <- rownames(obs)
# 存入data slot
seurat_obj[["RNA"]]$data <- normalized
}
添加降维数据
pca_file <- paste0(indir, "/pca.csv")
if(file.exists(pca_file)){
pca <- read.csv(pca_file, row.names = 1)
seurat_obj[["pca"]] <- CreateDimReducObject(
embeddings = as.matrix(pca),
key = "PC_",
assay = "RNA"
)
}
umap_file <- paste0(indir, "/umap.csv")
if(file.exists(umap_file)){
umap <- read.csv(umap_file, row.names = 1)
seurat_obj[["umap"]] <- CreateDimReducObject(
embeddings = as.matrix(umap),
key = "UMAP_",
assay = "RNA"
)
}
tsne_file <- paste0(indir, "/tsne.csv")
if(file.exists(tsne_file)){
tsne <- read.csv(tsne_file, row.names = 1)
seurat_obj[["tsne"]] <- CreateDimReducObject(
embeddings = as.matrix(tsne),
key = "tSNE_",
assay = "RNA"
)
}
查看结果
seurat_obj
An object of class Seurat
13714 features across 2638 samples within 1 assay
Active assay: RNA (13714 features, 0 variable features)
2 layers present: counts, data
2 dimensional reductions calculated: pca, umap
结果保存
saveRDS(seurat_obj, "object.rds",compress = T)
早期软件互转：h5ad与rds互转
