pyscenic做为scenic的python版本，相比R版本实现了分析速度上的巨大提升，尤其是最耗时的共表达网络构建步骤。实测时间如下，数据集1 (4257 cell x 36601 gene)，10 cores ：grn 网络构建2小时；ctx 网络纯化16分钟；aucell 细胞活性1分钟。另外，数据集2 (45000 cell x 33469 gene)，三个步骤的耗时分别为4小时、20分钟、5分钟。
从上面的流程示意图可知，整个转录因子分析步骤分成Co-expression、Motif-discovery、Cell-scoring，分别对应pyscenic软件的三个子命令grn、ctx、aucell。
RcisTarget数据库文件下载链接：
- TF：https://github.com/aertslab/pySCENIC/tree/master/resources
- Feather：motif-gene rank，https://resources.aertslab.org/cistarget/databases
- Tbl： motif-TF annotation，https://resources.aertslab.org/cistarget/motif2tf
测试数据来自《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》，从原始数据出发，降维聚类分群后，使用fibo细胞亚群做转录因子分析。
1、grn
输入表达矩阵和转录因子列表，软件基于共表达构建转录因子与潜在靶基因的调控网络。构建网络的算法可以选择GENIE3和GRNBoot2 (默认参数)，网络构建涉及随机过程，如果想要复盘时保持结果一致可以设置随机种子。
另外，表达矩阵的格式可以接受csv和loom，如果是csv (rows=cells x columns=genes) 则需要参数--transpose。表达矩阵在整个过程都有用到，为了方便可以使用loom格式。
pyscenic grn --seed 21 \
--num_workers 10 \
--method grnboost2 \
--output fibo_grn.tsv fibo_count.loom allTFs_hg38.txt
2、ctx
上一步得到了初始的调控网络，这一步基于motif和TF的关系以及motif对基因调控潜能的排序来修剪初始网络。首先，将motif注释到TF；接着，对motif调控的基因排序(feather格式的文件)，如果调控网络里面出现靶基因则折线上升，如此就会形成下图所示的曲线图，默认选择所有基因的前5%处 (标虚线的地方) 曲线下的面积计算AUC值，再接着计算一个标准化的AUC值做为NES，默认过滤阈值为3.0；最后，根据三种方式进一步修剪网络：1. thresholds 根据阈值选择motif调控基因排序的前百分比的靶基因，2. top_n_targets 每个TF保留前N个靶基因，3. top_n_regulators 每个靶基因保留前N个TF。此时，得到的TF与靶基因的调控网络叫做regulon (调控子)，其中保留的靶基因上游含有motif直接结合的位点。
其实，不难看出这一步包含的过程比较繁琐，经过各种修剪获得最终的网络。所以，这一步相当重要，可以适当调整阈值得到符合实际需求的网络。提示：这一步提到的实现过程为个人理解。
feather=hg38_refseq-r80_10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather
fname=motifs-v9-nr.hgnc-m0.001-o0.0.tbl
pyscenic ctx --annotations_fname $fname \
--expression_mtx_fname fibo_count.loom \
--mode dask_multiprocessing \
--mask_dropouts \
--num_workers 10 \
--output fibo_ctx.gmt fibo_grn.tsv $feather
3、aucell
最后，评估每个regulon在所有细胞里面的活性，富集方法与上面的motif富集相似。需要注意的是该步骤也有随机种子参数，如果想保证后面复现结果最好设置一下参数。
pyscenic aucell --seed 21 \
--num_workers 10 \
--output fibo_aucell.csv fibo_count.loom fibo_ctx.gmt
到此，pyscenic流程就结束了，后面的工作就是基于regulon活性矩阵探索数据了。当然，也可以将活性矩阵二值化，这样每个regulon在细胞中只有on和off两种状态，可以使用R包AUCell来完成。
library(AUCell)
library(reshape2)
auc <- read.table('fibo_aucell.csv',header=T,row.names=1,sep=',',stringsAsFactors=F,check.names=F)
auc <- t(auc)
rownames(auc) <- gsub(',','', rownames(auc))
auc[1:3, 1:5]
case1_AAACCTGGTAGTACCT-1 case1_AAAGCAAAGCCTTGAT-1
AHR(+) 0.0680935 0.07413331
ARID3A(+) 0.1266152 0.11878014
ARNT2(+) 0.0000000 0.00000000
case1_AAAGTAGGTCCATCCT-1 case1_AAATGCCGTTCCACTC-1
AHR(+) 0.07159951 0.07028004
ARID3A(+) 0.05228222 0.09272742
ARNT2(+) 0.00000000 0.00428051
case1_AACTCCCAGTAGTGCG-1
AHR(+) 0.0655522976
ARID3A(+) 0.0814689810
ARNT2(+) 0.0004553734
cells_assignment <- AUCell_exploreThresholds(auc, plotHist=F, assignCells=T)
thresholds <- getThresholdSelected(cells_assignment)
regulonsCells <- setNames(lapply(names(thresholds),
function(x) {
trh <- thresholds[x]
names(which(auc[x,]>trh))
}),names(thresholds))
regulonActivity <- melt(regulonsCells)
binaryRegulonActivity <- t(table(regulonActivity[,1], regulonActivity[,2]))
binaryRegulonActivity[1:3, 1:5]
case1_AAACCTGAGACAGAGA-1 case1_AAACCTGGTAGTACCT-1
ARID3A(+) 0 0
ARNT2(+) 0 0
ARNTL(+) 0 0
case1_AAAGCAAAGCCTTGAT-1 case1_AAAGTAGCACAACTGT-1
ARID3A(+) 0 0
ARNT2(+) 0 0
ARNTL(+) 0 0
case1_AAAGTAGGTCCATCCT-1
ARID3A(+) 0
ARNT2(+) 0
ARNTL(+) 0
最后的最后，展示一下重现文章的结果 (右边来自文章)，虽然缺失两个转录因子，可能是前期数据处理不同导致，但是其余的转录因子在两个细胞亚群里面的趋势与文章一致。
往期回顾
一网打尽scRNA矩阵格式读取和转化(h5 h5ad loom)
ggplot2 | 开发自己的画图函数
R包安装的4种姿势
clusterProfiler: No gene can be mapped | 怎么破？
R语言的碎碎念
