1. 簡介

SCENIC（single-cell regulatory network inference and clustering）是一個基于共表達(dá)和motif分析，計算單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)基因調(diào)控網(wǎng)絡(luò)重建以及細(xì)胞狀態(tài)鑒定的方法幔荒。

2017年發(fā)表在Nature Methods雜志上的SCENIC算法爹梁，利用單細(xì)胞RNA-seq數(shù)據(jù)，同時進(jìn)行基因調(diào)控網(wǎng)絡(luò)重建和細(xì)胞狀態(tài)鑒定念链，應(yīng)用于腫瘤和小鼠大腦單細(xì)胞圖譜數(shù)據(jù)掂墓，提出并證明了順式調(diào)控網(wǎng)絡(luò)分析能夠用于指導(dǎo)轉(zhuǎn)錄因子和細(xì)胞狀態(tài)的鑒定君编。

SCENIC通過使用生物學(xué)驅(qū)動的features自動清除腫瘤樣本特異性等批次效應(yīng)吃嘿。

有一些文章寫的挺好的兑燥，在這里匯總一下：

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）基因調(diào)控網(wǎng)絡(luò)包括TF（transcription factor轉(zhuǎn)錄因子）降瞳、cofactor（共調(diào)因子）與其調(diào)節(jié)的target gene 組成力崇，決定了某個狀態(tài)下的細(xì)胞的轉(zhuǎn)錄狀態(tài)亮靴。SCENIC流程包括三步驟：

（1）使用GENIE3或GRNBoost (Gradient Boosting) 基于共表達(dá)推斷轉(zhuǎn)錄因子與候選靶基因之間的共表達(dá)模塊于置。

（2）由于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）隆判；

（3）使用AUCell算法對每個細(xì)胞的每個regulon活性進(jìn)行打分侨嘀。對于一個regulon來說咬腕，比較細(xì)胞間的AUCell 得分可以鑒定出哪種細(xì)胞有顯著更高的subnetwork活性涨共。結(jié)果生成一個二進(jìn)制的regulon活性矩陣（binarized activity matrix）举反，這將確定Regulon在哪些細(xì)胞中處于“打開”狀態(tài)扒吁。

scenic分析原理

補(bǔ)充：

蛋白質(zhì)中功能的基本單元是domain,是一種特殊的三維結(jié)構(gòu)魁索，不同結(jié)構(gòu)的domain與其他分子特異性結(jié)合從而發(fā)揮功能粗蔚。與此類似鹏控，轉(zhuǎn)錄因子在于DNA序列結(jié)合時，其結(jié)合位點的序列也由于一定的特異性急前，不同轉(zhuǎn)錄因子結(jié)合的DNA序列的模式是不同的裆针。為了更好的描述結(jié)合位點序列的模式世吨，科學(xué)家們提出了motif的概念耘婚。

2.1 GENIE3

GENIE3是一種從基因表達(dá)量數(shù)據(jù)推斷基因調(diào)控網(wǎng)絡(luò)的方法沐祷。它訓(xùn)練預(yù)測數(shù)據(jù)集中每個基因表達(dá)的隨機(jī)森林模型赖临，并使用TF的表達(dá)作為輸入兢榨。然后使用不同的模型來得出TF的權(quán)重吵聪，并測量它們各自的相關(guān)性以預(yù)測每個靶基因的表達(dá)兼雄。

GENIE3的輸入為表達(dá)矩陣君旦，最好使用gene-summarized counts（可能是也可能不是UMIs）局蚀。其他單位琅绅，比如counts千扶，TPM和FPKM/RPKM也可以。但是要注意第一步的網(wǎng)絡(luò)相關(guān)分析基于共表達(dá)髓绽，一些作者建議也可以使用within-sample normalization比如TPM顺呕。

GENIE3的輸出是一個帶有基因株茶、潛在調(diào)節(jié)因子以及IM（importance measure）的表启盛。IM代表了TF（input gene）在預(yù)測靶標(biāo)時的權(quán)重僵闯。作者探索了幾種確定閾值的方法，最終選擇為每個TF建立多個潛在靶標(biāo)基因集：（1）設(shè)置幾個IM閾值（IM>0.001 and IM >0.005)

（2）選取每個TF的前50哥靶標(biāo)targets

（3）每個target gene保留top5堕阔，10棍厂，50個TFs颗味，然后按TF分開超陆。

在以上結(jié)果中，只有IM>0.001的links被算入浦马。

每個基因集接著被分為positive- and negetive-correlated targets來區(qū)分可能激活的和抑制的targets时呀。（TF和潛在靶標(biāo)的Spearman相關(guān)性計算）

最終谨娜，只有包含30個基因以上的基因集（TF共表達(dá)模型）被保留，用于下游分析坞靶。

2.2 RcisTarget

RcisTarget是i-cisTarget和iRegulon的motif富集框架的新R / Bioconductor實現(xiàn)。

RcisTarget從一個基因列表識別富集的TF-binding motifs和候選轉(zhuǎn)錄因子尿这。主要有兩步驟：

（1）選擇在基因集中基因TSS（transcription start site）附近顯著過表達(dá)的DNA motif碟摆。這一步通過在數(shù)據(jù)庫中應(yīng)用recovery-based method（基于恢復(fù)的方法）來實現(xiàn)的，該數(shù)據(jù)庫包含每個motif的全基因組跨物種排名嘉裤。保留注釋到對應(yīng)TF并且NES（normalized enrichment score）>3的motif。

（2）對于每一個motif和基因集典奉，RcisTarget預(yù)測候選靶標(biāo)基因踊淳，也就是基因集中排名領(lǐng)先的基因脱茉。這一步提供的結(jié)果跟i-cisTarget和iiRegulon相同。

為了構(gòu)建最終的regulon榜田，作者合并了每個TF module中預(yù)測的靶基因唱蒸，這些基因顯示了給定TF的任何motif的富集庆捺。但是在作者分析的數(shù)據(jù)中捉腥，這些modules數(shù)量很少而且motif富集很低坏匪。因此敦迄，最終決定從流程中去除對于直接表達(dá)的檢測，只使用positive-correlated targets進(jìn)行下游分析脾猛。

2.3 AUCell

對于一個給定的regulon，通過比較所有細(xì)胞間的AUCell（area under the recovery curve）打分值，我們可以識別哪些細(xì)胞具有更顯著高的regulon活性。

輸入為一個基因集，輸出為基因集每個細(xì)胞的‘a(chǎn)ctivity’糙箍。這些基因集即regulon，包含TFs和他們假定的的target雹拄≈式叮基于recovery analysis將根據(jù)表達(dá)水平將所有基因進(jìn)行排序禁悠。AUC代表了與細(xì)胞內(nèi)其他基因相比顾孽，特征基因中表達(dá)基因的比例及其相對表達(dá)值拦英。AUCell使用AUC來計算輸入基因集的關(guān)鍵子集是否在每個細(xì)胞的排名頂部都得到了富集霎冯。將輸出一個每個基因集在每個細(xì)胞的AUC score矩陣。

通過卡閾值得到的二元活性矩陣使矩陣維數(shù)減少（可理解為只有?0|1显晶，on|off）躏救，對于下游分析很有用崩掘。 例如抄瑟，基于regulon二元活性矩陣的聚類骂维，可以根據(jù)某個調(diào)控子網(wǎng)絡(luò)（regulon）的活性來識別細(xì)胞群類型和細(xì)胞狀態(tài)。由于regulon是整體評分的，而不是使用單個基因的表達(dá)，因此這種方法對于個別基因的dropouts很有效畏吓。

3.?分析方法-R版

3.1 輸入

單細(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é)果相差不大掂骏，但是不推薦TPM。

在進(jìn)行GENIE3分析前要對數(shù)據(jù)進(jìn)行過濾：

（1）過濾基因：對每個基因的總reads數(shù)進(jìn)行過濾，去除最可能不可信的只提供噪音的基因冬竟。具體的值取決于數(shù)據(jù)笑诅，文章中用到3 UMI counts × 30 (1% of cells) = minimum 90 counts per gene

（2）過濾基因：在多少細(xì)胞中被檢測到。去除只在一個活少量細(xì)胞中表達(dá)的基因，如果這些基因正好在一個細(xì)胞中集合砌梆，將獲得很大的權(quán)重烂瘫。推薦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ù)對應(yīng)物種下載相應(yīng)的數(shù)據(jù)庫文件到目錄下粟耻。下面示例為用R下載，也可以wget等其他方法。

### 準(zhǔn)備工作：下載評分?jǐn)?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 導(dǎo)入數(shù)據(jù)

根據(jù)不同的數(shù)據(jù)類型有不同的導(dǎo)入方式次哈，官方教程都有寫，這里針對seurat對象導(dǎo)入。

### 導(dǎo)入單細(xì)胞數(shù)據(jù)

## 準(zhǔn)備表達(dá)量矩陣

exprMat?<-?pt@assays$RNA@counts

exprMat?<-?as.matrix(exprMat)

##準(zhǔn)備細(xì)胞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è)置共啃，導(dǎo)入評分?jǐn)?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 為分析使用的線程挂滓，要考慮服務(wù)器配置情況設(shè)置

3.2.4 共表達(dá)網(wǎng)絡(luò)計算

使用GENIE3計算轉(zhuǎn)錄因子與每個基因的相關(guān)性贝椿，消耗資源非常大瑟蜈。我進(jìn)行測試時兩萬五細(xì)胞三萬基因的小鼠單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)跑了7天位迂，到現(xiàn)在還沒結(jié)束。處理大數(shù)據(jù)集的時候會非常非常慢，作者提出了兩種解決辦法：

（1）GENIE推斷共表達(dá)module時抽取少量細(xì)胞計算，在計算regulon活性（AUCell）時，使用所有細(xì)胞計算；

（2）采用python版本的SCENIC，具體步驟下面會說到。

我還是使用了最基礎(chǔ)的，而不是作者建議的方法。

##==轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)推斷==##

## 基因過濾，minCountsPerGene和minSamples采用默認(rèn)參數(shù)

# 過濾標(biāo)準(zhǔn)是基因表達(dá)量之和>細(xì)胞數(shù)*3%藕咏，且在1%的細(xì)胞中表達(dá)

genesKept?<-geneFiltering(exprMat, scenicOptions, minCountsPerGene?=?3?*?0.01?*ncol(exprMat), minSamples?=ncol(exprMat)?*?0.01)

exprMat_filtered?<-?exprMat[genesKept,]

## 計算相關(guān)性矩陣

runCorrelation(exprMat_filtered, scenicOptions)

## TF-Targets相關(guān)性回歸分析

exprMat_filtered_log?<-?log2(exprMat_filtered+1)

runGenie3(exprMat_filtered_log, scenicOptions)?## 這一步跑了七天

runGenie3會在int目錄產(chǎn)生不少結(jié)果文件：

·?1.2_corrMat.Rds 基因間相關(guān)性矩陣；

·?1.3_GENIE3_weightMatrix_part_1.Rds GENIE3中間結(jié)果坯沪，矩陣格式，有好多個，比如我這次分析，產(chǎn)生了10個；

·?1.4_GENIE3_linkList.Rds GENIE3的最終結(jié)果驻民，是把1.3那一堆文件合并在一起了虑乖。文件有三列：

TF為轉(zhuǎn)錄因子名稱，

Target為潛在靶基因的名稱，

weight時TF與Target之間的相關(guān)性權(quán)重。

3.2.5 Step1 - 推斷共表達(dá)模塊

GENIE3計算了轉(zhuǎn)錄因子與每個基因的相關(guān)性，接下來需要過濾掉低相關(guān)性的組合，形成共表達(dá)基因集。作者嘗試了多種過濾標(biāo)準(zhǔn)，沒有發(fā)現(xiàn)一種最佳策略，所以推薦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á)模塊绪钥；

# 推斷共表達(dá)模塊

runSCENIC_1_coexNetwork2modules(scenicOptions)

結(jié)果文件為1.6_tfModules_asDF.Rds，四列：

method為6種過濾方法

corr為runCorrelation(exprMat_filtered, scenicOptions)得到的結(jié)果，1代表激活，-1代表抑制，0代表中性。SCENIC只會采用值為1报亩，即激活的數(shù)據(jù)用于后續(xù)分析劲件。

3.2.6 Step2 - Motif驗證共表達(dá)模塊

上一步找到了每個轉(zhuǎn)錄因子強(qiáng)相關(guān)的靶基因牵辣，這一步要修建共表達(dá)模塊形成有生物學(xué)意義的調(diào)控單元（regulons）罢猪。

對每個共表達(dá)模塊進(jìn)行motif富集分析，保留顯著富集的motif；這一步依賴gene-motif評分?jǐn)?shù)據(jù)庫，行為基因泰鸡，列為motif余舶，值為排名，也就是我們下載的cisTarget數(shù)據(jù)庫。

使用數(shù)據(jù)庫對motif進(jìn)行TF注釋，得到高可信度铣焊、低可信度岛蚤。數(shù)據(jù)庫直接注釋和同源基因推斷的TF為高可信度她紫，使用motif序列相似性注釋的TF為低可信度結(jié)果朴乖。

用保留的motif對共表達(dá)模塊里的基因進(jìn)行打分（依據(jù)cisTarget數(shù)據(jù)庫），識別顯著高分的基因（理解為motif里這些基因的TSS很近）。

刪除共表達(dá)模塊中與motif評分不高的基因，剩下的基因集被稱為調(diào)控單元（regulon）梁剔。

## 推斷轉(zhuǎn)錄調(diào)控單元regulon

runSCENIC_2_createRegulons(scenicOptions)?# Toy run settings

# 可以通過coexMethod調(diào)整使用的過濾方法，默認(rèn)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：轉(zhuǎn)錄因子名稱

·?gene：TF靶基因名稱

·?nMotif：靶基因在數(shù)據(jù)庫的motif數(shù)量

·?bestMotif：最顯著富集的motif名稱

·?NES：標(biāo)準(zhǔn)富集分?jǐn)?shù)铃彰，分值越高越顯著

·?highConfAnnot：是不是高可信注釋

·?Genie3Weight：TF與靶基因的相關(guān)性權(quán)重

Step2_MotifEnrichment_preview.html

3.2.7 Step3 - Regulon活動評分與可視化

一個regulon=一個TF與其Target gene的基因集，regulon的名稱有兩種寫法：

·?TF名稱+extended+靶基因數(shù)目：轉(zhuǎn)錄因子與所有靶基因組成的基因調(diào)控網(wǎng)絡(luò)

·?TF名稱+靶基因數(shù)目：轉(zhuǎn)錄因子與高可信靶基因（即highConfAnnot=TRUE的基因）組成的基因調(diào)控網(wǎng)絡(luò)

AUCell對每個regulon在各個細(xì)胞中的活性進(jìn)行評分被饿。評分的基礎(chǔ)是基因表達(dá)值，分?jǐn)?shù)越高代表基因集的激活程度越高。這一步要用所有細(xì)胞做計算损俭。

##==regulon活性評分與可視化==##

##regulons計算AUC值并進(jìn)行下游分析

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 - 生成二進(jìn)制的regulon活性矩陣

runSCENIC_4_aucell_binarize(scenicOptions)

3.3 結(jié)果可視化

SCENIC通過幾個step即可繪制一些默認(rèn)的圖片

3.3.1 Regulon Activity heatmap

step3有所有regulon在細(xì)胞的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.

當(dāng)細(xì)胞種類比較多時可以用RSS來識別細(xì)胞類型特異性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評分恐似。

然后就可以挑選各個細(xì)胞類型代表性regulon繪制熱圖等等啦忧陪。

繪制過程中我發(fā)現(xiàn)有一個點就是，根據(jù)regulon list提取到norm的AUCscore之后蔗喂，再scale（計算z-score）比直接計算所有regulons的z-score再提取锈玉，繪制效果要好齐蔽。數(shù)據(jù)好的話勺美，可以得到下圖這樣的結(jié)果啦。

4. 分析方法-Python版

同樣不介紹安裝步驟，可查看官方說明蜈抓。

如果服務(wù)器配置了docker的話直接用docker比較簡便拾酝，拉一個鏡像就行了：

docker pull aertslab/pyscenic:0.10.0

這里主要記錄一下分析方法材诽。

4.1 輸入文件準(zhǔn)備

以下是官方教程用到的文件：

（1）表達(dá)量文件，cellinfo

### 導(dǎo)入單細(xì)胞數(shù)據(jù)

## 準(zhǔn)備表達(dá)量矩陣

exprMat?<-?epi@assays$RNA@counts

exprMat?<-?t(as.matrix(exprMat))

write.csv(exprMat,"expr.csv")

（2）cisTarget數(shù)據(jù)下載

從數(shù)據(jù)庫下載臀规，我做的是人的

wget?https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather .

wget?https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather .

（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)境的要求比較復(fù)雜吧斟赚。

4.2.1 grn

輸入文件為表達(dá)量文件公罕，這里要注意表達(dá)量矩陣expr.csv 為cell x gene张吉，即行名為基因，列名為細(xì)胞barcode。

輸出文件為expr_mat.adjacencies.tsv宦赠。

分析主要耗費的時間在這一步，我測的兩萬多細(xì)胞跑了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杂曲，和上一步的結(jié)果文件。

輸出文件為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

輸入文件為上一步的結(jié)果文件regulons.csv荚虚，和表達(dá)量文件expr.csv撬腾。

輸出為auc打分結(jié)果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. 結(jié)果展示方式

Construction of a human cell landscape at single-cell level