今天寫一篇關(guān)于單細(xì)胞轉(zhuǎn)錄因子的分析-scenic
在做分析之前你應(yīng)該看一下這篇文獻(xiàn)君纫，看看別人做了什么scenic分析
《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》

分析原理

對于這個包我想發(fā)表一下感想，這個包剛開始對于我來說比較難理解鲸阻，其實跟著官方教程running無非就是debug而已，這沒啥缝裤，重要的是這個分析做了些啥纵菌，最后得到啥，我能講啥溺忧，這是最難的咏连，反正我是不敢做那些我自己都不了解過程的分析，所以首先我們要了解這個包的分析流程：
GENIE3或GRNBoost→RcisTarget→AUCell（所以首先你得先了解這三個包做了啥鲁森，不過不太想了解也可以往后看祟滴，我會講講我的理解，希望對你們有所幫助）

1.使用GENIE3或GRNBoost(Gradient Boosting)基于共表達(dá)推斷轉(zhuǎn)錄因子與候選靶基因之間的共表達(dá)模塊
http://www.reibang.com/p/d71dcd4cff5a (可以學(xué)一下這個包歌溉，體驗一下垄懂，或許我理解的不對哈）
是不是看到這個就頭大 ，所有的教程都是這樣痛垛，其實簡單的說第一步推斷出轉(zhuǎn)錄因子(TF)與其作用的靶點/gene(一個轉(zhuǎn)錄因子可能調(diào)控很多的gene哈)

2.RcisTarget對每個共表達(dá)模塊進(jìn)行順式調(diào)控基序(cis-regulatory motif）分析
http://www.reibang.com/p/c9df97f9e6e0 (RcisTarget包）
由于GENIE3模型只是基于共表達(dá)草慧，會存在很多假陽性和間接靶標(biāo)，為了識別直接結(jié)合靶標(biāo)(direct-binding targets)匙头，使用RcisTarget對每個共表達(dá)模塊進(jìn)行順式調(diào)控基序(cis-regulatory motif）分析漫谷。進(jìn)行TF-motif富集分析，識別直接靶標(biāo)蹂析。(僅保留具有正確的上游調(diào)節(jié)子且顯著富集的motif modules舔示，并對它們進(jìn)行修剪以除去缺乏motif支持的間接靶標(biāo)。）這些處理后的每個TF及其潛在的直接targets gene被稱作一個調(diào)節(jié)子(regulon)电抚，所以第二步就得到了轉(zhuǎn)錄因子(TF)與其對應(yīng)的直接作用的靶點(targets genes),并且稱為regulon(所以每一個regulon都是1個TF和這個TF調(diào)控的的targets genes)

3.使用AUCell算法對每個細(xì)胞的每個regulon活性進(jìn)行打分
http://www.reibang.com/p/6a6705f12842（AUCell包）
對于一個regulon來說惕稻，比較細(xì)胞間的AUCell得分可以鑒定出該regulon在哪種細(xì)胞顯著更高的，這將確定regulon在哪些細(xì)胞中處于"打開"狀態(tài)蝙叛。簡單來說第三步得到了每一個細(xì)胞對每一個regulon的得分俺祠，通過得分我們可以知道每一個細(xì)胞的每一個regulon的激活情況，即TF的激活情況

所以綜上所得甥温，我的淺顯的認(rèn)知：SCENIC能讓我們知道每一個細(xì)胞他們不同的TF激活情況锻煌，不要跟RcisTarget包弄混了妓布，雖然都是得到都是TF的差異姻蚓，但是RcisTarget包是通過二者的差異基因來判斷是什么TF造成的，而SCENIC是能夠得出每一個細(xì)胞的TFs的激活情況匣沼。
ok 關(guān)于分析流程就將到這狰挡，當(dāng)然這只是我個人淺顯的理解，更標(biāo)準(zhǔn)释涛，正派解釋的還是要看原文獻(xiàn)和各個生物信息學(xué)大佬的解釋哈加叁。（歡迎大家討論與批評）

分析流程（R+linux）

我用過單純R跑，超級慢唇撬，并且跑不出來它匕，所以看別人的教程寫著用linux+R一起
1.在R語言中獲得單細(xì)胞的表達(dá)矩陣 scRNA是我自己的Seurat對象哈

write.csv(t(as.matrix(scRNA@assays$RNA@counts)), file = "scRNA.csv")
#提取表達(dá)矩陣，并命名為scRNA.csv窖认，注意矩陣一定要轉(zhuǎn)置豫柬，不然會報錯

2.從這里就要轉(zhuǎn)戰(zhàn)場去linux了告希，莫慌，學(xué)就是了
個人的linux學(xué)習(xí)流程：
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
建議先學(xué)第一個視頻再學(xué)jimmy老師的烧给，因為jimmy老師講的都是干貨燕偶，但是呢比較散，建議先把基礎(chǔ)視頻學(xué)一遍础嫡，在聽jimmy老師的指么，會更加有印象。
當(dāng)你學(xué)完了就可以run接下來的流程了榴鼎，
從現(xiàn)在開始都在linux運行
首先你要在Linux安裝conda伯诬，（這是第一道難關(guān)）conda是什么，怎么安裝檬贰，它能干什么：
conda是什么：我的理解姑廉，conda就像手機(jī)上的應(yīng)用超市，我們在手機(jī)上進(jìn)入應(yīng)用超市點擊下載就把a(bǔ)pp下載了翁涤，假如沒有應(yīng)用超市的話桥言，我們要安裝一個app則需要下載安裝包再解壓之類的，因此有了應(yīng)用超市葵礼，我們可以更便捷的安裝app号阿，那么conda就像服務(wù)器上的應(yīng)用超市，讓我們更方便的下載各種軟件和包鸳粉。
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

#設(shè)置國內(nèi)鏡像
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能創(chuàng)建一個環(huán)境
創(chuàng)建環(huán)境是什么意思呢挺智，首先要搞清楚創(chuàng)建環(huán)境跟mkdir不一樣，對于環(huán)境的理解就是艰山，類似于我們看書需要一個安靜的環(huán)境湖雹，所以我們可以conda一個環(huán)境，這個環(huán)境讓我們更好做這次分析曙搬，因此對于SCENIC分析來說 摔吏，由于r的分析過程很慢，而python的分析比較快纵装，所以我們在linux要創(chuàng)建一個python的環(huán)境
conda創(chuàng)建一個python環(huán)境

conda create -n pyscenic python=3.8  #創(chuàng)建一個名為pyscenic的python3.8環(huán)境征讲，這是一個難點，這里提醒一下橡娄，假如安裝failed的話诗箍，可能是.condarc文件配置的問題哈。
conda info -e    #查看全部的環(huán)境
conda activate pyscenic   #激活我們工作的環(huán)境挽唉，進(jìn)入到該環(huán)境
#在這個環(huán)境中配置一些依賴的東西
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.在服務(wù)器上創(chuàng)建一個工作文件夾氨距，并往文件夾中傳輸4個文件桑逝，進(jìn)入該文件夾，開始工作

mkdir pyscience  #創(chuàng)建名為pyscience的工作文件夾俏让，這樣我們就在一個名叫pyscenic 的環(huán)境下創(chuàng)建了一些pyscience的工作文件夾楞遏。
cd pyscience #進(jìn)入該文件夾

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這三個文件都是該上述那三個流程需要用到的參比數(shù)據(jù)集，假如分析不是人的就不是這三個哈首昔，并且必須確保這三個數(shù)據(jù)集是完好無損的
   4這個文件還記得嗎寡喝，是第一步在r語言操作獲得的單細(xì)胞表達(dá)矩陣
   假如大家需要123文件可以簡信我哈（佛系回復(fù)）

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包

當(dāng)我們run完這些代碼后预鬓，在我們工作的文件夾下會多一個sample.loom的文件，這樣就轉(zhuǎn)換成功了（這里跟jimmy大佬的代碼有一些不一樣赊颠，大佬直接用python腳本格二，而我怕出錯或者有一些依賴沒有安裝，所以用ipython包一行一行的run）
ipython運行示例圖

36872a6c3c3506f9cf0661770d55cc7.png

記得一定要exit退出ipython包

3.2 pyscenic 的3個步驟之 GRN(linux里是用這個算法竣蹦，R中是用GENIE3） 復(fù)制全部粘貼后運行即可 #這一步千萬不要改num_workers 設(shè)置20就行

pyscenic grn \
--num_workers 20 \
--output adj.sample.tsv \
--method grnboost2 \
sample.loom \
hs_hgnc_tfs.txt    #轉(zhuǎn)錄因子文件顶猜，1839  個基因的名字列表

3.3 pyscenic 的3個步驟之 cistarget 復(fù)制全部粘貼后運行即可 #這一步改num_workers 設(shè)置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 復(fù)制全部粘貼后運行即可 #這一步改num_workers 設(shè)置99也行

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

三個步驟run完后，用ls查看一下發(fā)現(xiàn)痘括，多了一個sample_SCENIC.loom文件长窄，這就是我們SCENIC分析的結(jié)果，將其下載到自己電腦并導(dǎo)入R中纲菌，進(jìn)行plot了挠日，by the way， GRN和cistarget這兩個步驟很慢翰舌，建議看文獻(xiàn)度過哈哈哈哈嚣潜。

4.R中plot

#在這里我們要理解SCENIC做了什么分析
#對每個細(xì)胞進(jìn)行轉(zhuǎn)錄因子富集分析  篩選出較為真實的轉(zhuǎn)錄因子  并以細(xì)胞為單位  即每一個細(xì)胞對每一個TF進(jìn)行打分（AUG)  
#所以我們可以按照細(xì)胞類型， 族灶芝，樣本來源郑原，將相應(yīng)的TFplot出來
#sample_SCENIC.loom導(dǎo)入R里面進(jìn)行可視化
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')                 #導(dǎo)入sample_SCENIC.loom文件

#首先我們要把導(dǎo)入的loom處理成R中的數(shù)據(jù)
#獲取regulon        regulon定義：TF與作用的genes
#1.提取每一個TF與每一個gene作用系數(shù)
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")   
regulons_incidMat[1:4,1:4]    #在這里就可以看出 每一個TF與每一個gene的作用數(shù)值
regulons <- regulonsToGeneLists(regulons_incidMat)      #做成一個list  TF與其作用gene的list  TF+genes  個人感覺這里假如后面想分析這個TF唉韭，則這里可以畫這個TF與其作用的gene的網(wǎng)絡(luò)圖                             
#2.獲得regulon的AUC  即TF在每一個細(xì)胞的激活程度
regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')
regulonAUC[1:4,1:4]  #regulonAUC這個文件含有每一個TF在各個細(xì)胞中的表達(dá)量  列名為細(xì)胞名   行名為TF
#3.找出在這單細(xì)胞數(shù)據(jù)中 高表達(dá)的TF
regulonAucThresholds <- get_regulon_thresholds(loom)     
tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))])   #這里可以看出哪一些TF是在這個單細(xì)胞數(shù)據(jù)中高表達(dá)的
#4.這兩步不知道是啥     
embeddings <- get_embeddings(loom)  #好像是什么嵌入   必須要做 否則后面會報錯
close_loom(loom)
#這樣就處理完成了

#接下來我們就可以挑選自己感興趣的TF夜涕，進(jìn)行個性化分析
#流程：1.查看富集到的全部的TF       2.從這些TF中挑選我們自己感興趣的或者與課題相關(guān)的


#1.查看富集到的全部的TF
#需要每一個細(xì)胞的每一個TF的表達(dá)量
#導(dǎo)入單細(xì)胞數(shù)據(jù)
scRNA <- readRDS("scRNA.rds")    #這個rds是我的seruat對象
#將每一個細(xì)胞的每一個TF的表達(dá)量與單細(xì)胞的每一個細(xì)胞匹配
sub_regulonAUC <- regulonAUC[,match(colnames(scRNA),colnames(regulonAUC))]
dim(sub_regulonAUC)
#確認(rèn)是否一致
identical(colnames(sub_regulonAUC), colnames(scRNA))
#[1] TRUE  一定要為TRUE，思考一下為啥属愤，假如你用過r跑的話女器，r是要求你除了表達(dá)矩陣還要細(xì)胞的meta信息，而我們用linux跑的話住诸，只用了細(xì)胞的表達(dá)矩陣驾胆，所以我們要判斷我們我們做的分析是這個單細(xì)胞數(shù)據(jù)的分析涣澡，并且必須一模一樣，否則你后面將分析的結(jié)果添加到seruat對象里就是不match的


#2.從這些TF中挑選我們自己感興趣的或者與課題相關(guān)的
#挑選感興趣的TF（所以要符合 1.在這個數(shù)據(jù)中有表達(dá) 2.不在所有細(xì)胞中都高表達(dá)丧诺，否則無法做差異分析）
names(regulons)        #我們可以通過這個函數(shù)來看在這個單細(xì)胞數(shù)據(jù)中 有哪些TF表達(dá)  在其中挑選感興趣的TF
#設(shè)置感興趣的TF   
regulonsToPlot = c("TWIST1(+)","NKX2-1(+)","ZNF831(+)","SIX4(+)","FOSB(+)","TBX21(+)")
#將感興趣的TF的表達(dá)量加入到seruat對象中
scRNA@meta.data = cbind(scRNA@meta.data,t(assay(sub_regulonAUC[regulonsToPlot,])))


#3.可視化
#設(shè)置畫圖的分組   用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哦入桂，終于長舒一口氣，這個包學(xué)了大概2個星期.........
寫在最后驳阎，大家還是多去學(xué)習(xí)jimmy大佬的教程抗愁，真的很有幫助，解決問題時面紅耳赤呵晚，解決完后心情舒坦蜘腌，這該死的成就感！６丁撮珠！

References:

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