今天寫(xiě)一篇關(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》

分析原理

對(duì)于這個(gè)包我想發(fā)表一下感想，這個(gè)包剛開(kāi)始對(duì)于我來(lái)說(shuō)比較難理解，其實(shí)跟著官方教程running無(wú)非就是debug而已，這沒(méi)啥，重要的是這個(gè)分析做了些啥，最后得到啥，我能講啥，這是最難的，反正我是不敢做那些我自己都不了解過(guò)程的分析，所以首先我們要了解這個(gè)包的分析流程：
GENIE3或GRNBoost→RcisTarget→AUCell（所以首先你得先了解這三個(gè)包做了啥，不過(guò)不太想了解也可以往后看，我會(huì)講講我的理解，希望對(duì)你們有所幫助）

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

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

3.使用AUCell算法對(duì)每個(gè)細(xì)胞的每個(gè)regulon活性進(jìn)行打分
http://www.itdecent.cn/p/6a6705f12842（AUCell包）
對(duì)于一個(gè)regulon來(lái)說(shuō)，比較細(xì)胞間的AUCell得分可以鑒定出該regulon在哪種細(xì)胞顯著更高的，這將確定regulon在哪些細(xì)胞中處于"打開(kāi)"狀態(tài)。簡(jiǎn)單來(lái)說(shuō)第三步得到了每一個(gè)細(xì)胞對(duì)每一個(gè)regulon的得分，通過(guò)得分我們可以知道每一個(gè)細(xì)胞的每一個(gè)regulon的激活情況，即TF的激活情況

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

分析流程（R+linux）

我用過(guò)單純R跑，超級(jí)慢，并且跑不出來(lái)，所以看別人的教程寫(xiě)著用linux+R一起
1.在R語(yǔ)言中獲得單細(xì)胞的表達(dá)矩陣 scRNA是我自己的Seurat對(duì)象哈

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

2.從這里就要轉(zhuǎn)戰(zhàn)場(chǎng)去linux了，莫慌，學(xué)就是了
個(gè)人的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é)第一個(gè)視頻再學(xué)jimmy老師的，因?yàn)閖immy老師講的都是干貨，但是呢比較散，建議先把基礎(chǔ)視頻學(xué)一遍，在聽(tīng)jimmy老師的，會(huì)更加有印象。
當(dāng)你學(xué)完了就可以run接下來(lái)的流程了，
從現(xiàn)在開(kāi)始都在linux運(yùn)行
首先你要在Linux安裝conda，（這是第一道難關(guān)）conda是什么，怎么安裝，它能干什么：
conda是什么：我的理解，conda就像手機(jī)上的應(yīng)用超市，我們?cè)谑謾C(jī)上進(jìn)入應(yīng)用超市點(diǎn)擊下載就把a(bǔ)pp下載了，假如沒(méi)有應(yīng)用超市的話，我們要安裝一個(gè)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è)置國(guó)內(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)建一個(gè)環(huán)境
創(chuàng)建環(huán)境是什么意思呢，首先要搞清楚創(chuàng)建環(huán)境跟mkdir不一樣，對(duì)于環(huán)境的理解就是，類似于我們看書(shū)需要一個(gè)安靜的環(huán)境，所以我們可以conda一個(gè)環(huán)境，這個(gè)環(huán)境讓我們更好做這次分析，因此對(duì)于SCENIC分析來(lái)說(shuō) ，由于r的分析過(guò)程很慢，而python的分析比較快，所以我們?cè)趌inux要?jiǎng)?chuàng)建一個(gè)python的環(huán)境
conda創(chuàng)建一個(gè)python環(huán)境

conda create -n pyscenic python=3.8  #創(chuàng)建一個(gè)名為pyscenic的python3.8環(huán)境，這是一個(gè)難點(diǎn)，這里提醒一下，假如安裝failed的話，可能是.condarc文件配置的問(wèn)題哈。
conda info -e    #查看全部的環(huán)境
conda activate pyscenic   #激活我們工作的環(huán)境，進(jìn)入到該環(huán)境
#在這個(gè)環(huán)境中配置一些依賴的東西
conda install pip
conda install -y numpy     
conda install -y -c anaconda cytoolz
conda install -y scanpy   #注意scanpy這個(gè)包，假如安裝不了用pip安裝：pip install -i https://mirrors.bfsu.edu.cn/pypi/web/simple scanpy 
#假如還是安裝不了，就去Google
pip install pyscenic=0.11.1   #必須安裝這個(gè)版本哈  因?yàn)樵?022年5月份更新了，所以要安裝以前的版本 否則后面運(yùn)行會(huì)報(bào)錯(cuò)

linux上分析流程
1.在服務(wù)器上創(chuàng)建一個(gè)工作文件夾，并往文件夾中傳輸4個(gè)文件，進(jìn)入該文件夾，開(kāi)始工作

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

2.使用文件傳輸軟件（Xftp）往我們pyscience文件夾中傳送4個(gè)文件

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

3.run Linux流程
3.1 將scRNA.csv這個(gè)文件變成一個(gè)loom文件

pip install ipython  #安裝ipython，便于運(yùn)行python代碼   #如果慢的話加其他鏡像
ipython   #啟動(dòng)ipython 為了后面用python語(yǔ)言 
#一行一行的輸入以下代碼，
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)我們r(jià)un完這些代碼后，在我們工作的文件夾下會(huì)多一個(gè)sample.loom的文件，這樣就轉(zhuǎn)換成功了（這里跟jimmy大佬的代碼有一些不一樣，大佬直接用python腳本，而我怕出錯(cuò)或者有一些依賴沒(méi)有安裝，所以用ipython包一行一行的run）
ipython運(yùn)行示例圖

36872a6c3c3506f9cf0661770d55cc7.png

記得一定要exit退出ipython包

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

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

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

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

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

4.R中plot

#在這里我們要理解SCENIC做了什么分析
#對(duì)每個(gè)細(xì)胞進(jìn)行轉(zhuǎn)錄因子富集分析  篩選出較為真實(shí)的轉(zhuǎn)錄因子  并以細(xì)胞為單位  即每一個(gè)細(xì)胞對(duì)每一個(gè)TF進(jìn)行打分（AUG)  
#所以我們可以按照細(xì)胞類型， 族，樣本來(lái)源，將相應(yīng)的TFplot出來(lái)
#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.提取每一個(gè)TF與每一個(gè)gene作用系數(shù)
regulons_incidMat <- get_regulons(loom, column.attr.name="Regulons")   
regulons_incidMat[1:4,1:4]    #在這里就可以看出 每一個(gè)TF與每一個(gè)gene的作用數(shù)值
regulons <- regulonsToGeneLists(regulons_incidMat)      #做成一個(gè)list  TF與其作用gene的list  TF+genes  個(gè)人感覺(jué)這里假如后面想分析這個(gè)TF，則這里可以畫(huà)這個(gè)TF與其作用的gene的網(wǎng)絡(luò)圖                             
#2.獲得regulon的AUC  即TF在每一個(gè)細(xì)胞的激活程度
regulonAUC <- get_regulons_AUC(loom,column.attr.name='RegulonsAUC')
regulonAUC[1:4,1:4]  #regulonAUC這個(gè)文件含有每一個(gè)TF在各個(gè)細(xì)胞中的表達(dá)量  列名為細(xì)胞名   行名為T(mén)F
#3.找出在這單細(xì)胞數(shù)據(jù)中 高表達(dá)的TF
regulonAucThresholds <- get_regulon_thresholds(loom)     
tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))])   #這里可以看出哪一些TF是在這個(gè)單細(xì)胞數(shù)據(jù)中高表達(dá)的
#4.這兩步不知道是啥     
embeddings <- get_embeddings(loom)  #好像是什么嵌入   必須要做 否則后面會(huì)報(bào)錯(cuò)
close_loom(loom)
#這樣就處理完成了

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


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


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


#3.可視化
#設(shè)置畫(huà)圖的分組   用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哦，終于長(zhǎng)舒一口氣，這個(gè)包學(xué)了大概2個(gè)星期.........
寫(xiě)在最后，大家還是多去學(xué)習(xí)jimmy大佬的教程，真的很有幫助，解決問(wèn)題時(shí)面紅耳赤，解決完后心情舒坦，這該死的成就感?。?！

References:

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