pySCENIC-linux

By Kaiyi Fu

安装+配置环境

#在linux上安装anaconda（略）

#设置镜像
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，并创建分析环境
conda create -n pyscenic python=3.7#创建一个pyscenic 的python环境，pyscenic要求python版本3.6及以上，目前python出到3.9了，我用3.9
conda activate pyscenic #激活pyscenic 环境
#安装依赖包
pip install pyscenic
pip install scanpy

#后面第一步grn时报错numpy与python版本不符，直接在这一步安装适配版本
#conda install -y numpy==1.21.3

下载数据

# 数据下载
#TF注释
#鼠
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.mgi-m0.001-o0.0.tbl
#人
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl

#转录组因子列表
#人
wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt

#鼠
wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_mm.txt

#整体基因组排名数据库feather格式的ranking排名数据库
#鼠的
wget https://resources.aertslab.org/cistarget/databases/mus_musculus/mm10/refseq_r80/mc_v10_clust/gene_based/mm10_10kbp_up_10kbp_down_full_tx_v10_clust.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 <- [e-mail][,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)

#按Esc，再输入:wq，保存并退出
:wq
# Vim 中保存编辑的操作步骤：
#进入编辑模式：默认是命令模式（Normal mode），切换到编辑模式（Insert mode）才能输入文本。按下 i 键进入插入模式（尾行显示INSERT），可以编辑文件内容。
#保存文件：按下 Esc 键回到命令模式（Normal mode），尾行不显示INSERT
#输入 :w，然后按回车键。这会将你的文件保存到当前的文件名（这里是 trans.py）。
#退出 Vim：输入 :q，然后按回车键。这会退出 Vim 编辑器。
#有修改过文件但不想保存，可以在命令模式下直接输入 :q!，然后按回车键。这会强制退出 Vim，丢弃所有未保存的修改。
#命令模式下使用 :wq 保存文件并退出 Vim。
#修改脚本：vim进入 修改后:wq保存退出


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

三段代码

#以鼠的为例，在linux终端运行

##开始分析
##第一步GRN---运行完得到sce.adj.csv文件
pyscenic grn --num_workers 10 \
--sparse \
--method grnboost2 \
--output sce.adj.csv \
sce.loom \
allTFs_mm.txt


#推断转录因子与提供的表达矩阵基因的共表达模块，基于grnboost2，R中时GENIE3
#线程数--num_workers 10
#输入转换好的表达矩阵sce.loom、转录因子列表mm_mgi_tfs.txt
#输出sce.adj.csv文件作为第二步输入
###erro：AttributeError: module 'numpy'has no attribute 'object'
#'np.object'was a deprecated alias for the builtin 'object', To avoid this error in existing code, use object by itself. Doing this will not modiy any behavior and is safe.
#numpy版本不符，卸载装低版本
#pip uninstall numpy
#pip install -U numpy==1.21.3
#查看帮助函数，看看各个参数的具体解释
pyscenic grn --help

##分析第二步：RcisTarget---运行完得到sce.regulons.csv文件
pyscenic ctx --num_workers 24 \
--output sce.regulons.csv \
--expression_mtx_fname sce.loom \
--all_modules \
--mask_dropouts \
--mode "dask_multiprocessing" \
--min_genes 10 \
--annotations_fname motifs-v10nr_clust-nr.mgi-m0.001-o0.0.tbl \
sce.adj.csv \
mm10_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather

#进行TF-motif富集分析，识别直接靶标
#得到转录因子(TF)与其对应的直接作用的靶点,称为regulon(每一个regulon是1个TF和其调控的靶基因)
#输出sce.regulons.csv作为第三步输入文件
#输入：转换好的表达矩阵sce.loom文件、motif注释文件、基因组rank文件、第一步得出的sce.adj.csv文件

##分析第三步：AUCell---运行完得到sce_SCENIC.loom文件，即分析结果
pyscenic aucell --num_workers 24 \
--output sce_SCENIC.loom \
sce.loom \
sce.regulons.csv

#使用AUCell对每个细胞的每个regulon活性进行评分。
#输入转换好的表达矩阵sce.loom和第二步得到的sce.regulons.csv
#输出sce_SCENIC.loom文件