SCENIC官網(wǎng):https://rawcdn.githack.com/aertslab/SCENIC/701cc7cc4ac762b91479b3bd2eaf5ad5661dd8c2/inst/doc/SCENIC_Setup.html
1. SCENIC簡介
組織內(nèi)細(xì)胞異質(zhì)性的基礎(chǔ)是細(xì)胞轉(zhuǎn)錄狀態(tài)的差異增蹭,轉(zhuǎn)錄狀態(tài)的特異性又是由轉(zhuǎn)錄因子主導(dǎo)的基因調(diào)控網(wǎng)絡(luò)(GRNs)決定并維持穩(wěn)定的闹蒜。因此分析單細(xì)胞的GRNs有助于深入挖掘細(xì)胞異質(zhì)性背后的生物學(xué)意義缩功,并為疾病的診斷害驹、治療以及發(fā)育分化的研究提供有價值的線索帅霜。然而單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)具有背景噪音高薄声、基因檢出率低和表達(dá)矩陣稀疏性的特點(diǎn)荒勇,給傳統(tǒng)統(tǒng)計學(xué)和生物信息學(xué)方法推斷高質(zhì)量的GRNs帶來了挑戰(zhàn)要糊。Single-cell regulatory network inference and clustering (SCENIC)是一種專為單細(xì)胞數(shù)據(jù)開發(fā)的GRNs算法萎胰,它的創(chuàng)新之處在于引入了轉(zhuǎn)錄因子motif序列驗(yàn)證統(tǒng)計學(xué)方法推斷的基因共表達(dá)網(wǎng)絡(luò)碾盟,從而識別高可靠性的由轉(zhuǎn)錄因子主導(dǎo)的GRNs。
2. 輸入
輸入:SCENIC需要輸入的是單細(xì)胞RNA-seq表達(dá)矩陣—— 每列對應(yīng)于樣品(細(xì)胞)技竟,每行對應(yīng)一個基因冰肴。基因ID應(yīng)該是gene-symbol并存儲為rownames (尤其是基因名字部分是為了與RcisTarget數(shù)據(jù)庫兼容)榔组;表達(dá)數(shù)據(jù)是Gene的reads count熙尉。根據(jù)作者的測試,提供原始的或Normalized UMI count搓扯,無論是否log轉(zhuǎn)換检痰,或使用TPM值,結(jié)果相差不大锨推。(Overall, SCENIC is quite robust to this choice, we have applied SCENIC to datasets using raw (logged) UMI counts, normalized UMI counts, and TPM and they all provided reliable results (see Aibar et al. (2017)).)
另外铅歼,運(yùn)行單細(xì)胞轉(zhuǎn)錄因子分析之SCENIC流程還需要下載配套數(shù)據(jù)庫公壤,不同物種不一樣, 在 https://resources.aertslab.org/cistarget/ 查看自己的物種谭贪,按需下載:
# https://resources.aertslab.org/cistarget/
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-500bp-upstream-7species.mc9nr.feather",
"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg19/refseq_r45/mc9nr/gene_based/hg19-tss-centered-10kb-7species.mc9nr.feather")
# dir.create("cisTarget_databases"); setwd("cisTarget_databases") # if needed
dbFiles <- c("https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-500bp-upstream-7species.mc9nr.feather",
"https://resources.aertslab.org/cistarget/databases/mus_musculus/mm9/refseq_r45/mc9nr/gene_based/mm9-tss-centered-10kb-7species.mc9nr.feather")
# mc9nr: Motif collection version 9: 24k motifs
dbFiles
for(featherURL in dbFiles)
{
download.file(featherURL, destfile=basename(featherURL)) # saved in current dir
# (1041.7 MB)
#
}
3. SCENIC分析流程
SCENIC主要包含三個步驟:
-
GENIE3
(隨機(jī)森林)/GRNBoost
(Gradient Boosting):基于共表達(dá)情況鑒定每個TF的潛在靶點(diǎn)境钟,推斷轉(zhuǎn)錄因子與候選靶基因之間的共表達(dá)模塊。每個模塊包含一個轉(zhuǎn)錄因子及其靶基因俭识。(R語言版SCENIC這一步用的隨機(jī)森林的算法慨削,Python版用的是Gradient Boosting算法,兩種算法不太一樣套媚,但算出來的結(jié)果是類似的) -
RcisTarget
:由于GENIE3只是推斷共表達(dá)缚态,因此會有假陽性和間接targets。使用RcisTarget基于DNA-motif分析識別具有正確上游調(diào)控子且顯著富集的motif(轉(zhuǎn)錄因子直接結(jié)合的motif)堤瘤,修剪掉缺乏motif支持的間接靶標(biāo)玫芦。修剪后的每個TF和其潛在的直接作用的target genes被稱為為一個regulon。(這一步是SCENIC和其他大多數(shù)共表達(dá)算法的重要區(qū)別) -
AUCell
:分析每個細(xì)胞的regulons活性并進(jìn)行打分本辐,打分的基礎(chǔ)是基因的表達(dá)值桥帆。對于regulon來說,比較細(xì)胞間的AUCell得分可以鑒定出哪種細(xì)胞具有顯著更高的regulon活性慎皱。打分值可以進(jìn)一步轉(zhuǎn)化為二進(jìn)制regulon活性矩陣(binarized activity matrix)老虫,這將最大化細(xì)胞類型的差異,確定regulon在哪些細(xì)胞中處于“開放”狀態(tài)茫多。
3.1 上游分析
在分析前首先要對initializeScenic()進(jìn)行設(shè)置祈匙,后續(xù)的分析基于這個設(shè)定。
1. runCorrelation()
共表達(dá)分析的結(jié)果中既有正向調(diào)控也有負(fù)向調(diào)控天揖,GENIE3無法區(qū)分夺欲。因此需要相關(guān)性矩陣幫助篩選共表達(dá)模塊中和TF正相關(guān)的基因。
2. runGenie3()
參考數(shù)據(jù)庫找出輸入基因中哪些是TF今膊,計算每個TF和各個基因之間的相關(guān)性權(quán)重些阅。權(quán)重其實(shí)也就是TF對gene表達(dá)量的貢獻(xiàn)。
3. runSCENIC_1_coexNetwork2modules()
得到上面的TF和gene的權(quán)重矩陣以后斑唬,就可以生成以TF為核心的geneset扑眉。(比如TF1它可能對gene1、gene5赖钞、gene12...的預(yù)測都比較好,就可以得到以TF1為核心的geneset)聘裁。隨后使用以下6種方法過濾掉低相關(guān)性的TF-genes共表達(dá)基因集雪营,得到以TF為核心的共表達(dá)基因集。這一步運(yùn)行得到的結(jié)果中包含了一列corr衡便,是runCorrelation()得到的結(jié)果献起。1代表激活洋访,-1代表抑制,0代表中性谴餐,SCENIC只會采用corr值為1的數(shù)據(jù)用于后續(xù)分析姻政,以得到正調(diào)控共表達(dá)模塊。
6種方法 | 含義 |
---|---|
w001 | 以每個TF為核心保留weight>0.001的基因形成共表達(dá)模塊岂嗓; |
w005 | 以每個TF為核心保留weight>0.005的基因形成共表達(dá)模塊汁展; |
top50 | 以每個TF為核心保留weight值top50的基因形成共表達(dá)模塊; |
top5perTarget | 每個基因保留weight值top5的TF得到精簡的TF-Target關(guān)聯(lián)表厌殉,然后把基因分配給TF構(gòu)建共表達(dá)模塊食绿; |
top10perTarget | 每個基因保留weight值top10的TF得到精簡的TF-Target關(guān)聯(lián)表,然后把基因分配給TF構(gòu)建共表達(dá)模塊公罕; |
top50perTarget | 每個基因保留weight值top50的TF得到精簡的TF-Target關(guān)聯(lián)表器紧,然后把基因分配給TF構(gòu)建共表達(dá)模塊; |
4. runSCENIC_2_createRegulons()
通過RcisTarget數(shù)據(jù)庫對得到的共表達(dá)模塊進(jìn)行修剪楼眷。通俗的講就是在計算出的TF-gene對中铲汪,結(jié)合數(shù)據(jù)庫查看該gene上游序列是否存在該TF結(jié)合的motif。從而剔除TF非直接調(diào)控基因的共表達(dá)模塊罐柳,保留Motif分析共表達(dá)模塊內(nèi)與TF有直接調(diào)控關(guān)系的基因掌腰,得到Regulon。(RcisTarget數(shù)據(jù)庫中的數(shù)據(jù)目前只支持人硝清,鼠辅斟,果蠅三個物種。)
5. runSCENIC_3_scoreCells()
每個Regulon就是一個轉(zhuǎn)錄因子及其直接調(diào)控靶基因的基因集芦拿,SCENIC接下來的工作就是對每個regulon在各個細(xì)胞中的活性評分士飒,得到每個基因集在每個細(xì)胞的AUC score矩陣(AUC代表來與細(xì)胞內(nèi)其他基因相比,特征基因中代表基因的比例及其相對表達(dá)值)蔗崎。評分是基于recovery analysis酵幕,根據(jù)基因的表達(dá)值進(jìn)行,分?jǐn)?shù)越高代表基因集的激活程度越高缓苛。
6. runSCENIC_4_aucell_binarize()
對于細(xì)胞類型清晰的數(shù)據(jù)集而言芳撒,構(gòu)建regulons并對其活性打分足夠后續(xù)分析。但是未桥,在很多情況下將評分轉(zhuǎn)換為二進(jìn)制的“開/關(guān)”(on|off)笔刹,則既方便解釋,又能最大化體現(xiàn)細(xì)胞類型的差異冬耿。將特定的regulon轉(zhuǎn)換為“0/1”有利于探索和解釋關(guān)鍵轉(zhuǎn)錄因子舌菜。將所有的regulons轉(zhuǎn)換為“0/1”后創(chuàng)建二進(jìn)制的活性矩陣,則可以用于細(xì)胞聚類亦镶,由于regulon是基于整體評分的日月,對消除技術(shù)偏倚如個別基因的dropout特別有用袱瓮。AUCell會自動計算可能的閾值進(jìn)行“0/1”轉(zhuǎn)換,作者建議在轉(zhuǎn)換之前手工檢查并調(diào)整這些閾值爱咬。
??在內(nèi)存不足的情況下尺借,可以使用抽樣的方法,選取部分基因來得到Regulons精拟,在runSCENIC_3_scoreCells()和runSCENIC_4_aucell_binarize()這兩步時燎斩,再把表達(dá)矩陣替換成全部基因的表達(dá)矩陣來計算評分。
3.2 下游分析
三個方向:
1串前、降維聚類發(fā)現(xiàn)新亞群(cell type/state由轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的差異決定)
2瘫里、case-control之間的regulons差異分析
3、尋找cell type/state特異性的regulon/TF
延伸參考:
http://www.reibang.com/p/cd967c449177
https://cloud.tencent.com/developer/article/1692240