# pyscenic **Repository Path**: mo_xinwu/pyscenic ## Basic Information - **Project Name**: pyscenic - **Description**: 基于 pySCENIC 的单细胞转录因子分析 和 基于 R 的 SCENIC 分析可视化 - **Primary Language**: Unknown - **License**: GPL-3.0 - **Default Branch**: master - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 10 - **Forks**: 0 - **Created**: 2023-03-29 - **Last Updated**: 2025-04-08 ## Categories & Tags **Categories**: Uncategorized **Tags**: None ## README # README ## SCENIC 分析 转录因子 (transcription factors, TFs) 是直接作用于基因组，与特定DNA序列结合 (TFBS/motif) ，调控DNA转录过程的一类蛋白质。转录因子可以调节基因组DNA开放性、募集RNA聚合酶进行转录过程、募集辅助因子调节特定的转录阶段，调控诸多生命进程，诸如免疫反应、发育模式等。因此，分析转录因子表达及其调控活性对于解析复杂生命活动具有重要意义。 SCENIC（single-cell regulatory network inference and clustering）是一个基于共表达和motif分析，计算单细胞转录组数据基因调控网络重建以及细胞状态鉴定的方法。在输入单细胞基因表达量矩阵后，SCENIC经过以下三个步骤完成转录因子分析： - 第一步，GENIE3（随机森林)/GRNBoost (Gradient Boosting) 推断转录因子与候选靶基因之间的共表达模块，每个模块包含一个转录因子及其靶基因，纯粹基于共表达； - 第二步，RcisTatget分析每个共表达模块中的基因，以鉴定enriched motifs，仅保留TF motif富集的模块和targets，构建TF-targets网络，每个TF及其潜在的直接targets gene被称作一个调节因子（Regulons）； - 第三步，AUCelll计算调节因子（Regulons）的活性，这将确定Regulon在哪些细胞中处于“打开”状态。 ## 数据库配置 目前 scenic 分析支持的物种仅为 人、小鼠 和 果蝇，其他物种暂时还没有相应的数据库，无法分析。题外话：如果是这三个物种之外的物种想要做 scenic 分析的话, 可以使用同源基因进行分析, 将要分析的物种的 gene Symbol 转换到 人、小鼠 或 果蝇 的 gene Symbol 来进行 scenic 分析。 1. 转录因子列表文件：人和小鼠 ```shell wget https://github.com/aertslab/pySCENIC/blob/master/resources/hs_hgnc_tfs.txt wget https://github.com/aertslab/pySCENIC/blob/master/resources/mm_mgi_tfs.txt ``` 2. Gene-motif 排名文件：人和小鼠 ```shell # 500bp curl -O https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather curl -O https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather # 10k curl -O https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather curl -O https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather ``` 3. 转录因子的 Motif annotation 文件：人和小鼠 ```shell curl -O https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl curl -O https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.mgi-m0.001-o0.0.tbl ``` 踩过一个坑 `pyarrow.lib.ArrowInvalid: Not a feather file`原因是：脚本写的是 `pyscenic 0.12.1`, 但是我原有的环境是 `pyscenic 0.11.2` 两个版本不一致, 它们的`feather`数据是不通用的, 这一点一定要注意。 ## 软件安装 在 `conda_env/environment.yml` 文件中已经安装了部分软件，可以满足 pySCENIC 分析的基本要求了，但是为了能够兼容 R 版本的 SCENIC 分析及更好的可视化，还需要安装些软件。 ```bash conda env create -f environment.yml ``` ```R # conda activate pyscenic BiocManager::install(c("AUCell", "RcisTarget","GENIE3")) BiocManager::install(c("zoo", "mixtools", "rbokeh")) BiocManager::install(c("DT", "NMF", "ComplexHeatmap", "R2HTML", "Rtsne")) BiocManager::install(c("doMC", "doRNG")) devtools::install_github("aertslab/SCopeLoomR", build_vignettes = TRUE) devtools::install_github("aertslab/SCENIC") ``` ## SCENIC 分析 ### 提取并过滤 gene x cell 矩阵 以 Seurat 对象为中心, 进行操作, 提取与过滤操作并不会对细胞进行过滤, 注意矩阵的行名必须为 symbol ID(gene_short_name), 不能是 gene_id(ENS***). ```R library(Seurat) library(dplyr) source("function.R") matrix_for_pyscenic <- filter_matrix(project = project,minCountsPerGene = 1,minPercent = 0.01) %>% as.data.frame() %>% tibble::rownames_to_column(var = "Symbol") write.csv(matrix_for_pyscenic,file = 'matrix_for_pyscenic.csv', quote=FALSE, row.names=FALSE) ``` |Symbol |CAGCTGGCAAGCGTAG.2|TAGACCATCTATCCCG.2|GTAGTCATCTCCCTGA.1|GACCTGGTCCGCATAA.2| |---------|------------------|------------------|------------------|------------------| |LINC01128| 0| 0| 0| 0| |HES4 | 0| 0| 0| 0| |ISG15 | 0| 0| 0| 0| |TNFRSF4 | 0| 0| 0| 0| ### 运行 pyscenic `pyscenic` 的三步运行过程已经封装在 `run_pyscenic.pl` 中, 按照脚本需求准备好 数据库 和 表达矩阵(csv), 进行运行即可完成, 注意：表达矩阵行为细胞, 列为基因; 如果行为基因, 列为细胞要加参数 `--transpose` 运行。 gene counts 矩阵应该是经过一定的过滤后的矩阵, 运行时间与矩阵大小是直接相关的，测试发现 500MB 的矩阵 `num_workers = 20` 的情况下运行时间为 `2h52min`, 也不算太慢。 ```shell perl run_pyscenic.pl -h Usage: perl run_pyscenic.pl --mtx --outdir --tax [--db_dir db_dir] [--num_workers num_workers] [--transpose ] # Options: # --mtx expression matrix for the single cell experiment, formats are supported: csv(rows=cells x columns=genes). # --transpose transpose the expression matrix (rows=genes x columns=cells). # --outdir result output dir. # --tax the conventional species such as human,mouse. # --num_workers the number of workers to use(default: 16). # --db_dir where tf_file feather_file and tbl_file saved. # -h print help messages. # # Database directory structure # ├── hg19-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather # ├── hg19-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather # ├── hs_hgnc_tfs.txt # ├── mm9-500bp-upstream-10species.mc9nr.genes_vs_motifs.rankings.feather # ├── mm9-tss-centered-10kb-10species.mc9nr.genes_vs_motifs.rankings.feather # ├── mm_mgi_tfs.txt # ├── motifs-v9-nr.hgnc-m0.001-o0.0.tbl # └── motifs-v9-nr.mgi-m0.001-o0.0.tbl ``` ### 可视化 对使用 pyscenic 跑完后的结果, 我们使用 R 进行可视化. ```R library(Seurat) library(tidyverse) library(ComplexHeatmap) library(RColorBrewer) source("function.R") load("project.rdata") AUC_score <- data.table::fread("finalout_auc_score_mtx.xls",sep="\t", header=TRUE,data.table=F) %>% tibble::column_to_rownames(.,var = "Regulon") TF_list <- c("ACAA1(20g)","AHR(250g)","ANXA1(5g)") # 使用 TF_list 的 AUC_score 进行 tsne 和 umap 着色 plot_reductions(project = project, auc_matrix = AUC_score, TF_list = TF_list) # 使用 TF_list 的 AUC_score 绘制 小提琴图 plot_violinplot(project = project, auc_matrix = AUC_score, TF_list = TF_list, group_by = "CellType") # 使用 TF_list 的 AUC_score 绘制 山脊图 plot_ggridgesplot(project = project, auc_matrix = AUC_score, TF_list = TF_list, group_by = "CellType") # 绘制单细胞水平 AUC_score 热图, 可以自行调整颜色范围, 如(-10,0), 具体调整见 function.R 代码 163、164、167 行 plot_sc_heatmap(project = project, auc_matrix = AUC_score, group_by = "CellType") # 按组绘制 AUC_score 热图, 组内取平均值, 可以自行调整颜色范围, 如(-10,0), 具体调整见 function.R 代码 196、197、200 行 plot_group_heatmap(project = project, auc_matrix = AUC_score, group_by = "CellType") ``` ## pyscenic 结果说明 1. 转录因子与细胞的调控得分矩阵(finalout_auc_score_mtx.xls) 该文件为每个转录因子在每个细胞的调控得分，TF/Regulon表示转录因子，括号里面的是该转录因子对应的靶基因数量(g:gene的首字母)；barcode序列为细胞，数值代表每个转录因子在每个细胞的调控得分。 |Regulon |CAGCTGGCAAGCGTAG-2 |TAGACCATCTATCCCG-2 |GTAGTCATCTCCCTGA-1 |GACCTGGTCCGCATAA-2 | |--------------|-------------------|--------------------|---------------------|---------------------| |ARID3A(63g) |0.0 |0.00823045267489712 |0.0036743092298647854|0.0019106407995296885| |BACH1(625g) |0.016 |0.021777777777777778|0.018548148148148147 |0.01405925925925926 | |BACH2(5g) |0.06481481481481481|0.07037037037037037 |0.08333333333333333 |0.12777777777777777 | |BATF(369g) |0.04509183980728696|0.0725685034628124 |0.04815316671685235 |0.039345578640971594 | 2. 可视化 - 单细胞层面热图 可以自行调整颜色范围, 如(-10,0), 具体调整见 function.R 代码 163、164、167 行 - 以 CellType 为分组依据的热图(组内取平均) 可以自行调整颜色范围, 如(-10,0), 具体调整见 function.R 代码 196、197、200 行 - AUC sore 的 tsne 和 umap 着色图, 以 ACAA1(20g) 为例 - 以 CellType 为分组依据的小提琴图 和 山脊图(组内取平均), 以 ACAA1(20g) 为例 3. 分组求平均的 AUC score , 文件名：finalout_auc_group.xls ## 其他 除了基于 AUC score 进行可视化展示外, 还可以使用 SCENIC R package 进行计算 Regulon Specificity Score (RSS) , 进行可视化等等, 但是安装 SCENIC 一直装不上, 这一块就没做了。 具体可以参考官方文档：[SCENIC](http://htmlpreview.github.io/?https://github.com/aertslab/SCENIC/blob/master/inst/doc/SCENIC_Running.html) ## 非人、小鼠、果蝇物种 如果是这三个物种之外的物种想要做 scenic 分析的话, 可以使用同源基因进行分析, 将要分析的物种的 gene Symbol 转换到 人、小鼠 或 果蝇 的 gene Symbol 来进行 scenic 分析。进行同源转换的方式可以使用：R package `homologene` 或者 [ensemble的BioMart数据库](http://asia.ensembl.org/index.html)的网页工具 或者 ensembl 推出的 biomaRt 包。 其实, 从 ensemble 下载的 GTF 文件, 里面就有 ensembl_gene_id 和 symble_name, 偏门物种可能没有, 给的 symble_name 好像就是对应到 人 的 hgnc_symbol。 ```R library(biomaRt) listMarts() # biomart version # 1 ENSEMBL_MART_ENSEMBL Ensembl Genes 109 # 2 ENSEMBL_MART_MOUSE Mouse strains 109 # 3 ENSEMBL_MART_SNP Ensembl Variation 109 # 4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 109 mart <- useMart('ENSEMBL_MART_ENSEMBL') listDatasets(mart) searchDatasets(mart = mart, pattern = "sscrofa") # 一个一个的查找很麻烦, 可以使用 searchDatasets 函数来查找需要的物种 # listMarts() -> useMart('ENSEMBL_MART_ENSEMBL') -> listDatasets(mart) 是用来 查看biomaRt支持ID转换的物种 pig <- useMart("ensembl", dataset = "sscrofa_gene_ensembl") gene <- readLines("ensemble.txt") # ENSSSCG00000002633 # ENSSSCG00000002635 # ENSSSCG00000002637 listAttributes(pig) # 获取indica的mart的描述信息, 重点查看name列, 这里包含了所有类型的基因ID编号类型，如 ensembl_gene_id、bgi_gene、rap_gene 等, filters、attributes 需要填写这里面的name searchFilters(mart = pig, 'ensembl|hgnc') # listAttributes(pig) 的信息太多了, 使用 searchFilters(mart = pig, 'ensembl|hgnc') 来进行查找 name 中包含 ensembl 或 hgnc gene_name <- getBM(attributes=c("hgnc_symbol","ensembl_gene_id"),filters = "ensembl_gene_id",values = gene, mart = pig) # hgnc_symbol ensembl_gene_id # GAS8 ENSSSCG00000002633 # SPIRE2 ENSSSCG00000002635 # DBNDD1 ENSSSCG00000002637 # 这样就拿到了 GTF 文件中, ensembl_gene_id 对应到物种人的 symbol name了, 它其实与 GTF 文件中的 gene_name 字段的内容是一致的, 如有的话. # gene_id "ENSSSCG00000002633"; gene_version "3"; gene_name "GAS8"; gene_source "ensembl"; gene_biotype "protein_coding" # gene_id "ENSSSCG00000002635"; gene_version "3"; gene_name "SPIRE2"; gene_source "ensembl"; gene_biotype "protein_coding" # gene_id "ENSSSCG00000002637"; gene_version "3"; gene_source "ensembl"; gene_biotype "protein_coding" # 没有 gene_name 字段 ``` 综上, 要进行非人、小鼠、果蝇的 scenic 分析, 是不推荐的, 与人、小鼠、果蝇亲缘关系近的还是可以的, 如 猴、猩猩、大鼠, 一定要进行 scenic 分析的话, 可以使用上述方法进行同源转换到人或小鼠, 再把表达矩阵中非同源基因去除, 一方面可以减小计算量, 另一方面数据库中都没有这些基因, 也不会影响结果. ## 2023-06 更新 1. 更新 `conda_env/environment.yml` 文件，更新后的 conda 环境文件可以完成大多数软件的安装，`SCENIC` 除外，需要进入 R 后进行安装。需要注意的是，安装` pyscenic` 一定要注意 `numpy 的版本，不能太高`，在 `numpy 1.24` 后，删除了 `np.object` 数据类型，会导致 pyscenic 报错。 ```R # 安装 SCENIC 软件包 devtools::install_github("aertslab/SCENIC") ``` 2. 关于 R 版本的 SCENIC 还是 python 版本的 SCENIC 的选择，R 版本的中间结果多，运行速度慢，有些结果可能文章中能用；python 版本运行速度快，但没有中间结果，只有最终的 AUCell 的打分矩阵，后续的可视化倒是没问题。就看个人选择了。 3. 其实可以有第三种解决方案，可以结合 R 版本 和 pyscenic 一起使用，在 R 版本中，只有 `runGenie3` 运行时间很长，我们可以使用 `pyscenic grn` 来替换 `runGenie3` 其他步骤耗时相对较少，还是能够接受的。 4. R 版本的 SCENIC(1.3.1) 只能使用以下链接的 RcisTarget 数据库的数据，与 pyscenic(0.12.1)使用的数据是不太一样的，要根据使用的 pyscenic 版本去进行选择 RcisTarget 数据库的数据 ```bash # 500bp curl -O https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather curl -O https://resources.aertslab.org/cistarget/databases/old/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-7species.mc9nr.feather # 10k curl -O https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather curl -O https://resources.aertslab.org/cistarget/databases/old/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-7species.mc9nr.feather ``` 5. 代码更新，新增 `run_scenic.r` 的 R 脚本。 6. 官方教程： - R: https://rdrr.io/github/aertslab/SCENIC/f/vignettes/SCENIC_Running.Rmd - python: https://pyscenic.readthedocs.io/en/latest/tutorial.html