今天写一篇关于单细胞转录因子的分析-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运行示例图

36872a6c3c3506f9cf0661770d55cc7.png

记得一定要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