最近可算拿到了黑神话悟空的全成就，连续两周目的游戏实在是有点肝疼了。相信不少热衷3A游戏的朋友为了玩黑神话也入手了几款中高端显卡，游戏通关了，显卡可不能闲着，必须拿它盘点单细胞数据。
简单介绍下rapids-singlecell：
Rapids-singlecell 作为 scanpy 的替代品，提供增强的单细胞数据分析，并且同时整合了 squidpy 和 decoupler 的部分功能。它通过 cupy 和 Nvidia 的 RAPIDS 实现 GPU 计算，强调高效计算，并且最终作为 scverse 生态系统的一部分，通过社区协作不断保持兼容性和发展。此外rapids-singlecell还有几个独有的高变基因筛选方法：person_residual和poisson_gene_selection。
简而言之rapids-singlecell具有以下特性：
GPU 优化：利用 GPU 加速大数据集的处理。
集成 scverse 库：整合 scanpy 的主要功能，并增加 squidpy 和 decoupler 的特性。
简单的安装过程：可通过 Conda 和 PyPI 获取。
整体而言rapids-singlecell使用方法和scanpy几乎一致，可视化部分完全是调用scanpy函数完成，个别大计算量的部分转存到GPU的 RAM完成。差不多20w细胞数量的话，16G 显存的显卡刚好能够完成，时间差不多在10min左右的，还是很推荐高细胞数量的同学尝试使用下rsc的。
部署rapids-singlecell
git clone https://github.com/scverse/rapids_singlecell.git#mamba 创建虚拟环境#推荐24.06版本，方便部署pytorch使用poisson_gene_selection筛选高边基因。mamba env create -f conda/rsc_rapids_24.06.yml #default CUDA-11.8# ormamba env create -f conda/rsc_rapids_24.08.yml #default CUDA-12.5
激活环境
conda activate rapids_singlecel加载包
# 导入 scanpy 库，用于单细胞基因组数据的处理import scanpy as sc# 导入 cupy 库，用于在 GPU 上进行加速计算import cupy as cp# 导入 time 库，用于计算运行时间import time# 导入 rapids_singlecell 库，用于增强的单细胞数据分析import rapids_singlecell as rsc# 导入 warnings 库并忽略警告信息，清理输出import warningswarnings.filterwarnings("ignore")# 导入 rmm 库（RAPIDS 内存管理），设置自定义内存分配器import rmmfrom rmm.allocators.cupy import rmm_cupy_allocator# 初始化 RMM（RAPIDS Memory Manager）进行 GPU 内存管理rmm.reinitialize(managed_memory=False, # 禁用托管内存，避免超限pool_allocator=False, # 禁用内存池分配器（默认值为 False）devices=0, # 注册 GPU 设备 ID，这里仅注册 GPU 0)# 将 cupy 的内存分配器设置为 RMM 的分配器，以便更高效的 GPU 内存管理cp.cuda.set_allocator(rmm_cupy_allocator)
加载预处理数据
# 计时代码执行时间%%time# 读取单细胞数据文件（.h5ad 格式），并存储为 AnnData 对象adata = sc.read("h5/adata.raw.h5ad")# 将 AnnData 对象转移到 GPU 上以加速计算%%timersc.get.anndata_to_GPU(adata)# 标记基因家族，找到以 "MT-" 为前缀的基因并归类到 "MT" 家族%%timersc.pp.flag_gene_family(adata, gene_family_name="MT", gene_family_prefix="MT-")# 标记基因家族，找到以 "RPS" 为前缀的基因并归类到 "RIBO" 家族%%timersc.pp.flag_gene_family(adata, gene_family_name="RIBO", gene_family_prefix="RPS")# 计算质量控制指标，基于已标记的 "MT" 和 "RIBO" 基因家族%%timersc.pp.calculate_qc_metrics(adata, qc_vars=["MT", "RIBO"])# 绘制散点图，展示 total_counts 和 pct_counts_MT 的关系sc.pl.scatter(adata, x="total_counts", y="pct_counts_MT")# 绘制散点图，展示 total_counts 和 n_genes_by_counts 的关系sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts")# 绘制小提琴图，显示各患者（PatientNumber）中的 n_genes_by_counts 分布，添加抖动效果sc.pl.violin(adata, "n_genes_by_counts", jitter=0.4, groupby="PatientNumber")# 绘制小提琴图，显示各患者中的 total_counts 分布，添加抖动效果sc.pl.violin(adata, "total_counts", jitter=0.4, groupby="PatientNumber")# 绘制小提琴图，显示各患者中的 pct_counts_MT 分布，添加抖动效果sc.pl.violin(adata, "pct_counts_MT", jitter=0.4, groupby="PatientNumber")# 根据基因数进行数据筛选，保留基因数在 200 到 5000 之间的细胞%%timeadata = adata[adata.obs["n_genes_by_counts"] > 200]adata = adata[adata.obs["n_genes_by_counts"] < 5000]adata.shape# 根据线粒体基因比例筛选，保留线粒体基因比例小于 20% 的细胞%%timeadata = adata[adata.obs["pct_counts_MT"] < 20]adata.shape# 对基因进行过滤，保留至少出现 3 次的基因%%timersc.pp.filter_genes(adata, min_count=3)# 将原始计数数据复制到 AnnData 对象的 "counts" 层中，便于进一步分析adata.layers["counts"] = adata.X.copy()
标准化数据
归一化总计数，使每个细胞的总基因表达量达到指定的目标值（1e4），以便进行比较%timersc.pp.normalize_total(adata, target_sum=1e4)# 对数据进行对数转换（log1p），即 log(1 + x)，使数据更适合下游分析%timersc.pp.log1p(adata)
挑选高变基因
# 标记高度变异的基因，使用 Seurat v3 方法，每批次提取前 5000 个变异较大的基因（按患者编号分批次）%%timersc.pp.highly_variable_genes(adata,n_top_genes=5000,flavor="seurat_v3",batch_key="PatientNumber",layer="counts",)# 将原始数据存储到 adata.raw 中，便于后续回溯到未处理的基因表达数据%%timeadata.raw = adata# 根据高度变异的基因筛选数据，仅保留这些基因用于下游分析%%timeadata = adata[:, adata.var["highly_variable"]]# 回归处理，去除 total_counts 和 pct_counts_MT 的影响，以减少潜在的批次效应%%timersc.pp.regress_out(adata, keys=["total_counts", "pct_counts_MT"])
数据缩放
标准化数据，使每个基因的平均值为 0，标准差为 1，并限制最大值为 10，以减少异常值的影响%timersc.pp.scale(adata, max_value=10)
主成分分析
# 进行主成分分析（PCA），提取前 100 个主成分，用于降维和特征提取%%timersc.pp.pca(adata, n_comps=100)# 绘制主成分方差比图（PCA 解释方差），展示前 100 个主成分的对数方差贡献sc.pl.pca_variance_ratio(adata, log=True, n_pcs=100)# 将 AnnData 对象从 GPU 转移到 CPU，指定 convert_all=True 以确保所有数据转移%%timersc.get.anndata_to_CPU(adata, convert_all=True)
批次效应矫正
使用 Harmony 算法进行数据整合，以 "PatientNumber" 为批次键，调整不同患者间的批次效应%timersc.pp.harmony_integrate(adata, key="PatientNumber")
聚类可视化
# 计算邻居关系，基于 40 个主成分，并使用 15 个邻居，用于后续的降维和聚类分析%%timersc.pp.neighbors(adata, n_neighbors=15, n_pcs=40)# 进行 UMAP 降维，得到二维表示，用于可视化%%timersc.tl.umap(adata)# 使用 Louvain 算法进行聚类，设置分辨率为 0.6%%timersc.tl.louvain(adata, resolution=0.6)# 使用 Leiden 算法进行聚类，设置分辨率为 0.6%%timersc.tl.leiden(adata, resolution=0.6)# 绘制 UMAP 图，显示 Louvain 和 Leiden 聚类结果%%timesc.pl.umap(adata, color=["louvain", "leiden"], legend_loc="on data")# 绘制 UMAP 图，显示细胞类型（CellType）和患者编号（PatientNumber）信息sc.pl.umap(adata, color=["CellType", "PatientNumber"])# 计算基于患者编号（PatientNumber）的嵌入密度，并进行可视化%%timersc.tl.embedding_density(adata, groupby="PatientNumber")sc.pl.embedding_density(adata, groupby="PatientNumber")# 进行 t-SNE 降维，使用前 40 个主成分进行可视化%%timersc.tl.tsne(adata, n_pcs=40)# 使用 K-means 算法进行聚类，设置聚类数为 8rsc.tl.kmeans(adata, n_clusters=8)# 绘制 t-SNE 图，显示 K-means 聚类结果%%timesc.pl.tsne(adata, color=["kmeans"], legend_loc="on data")# 绘制 t-SNE 图，显示细胞类型（CellType）和患者编号（PatientNumber）信息sc.pl.tsne(adata, color=["CellType", "PatientNumber"])
差异基因
# 使用逻辑回归（logreg）方法对不同细胞类型（CellType）进行差异基因分析，使用非原始数据（use_raw=False）%%timersc.tl.rank_genes_groups_logreg(adata, groupby="CellType", use_raw=False)# 绘制前 20 个差异基因的图表%%timesc.pl.rank_genes_groups(adata, n_genes=20)
双细胞过滤
rsc.pp.scrublet(adata,batch_key="PatientNumber")rsc的单细胞的部分就介绍到这边了，主要的功能和模块就这些部分了。当然后续有时间会和大家分享rsc包的空转部分和富集分析模块，其中的富集分析还是很值得学习一下的。最后祝愿大家在这条探索单细胞世界的道路上，始终保持初心，勇往直前，最终实现知识的圆满，完成学术的理想，全始全终！
