PySCENIC（一）：python版单细胞转录组转录因子分析

关于单细胞转录组转录因子的分析我们之前在单细胞系列讲过R语言版本的，参考：跟着Cell学单细胞转录组分析(十二):转录组因子分析，但是R语言分析起来速度非常慢，如果你动辄上万的单细胞可能要运行好几周，这显然不现实。pySCENIC则很好的解决了这个问题，分析速度很快。官方教程参考：

https://pyscenic.readthedocs.io/en/latest/

一、软件安装

老样子，还是先说一下安装和分析文件的准备，前面环境的配置和之前cellphonedb一样，如果已经操作过的，可以跳过：

#安装下载及环境设置
# 安装一个conda，为什么安装他可以理解为Rstuido之于R,后期在环境设置、软件安装上很方便。
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# bash安装，按照指引，都选yes，这样一些依赖的python包都安装了。
bash Miniconda3-latest-Linux-x86_64.sh 
#  激活conda环境
source ~/.bashrc

#设置镜像
conda config --add channels r 
conda config --add channels conda-forge 
conda config --add channels bioconda
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes

然后就是安装pySCENIC及环境设置：

#安装pyscenic，并创建分析环境
conda create -n pyscenic python=3.9#创建一个pyscenic 的python环境，pyscenic要求python版本3.6及以上，目前python出到3.9了，我用3.9
conda activate pyscenic #激活pyscenic 环境
#安装依赖包
conda install -y numpy
conda install -y -c anaconda cytoolz
conda install -y scanpy
#安装pyscenic
pip install pyscenic

pip安装即可，如果安装失败或者觉得下载速度太慢，可以使用下面代码，使用与其他软件的pip安装。

pip install cellphonedb -i http://pypi.douban.com/simple/ --trusted-host pypi.douban.com

二、数据文件准备

一些分析用的依赖数据库及文件：

#TF注释
#鼠的
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.mgi-m0.001-o0.0.tbl
#人的
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl

#转录组因子列表
下载地址
https://github.com/aertslab/pySCENIC/tree/master/resources
#人的文件名：hs_hgnc_tfs.txt，复制为txt文件即可
#鼠的文件名：mm_mgi_tfs.txt，复制为txt文件即可

#reference数据库，之前一些网上教程的链接文件已经不行了，因为做了更新，跑的时候会出错，我是根据报错选择了下面的文件
#鼠的
wget https://resources.aertslab.org/cistarget/databases/mus_musculus/mm10/refseq_r80/mc_v10_clust/gene_based/mm10_10kbp_up_10kbp_down_full_tx_clustered.genes_vs_motifs.rankings.feather
#人的
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather

三、分析文件准备

这一步是从我们的R种seurat对象提取表达矩阵，并转化为loom文件。这里使用一个python脚本，转化一下即可。

#文件准备
#pyscenic的输入文件是行为基因名，列为细胞ID的矩阵，所以在seurat对象中导出矩阵的时候需要转置一下，可以用标准化矩阵，也可以用counts矩阵，影响不大！
#表达矩阵、meta----R中进行
write.csv(t(as.matrix(sce@assays$RNA@counts)),file = "sce_exp.csv")
#cellInfo <- sce@meta.data[,c("celltype","nCount_RNA","nFeature_RNA")]
#colnames(cellInfo) <- c('CellType', 'nGene' ,'nUMI')
#head(cellInfo)
#write.csv(cellInfo, file = "cellInfo.csv")

#转化为loom文件，Linux下的python脚本
#编辑脚本
vim trans.py
#输入以下内容
import os, sys
os.getcwd()
os.listdir(os.getcwd())
import loompy as lp;
import numpy as np;
import scanpy as sc;
x=sc.read_csv("sce_exp.csv");#R中导出的表达矩阵
row_attrs = {"Gene": np.array(x.var_names),};
col_attrs = {"CellID": np.array(x.obs_names)};
lp.create("sce.loom",x.X.transpose(),row_attrs,col_attrs)

#保存并退出
#运行trans.py
python trans.py
ls
#这样在文件夹中会出现sce.loom文件，就是接下来输入pyscenic的文件。

这样整个准备工作就完成了，后续的分析很简洁，只需要三步。然后就是可视化了，这个下回分解，觉得分享有用对你有帮助的，点个赞、分享下再走呗！记得关注下，之后不迷路！更多精彩请访问我的公众号《KS科研分享与服务》