本文记录了在Win10平台通过Rstudio使用reticulate为 Seurat::FindClusters
链接Python环境下的Leidenalg算法进行聚类的实现过程Leidenalg error · Issue #2512 · satijalab/seurat (github.com)。同时,在Seurat和Scanpy流程框架下,对Louvain和Leiden算法在处理10万个细胞样本量的虚拟表达矩阵时的速度表现进行了测试。测试结果表明,声称比Louvain算法有计算速度提升的Leiden算法并没有表现出宣称的速度优势。在大型数据集的搜索速度上,Louvain算法会领先Leiden算法几倍的时间。
:
1、配置win10 Python环境 + Seurat::FindClusters 中链接 Leidenalg
背景: Leiden 在2018年被提出,算法相比 Louvain 优化了社区分割策略,可以产生更高质量的连通社区。在计算速度上也比较 Louvain 更快。Louvain 算法作者也曾表示推荐使用 Leiden algorithm 替代Louvain 进行大规模的社区识别。
$\,\,\,\,\,$安装好conda之后,可以根据需求配置conda channel, 参考 :python版本管理神器miniconda使用指南;
第二步,配置 Anaconda Prompt (anaconda3) 到 Windows terminal 为了方便管理windows平台下的conda env,主要参考教程 【Incorporating an Anaconda Prompt to a Windows Terminal】 和 【将Anaconda Prompt配置进Windows Terminal】
**参考 config.json**
{
"colorScheme": "Tango Dark",
"commandline": "cmd.exe /K C:\Users\<user_name>\anaconda3\Scripts\activate.bat,
"cursorShape": "filledBox",
"experimental.retroTerminalEffect": false,
"guid": "{XXXX-XXXX}",
"hidden": false,
"icon": "C:\Users\<user_name>\anaconda3\Menu\anaconda-navigator.ico",
"name": "Anaconda Prompt",
"scrollbarState": "hidden",
"startingDirectory": "%USERPROFILE%",
"useAcrylic": true
},
第三步,安装好Anaconda之后,在Anaconda Prompt 中创建conda环境和安装scanpy系列软件:
conda create -n scRNA python=3.9
conda activate scRAN
conda install -c conda-forge scanpy python-igraph leidenalg
pip install scanpy
第四步,在Rstudio中调用
reticulate
配置python解释器:
参考:Seurat 4 R包源码解析 22: step10 细胞聚类 FindClusters() | 社群发现
reticulate::use_python("C:\\Users\\<user>\\anaconda3\\envs\\scRNA\\python.exe") ### Windows中的 python为.exe
reticulate::py_module_available("leidenalg") ### 检查leiden是否可用,显示为 True 即代表正确链接到刚刚配置的 scRNA conda环境
Linux conda环境下配置python 代码示例: reticulate::use_python(system("which python", intern = T))
NOTE: Windows系统下每次打开R会话加载 Seurat
之前,必须首先通过以上代码配置 python解释器,然后再跑 Seurat::FindClusters(),不然会无法正确链接到环境。
Linux 系统不存在此问题
第五步,在Rstudio中使用
Seurat::FindClusters()
调用 leiden 算法:
pacman::p_load(Seurat,dplyr)
pbmc_small %>% FindClusters(algorithm = 4,resolution = 0.5)
r-miniconda 题外话:
关于r-miniconda 环境的存在,如果系统不存在conda环境,reticulate 将提供安装基于Miniconda的轻量级R语言环境 r-miniconda环境服务,存在于用户的 "C:/Users//AppData/Local/r-miniconda" 目录下。如果不需要使用,可以简单通过以下命令快速卸载:reticulate::miniconda_uninstall(path = reticulate::miniconda_path())
, 这样在 Anaconda Prompt 通过 "conda env list"查询环境的时候就不会显示 r-miniconda内建立的环境了。
2、Louvain vs Leiden 算法速度表现测试
2.1 Seurat 测试5w随机细胞数下 Louvain vs Leiden 聚类速度
虽然 Leiden 是 Louvain的优化版本,但 不建议在R平台下链接python的 Leiden ,实测使用 5w细胞数时该算法会抛出异常 long vectors not supported yet: ../../src/include/Rinlinedfuns.h:537 。
测试代码
### Prepare env
pacman::p_load(Seurat,dplyr)
reticulate::use_python(system("which python", intern = T))
reticulate::py_module_available("leidenalg")
### PrePare Test OBT
C_50k <- Matrix::Matrix(ncol = 5e4, data = rpois((1e4) * (5e4),lambda = 0.3), sparse = TRUE)
rownames(C_50k) <- paste0("g.",1:nrow(C_50k)); colnames(C_50k) <- paste0("c.",1:ncol(C_50k))
C_50k <- CreateSeuratObject(counts = C_50k) %>%
NormalizeData(verbose = 0) %>% FindVariableFeatures(verbose = 0) %>%
ScaleData(verbose = 0) %>% RunPCA(verbose = 0) %>%
FindNeighbors(dims = 1:30)
### Test Clustering Algorithm
###############################################################################################################
Louvain_.time <- system.time({cur_seu %>% FindClusters(resolution = 1,algorithm = 1,verbose = 0)})[["elapsed"]] # 1 = original Louvain algorithm;
LouvainM.time <- system.time({cur_seu %>% FindClusters(resolution = 1,algorithm = 2,verbose = 0)})[["elapsed"]] # 2 = Louvain algorithm with multilevel refinement;
SLM_____.time <- system.time({cur_seu %>% FindClusters(resolution = 1,algorithm = 3,verbose = 0)})[["elapsed"]] # 3 = SLM algorithm;
Leiden__.time <- system.time({cur_seu %>% FindClusters(resolution = 1,algorithm = 4,verbose = 0)})[["elapsed"]] # 4 = Leiden algorithm
2.2 Scanpy 测试10w随机细胞数下 Louvain vs Leiden 聚类速度
在面对10万个细胞的数据集时,本测试发现在Sanpy中使用Louvain聚类的耗时最短,仅为8秒。相比之下,同样在该平台下使用Leiden算法的速度表现略显逊色,用时接近1分钟,并没有宣传中提到的速度优化表现。受到R语言本身特性的影响,Seurat中同样使用Louvain进行聚类的耗时高达30分钟,相比Python平台差了100余倍。
测试代码
{r}
pacman::p_load(Seurat,dplyr)
#### PrePare Test OBT
C_100k <- Matrix::Matrix(ncol = 1e5, data = rpois((1e4) * (1e5),lambda = 0.3), sparse = TRUE)
rownames(C_100k) <- paste0("g.",1:nrow(C_100k)); colnames(C_100k) <- paste0("c.",1:ncol(C_100k))
Seurat_100k <- CreateSeuratObject(counts = C_100k) %>%
NormalizeData(verbose = 0) %>% FindVariableFeatures(verbose = 0) %>%
ScaleData(verbose = 0) %>% RunPCA(verbose = 0) %>%
FindNeighbors(dims = 1:30)
system.time({Seurat_100k %>% FindClusters(resolution = 1,algorithm = 1,verbose = 0)})[["elapsed"]]
#### Format Seurat to H5
dior::write_h5(Seurat_100k,file = "Seurat_100k.h5")
Switch to Python -------------------------------------------------
{py}
import scanpy as sc
import diopy
import h5py
### Format H5 to Adata
h5 = h5py.File('Seurat_100k.h5', 'r')
adata = diopy.input.h5_to_adata(h5 = h5, assay_name='RNA')
h5.close()
### Add SNN key_dict that was missing when dior was written out h5
adata.uns['neighbors'] = {'connectivities_key': 'connectivities','distances_key': 'distances','params': {'n_neighbors': 20, 'method': 'umap', 'random_state': 0, 'metric': 'euclidean'}}
#__________________________________________
from timeit import default_timer as timer
start = timer()
sc.tl.leiden(adata, resolution = 1, key_added = "leiden_cluster")
print("leiden Time Elaped: ",timer() - start)
#__________________________________________
from timeit import default_timer as timer
start = timer()
sc.tl.louvain(adata, resolution = 1, key_added = "louvain_cluster")
print("louvain Time Elaped: ",timer() - start)
Reference
Preprocessing and clustering 3k PBMCs — Scanpy documentation (scanpy-tutorials.readthedocs.io)
Installing Python Packages • reticulate (rstudio.github.io)
单细胞转录组经典聚类算法:Louvain与Leiden
Scanpy 与 Seurat 数据互通