单细胞h5ad转成seurat对象后metadata信息丢失怎么办？(GSE156625)

生信技能树

发布于 2025-03-20 18:47:04

4900

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

我们的生信马拉松授课群里有一个学员学习练手的单细胞数据为h5ad格式，他想转成seurat对象后发现metadata信息丢失了，下面来看一看~

数据情况

数据编号为：GSE156625，

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE156625，下载下面的几个文件：

GSE156625_HCCbarcodes.tsv.gz 128.7 Mb (ftp)(http) TSV
GSE156625_HCCgenes.tsv.gz 258.6 Kb (ftp)(http) TSV
GSE156625_HCCmatrix.mtx.gz 561.9 Mb (ftp)(http) MTX
GSE156625_HCCscanpyobj.h5ad.gz 800.2 Mb (ftp)(http) H5AD

然后整理成下面的格式：genes.tsv.gz 改名为 features.tsv.gz

├── GSE156625_HCC
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
└── GSE156625_HCCscanpyobj.h5ad

上一次我写了一个h5ad转为seurat对象的稿子：一步一个坑：单细胞数据的h5ad格式转换成R可读取对象

底下有个留言说可以使用zellkonverter，然后我去试了一下，发现这个软件在转换的时候也需要用到conda环境，并且默认安装，网速慢，然后还不能指定已有的conda环境，果断放弃了！

https://www.bioconductor.org/packages/release/bioc/vignettes/zellkonverter/inst/doc/zellkonverter.html

转换实战

我们还是使用 一步一个坑：单细胞数据的h5ad格式转换成R可读取对象 这个帖子中的办法吧！

终端/linux：首先是创建一个 conda的小环境 sc

# 创建小环境sc
conda create -n sc -y python=3.10
conda activate sc

# 安装软件
conda install -y loompy
conda install anndata -y
conda install scipy -y

然后是在R语言中：

rm(list=ls())
library(sceasy)
library(reticulate)
library(Seurat)
packageVersion("Seurat")
# [1] ‘5.1.0’

library(BiocParallel)
register(MulticoreParam(workers = 4, progressbar = TRUE)) 
# 这个conda环境就是刚刚创建的，换成自己的目录即可
use_condaenv(condaenv = '/nas2/zhangj/biosoft/miniconda3/envs/sc')
# h5ad转为Seurat
sceasy::convertFormat(obj = "GSE156625/GSE156625_HCCscanpyobj.h5ad", from="anndata",to="seurat",outFile = 'GSE156625_HCCscanpyobj.rds')

# 读取刚刚转换成功的rds对象
sce.all <- readRDS("GSE156625_HCCscanpyobj.rds")
sce.all
dim(sce.all)
rownames(sce.all)

# 查看特征
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
dim(sce.all@meta.data)
phe <- sce.all@meta.data
save(phe, file = "phe.RData")

这里读进来，发现数据只有2000多个基因了，所以我们就把metadata信息保存出去，然后读取上面的三个 10x 标准文件吧。

下面还是R语言：

## 读取表达矩阵
counts <- Read10X("GSE156625/GSE156625_HCC/", gene.column = 2)
dim(counts)
as.data.frame(counts[1:5,1:5])
# 创建seurat对象
sce <- CreateSeuratObject(counts, meta.data = phe, min.cells = 3)
sce

# > sce
# An object of class Seurat 
# 26055 features across 42762240 samples within 1 assay 
# Active assay: RNA (26055 features, 0 variable features)
# 1 layer present: counts

# 两者取交集
head(sce@meta.data)
a <- intersect(colnames(counts),rownames(phe))
length(a)
head(a)
sce_sub <- subset(sce, cells = a)
sce_sub
table(sce_sub$orig.ident)
table(sce_sub$patient_id)
table(sce_sub$patientno)

需要注意的是，上面读取的三个文件的表达矩阵有42762240，4千多万个细胞，而GSE156625_HCCscanpyobj.h5ad中是73589多个，需要取个交集~

后面就可以走标准流程分析了~