pyscenic做為scenic的python版本乾巧，相比R版本實(shí)現(xiàn)了分析速度上的巨大提升句喜，尤其是最耗時(shí)的共表達(dá)網(wǎng)絡(luò)構(gòu)建步驟。實(shí)測(cè)時(shí)間如下沟于，數(shù)據(jù)集1 (4257 cell x 36601 gene)咳胃，10 cores ：grn 網(wǎng)絡(luò)構(gòu)建2小時(shí)；ctx 網(wǎng)絡(luò)純化16分鐘旷太；aucell 細(xì)胞活性1分鐘拙绊。另外，數(shù)據(jù)集2 (45000 cell x 33469 gene)泳秀，三個(gè)步驟的耗時(shí)分別為4小時(shí)标沪、20分鐘、5分鐘嗜傅。
??從上面的流程示意圖可知金句，整個(gè)轉(zhuǎn)錄因子分析步驟分成Co-expression、Motif-discovery吕嘀、Cell-scoring违寞，分別對(duì)應(yīng)pyscenic軟件的三個(gè)子命令grn、ctx偶房、aucell趁曼。
RcisTarget數(shù)據(jù)庫(kù)文件下載鏈接：
- TF：https://github.com/aertslab/pySCENIC/tree/master/resources
- Feather：motif-gene rank，https://resources.aertslab.org/cistarget/databases
- Tbl： motif-TF annotation棕洋，https://resources.aertslab.org/cistarget/motif2tf
測(cè)試數(shù)據(jù)來自《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》挡闰，從原始數(shù)據(jù)出發(fā)，降維聚類分群后掰盘，使用fibo細(xì)胞亞群做轉(zhuǎn)錄因子分析摄悯。
1、grn
??輸入表達(dá)矩陣和轉(zhuǎn)錄因子列表愧捕，軟件基于共表達(dá)構(gòu)建轉(zhuǎn)錄因子與潛在靶基因的調(diào)控網(wǎng)絡(luò)奢驯。構(gòu)建網(wǎng)絡(luò)的算法可以選擇GENIE3和GRNBoot2 (默認(rèn)參數(shù))，網(wǎng)絡(luò)構(gòu)建涉及隨機(jī)過程次绘，如果想要復(fù)盤時(shí)保持結(jié)果一致可以設(shè)置隨機(jī)種子瘪阁。
??另外撒遣，表達(dá)矩陣的格式可以接受csv和loom，如果是csv (rows=cells x columns=genes) 則需要參數(shù)--transpose管跺。表達(dá)矩陣在整個(gè)過程都有用到义黎，為了方便可以使用loom格式。
pyscenic grn --seed 21 \
--num_workers 10 \
--method grnboost2 \
--output fibo_grn.tsv fibo_count.loom allTFs_hg38.txt
2伙菜、ctx
??上一步得到了初始的調(diào)控網(wǎng)絡(luò)，這一步基于motif和TF的關(guān)系以及motif對(duì)基因調(diào)控潛能的排序來修剪初始網(wǎng)絡(luò)命迈。首先贩绕，將motif注釋到TF；接著壶愤，對(duì)motif調(diào)控的基因排序(feather格式的文件)淑倾，如果調(diào)控網(wǎng)絡(luò)里面出現(xiàn)靶基因則折線上升，如此就會(huì)形成下圖所示的曲線圖征椒，默認(rèn)選擇所有基因的前5%處 (標(biāo)虛線的地方) 曲線下的面積計(jì)算AUC值娇哆，再接著計(jì)算一個(gè)標(biāo)準(zhǔn)化的AUC值做為NES，默認(rèn)過濾閾值為3.0勃救；最后碍讨，根據(jù)三種方式進(jìn)一步修剪網(wǎng)絡(luò)：1. thresholds 根據(jù)閾值選擇motif調(diào)控基因排序的前百分比的靶基因，2. top_n_targets 每個(gè)TF保留前N個(gè)靶基因蒙秒，3. top_n_regulators 每個(gè)靶基因保留前N個(gè)TF勃黍。此時(shí)，得到的TF與靶基因的調(diào)控網(wǎng)絡(luò)叫做regulon (調(diào)控子)晕讲，其中保留的靶基因上游含有motif直接結(jié)合的位點(diǎn)覆获。
??其實(shí)，不難看出這一步包含的過程比較繁瑣瓢省，經(jīng)過各種修剪獲得最終的網(wǎng)絡(luò)弄息。所以，這一步相當(dāng)重要勤婚，可以適當(dāng)調(diào)整閾值得到符合實(shí)際需求的網(wǎng)絡(luò)摹量。提示：這一步提到的實(shí)現(xiàn)過程為個(gè)人理解。
feather=hg38_refseq-r80_10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather
fname=motifs-v9-nr.hgnc-m0.001-o0.0.tbl
pyscenic ctx --annotations_fname $fname \
--expression_mtx_fname fibo_count.loom \
--mode dask_multiprocessing \
--mask_dropouts \
--num_workers 10 \
--output fibo_ctx.gmt fibo_grn.tsv $feather
3馒胆、aucell
??最后荆永，評(píng)估每個(gè)regulon在所有細(xì)胞里面的活性，富集方法與上面的motif富集相似国章。需要注意的是該步驟也有隨機(jī)種子參數(shù)具钥，如果想保證后面復(fù)現(xiàn)結(jié)果最好設(shè)置一下參數(shù)。
pyscenic aucell --seed 21 \
--num_workers 10 \
--output fibo_aucell.csv fibo_count.loom fibo_ctx.gmt
??到此液兽，pyscenic流程就結(jié)束了骂删，后面的工作就是基于regulon活性矩陣探索數(shù)據(jù)了掌动。當(dāng)然，也可以將活性矩陣二值化宁玫，這樣每個(gè)regulon在細(xì)胞中只有on和off兩種狀態(tài)粗恢，可以使用R包AUCell來完成。
library(AUCell)
library(reshape2)
auc <- read.table('fibo_aucell.csv',header=T,row.names=1,sep=',',stringsAsFactors=F,check.names=F)
auc <- t(auc)
rownames(auc) <- gsub(',','', rownames(auc))
auc[1:3, 1:5]
case1_AAACCTGGTAGTACCT-1 case1_AAAGCAAAGCCTTGAT-1
AHR(+) 0.0680935 0.07413331
ARID3A(+) 0.1266152 0.11878014
ARNT2(+) 0.0000000 0.00000000
case1_AAAGTAGGTCCATCCT-1 case1_AAATGCCGTTCCACTC-1
AHR(+) 0.07159951 0.07028004
ARID3A(+) 0.05228222 0.09272742
ARNT2(+) 0.00000000 0.00428051
case1_AACTCCCAGTAGTGCG-1
AHR(+) 0.0655522976
ARID3A(+) 0.0814689810
ARNT2(+) 0.0004553734
cells_assignment <- AUCell_exploreThresholds(auc, plotHist=F, assignCells=T)
thresholds <- getThresholdSelected(cells_assignment)
regulonsCells <- setNames(lapply(names(thresholds),
function(x) {
trh <- thresholds[x]
names(which(auc[x,]>trh))
}),names(thresholds))
regulonActivity <- melt(regulonsCells)
binaryRegulonActivity <- t(table(regulonActivity[,1], regulonActivity[,2]))
binaryRegulonActivity[1:3, 1:5]
case1_AAACCTGAGACAGAGA-1 case1_AAACCTGGTAGTACCT-1
ARID3A(+) 0 0
ARNT2(+) 0 0
ARNTL(+) 0 0
case1_AAAGCAAAGCCTTGAT-1 case1_AAAGTAGCACAACTGT-1
ARID3A(+) 0 0
ARNT2(+) 0 0
ARNTL(+) 0 0
case1_AAAGTAGGTCCATCCT-1
ARID3A(+) 0
ARNT2(+) 0
ARNTL(+) 0
??最后的最后欧瘪，展示一下重現(xiàn)文章的結(jié)果 (右邊來自文章)眷射，雖然缺失兩個(gè)轉(zhuǎn)錄因子，可能是前期數(shù)據(jù)處理不同導(dǎo)致佛掖，但是其余的轉(zhuǎn)錄因子在兩個(gè)細(xì)胞亞群里面的趨勢(shì)與文章一致妖碉。
往期回顧
一網(wǎng)打盡scRNA矩陣格式讀取和轉(zhuǎn)化(h5 h5ad loom)
ggplot2 | 開發(fā)自己的畫圖函數(shù)
R包安裝的4種姿勢(shì)
clusterProfiler: No gene can be mapped | 怎么破？
R語(yǔ)言的碎碎念
