Python中Seurat与h5ad数据双向转换及10X多样本数据整理读取
当下已完成基于R语言的单细胞分析实战相关系列，涵盖初级、中级和高级三个阶段，每个阶段目前各有4篇内容，后续会陆续更新。接下来，打算开展基于Python的单细胞分析实战系列。不过在该系列内容正式开启之前，笔者觉得有必要先学习数据转化与读取的流程。
本次涉及的工程文件可通过网盘获取：中级篇2，链接: https://pan.baidu.com/s/1y-HHLXoXsJbgWKCdz26-gQ，提取码: yx93 。另外，可向“生信技能树”公众号发送关键词‘单细胞’，直接获取Seurat V5版本的完整代码。
首先推荐python官网：https://www.python.org/doc/
其教程有多种语言版本，能选择中文。
快速检索就能查到所需信息，比如查询os库。
此网站可当作工具手册助力使用者快速上手，当然大模型辅助也是不可或缺的。
Seurat和h5ad数据相互转化
后续肯定会用到数据的相互转化，所以得掌握这个步骤。
1.导入(Seurat V5对象转换为h5ad数据)
R语言部分
rm(list = ls())
library(Seurat)
library(qs)
library(reticulate)
library(hdf5r)
library(sceasy)
library(BiocParallel)
register(MulticoreParam(workers = 4, progressbar = TRUE))
scRNA_V5 <- qread("CD4+T_final.qs")
scRNA_V5
# 一个Seurat类的对象
# 21308个特征分布在5912个样本中，属于1个分析测定
# 活跃的测定: RNA (21308个特征，2000个可变特征)
# 存在3个层: counts, data, scale.data
# 计算了3个维度约简: pca, harmony, umap
2.配置python环境
在linux/终端中，numpy不能用最新版本。
conda create -n sceasy python=3.9
conda activate sceasy
conda install anndata
conda install scipy
conda install loompy
conda install numpy==1.24.4
3.R语言转换
再回到R语言，在R语言中调用Python环境。
# 在R语言中加载python环境
use_condaenv('sceasy')
loompy <- reticulate::import('loompy')
py_config()
# python: /opt/anaconda3/envs/sceasy/bin/python
# libpython: /opt/anaconda3/envs/sceasy/lib/libpython3.9.dylib
# pythonhome: /opt/anaconda3/envs/sceasy:/opt/anaconda3/envs/sceasy
# version: 3.9.20 | packaged by conda-forge | (main, Sep 22 2024, 14:06:09) [Clang 17.0.6 ]
# numpy: /opt/anaconda3/envs/sceasy/lib/python3.9/site-packages/numpy
# numpy_version: 1.24.4
# loompy: /opt/anaconda3/envs/sceasy/lib/python3.9/site-packages/loompy
开始转换
# Seurat 转 AnnData
scRNA_V5[["RNA"]] <- as(scRNA_V5[["RNA"]], "Assay")
sceasy::convertFormat(scRNA_V5, from="seurat", to="anndata",
outFile='scRNA_V5.h5ad')
# AnnData对象，n_obs × n_vars = 5042 × 20124
# obs: 'nCount_RNA', 'nFeature_RNA', 'Sample', 'Cell.Barcode', 'Type', 'RNA_snn_res.0.1', 'RNA_snn_res.0.2', 'RNA_snn_res.0.3', 'RNA_snn_res.0.4', 'RNA_snn_res.0.5', 'RNA_snn_res.0.6', 'RNA_snn_res.0.7', 'RNA_snn_res.0.8', 'RNA_snn_res.0.9', 'RNA_snn_res.1', 'RNA_snn_res.1.1', 'RNA_snn_res.1.2', 'seurat_clusters', 'celltype', 'seurat_annotation'
# var: 'vf_vst_counts_mean', 'vf_vst_counts_variance', 'vf_vst_counts_variance.expected', 'vf_vst_counts_variance.standardized', 'vf_vst_counts_variable', 'vf_vst_counts_rank', 'var.features', 'var.features.rank'
# obsm: 'X_pca', 'X_harmony', 'X_umap'
# 警告信息:
# 在.regularise_df(obj@meta.data, drop_single_values = drop_single_values) :
# 丢弃单个类别的变量:orig.ident
4.在jupyter lab中确认一下
打开jupyter lab
首先得安装anaconda/miniconda等；
- 直接点击电脑上的anaconda软件（笔者安装的是anaconda，所以这么写），进入界面后点击jupyter lab
- 在终端中输入jupyter lab，电脑会自动打开该界面，如果直接以家为自己的工作目录，就直接输入jupyter lab，如果想要切换工作目录，那就先cd进去，再输入jupyter lab
# mkdir ~/xxxx
cd ~/xxxx
jupyter lab
注册python内核
激活环境之后，注册小环境对应的内核，在Jupyter Notebook中运行Python代码时，真正执行代码的“引擎”是ipykernel，Jupyter提供前端界面（网页Notebook），ipykernel是后端：它运行代码并把结果返回给前端显示。所以先要安装一下ipykernel。
conda activate sceasy
pip install ipykernel
python -m ipykernel install --user --name sceasy --display-name "Python(sceasy)"
增加了一个新的sceasy内核，之后涉及数据转化就可以连接该内核。
如何删除
# 查看当前列表
jupyter kernelspec list
# 删除
jupyter kernelspec uninstall sceasy
# 删除conda环境
conda remove --name sceasy --all
为何创建环境后还需注册内核？
在回答这个问题前，需要了解Conda、Anaconda、Miniconda
-
Conda: 一个开源的软件包管理和环境管理系统，最初为支持Python语言设计，但现在也能管理和部署任何语言的软件包与依赖项。它能帮用户轻松安装不同软件版本，还能在相互隔离的环境中管理这些软件包，避免版本冲突。Conda可用命令行使用，支持多个平台（Windows、macOS和Linux）。
-
Anaconda: 基于Conda的发行版，专为数据科学和机器学习等任务设计。它包含Conda以及大量预装软件包、库和工具，比如Jupyter Notebook、Spyder IDE等。Anaconda发行版很适合刚开始从事数据分析、科学计算、机器学习等领域的新手，因为它几乎包含了所有常用的库和工具，减少了单独安装配置的时间。不过，由于其庞大的体积（通常超过数GB），对只需要少量特定工具或库的用户来说，可能过于臃肿。
-
Miniconda: 是轻量级的Conda安装器，只包含Conda包管理器和Python基础环境。相比完整的Anaconda发行版，Miniconda更精简，允许用户根据自身需求自定义安装所需的软件包和库。这种方式对已有特定要求或希望从头构建自己环境的高级用户很合适。Miniconda提供灵活基础，用户可在此基础上通过Conda安装任何需要的包，不必携带不必要的预装软件。所以Anaconda和miniconda都属于conda，且conda更像管理平台（工具），管理着软件和环境。
而Conda环境是独立的Python环境，包含特定版本的Python和安装的库，旨在隔离项目依赖并轻松管理依赖。不同项目可能需要不同版本的Python或第三方库，通过创建独立的Conda环境可避免依赖冲突，并允许为每个环境单独安装、更新或删除库而不影响其他环境。Jupyter Notebook或Jupyter Lab是交互式编程工具，本身不知道Conda环境的存在。为让Jupyter能使用某个Conda环境中的Python和库，需将该环境注册为Jupyter内核。内核是Jupyter的执行引擎，每个Jupyter实例会绑定到一个内核，负责运行代码并返回结果，个人理解这里的内核相当于特定的门或钥匙，将jupyter lab和conda环境打通。即使创建了Conda环境，Jupyter也不会自动识别，所以需手动将Conda环境注册为Jupyter内核，以便在Jupyter中选择特定环境对应的内核，确保代码运行在正确环境中，并能访问该环境中安装的库（如pandas、numpy等）。若觉得迷糊也没关系，简单来说conda环境和jupyter内核相互独立，使用者只需知道创建环境后注册内核就能让它们“协同工作”（不然会调用base中的库）。
在python中确认h5ad文件
安装必要的库
# 安装scanpy
pip install scanpy
加载库
import scanpy as sc
import os
# 确认路径
os.getcwd()
# '/Users/zaneflying/Desktop/PythonGSE188711'
读取数据
# 读取数据
adata = sc.read_h5ad('scRNA_V5.h5ad')
adata
# AnnData对象，n_obs × n_vars = 5912 × 21308
# obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'percent_mito', 'percent_ribo', 'percent_hb', 'seurat_clusters', 'RNA_snn_res.0.1', 'RNA_snn_res.0.2', 'RNA_snn_res.0.3', 'RNA_snn_res.0.5', 'RNA_snn_res.0.8', 'RNA_snn_res.1', 'celltype', 'location', 'pANN_0.25', 'contamination', 'RNA_snn_res.0.4', 'RNA_snn_res.0.6', 'RNA_snn_res.0.7', 'RNA_snn_res.0.9', 'RNA_snn_res.1.1', 'RNA_snn_res.1.2', 'RNA_snn_res.1.3', 'RNA_snn_res.1.4', 'RNA_snn_res.1.5', 'RNA_snn_res.1.6', 'RNA_snn_res.1.7', 'RNA_snn_res.1.8', 'RNA_snn_res.1.9', 'RNA_snn_res.2'
# var: 'vf_vst_counts_mean', 'vf_vst_counts_variance', 'vf_vst_counts_variance.expected', 'vf_vst_counts_variance.standardized', 'vf_vst_counts_variable', 'vf_vst_counts_rank', 'var.features', 'var.features.rank'
# obsm: 'X_harmony', 'X_pca', 'X_umap'
5.导入(h5ad数据转换为Seurat对象)
R语言部分
rm(list = ls())
library(sceasy)
library(qs)
library(reticulate)
library(Seurat)
library(BiocParallel)
register(MulticoreParam(workers = 4, progressbar = TRUE))
6.数据转换
# h5ad转为Seurat
sceasy::convertFormat(obj = "scRNA_V5.h5ad",
from="anndata",to="seurat",
outFile = 'scRNA_CD4.rds')
# X -> counts
# 一个Seurat类的对象
# 21308个特征分布在5912个样本中，属于1个分析测定
# 活跃的测定: RNA (21308个特征，0个可变特征)
# 存在2个层: counts, data
# 计算了3个维度约简: harmony, pca, umap
7.数据确认
scRNA <- readRDS("./scRNA_CD4.rds")
rownames(scRNA)
colnames(scRNA)
数据读取流程(Python)
1.创建和激活环境
# 创建名为sc的python版本为3.9的环境
# -n是-name的简称
conda create -n sc python=3.9
# 激活
conda activate sc
2.在Jupyter lab中创建内核
python -m ipykernel install --user --name sc --display-name "Python(sc)"
3.安装库
conda install scanpy
正式移动和读取文件(Python)
1.导入所需的库
import os
import shutil
import re
2.查询当前的路径
os.getcwd()
#'/Users/zaneflying/Desktop/PythonGSE188711'
3.列出目录下所有包含 'features' 的文件
fs = [
os.path.join('GSE188711_RAW', f) # 拼接出完整路径，如 GSE188711_RAW/GSM5688707_features.tsv.gz
for f in os.listdir('GSE188711_RAW') # 遍历目录下的所有文件名
if 'features' in f # 只选取文件名中包含 'features' 的文件
]
print(fs)
samples1 = fs
# ['GSE188711_RAW/GSM56887********_features_JCA.tsv.gz', 'GSE188711_RAW/GSM56887********_features_WGC.tsv.gz',
# 'GSE188711_RAW/GSM56887********_features_R_CRC4.tsv.gz', 'GSE188711_RAW/GSM56887********_features_RS-CRC1.tsv.gz',
# 'GSE188711_RAW/GSM56887********_features_R_CRC3.tsv.gz', 'GSE188711_RAW/GSM56887********_features_LS-CRC3.tsv.gz']
4.用正则表达式替换字符串，得到样本名
samples2 = [re.sub(r'_features', '', f) for f in fs]
samples2 = [re.sub(r'\.tsv\.gz$', '', f) for f in samples2]
print(samples2)
其实这里的样本名可以清洗得更简洁，比如只保留GSM号码，简洁有诸多好处，建议保持简洁。
samples2 = [re.findall(r'GSM\d+', f)[0] for f in fs]
print(samples2)
# ['GSM5688707', 'GSM5688706', 'GSM5688711', 'GSM5688709', 'GSM5688710', 'GSM5688708']
5.创建文件夹并重命名和移动文件
# 3. 创建 outputs 文件夹（如果不存在）
os.makedirs('outputs', exist_ok=True)
# 4. 遍历样本，创建子文件夹并复制文件
for x, y in zip(samples2, samples1):
output_dir = os.path.join('outputs', x)
os.makedirs(output_dir, exist_ok=True)
# 复制并重命名 features 文件
shutil.copy(y, os.path.join(output_dir, 'features.tsv.gz'))
# 对应的 matrix 文件名替换
matrix_file = y.replace('features', 'matrix').replace('.tsv', '.mtx')
shutil.copy(matrix_file, os.path.join(output_dir, 'matrix.mtx.gz'))
# 对应的 barcodes 文件名替换
barcodes_file = y.replace('features', 'barcodes')
shutil.copy(barcodes_file, os.path.join(output_dir, 'barcodes.tsv.gz'))
相关文章
