RSS特异性指数
继前一篇帖子[网络CSI评估基因关联性及regulon聚类模块化]之后，今天咱们接着聊点单细胞regulon方面的内容。
图注描述：
(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在每个细胞中的活性值RAS(regulon activity scores) ，这已经是想要的数据，可以从里面找到细胞类型特异性的regulon，以便后续生物学意义的探究。那么，如何从中找到细胞类型特异性的regulon呢？
显然，最简单粗暴且有效的方法是通过图形展示的方法 (如上图Fig C) 人为来筛选，只需将所有regulon在关注的细胞类型中的活性值通过图形过目一遍即可，这种方式除了有点费人和时间外其他都挺好。
或者，也可以通过简单的数据计算来筛选，比如可以计算regulon在每种细胞类型中的平均活性值并在不同细胞类型间比较来获取特异性的regulon。
不过，现在已经有人通过统计模型的方式造出了更加全面系统的轮子，定义了一个概念regulon specificity score (RSS)，可以方便快速地寻找到细胞类型特异性的regulon，你确定不试一试还要继续人肉查找？
那么，下面咱们就简单来说一说RSS的相关内容。Jensen-Shannon Divergence (JSD)可以用来评估两个随机变量分布间的差异，RSS基于此来计算特异性，最终的计算公式如下：
JSD取值范围在0-1之间，0表示分布一致，1表示分布完全不同。那么，从上面的公式咱们可以知道，RSS值越大则regulon的细胞类型特异性越高。
当然，RSS的计算无需自己实现，可以直接使用R包SCENIC中的函数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不能在细胞类型间直接比较！
为什么在不同细胞类型间不能直接比较RSS呢？回答这个问题，咱们得先搞清楚RSS的计算过程。RSS的值是在特定细胞类型中计算出得每一个regulon的特异性指数，每一种细胞类型都循环进行一次这样的计算。所以，不同regulon的RSS值在同种细胞类型中可以直接比较，而相同regulon的RSS值在不同细胞类型间没法直接比较。
为了给出更直观的解释，咱打个比方来说更通俗易懂些。就像一个班级进行了多科考试，如语文和数学，语文考了90，数学考了80，是不是说语文成绩就比数据好呢？答案是否定的。因为不同学科之间的内容不一样，所以成绩的好坏没法直接横向比较。如果想要比较两科间的成绩，可以通过数据转换的方式将成绩变为z-score，这样就可以反映出各科成绩偏离班级平均值的程度，从而反映出成绩的好坏。
同理，比较相同regulon在不同细胞类型间的特异性时也可以将RSS转换为z-score，从而反映出regulon在哪种细胞类型更具特异性。R包SCENIC中的函数plotRSS可以完成RSS的转换，并且可以用热图来展示regulon在各细胞类型间的特异性。
查看特定细胞类型：
plotRSS_oneSet(rss,'iCAF')
结果如下：
查看所有细胞类型：
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
结果如下：
plotRSS虽然好用，但其中的参数缺乏注释，使用时需要留心，尤其是zThreshold、thr两个参数。或者，也可以自行使用scale函数将细胞类型内的RSS转换为z-score。
注意了，细心的朋友也许发现了，标准化后的regulon个数由原先的387变成了142。这是因为plotRSS对结果进行了过滤，zThreshold参数可以控制过滤的阈值。除此之外，热图所展示的内容是用参数thr进一步过滤的结果。
