目前已完成基于R语言的单细胞分析实战系列,涵盖初级、中级和高级三个阶段(每个阶段暂时各包含4篇内容,后续将陆续更新)。接下来,将尝试开展基于Python的单细胞分析实战系列。但在系列内容正式开始之前,笔者认为有必要先学习一下数据转化和读取流程。
本次内容涉及到的工程文件可通过网盘获得:中级篇2,链接: https://pan.baidu.com/s/1y-HHLXoXsJbgWKCdz26-gQ 提取码: yx93 。
此外,可以向“生信技能树”公众号发送关键词‘单细胞’,直接获取Seurat V5版本的完整代码。
首先推荐python官网:https://www.python.org/doc/
里面的教程有多语言版本,可以选择中文
快速检索就可以查询想要的信息,比如查了一下os库
这个网站可以作为工具手册帮助使用者快速上手,当然大模型辅助也一定是必不可少的。
后续一定会涉及到数据的相互转化,因此必须学会这个步骤。
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
# An object of class Seurat
# 21308 features across 5912 samples within 1 assay
# Active assay: RNA (21308 features, 2000 variable features)
# 3 layers present: counts, data, scale.data
# 3 dimensional reductions calculated: pca, harmony, umap
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
再次回到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 to AnnData
scRNA_V5[["RNA"]] <- as(scRNA_V5[["RNA"]], "Assay")
sceasy::convertFormat(scRNA_V5, from="seurat", to="anndata",
outFile='scRNA_V5.h5ad')
# AnnData object with 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'
# Warning message:
# In .regularise_df(obj@meta.data, drop_single_values = drop_single_values) :
# Dropping single category variables:orig.ident
首先必须要安装anaconda/miniconda等;
# mkdir ~/xxxx
cd ~/xxxx
jupyter lab
激活环境之后,注册小环境对应的内核, 在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环境,包含特定版本的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 object with 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'
R语言部分
rm(list = ls())
library(sceasy)
library(qs)
library(reticulate)
library(Seurat)
library(BiocParallel)
register(MulticoreParam(workers = 4, progressbar = TRUE))
# h5ad转为Seurat
sceasy::convertFormat(obj = "scRNA_V5.h5ad",
from="anndata",to="seurat",
outFile = 'scRNA_CD4.rds')
# X -> counts
# An object of class Seurat
# 21308 features across 5912 samples within 1 assay
# Active assay: RNA (21308 features, 0 variable features)
# 2 layers present: counts, data
# 3 dimensional reductions calculated: harmony, pca, umap
scRNA <- readRDS("./scRNA_CD4.rds")
rownames(scRNA)
colnames(scRNA)
# 创建一名为sc的python版本为3.9的环境
# -n是 -name的简称
conda create -n sc python=3.9
# 激活
conda activate sc
python -m ipykernel install --user --name sc --display-name "Python(sc)"
conda install scanpy
import os
import shutil
import re
os.getcwd()
#'/Users/zaneflying/Desktop/PythonGSE188711'
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/GSM5688707_features_JCA.tsv.gz', 'GSE188711_RAW/GSM5688706_features_WGC.tsv.gz',
# 'GSE188711_RAW/GSM5688711_features_R_CRC4.tsv.gz', 'GSE188711_RAW/GSM5688709_features_RS-CRC1.tsv.gz',
# 'GSE188711_RAW/GSM5688710_features_R_CRC3.tsv.gz', 'GSE188711_RAW/GSM5688708_features_LS-CRC3.tsv.gz']
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']
# 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'))
6.0 安装库,linux/终端
conda install loompy
pip install numpy==1.24.4
6.1 导入
import os
import scanpy as sc
6.2 检查路径
os.getcwd()
# '/Users/zaneflying/Desktop/PythonGSE188711'
6.3 读取样本路径和信息(排除了“.”开头的隐藏文件)
dir_path = '/Users/zaneflying/Desktop/PythonGSE188711/outputs/'
samples = [s for s in os.listdir(dir_path) if not s.startswith('.')]
print(samples)
# ['GSM5688711', 'GSM5688710', 'GSM5688709', 'GSM5688707', 'GSM5688706', 'GSM5688708']
6.4 读取样本构建sce_list
sce_list = []
for sample in samples:
sample_path = os.path.join(dir_path, sample)
# 读取10X格式数据
adata = sc.read_10x_mtx(sample_path, var_names='gene_symbols', cache=True)
# 过滤:至少5个细胞的基因,至少300个基因的细胞 (对应min.cells和min.features)
sc.pp.filter_genes(adata, min_cells=5)
sc.pp.filter_cells(adata, min_genes=300)
# 给AnnData添加项目名称
adata.obs['project'] = sample
sce_list.append(adata)
# sce_list 是包含多个AnnData对象的列表
sce_list
6.5 合并
adata = sc.concat(sce_list,label='sampleid',index_unique='_')
adata.obs_names_make_unique()
adata
参数index_unique='-'帮助使用者把样本名和细胞名进行连接。
6.6 信息查看
adata.obs_names
adata.obs_keys
adata.var_names
保存h5ad文件
adata.write('GSE188711.h5ad')
保存loom文件
import loompy
adata.write_loom('GSE188711.loom')
致谢:感谢曾老师以及生信技能树团队全体成员。更多精彩内容可关注公众号:生信技能树,单细胞天地,生信菜鸟团等公众号。
注:若对内容有疑惑或者有发现明确错误的朋友,请联系后台(欢迎交流)。更多相关内容可关注公众号:生信方舟 。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。