|
1. 簡介
SCENIC（single-cell regulatory network inference and clustering）是一個基于共表達和motif分析，計算單細胞轉錄組數(shù)據(jù)基因調控網絡重建以及細胞狀態(tài)鑒定的方法季率。
2017年發(fā)表在Nature Methods雜志上的SCENIC算法仙蚜，利用單細胞RNA-seq數(shù)據(jù)，同時進行基因調控網絡重建和細胞狀態(tài)鑒定，應用于腫瘤和小鼠大腦單細胞圖譜數(shù)據(jù)患亿，提出并證明了順式調控網絡分析能夠用于指導轉錄因子和細胞狀態(tài)的鑒定慌核。
SCENIC通過使用生物學驅動的features自動清除腫瘤樣本特異性等批次效應。
有一些文章寫的挺好的央拖，在這里匯總一下：
https://www.cnblogs.com/raisok/p/12425225.html
http://www.reibang.com/p/cd967c449177
https://cloud.tencent.com/developer/article/1692240
2. 原理
GRN（gene regulatory network）基因調控網絡包括TF（transcription factor轉錄因子）祭阀、cofactor（共調因子）與其調節(jié)的target gene 組成，決定了某個狀態(tài)下的細胞的轉錄狀態(tài)鲜戒。SCENIC流程包括三步驟：
（1）使用GENIE3或GRNBoost (Gradient Boosting) 基于共表達推斷轉錄因子與候選靶基因之間的共表達模塊专控。
（2）由于GENIE3模型只是基于共表達，會存在很多假陽性和間接靶標遏餐，為了識別直接結合靶標（direct-binding targets）伦腐，使用RcisTarget對每個共表達模塊進行順式調控基序（cis-regulatory motif）分析。進行TF-motif富集分析失都，識別直接靶標柏蘑。（僅保留具有正確的上游調節(jié)子且顯著富集的motif modules，并對它們進行修剪以除去缺乏motif支持的間接靶標粹庞。）這些處理后的每個TF及其潛在的直接targets gene被稱作一個調節(jié)子（regulon）咳焚；
（3）使用AUCell算法對每個細胞的每個regulon活性進行打分。對于一個regulon來說信粮，比較細胞間的AUCell 得分可以鑒定出哪種細胞有顯著更高的subnetwork活性黔攒。結果生成一個二進制的regulon活性矩陣（binarized activity matrix），這將確定Regulon在哪些細胞中處于“打開”狀態(tài)强缘。
scenic分析原理
補充：
蛋白質中功能的基本單元是domain,是一種特殊的三維結構督惰，不同結構的domain與其他分子特異性結合從而發(fā)揮功能。與此類似旅掂，轉錄因子在于DNA序列結合時赏胚，其結合位點的序列也由于一定的特異性，不同轉錄因子結合的DNA序列的模式是不同的商虐。為了更好的描述結合位點序列的模式觉阅，科學家們提出了motif的概念崖疤。
2.1 GENIE3
GENIE3是一種從基因表達量數(shù)據(jù)推斷基因調控網絡的方法。它訓練預測數(shù)據(jù)集中每個基因表達的隨機森林模型典勇，并使用TF的表達作為輸入劫哼。然后使用不同的模型來得出TF的權重，并測量它們各自的相關性以預測每個靶基因的表達割笙。
GENIE3的輸入為表達矩陣权烧，最好使用gene-summarized counts（可能是也可能不是UMIs）。其他單位伤溉，比如counts般码，TPM和FPKM/RPKM也可以。但是要注意第一步的網絡相關分析基于共表達乱顾，一些作者建議也可以使用within-sample normalization比如TPM板祝。
GENIE3的輸出是一個帶有基因、潛在調節(jié)因子以及IM（importance measure）的表走净。IM代表了TF（input gene）在預測靶標時的權重券时。作者探索了幾種確定閾值的方法，最終選擇為每個TF建立多個潛在靶標基因集：（1）設置幾個IM閾值（IM>0.001 and IM >0.005)
（2）選取每個TF的前50哥靶標targets
（3）每個target gene保留top5温技，10革为，50個TFs，然后按TF分開舵鳞。
在以上結果中震檩，只有IM>0.001的links被算入。
每個基因集接著被分為positive- and negetive-correlated targets來區(qū)分可能激活的和抑制的targets蜓堕。（TF和潛在靶標的Spearman相關性計算）
最終抛虏，只有包含30個基因以上的基因集（TF共表達模型）被保留，用于下游分析套才。
2.2 RcisTarget
RcisTarget是i-cisTarget和iRegulon的motif富集框架的新R / Bioconductor實現(xiàn)迂猴。
RcisTarget從一個基因列表識別富集的TF-binding motifs和候選轉錄因子。主要有兩步驟：
（1）選擇在基因集中基因TSS（transcription start site）附近顯著過表達的DNA motif背伴。這一步通過在數(shù)據(jù)庫中應用recovery-based method（基于恢復的方法）來實現(xiàn)的沸毁，該數(shù)據(jù)庫包含每個motif的全基因組跨物種排名。保留注釋到對應TF并且NES（normalized enrichment score）>3的motif傻寂。
（2）對于每一個motif和基因集息尺，RcisTarget預測候選靶標基因，也就是基因集中排名領先的基因疾掰。這一步提供的結果跟i-cisTarget和iiRegulon相同搂誉。
為了構建最終的regulon，作者合并了每個TF module中預測的靶基因静檬，這些基因顯示了給定TF的任何motif的富集炭懊。但是在作者分析的數(shù)據(jù)中并级，這些modules數(shù)量很少而且motif富集很低。因此侮腹，最終決定從流程中去除對于直接表達的檢測嘲碧，只使用positive-correlated targets進行下游分析。
2.3 AUCell
對于一個給定的regulon凯旋，通過比較所有細胞間的AUCell（area under the recovery curve）打分值呀潭，我們可以識別哪些細胞具有更顯著高的regulon活性钉迷。
輸入為一個基因集至非，輸出為基因集每個細胞的'activity’。這些基因集即regulon糠聪，包含TFs和他們假定的的target荒椭。基于recovery analysis將根據(jù)表達水平將所有基因進行排序舰蟆。AUC代表了與細胞內其他基因相比趣惠，特征基因中表達基因的比例及其相對表達值。AUCell使用AUC來計算輸入基因集的關鍵子集是否在每個細胞的排名頂部都得到了富集身害。將輸出一個每個基因集在每個細胞的AUC score矩陣味悄。
通過卡閾值得到的二元活性矩陣使矩陣維數(shù)減少（可理解為只有 0|1，on|off）塌鸯，對于下游分析很有用侍瑟。 例如，基于regulon二元活性矩陣的聚類丙猬，可以根據(jù)某個調控子網絡（regulon）的活性來識別細胞群類型和細胞狀態(tài)涨颜。由于regulon是整體評分的，而不是使用單個基因的表達茧球，因此這種方法對于個別基因的dropouts很有效庭瑰。
3. 分析方法-R版
3.1 輸入
單細胞RNA-seq表達矩陣
每列對應于樣品（細胞），每行對應一個基因抢埋〉穑基因ID應該是gene-symbol并存儲為rownames （尤其是基因名字部分是為了與RcisTarget數(shù)據(jù)庫兼容）；
表達數(shù)據(jù)是Gene的reads count揪垄。根據(jù)作者的測試饲趋，提供原始的或Normalized UMI count絮重，無論是否log轉換，或使用TPM值，結果相差不大碘裕，但是不推薦TPM。
在進行GENIE3分析前要對數(shù)據(jù)進行過濾：
（1）過濾基因：對每個基因的總reads數(shù)進行過濾，去除最可能不可信的只提供噪音的基因。具體的值取決于數(shù)據(jù)辽社，文章中用到3 UMI counts × 30 (1% of cells) = minimum 90 counts per gene
（2）過濾基因：在多少細胞中被檢測到。去除只在一個活少量細胞中表達的基因翘鸭，如果這些基因正好在一個細胞中集合滴铅，將獲得很大的權重。推薦1%就乓。
3.2 分析步驟
官方github教程
3.2.1 加載R包
library(foreach)
library(SCENIC)
library(Seurat)
library(GENIE3)
library(AUCell)
library(RcisTarget)
library(loomR)
3.2.2 下載數(shù)據(jù)庫
比如說創(chuàng)建一個cisTarget_databases目錄汉匙，然后根據(jù)對應物種下載相應的數(shù)據(jù)庫文件到目錄下。下面示例為用R下載生蚁，也可以wget等其他方法噩翠。
準備工作：下載評分數(shù)據(jù)庫
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
download.file(dbFiles[1],'SCENIC/cisTarget_databases/mm9-500bp-upstream-7species.mc9nr.feather')
download.file(dbFiles[2],'SCENIC/cisTarget_databases/mm9-tss-centered-10kb-7species.mc9nr.feather')
3.2.3 導入數(shù)據(jù)
根據(jù)不同的數(shù)據(jù)類型有不同的導入方式，官方教程都有寫邦投，這里針對seurat對象導入伤锚。
導入單細胞數(shù)據(jù)
準備表達量矩陣
exprMat <- pt@assays$RNA@counts
exprMat <- as.matrix(exprMat)
準備細胞meta信息
cellInfo <- data.frame(pt@meta.data)
colnames(cellInfo)[which(colnames(cellInfo)=='orig.ident')] <- 'sample'
colnames(cellInfo)[which(colnames(cellInfo)=='seurat_clusters')] <- 'cluster'
colnames(cellInfo)[which(colnames(cellInfo)=='sub_type')] <- 'celltype'
cellInfo <- cellInfo[,c('sample','cluster','celltype','Group')]
saveRDS(cellInfo, file='int/cellInfo.Rds')
Initialize settings 初始設置，導入評分數(shù)據(jù)庫
scenicOptions <- initializeScenic(org='mgi',
dbDir='cisTarget_databases', nCores=10)scenicOptions@inputDatasetInfo$cellInfo <- 'int/cellInfo.Rds'
saveRDS(scenicOptions, file='int/scenicOptions.Rds')
scenicOptions <- initializeScenic()這一步的參數(shù)解釋：
dbDir 為存放數(shù)據(jù)庫的目錄
nCores 為分析使用的線程志衣，要考慮服務器配置情況設置
3.2.4 共表達網絡計算
使用GENIE3計算轉錄因子與每個基因的相關性屯援，消耗資源非常大。我進行測試時兩萬五細胞三萬基因的小鼠單細胞轉錄組數(shù)據(jù)跑了7天念脯，到現(xiàn)在還沒結束狞洋。處理大數(shù)據(jù)集的時候會非常非常慢，作者提出了兩種解決辦法：
（1）GENIE推斷共表達module時抽取少量細胞計算绿店，在計算regulon活性（AUCell）時吉懊，使用所有細胞計算；
（2）采用python版本的SCENIC惯吕，具體步驟下面會說到惕它。
我還是使用了最基礎的，而不是作者建議的方法废登。
==轉錄調控網絡推斷==##
基因過濾淹魄，minCountsPerGene和minSamples采用默認參數(shù)
過濾標準是基因表達量之和>細胞數(shù)*3%，且在1%的細胞中表達
genesKept <-geneFiltering(exprMat, scenicOptions, minCountsPerGene = 3 * 0.01 *ncol(exprMat), minSamples =ncol(exprMat) * 0.01)
exprMat_filtered <- exprMat[genesKept,]
計算相關性矩陣
runCorrelation(exprMat_filtered, scenicOptions)
TF-Targets相關性回歸分析
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions) ## 這一步跑了七天
runGenie3會在int目錄產生不少結果文件：
· 1.2_corrMat.Rds 基因間相關性矩陣堡距；
· 1.3_GENIE3_weightMatrix_part_1.Rds GENIE3中間結果甲锡，矩陣格式，有好多個羽戒，比如我這次分析缤沦，產生了10個；
· 1.4_GENIE3_linkList.Rds GENIE3的最終結果易稠，是把1.3那一堆文件合并在一起了缸废。文件有三列：
TF為轉錄因子名稱，
Target為潛在靶基因的名稱，
weight時TF與Target之間的相關性權重企量。
3.2.5 Step1 - 推斷共表達模塊
GENIE3計算了轉錄因子與每個基因的相關性测萎，接下來需要過濾掉低相關性的組合，形成共表達基因集届巩。作者嘗試了多種過濾標準硅瞧，沒有發(fā)現(xiàn)一種最佳策略，所以推薦6種都用恕汇，這六種分別為：
· w001：以每個TF為核心保留weight>0.001的基因形成共表達模塊腕唧；
· w005：以每個TF為核心保留weight>0.005的基因形成共表達模塊；
· top50：以每個TF為核心保留weight值top50的基因形成共表達模塊瘾英；
· top5perTarget：每個基因保留weight值top5的TF得到精簡的TF-Target關聯(lián)表枣接，然后把基因分配給TF構建共表達模塊；
· top10perTarget：每個基因保留weight值top10的TF得到精簡的TF-Target關聯(lián)表方咆，然后把基因分配給TF構建共表達模塊月腋；
· top50perTarget：每個基因保留weight值top50的TF得到精簡的TF-Target關聯(lián)表，然后把基因分配給TF構建共表達模塊瓣赂；
推斷共表達模塊
runSCENIC_1_coexNetwork2modules(scenicOptions)
結果文件為1.6_tfModules_asDF.Rds，四列：
method為6種過濾方法
corr為runCorrelation(exprMat_filtered, scenicOptions)得到的結果片拍，1代表激活煌集，-1代表抑制，0代表中性捌省。SCENIC只會采用值為1苫纤，即激活的數(shù)據(jù)用于后續(xù)分析。
3.2.6 Step2 - Motif驗證共表達模塊
上一步找到了每個轉錄因子強相關的靶基因纲缓，這一步要修建共表達模塊形成有生物學意義的調控單元（regulons）卷拘。
對每個共表達模塊進行motif富集分析，保留顯著富集的motif祝高；這一步依賴gene-motif評分數(shù)據(jù)庫栗弟，行為基因，列為motif工闺，值為排名乍赫，也就是我們下載的cisTarget數(shù)據(jù)庫。
使用數(shù)據(jù)庫對motif進行TF注釋陆蟆，得到高可信度雷厂、低可信度。數(shù)據(jù)庫直接注釋和同源基因推斷的TF為高可信度叠殷，使用motif序列相似性注釋的TF為低可信度結果改鲫。
用保留的motif對共表達模塊里的基因進行打分（依據(jù)cisTarget數(shù)據(jù)庫），識別顯著高分的基因（理解為motif里這些基因的TSS很近）。
刪除共表達模塊中與motif評分不高的基因像棘，剩下的基因集被稱為調控單元（regulon）纫塌。
推斷轉錄調控單元regulon
runSCENIC_2_createRegulons(scenicOptions) # Toy run settings
可以通過coexMethod調整使用的過濾方法，默認6中方法都做
如果做測試可以只選擇幾種方法減少計算量
output/Step2_MotifEnrichment.tsv讲弄，有motifDb geneSet motif NES AUC highlightedTFs TFinDB TF_highConf TF_lowConf nEnrGenes rankAtMax enrichedGenes措左，這11列。
· geneSet: Name of the gene set
· motif: ID of the motif
· NES: Normalized enrichment score of the motif in the gene-set
· AUC: Area Under the Curve (used to calculate the NES)
· TFinDB: Indicates whether the highlightedTFs are included within the high confidence annotation (two asterisks) or low confidence annotation (one asterisk).
· TF_highConf: Transcription factors annotated to the motif according to 'motifAnnot_highConfCat’.
· TF_lowConf: Transcription factors annotated to the motif according to 'motifAnnot_lowConfCat’.
· enrichedGenes: Genes that are highly ranked for the given motif.
· nErnGenes: Number of genes highly ranked
· rankAtMax: Ranking at the maximum enrichment, used to determine the number of enriched genes.
output/Step2_regulonTargetsInfo.tsv
最終的regulon文件避除，將第一個文件的數(shù)據(jù)整合在了一起怎披，表頭含義如下：
· TF：轉錄因子名稱
· gene：TF靶基因名稱
· nMotif：靶基因在數(shù)據(jù)庫的motif數(shù)量
· bestMotif：最顯著富集的motif名稱
· NES：標準富集分數(shù)，分值越高越顯著
· highConfAnnot：是不是高可信注釋
· Genie3Weight：TF與靶基因的相關性權重
Step2_MotifEnrichment_preview.html
3.2.7 Step3 - Regulon活動評分與可視化
一個regulon=一個TF與其Target gene的基因集瓶摆，regulon的名稱有兩種寫法：
· TF名稱+extended+靶基因數(shù)目：轉錄因子與所有靶基因組成的基因調控網絡
· TF名稱+靶基因數(shù)目：轉錄因子與高可信靶基因（即highConfAnnot=TRUE的基因）組成的基因調控網絡
AUCell對每個regulon在各個細胞中的活性進行評分凉逛。評分的基礎是基因表達值，分數(shù)越高代表基因集的激活程度越高群井。這一步要用所有細胞做計算状飞。
==regulon活性評分與可視化==##
regulons計算AUC值并進行下游分析
exprMat_all <- as.matrix(scRNA@assays$RNA@counts)
exprMat_all <- log2(exprMat_all+1)
runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)
3.2.8 Step4 - 生成二進制的regulon活性矩陣
runSCENIC_4_aucell_binarize(scenicOptions)
3.3 結果可視化
SCENIC通過幾個step即可繪制一些默認的圖片
3.3.1 Regulon Activity heatmap
step3有所有regulon在細胞的AUCscore熱圖Step3_RegulonActivity_heatmap.pdf
3.3.2 Binary Regulon Activity Heatmap
step4則生成多個熱圖Step4_BinaryRegulonActivity_Heatmap*
3.3.3 ****Cell-type specific regulators (RSS)
****Cell-type specific****regulators (based on the ****Regulon Specificity Score (RSS)****proposed by Suo et al. for the Mouse Cell Atlas in 2018).Useful for big analysis with many cell types, to identify the cell-type specific regulons.
當細胞種類比較多時可以用RSS來識別細胞類型特異性regulons
regulonAUC <-loadInt(scenicOptions, 'aucell_regulonAUC')
cellInfo1 <- readRDS('cellInfo1.rds')
rss <- calcRSS(AUC=getAUC(regulonAUC), cellAnnotation=cellInfo1[colnames(regulonAUC), 'type'])
rssPlot <-plotRSS(rss) #大小 rss評分，顏色 Z-score
rssPlot$plot
得到熱點圖书斜，顏色深淺代表z-score值诬辈，點的大小代表RSS評分。
然后就可以挑選各個細胞類型代表性regulon繪制熱圖等等啦荐吉。
繪制過程中我發(fā)現(xiàn)有一個點就是焙糟，根據(jù)regulon list提取到norm的AUCscore之后，再scale（計算z-score）比直接計算所有regulons的z-score再提取样屠，繪制效果要好穿撮。數(shù)據(jù)好的話，可以得到下圖這樣的結果啦痪欲。
4. 分析方法-Python版
同樣不介紹安裝步驟悦穿，可查看官方說明。
如果服務器配置了docker的話直接用docker比較簡便业踢，拉一個鏡像就行了：
docker pull aertslab/pyscenic:0.10.0
這里主要記錄一下分析方法栗柒。
4.1 輸入文件準備
以下是官方教程用到的文件：
（1）表達量文件，cellinfo
導入單細胞數(shù)據(jù)
準備表達量矩陣
exprMat <- epi@assays$RNA@counts
exprMat <- t(as.matrix(exprMat))
write.csv(exprMat,'expr.csv')
（2）cisTarget數(shù)據(jù)下載
從數(shù)據(jù)庫下載陨亡，我做的是人的
（3）motif數(shù)據(jù)下載
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
（4）需要手動處理的tf名稱列表
需要自己處理一下數(shù)據(jù)庫文件來生成傍衡，作者也提供了處理方法；官方也提供了人和小鼠的所有TFs列表
BASEFOLDER_NAME = '../resources/'
MOTIFS_HGNC_FNAME = os.path.join(BASEFOLDER_NAME, 'motifs-v9-nr.hgnc-m0.001-o0.0.tbl')
df_motifs_hgnc = pd.read_csv(MOTIFS_HGNC_FNAME, sep='\t')
hs_tfs = df_motifs_hgnc.gene_name.unique()
with open(OUT_TFS_HGNC_FNAME, 'wt') as f:
f.write('\n'.join(hs_tfs) + '\n')len(hs_tfs)
4.2 分析步驟
嘗試了jupyter教程负蠕，很奇怪同樣的環(huán)境命令行都可以加載蛙埂，但是jupyter無法加載scenic下的一個包。
嘗試了命令行方式運行遮糖，也很奇怪有報錯绣的。
最后用docker運行成功，可能對環(huán)境的要求比較復雜吧。
4.2.1 grn
輸入文件為表達量文件屡江，這里要注意表達量矩陣expr.csv 為cell x gene芭概，即行名為基因，列名為細胞barcode惩嘉。
輸出文件為expr_mat.adjacencies.tsv罢洲。
分析主要耗費的時間在這一步，我測的兩萬多細胞跑了15h左右文黎。
docker run -it -d --rm -v/xxxxx/SCENICanalysis/:/data aertslab/pyscenic pyscenic grn --num_workers 30 -o /data/malignant.out/expr_mat.adjacencies.tsv /data/expr.csv /data/hg38/hs_hgnc_tfs.txt
4.2.2 ctx
這里用到了數(shù)據(jù)庫文件惹苗，expr.csv，和上一步的結果文件耸峭。
輸出文件為regulons.csv桩蓉。
docker run -it -d --rm -v/xxxxx/SCENICanalysis/:/data aertslab/pyscenic pyscenic ctx --num_workers 30/data/epi.out/expr_mat.adjacencies.tsv /data/hg38/hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather /data/hg38/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather --annotations_fname/data/hg38/motifs-v9-nr.hgnc-m0.001-o0.0.tbl --expression_mtx_fname/data/expr.csv --modedask_multiprocessing --output /data/epi.out/regulons.csv
4.2.3 aucell
輸入文件為上一步的結果文件regulons.csv，和表達量文件expr.csv劳闹。
輸出為auc打分結果auc_mtx.csv院究。
docker run -it -d --rm -v/xxxxx/SCENICanalysis/:/data aertslab/pyscenic pyscenic aucell /data/expr.csv /data/epi.out/regulons.csv -o/data/epi.out/auc_mtx.csv --num_workers 30
5. 結果展示方式
Construction of a human cell landscape at single-cell level
原文：SCENIC單細胞轉錄因子分析
