RSS特異性指數(shù)
??繼前一篇帖子[網(wǎng)絡(luò)CSI評估基因關(guān)聯(lián)性及regulon聚類模塊化]之后，今天咱們接著聊點單細(xì)胞regulon方面的內(nèi)容。
圖注描述：
(A) Rank for regulons in erythroblast cell based on regulon specificity score (RSS).
(B) Erythroblast cells are highlighted in the t-SNE map (red dots).
(C) Binarized regulon activity scores (RAS) (do Z score normalization across all samples, and set 2.5 as cutoff to convert to 0 and 1) for top regulon Lmo2 on t-SNE map (dark green dots).
??SCENIC流程會得到regulon在每個細(xì)胞中的活性值RAS(regulon activity scores) ，這已經(jīng)是想要的數(shù)據(jù)，可以從里面找到細(xì)胞類型特異性的regulon，以便后續(xù)生物學(xué)意義的探究。那么，如何從中找到細(xì)胞類型特異性的regulon呢？
??顯然，最簡單粗暴且有效的方法是通過圖形展示的方法 (如上圖Fig C) 人為來篩選，只需將所有regulon在關(guān)注的細(xì)胞類型中的活性值通過圖形過目一遍即可，這種方式除了有點費人和時間外其他都挺好。
??或者，也可以通過簡單的數(shù)據(jù)計算來篩選，比如可以計算regulon在每種細(xì)胞類型中的平均活性值并在不同細(xì)胞類型間比較來獲取特異性的regulon。
??不過，現(xiàn)在已經(jīng)有人通過統(tǒng)計模型的方式造出了更加全面系統(tǒng)的輪子，定義了一個概念regulon specificity score (RSS)，可以方便快速地尋找到細(xì)胞類型特異性的regulon，你確定不試一試還要繼續(xù)人肉查找？
??那么，下面咱們就簡單來說一說RSS的相關(guān)內(nèi)容。Jensen-Shannon Divergence (JSD)可以用來評估兩個隨機變量分布間的差異，RSS基于此來計算特異性，最終的計算公式如下：
??JSD取值范圍在0-1之間，0表示分布一致，1表示分布完全不同。那么，從上面的公式咱們可以知道，RSS值越大則regulon的細(xì)胞類型特異性越高。
??當(dāng)然，RSS的計算無需自己實現(xiàn)，可以直接使用R包SCENIC中的函數(shù)calcRSS來完成，代碼示例如下：
library(SCENIC)
auc[1:5,1:2]
case1_AAACCTGGTAGTACCT-1 case1_AAAGCAAAGCCTTGAT-1
AHR(+) 0.06809350 0.07413331
ARID3A(+) 0.12661524 0.11878014
ARNT2(+) 0.00000000 0.00000000
ARNTL(+) 0.06979274 0.04075715
ATF1(+)
head(anno_col)
celltype
case1_AAACCTGGTAGTACCT-1 iCAF
case1_AAAGCAAAGCCTTGAT-1 iCAF
case1_AAAGTAGGTCCATCCT-1 iCAF
case1_AAATGCCGTTCCACTC-1 iCAF
case1_AACTCCCAGTAGTGCG-1 iCAF
case1_AACTCTTCAGGTTTCA-1 iCAF
rss <- calcRSS(AUC=auc, cellAnnotation=anno_col$celltype)
rss[1:5,]
iCAF mCAF
AHR(+) 0.4074940 0.46928566
ARID3A(+) 0.3908104 0.47588790
ARNT2(+) 0.1097911 0.08567146
ARNTL(+) 0.4199884 0.42562778
ATF1(+) 0.4348941 0.43056537
dim(rss)
[1] 387 2
RSS不能在細(xì)胞類型間直接比較！
??為什么在不同細(xì)胞類型間不能直接比較RSS呢？回答這個問題，咱們得先搞清楚RSS的計算過程。RSS的值是在特定細(xì)胞類型中計算出得每一個regulon的特異性指數(shù)，每一種細(xì)胞類型都循環(huán)進(jìn)行一次這樣的計算。所以，不同regulon的RSS值在同種細(xì)胞類型中可以直接比較，而相同regulon的RSS值在不同細(xì)胞類型間沒法直接比較。
??為了給出更直觀的解釋，咱打個比方來說更通俗易懂些。就像一個班級進(jìn)行了多科考試，如語文和數(shù)學(xué)，語文考了90，數(shù)學(xué)考了80，是不是說語文成績就比數(shù)據(jù)好呢？答案是否定的。因為不同學(xué)科之間的內(nèi)容不一樣，所以成績的好壞沒法直接橫向比較。如果想要比較兩科間的成績，可以通過數(shù)據(jù)轉(zhuǎn)換的方式將成績變?yōu)?code>z-score，這樣就可以反映出各科成績偏離班級平均值的程度，從而反映出成績的好壞。
??同理，比較相同regulon在不同細(xì)胞類型間的特異性時也可以將RSS轉(zhuǎn)換為z-score，從而反映出regulon在哪種細(xì)胞類型更具特異性。R包SCENIC中的函數(shù)plotRSS可以完成RSS的轉(zhuǎn)換，并且可以用熱圖來展示regulon在各細(xì)胞類型間的特異性。
查看特定細(xì)胞類型：
plotRSS_oneSet(rss,'iCAF')
結(jié)果如下：
查看所有細(xì)胞類型：
rssPlot <- plotRSS(rss, col.low="#473172", col.mid="#20988b", col.high="#f9e920", verbose=F)
head(rssPlot$df)
Topic cellType RSS Z
9 ATF1(+) iCAF 0.4348941 1.0295290
10 ATF1(+) mCAF 0.4305654 0.1494359
51 CD59(+) iCAF 0.2870382 0.0000000
52 CD59(+) mCAF 0.5512076 1.1032444
53 CDX1(+) iCAF 0.4592326 1.2562512
54 CDX1(+) mCAF 0.3015446 0.0000000
dim(rssPlot$df)
[1] 142 4
rssPlot$plot
結(jié)果如下：
??plotRSS雖然好用，但其中的參數(shù)缺乏注釋，使用時需要留心，尤其是zThreshold、thr兩個參數(shù)。或者，也可以自行使用scale函數(shù)將細(xì)胞類型內(nèi)的RSS轉(zhuǎn)換為z-score。
??注意了，細(xì)心的朋友也許發(fā)現(xiàn)了，標(biāo)準(zhǔn)化后的regulon個數(shù)由原先的387變成了142。這是因為plotRSS對結(jié)果進(jìn)行了過濾，zThreshold參數(shù)可以控制過濾的閾值。除此之外，熱圖所展示的內(nèi)容是用參數(shù)thr進(jìn)一步過濾的結(jié)果。
