今天写一篇关于单细胞转录因子的分析-scenic
在做分析之前你应该看一下这篇文献，看看别人做了什么scenic分析
《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》

分析原理

对于这个包我想发表一下感想，这个包刚开始对于我来说比较难理解，其实跟着官方教程running无非就是debug而已，这没啥，重要的是这个分析做了些啥，最后得到啥，我能讲啥，这是最难的，反正我是不敢做那些我自己都不了解过程的分析，所以首先我们要了解这个包的分析流程：
GENIE3或GRNBoost→RcisTarget→AUCell（所以首先你得先了解这三个包做了啥，不过不太想了解也可以往后看，我会讲讲我的理解，希望对你们有所帮助）

1.使用GENIE3或GRNBoost(Gradient Boosting)基于共表达推断转录因子与候选靶基因之间的共表达模块
https://www.jianshu.com/p/d71dcd4cff5a (可以学一下这个包，体验一下，或许我理解的不对哈）
是不是看到这个就头大 ，所有的教程都是这样，其实简单的说第一步推断出转录因子(TF)与其作用的靶点/gene(一个转录因子可能调控很多的gene哈)

2.RcisTarget对每个共表达模块进行顺式调控基序(cis-regulatory motif）分析
https://www.jianshu.com/p/c9df97f9e6e0 (RcisTarget包）
由于GENIE3模型只是基于共表达，会存在很多假阳性和间接靶标，为了识别直接结合靶标(direct-binding targets)，使用RcisTarget对每个共表达模块进行顺式调控基序(cis-regulatory motif）分析。进行TF-motif富集分析，识别直接靶标。(仅保留具有正确的上游调节子且显著富集的motif modules，并对它们进行修剪以除去缺乏motif支持的间接靶标。）这些处理后的每个TF及其潜在的直接targets gene被称作一个调节子(regulon)，所以第二步就得到了转录因子(TF)与其对应的直接作用的靶点(targets genes),并且称为regulon(所以每一个regulon都是1个TF和这个TF调控的的targets genes)

3.使用AUCell算法对每个细胞的每个regulon活性进行打分
https://www.jianshu.com/p/6a6705f12842（AUCell包）
对于一个regulon来说，比较细胞间的AUCell得分可以鉴定出该regulon在哪种细胞显著更高的，这将确定regulon在哪些细胞中处于"打开"状态。简单来说第三步得到了每一个细胞对每一个regulon的得分，通过得分我们可以知道每一个细胞的每一个regulon的激活情况，即TF的激活情况

所以综上所得，我的浅显的认知：SCENIC能让我们知道每一个细胞他们不同的TF激活情况，不要跟RcisTarget包弄混了，虽然都是得到都是TF的差异，但是RcisTarget包是通过二者的差异基因来判断是什么TF造成的，而SCENIC是能够得出每一个细胞的TFs的激活情况。
ok 关于分析流程就将到这，当然这只是我个人浅显的理解，更标准，正派解释的还是要看原文献和各个生物信息学大佬的解释哈。（欢迎大家讨论与批评）

分析流程（R+linux）

我用过单纯R跑，超级慢，并且跑不出来，所以看别人的教程写着用linux+R一起
1.在R语言中获得单细胞的表达矩阵 scRNA是我自己的Seurat对象哈

write.csv(t(as.matrix(scRNA@assays$RNA@counts)), file = "scRNA.csv") #提取表达矩阵，并命名为scRNA.csv，注意矩阵一定要转置，不然会报错

2.从这里就要转战场去linux了，莫慌，学就是了
个人的linux学习流程：
1.https://www.bilibili.com/video/BV1pE411C7ho?p=42&spm_id_from=333.880.my_history.page.click
2.https://www.bilibili.com/video/BV1ds411g7eg?p=7&spm_id_from=333.880.my_history.page.click
建议先学第一个视频再学jimmy老师的，因为jimmy老师讲的都是干货，但是呢比较散，建议先把基础视频学一遍，在听jimmy老师的，会更加有印象。
当你学完了就可以run接下来的流程了，
从现在开始都在linux运行
首先你要在Linux安装conda，（这是第一道难关）conda是什么，怎么安装，它能干什么：
conda是什么：我的理解，conda就像手机上的应用超市，我们在手机上进入应用超市点击下载就把app下载了，假如没有应用超市的话，我们要安装一个app则需要下载安装包再解压之类的，因此有了应用超市，我们可以更便捷的安装app，那么conda就像服务器上的应用超市，让我们更方便的下载各种软件和包。
conda安装：

#下载安装包，解压，激活 wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh #下载安装包 bash Miniconda3-latest-Linux-x86_64.sh #使用bash命令解压安装，记得是一路yes下去 source ~/.bashrc #激活安装的conda #设置国内镜像 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

conda能干什么：1.conda能够更好的安装软件与r包 2.conda能创建一个环境
创建环境是什么意思呢，首先要搞清楚创建环境跟mkdir不一样，对于环境的理解就是，类似于我们看书需要一个安静的环境，所以我们可以conda一个环境，这个环境让我们更好做这次分析，因此对于SCENIC分析来说 ，由于r的分析过程很慢，而python的分析比较快，所以我们在linux要创建一个python的环境
conda创建一个python环境

conda create -n pyscenic python=3.8 #创建一个名为pyscenic的python3.8环境，这是一个难点，这里提醒一下，假如安装failed的话，可能是.condarc文件配置的问题哈。 conda info -e #查看全部的环境 conda activate pyscenic #激活我们工作的环境，进入到该环境 #在这个环境中配置一些依赖的东西 conda install pip conda install -y numpy conda install -y -c anaconda cytoolz conda install -y scanpy #注意scanpy这个包，假如安装不了用pip安装：pip install -i https://mirrors.bfsu.edu.cn/pypi/web/simple scanpy #假如还是安装不了，就去Google pip install pyscenic=0.11.1 #必须安装这个版本哈 因为在2022年5月份更新了，所以要安装以前的版本 否则后面运行会报错

linux上分析流程
1.在服务器上创建一个工作文件夹，并往文件夹中传输4个文件，进入该文件夹，开始工作

mkdir pyscience #创建名为pyscience的工作文件夹，这样我们就在一个名叫pyscenic 的环境下创建了一些pyscience的工作文件夹。 cd pyscience #进入该文件夹

2.使用文件传输软件（Xftp）往我们pyscience文件夹中传送4个文件

1. hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather
2. hs_hgnc_tfs.txt
3. motifs-v9-nr.hgnc-m0.001-o0.0.tbl
4. scRNA.csv
   解释这些文件
   123这三个文件都是该上述那三个流程需要用到的参比数据集，假如分析不是人的就不是这三个哈，并且必须确保这三个数据集是完好无损的
   4这个文件还记得吗，是第一步在r语言操作获得的单细胞表达矩阵
   假如大家需要123文件可以简信我哈（佛系回复）

3.run Linux流程
3.1 将scRNA.csv这个文件变成一个loom文件

pip install ipython #安装ipython，便于运行python代码 #如果慢的话加其他镜像 ipython #启动ipython 为了后面用python语言 #一行一行的输入以下代码， 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("scRNA.csv") row_attrs = {"Gene": np.array(x.var_names),} col_attrs = {"CellID": np.array(x.obs_names)} lp.create("sample.loom",x.X.transpose(),row_attrs,col_attrs) exit #退出ipython包

当我们run完这些代码后，在我们工作的文件夹下会多一个sample.loom的文件，这样就转换成功了（这里跟jimmy大佬的代码有一些不一样，大佬直接用python脚本，而我怕出错或者有一些依赖没有安装，所以用ipython包一行一行的run）
ipython运行示例图

记得一定要exit退出ipython包

3.2 pyscenic 的3个步骤之 GRN(linux里是用这个算法，R中是用GENIE3） 复制全部粘贴后运行即可 #这一步千万不要改num_workers 设置20就行

pyscenic grn \ --num_workers 20 \ --output adj.sample.tsv \ --method grnboost2 \ sample.loom \ hs_hgnc_tfs.txt #转录因子文件，1839 个基因的名字列表

3.3 pyscenic 的3个步骤之 cistarget 复制全部粘贴后运行即可 #这一步改num_workers 设置99也行

pyscenic ctx \ adj.sample.tsv \ hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather \ --annotations_fname motifs-v9-nr.hgnc-m0.001-o0.0.tbl \ --expression_mtx_fname sample.loom \ --mode "dask_multiprocessing" \ --output reg.csv \ --num_workers 3 \ --mask_dropouts

3.4 pyscenic 的3个步骤之 AUCell 复制全部粘贴后运行即可 #这一步改num_workers 设置99也行

pyscenic aucell \ sample.loom \ reg.csv \ --output sample_SCENIC.loom \ --num_workers 3

三个步骤run完后，用ls查看一下发现，多了一个sample_SCENIC.loom文件，这就是我们SCENIC分析的结果，将其下载到自己电脑并导入R中，进行plot了，by the way， GRN和cistarget这两个步骤很慢，建议看文献度过哈哈哈哈。

4.R中plot

#在这里我们要理解SCENIC做了什么分析 #对每个细胞进行转录因子富集分析 筛选出较为真实的转录因子 并以细胞为单位 即每一个细胞对每一个TF进行打分（AUG) #所以我们可以按照细胞类型， 族，样本来源，将相应的TFplot出来 #sample_SCENIC.loom导入R里面进行可视化 library(SCENIC) packageVersion("SCENIC") library(Seurat) library(SCopeLoomR) library(AUCell) library(SCENIC) library(dplyr) library(KernSmooth) library(RColorBrewer) library(plotly) library(BiocParallel) library(grid) library(ComplexHeatmap) library(data.table) library(scRNAseq) rm(list=ls()) # 提取sample_SCENIC.loom 信息 loom <- open_loom('sample_SCENIC.loom') #导入sample_SCENIC.loom文件 #首先我们要把导入的loom处理成R中的数据 #获取regulon regulon定义：TF与作用的genes #1.提取每一个TF与每一个gene作用系数 regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons") regulons_incidMat[1:4,1:4] #在这里就可以看出 每一个TF与每一个gene的作用数值 regulons <- regulonsToGeneLists(regulons_incidMat) #做成一个list TF与其作用gene的list TF+genes 个人感觉这里假如后面想分析这个TF，则这里可以画这个TF与其作用的gene的网络图 #2.获得regulon的AUC 即TF在每一个细胞的激活程度 regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC') regulonAUC[1:4,1:4] #regulonAUC这个文件含有每一个TF在各个细胞中的表达量 列名为细胞名 行名为TF #3.找出在这单细胞数据中 高表达的TF regulonAucThresholds <- get_regulon_thresholds(loom) tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))]) #这里可以看出哪一些TF是在这个单细胞数据中高表达的 #4.这两步不知道是啥 embeddings <- get_embeddings(loom) #好像是什么嵌入 必须要做 否则后面会报错 close_loom(loom) #这样就处理完成了 #接下来我们就可以挑选自己感兴趣的TF，进行个性化分析 #流程：1.查看富集到的全部的TF 2.从这些TF中挑选我们自己感兴趣的或者与课题相关的 #1.查看富集到的全部的TF #需要每一个细胞的每一个TF的表达量 #导入单细胞数据 scRNA <- readRDS("scRNA.rds") #这个rds是我的seruat对象 #将每一个细胞的每一个TF的表达量与单细胞的每一个细胞匹配 sub_regulonAUC <- regulonAUC[,match(colnames(scRNA),colnames(regulonAUC))] dim(sub_regulonAUC) #确认是否一致 identical(colnames(sub_regulonAUC), colnames(scRNA)) #[1] TRUE 一定要为TRUE，思考一下为啥，假如你用过r跑的话，r是要求你除了表达矩阵还要细胞的meta信息，而我们用linux跑的话，只用了细胞的表达矩阵，所以我们要判断我们我们做的分析是这个单细胞数据的分析，并且必须一模一样，否则你后面将分析的结果添加到seruat对象里就是不match的 #2.从这些TF中挑选我们自己感兴趣的或者与课题相关的 #挑选感兴趣的TF（所以要符合 1.在这个数据中有表达 2.不在所有细胞中都高表达，否则无法做差异分析） names(regulons) #我们可以通过这个函数来看在这个单细胞数据中 有哪些TF表达 在其中挑选感兴趣的TF #设置感兴趣的TF regulonsToPlot = c("TWIST1(+)","NKX2-1(+)","ZNF831(+)","SIX4(+)","FOSB(+)","TBX21(+)") #将感兴趣的TF的表达量加入到seruat对象中 scRNA@meta.data = cbind(scRNA@meta.data,t(assay(sub_regulonAUC[regulonsToPlot,]))) #3.可视化 #设置画图的分组 用idents 也可以用其他的分组 Idents(scRNA) <- sce$celltype table(Idents(scRNA)) DotPlot(scRNA, features = unique(regulonsToPlot)) + RotatedAxis() RidgePlot(scRNA, features = regulonsToPlot , ncol = 1) VlnPlot(scRNA, features = regulonsToPlot,pt.size = 0 ) FeaturePlot(scRNA, features = regulonsToPlot )

更多的plot形式大家去Google哦，终于长舒一口气，这个包学了大概2个星期.........
写在最后，大家还是多去学习jimmy大佬的教程，真的很有帮助，解决问题时面红耳赤，解决完后心情舒坦，这该死的成就感！！！

References:

https://cn.bing.com/search?q=pyscenic&form=QBLHCN&sp=-1&pq=pyscenic&sc=8-8&qs=n&sk=&cvid=875AE9ECDA3D4EE5BEE138BA629C438E
https://www.jianshu.com/p/0a29ecfaf21e
https://blog.csdn.net/qq_45688354/article/details/108014189