该文首发于生信技能树,推文链接:https://mp.weixin.qq.com/s/W23Reg6Hi4XWxpMvfctP8g
明白了基因调控网络的基础知识之后,就可以尝试实际操作一下SCENIC分析。
基因调控网络的基础知识可见推文:https://mp.weixin.qq.com/s/sL_8YFulHsZ42L8G5DyY8w
由于R语言跑SCENIC流程实在是慢,还是决定尝试使用python版的SCENIC进行分析。
这次pySCENIC流程实践中出现了很多报错,花了很大时间和精力去解决。希望这个推文能够帮助到遇到同样问题的你。
同时十分感谢曾老师的鼎力相助!~
如果没有耐心读完报错和解决的过程,可以跳到文末。在文末处提供了完整流程和一些解决报错的小思考。
首先展示一下战友的情况:Macbook Pro, M2芯片,16GB内存,老爷车,但也够用。
分析前数据准备-R studio
1、导入数据
load("sce_epi.Rdata")
table(sce$location)
table(sce$seurat_clusters)
table(sce$celltype)
dim(sce)
# [1] 20930 2899
2、R语言中提取exp矩阵文件
library(Seurat)
library(SCopeLoomR)
# 将sce对象中的数据转置并保存为CSV
counts <- GetAssayData(sce, slot = "counts")
counts <- counts[,sample(1:ncol(counts), 100)]
# [1] 20930 100
## 保存为loom文件
library(loomR)
loom <- create(filename = "data/input_new.loom",
data = counts, overwrite = TRUE)
loom$close()
rm(loom)
后续都会用这个input_new.loom文件数据,如果loom经常出现被锁定的情况就用python语言去转换。
GRN分析报错情况及解决流程
环境配置—终端/linux (不在R studio了哈!!)
先跟着众多公众号的教程进行环境配置, 本次采用python 3.9
conda create -y -n pyscenic_3.9 python=3.9
conda activate pyscenic_3.9
pip install pyscenic
conda install numpy
conda install pandas
GRN分析—终端/linux (不在R studio了哈!!)
step1-基因调控网络的构建
输入:表达矩阵input_new.loom文件,转录因子文件 请注意需要把分析文件和工作路径均调整到同一个文件夹中,当然熟练之后也可以自己设定路径啥的。
# pyscenic grn方式
pyscenic grn \
--num_workers 10 \
--output step1out_grn.tsv \
--method grnboost2 \
input_new.loom \
allTFs_hg38.txt
由于对python是真的非常不熟悉,就尝试依靠大模型。
尝试了这三种方法,都失败了! 具体细节就不放出来了,真是"又臭又长"。
接下来找了一个新的环境配置教程,把python提升到3.10并对软件工具进行版本设定
conda create -y -n pyscenic_env python=3.10
conda activate pyscenic_env
pip install pyscenic==0.12.1
conda install numpy==1.23.5
conda install pandas==1.5.3
运行上述的同样代码,依旧报错。
尝试把pyscenic grn改成arboreto_with_multiprocessing.py。
arboreto_with_multiprocessing.py \
--num_workers 10 \
--output step1out_grn.tsv \
--method grnboost2 \
input_new.loom \
allTFs_hg38.txt
可以运行了!
那为什么pyscenic grn改成arboreto_with_multiprocessing.py就可以运行了呢?
好吧,大模型给的解释其实是浅尝辄止。
step2-调控模块的识别
输入:起始位点的文件,motif和TFs映射关系文件和刚开始构建的loom文件
pyscenic ctx \
step1out_grn.tsv \
hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
--annotations_fname motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
--expression_mtx_fname input.loom \
--mode "dask_multiprocessing" \
--output step2out_ctx.tsv \
--num_workers 10
第二步可以运行,唯一的小问题就是如果电脑不够给力的话,把--mode "dask_multiprocessing" \ 去掉。
step3-AUCell打分 输入: 一开始的loom文件和第二步得到的ctx文件。
pyscenic aucell \
input_new.loom \
step2out_ctx.tsv \
--output out_SCENIC.loom \
--num_workers 10
上述的方法都尝试了,然而还是没能成功!同时还尝试了更换python版本环境等,也未能成功。
突然想到最开始在配置环境的时候尝试的配置是python=3.7的环境(基于曾老师最初的推文),但电脑会提醒由于是基于osx-arm64的构架,所以不能配置python=3.7的环境。因此我在猜想会不会跟电脑的硬件构架的关?
查阅了osx-64 和 osx-arm64 相关的知识,发现确实有一些差异!
1. osx-64:
架构: 针对 x86_64 架构,也就是传统的 Intel 和 AMD 处理器架构。这是 Mac 电脑在 Apple Silicon (M1, M2) 之前使用的主流架构。
适用设备: 通常用于基于 Intel 处理器的 Mac 电脑(比如 MacBook、iMac、Mac Pro 等)。
软件兼容性: 针对 x86_64 架构编译的软件包,使用的是 Intel 指令集。
2. osx-arm64:
架构: 针对 ARM64 架构,也就是 Apple Silicon 芯片(如 M1, M2)。这些芯片基于 ARM 架构,与传统的 Intel x86_64 架构不同。
适用设备: 专为 Apple Silicon 处理器的 Mac 设备设计(如 M1/M2 MacBook Air, MacBook Pro, Mac mini, iMac 等)。
软件兼容性: 针对 ARM64 架构编译的软件包,使用的是 ARM 指令集。Apple Silicon 设备通常可以通过 Rosetta 2 运行 osx-64 编译的软件包,但 osx-arm64 的软件包会有更好的性能和优化。
接下来在构建环境的时候尝试使用osx-64的硬件构架。并且重新配置python3.7版本。
CONDA_SUBDIR=osx-64 conda create -n pyscenic_3.7 python=3.7
conda activate pyscenic_3.7
pip install pyscenic
conda install numpy
conda install pandas
pyscenic grn \
--num_workers 10 \
--output step1out_grn.tsv \
--method grnboost2 \
input_new.loom \
allTFs_hg38.txt
成功运行了!
pyscenic ctx \
step1out_grn.tsv \
hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
--annotations_fname motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
--expression_mtx_fname input_new.loom \
--output step2out_ctx.tsv \
--mode "dask_multiprocessing" \
--num_workers 10 \
--mask_dropouts
pyscenic aucell \
input_new.loom \
step2out_ctx.tsv \
--output out_SCENIC.loom \
--num_workers 10
第三步的时候出现了loompy和pyscenic版本不兼容的问题,因此把loompy降级即可!
pip install loompy==2.0.17
然后重新运行第三步代码
在R studio中稍微试了一下
library(SCopeLoomR)
library(SCENIC)
scenicLoomPath='out_SCENIC.loom'
loom <- open_loom(scenicLoomPath)regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")regulons <- regulonsToGeneLists(regulons_incidMat)7 regulonAUC <- get_regulons_AUC(loom, column.attr.name="RegulonsAUC")
运行成功~
完整流程:
1、数据准备-R studio
load("sce_epi.Rdata")
table(sce$location)
table(sce$seurat_clusters)
table(sce$celltype)
dim(sce)
# [1] 20930 2899
library(Seurat)
library(SCopeLoomR)
# 将sce对象中的数据转置并保存为CSV
counts <- GetAssayData(sce, slot = "counts")
counts <- counts[,sample(1:ncol(counts), 100)]
# [1] 20930 100
## 保存为loom文件
library(loomR)
loom <- create(filename = "data/input_new.loom",
data = counts, overwrite = TRUE)
loom$close()
rm(loom)
2、环境配置及GRN运行-终端/linux
请注意把需要的文件放到同一个文件夹中,或者设置文件路径。
# 环境配置
# 当然如果没有遇到这样问题的,只需要按照常规方法配置环境!!
# conda create -n pyscenic_3.7 python=3.7
CONDA_SUBDIR=osx-64 conda create -n pyscenic_3.7 python=3.7
conda activate pyscenic_3.7
pip install pyscenic
conda install numpy
conda install pandas
# GRN-step1 (两种方式)
pyscenic grn \
--num_workers 10 \
--output step1out_grn.tsv \
--method grnboost2 \
input_new.loom \
allTFs_hg38.txt
# arboreto_with_multiprocessing.py \
# --num_workers 10 \
# --output step1out_grn.tsv \
# --method grnboost2 \
# input_new.loom \
# allTFs_hg38.txt
# step2
pyscenic ctx \
step1out_grn.tsv \
hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
--annotations_fname motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
--expression_mtx_fname input_new.loom \
--output step2out_ctx.tsv \
--mode "dask_multiprocessing" \
--num_workers 10 \
--mask_dropouts
# step3
pip install loompy==2.0.17
pyscenic aucell \
input_new.loom \
step2out_ctx.tsv \
--output out_SCENIC.loom \
--num_workers 10
得到了out_SCENIC.loom文件之后就可以放到R stduio中进行处理了!
小思考:
1、计算机是极其"理性"的,它不会以人或者物的意志而随意变化。也就是说如果它出现了报错,一定是能够追根溯源的,之前能运行的代码,只要按照相应条件一定也是能够运行的。因此遇到报错不要害怕,追根溯源即可!
2、尽量不要遗漏每一个报错并努力去理解每一个报错,解决之道就在“其”中。
3、熟悉基础知识和善用辅助工具这两者是并重的,大模型有它独到的优势,但始终不是万能,可能更多的时候我们还是需要靠自己去学习一些计算机/计算机语言基础知识,在这个基础之上更好的使用辅助工具。如果笔者自己能够更加熟悉硬件/计算机语言/分析框架等基础知识的话,应该能更快的做出调整吧。
参考资料:
1、生信技能树/生信随笔:
https://mp.weixin.qq.com/s/IqtNTMb4Jet0VF2-CzgrfA https://mp.weixin.qq.com/s/pN4qWdUszuGqr8nOJstn8w https://mp.weixin.qq.com/s/mBR3IwWvQDcTOXNwM_YCEg https://mp.weixin.qq.com/s/eRFUFPvNDfcU2kW6RkuLSw https://mp.weixin.qq.com/s/ncSW8EXrpzqD-3b7uXy5Mg
2、生信菜鸟团:
https://mp.weixin.qq.com/s/2czUaSzDWrkMF-p8FKDKqg
3、Ks科研分享与服务:
https://mp.weixin.qq.com/s/PS-iI3gtLhMpiM_09pk6sQ https://www.bilibili.com/video/BV1Yj411r7dE/?spm_id_from=333.999.0.0&vd_source=05b5479545ba945a8f5d7b2e7160ea34
4、科研小徐:
https://www.jianshu.com/p/f92e5a063e7a
5、arboreto_with_multiprocessing.py
https://github.com/aertslab/pySCENIC/blob/master/src/pyscenic/cli/arboreto_with_multiprocessing.py
致谢:感谢曾老师以及生信技能树团队全体成员。
注:若对内容有疑惑或者有发现明确错误的朋友,请联系后台(欢迎交流)。更多内容可关注公众号:生信方舟
- END -