使用cytoTRACE评估不同单细胞亚群的分化潜能
之前我在单细胞天地公众号看到了：单细胞水平看骨髓瘤的细胞状态和基因调控 ，文章使用cytoTRACE得到一个由低分化向高分化的轨迹，正常样本的浆细胞聚集在尾部（标志着分化的结束），而正在分化的细胞聚在头部。使用CytoTRACE看到骨髓瘤细胞具有比正常浆细胞更高的CytoTRACE score，标志着更好的分化潜能。如下所示：
骨髓瘤细胞具有比正常浆细胞更高的CytoTRACE score
感兴趣的可以自己去阅读该文章：《Dynamic transcriptional reprogramming leads to immunotherapeutic vulnerabilities in myeloma》
下面，我们来演示一下如何使用 cytoTRACE包 处理自己的单细胞表达量矩阵，并且得到 CytoTRACE score来辅助解释生物学现象。
在R里面安装cytoTRACE包
cytoTRACE目前有网页版，见：https://cytotrace.stanford.edu/，但是我们仍然是喜欢自己安装好后在R里面操作，需要去它的主页简单的同意后就可以看到下载链接：Download the CytoTRACE R package v0.3.3 using the following link: Download CytoTRACE v0.3.3.
下载这个 https://cytotrace.stanford.edu/CytoTRACE_0.3.3.tar.gz 后，就本地使用它这个文件安装，代码如下所示：
devtools::install_local("./CytoTRACE_0.3.3.tar.gz")
因为它的使用，同时依赖于操作系统里面的Python环境，并且需要安装两个依赖的Python模块，我们仍然是使用 reticulate 代替我们操作自己的conda来安装这两个模块，全部的R代码如下所示：
library(reticulate)
conda_create("cytoTRACE",python_version = '3.7')
use_condaenv("cytoTRACE")
conda_install("cytoTRACE", "numpy")
本来以为是无需去Linux终端命令行里面直接操作Python或者conda的，全程使用R语言，但是scanoramaCT这个Python模块一直失败，如下所示的代码 ：
conda_install("cytoTRACE", "scanoramaCT")
reticulate::py_install(packages = c ( "scanoramaCT" ))
所以仍然是不得不去Linux终端命令行里面直接操作Python，代码如下所示；
MacBook-Pro:CytoTRACE jmzeng$ which pip
/Users/jmzeng/Library/r-miniconda/envs/r-reticulate/bin/pip
MacBook-Pro:CytoTRACE jmzeng$ pip install ScanoramaCT
如果这个ScanoramaCT的Python模块安装失败，主要是影响 iCytoTRACE 功能，并不会影响 CytoTRACE功能。
认识cytoTRACE包
主要是需要认识3个函数：
- CytoTRACE: 对自定义 scRNA-seq 数据集进行CytoTRACE 分析的函数
- iCytoTRACE：跨多个 scRNA-seq 批次/数据集运行 CytoTRACE 的功能
- plotCytoTRACE: 用于生成 CytoTRACE、表型和基因表达的 2D 可视化的函数
然后，cytoTRACE包自带两个数据集：两个不同平台的骨髓分化 单细胞数据集（marrow_10x_expr和marrow_plate_expr），而且两个单细胞数据集都是带有相应的表型标签（marrow_10x_pheno和marrow_plate_pheno）
首先带领大家认识一下自带的数据集：
library(CytoTRACE)
dim(marrow_10x_expr)
marrow_10x_expr[1:4,1:2]
table(marrow_10x_pheno)
dim(marrow_plate_expr)
marrow_plate_expr[1:4,1:2]
table(marrow_plate_pheno)
可以看到， 里面的 细胞亚群数量都很多，但是每个单细胞亚群里面的单细胞数量就几百个，而且可以看到都是小鼠的单细胞数据集哦 ：
> dim(marrow_10x_expr)
[1] 13526 3427
> marrow_10x_expr[1:4,1:2]
10X_P7_2_AAACCTGCAGTAACGG 10X_P7_2_AAACGGGAGGACGAAA
Mrpl15 1 1
Lypla1 1 1
Tcea1 0 1
Atp6v1h 0 0
> table(marrow_10x_pheno)
marrow_10x_pheno
Erythrocytes Erythroid progenitors and erythroblasts
142 268
Granulocyte progenitors Granulocytes
330 770
Immature B cells Macrophages
299 222
Mature B cells Megakaryocyte progenitors
91 62
Monocyte progenitors Monocytes
322 522
Stem and progenitors
399
> dim(marrow_plate_expr)
[1] 17479 4442
> marrow_plate_expr[1:4,1:2]
A22.D042044.3_9_M.1.1 C5.D042044.3_9_M.1.1
0610005C13Rik 0 0
0610007C21Rik 0 112
0610007L01Rik 0 0
0610007N19Rik 0 0
> table(marrow_plate_pheno)
marrow_plate_pheno
Early immature B cells Early immature granulocytes
285 209
Late immature B cells Late immature granulocytes
863 390
Mature B cells Mature granulocytes
700 380
Monocyte progenitors and monocytes Stem and progenitors
324 1291
>
首先在单个数据集运行CytoTRACE
前面我们提到了 CytoTRACE函数: 对自定义 scRNA-seq 数据集进行CytoTRACE 分析的函数，而且它有两个自带的单细胞数据集，我们随意选择一个运行CytoTRACE函数即可：
library(CytoTRACE)
marrow_10x_expr_results <- CytoTRACE(mat = marrow_10x_expr)
#phenotype为每个barcodes的
plotCytoGenes(marrow_10x_expr_results, numOfGenes = 10)
前面我们看到了 marrow_10x_expr是 三千多个细胞的一万多个基因的表达量矩阵，在我的Mac电脑，一分钟不到，就出了结果：
CytoTRACE will be run on 3 sub-sample(s) of approximately 1142 cells each using 1 / 1 core(s)
Pre-processing data and generating similarity matrix...
Calculating gene counts signature...
Smoothing values with NNLS regression and diffusion...
Calculating genes associated with CytoTRACE...
Done
然后是针对结果的一些可视化：
plotCytoTRACE(marrow_10x_expr_results, phenotype = marrow_10x_pheno, gene = "Top2a")
ggplot2::ggsave('marrow_10x_expr_results_plotCytoTRACE.pdf',width = 12)
蛮奇怪的，上面的默认出图占画布窗口很大，而ggplot2::ggsave仅仅是能保存下面的图表里面的左边部分：
默认出图
所以我仔细看了看 plotCytoTRACE 代码，发现它自己有保存图表的功能，全部的完整的代码是：
plotCytoTRACE(marrow_10x_expr_results, phenotype = marrow_10x_pheno,
gene = "Top2a",outputDir = "marrow_10x")
一次性输出了很多图表哦，大家可以慢慢看，读懂它们。
可以看到这个时候每个单细胞都有一个CytoTRACE的打分，其实就可以把它当做是gsea或者gsva对基因集的打分，或者转录因子活性指数，一样的可视化啦。
比如不同单细胞亚群的CytoTRACE的打分的差异箱线图展示：
CytoTRACE的打分的差异箱线图
cytoTRACE还可以整合不同平台和数据集
前面我们提到了iCytoTRACE：跨多个 scRNA-seq 批次/数据集运行 CytoTRACE 的功能，然后，cytoTRACE包自带两个数据集：两个不同平台的骨髓分化 单细胞数据集（marrow_10x_expr和marrow_plate_expr）。现在我们也来试试看：
##整合不同平台不同数据集。
datasets <- list(marrow_10x_expr, marrow_plate_expr)
results <- iCytoTRACE(datasets)
###可视化
plotCytoTRACE(results,
phenotype = c(marrow_10x_pheno,marrow_plate_pheno),
gene = "Top2a",outputDir = "merge")
这个运行时间稍微长一点哦，三五分钟吧。可以看到了这两个数据集，整合的蛮好的。
整合的蛮好的
但是有一个疑惑，似乎这个cytoTRACE打分，其实就是细胞增殖这个生物学现象的反映，比如top2a和mki67这样的基因的表达量高低。
自己的单细胞数据集如何运行cytoTRACE呢
前面我们演示的都是cytoTRACE包自带两个数据集：两个不同平台的骨髓分化 单细胞数据集（marrow_10x_expr和marrow_plate_expr）。但是真实情况下，我们往往是需要量化自己的单细胞数据集，以大家熟知的pbmc3k数据集为例。大家先安装这个数据集对应的包，并且对它进行降维聚类分群，参考前面的例子：人人都能学会的单细胞聚类分群注释 ，而且每个亚群找高表达量基因，都存储为Rdata文件。
# 0.安装R包 ----
# devtools::install_github('caleblareau/BuenColors')
# utils::install.packages(pkgs = "ggstatsplot")
# InstallData("pbmc3k")
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
p1=DimPlot(sce,label = T)
因为这个数据集已经进行了不同单细胞亚群的标记，所以我们很容易拿到单细胞表达量矩阵和细胞对应的表型信息。全部的代码 如下所示：
table(is.na(pbmc3k.final$seurat_annotations))
####提取表型文件
table(pbmc3k.final$seurat_annotations)
phe <- pbmc3k.final$seurat_annotations
phe = as.character(phe)
names(phe) <- rownames(pbmc3k.final@meta.data)
####提取表达矩阵
mat_3k <- as.matrix(pbmc3k.final@assays$RNA@counts)
mat_3k[1:4,1:4]
results <- CytoTRACE(mat = mat_3k)
plotCytoGenes(results, numOfGenes = 10,outputDir = "pbmc_3k")
plotCytoTRACE(results, phenotype = phe,gene = "CCR7",outputDir = "pbmc_3k")
得到如下所示的结果：
终于解答我的一个疑问，就是我们前面的拟时序分析的时候，见：拟时序分析就是差异分析的细节剖析，提取了 CD14+ Mono 和 FCGR3A+ Mono这两个单细胞亚群，看其细节的变化，但是因为monocle其实得到的pseudotime值是没有方向，大小的顺序是可以调整的， 所以我们仅仅是凭借下面的拟时序图表没办法判断CD14+ Mono 和 FCGR3A+ Mono的分化关系：
但是这个时候，我们的CytoTRACE打分说明了CD14+ Mono 是远低于 FCGR3A+ Mono的，所以拟时序应该是FCGR3A+ Mono朝着CD14+ Mono 的方向。
当然了，后续我们还会使用RNA速率分析来辅助说明这一点，敬请期待哈， 持续关注我们的这个公众号专辑。本文开头介绍的文章：《Dynamic transcriptional reprogramming leads to immunotherapeutic vulnerabilities in myeloma》也是结合了 RNA速率分析来辅助说明CytoTRACE打分的生物学意义。
CytoTRACE 运行原理和步骤（摘抄）
下面的长篇大论，可以不看：
（1）基因计数：第一步是计算每个细胞中可检测表达的基因数量。这是通过对每个单细胞表达大于零的基因总数求和来完成的。
（2）基因计数特征（GCS）：第二步是捕捉表达模式与基因计数相关的基因。这是通过以下步骤完成的：
- 输入基因表达表被重新调整为每百万转录本 (TPM) 或每百万计数 (CPM)。
- 将每个单个细胞的转录本总和设置为该细胞中可检测表达的基因总数。这样做是为了将基因表达矩阵转换为相对转录物计数，或细胞裂解物中 mRNA 分子的估计丰度，我们和其他人已经证明这可以改进差异表达基因的检测（Gulati 等人，2020 年；Qiu 等人等，2017）。
- 生成的表达式矩阵是 log 2归一化的，伪计数为 1。
- 为了测量每个基因与基因计数的关系，计算每个基因的标准化表达和基因计数之间的 Pearson 相关性。
- 与基因计数最正相关的前 200 个基因的几何平均表达是基因计数特征 (GCS)。
(3) CytoTRACE：最后一步是通过利用细胞之间的局部相似性并应用两步平滑程序来迭代改进我们对 GCS 向量的估计：
- 创建我们的最近邻图，我们将归一化的表达矩阵（见上文）转换为马尔可夫过程，以捕捉单元格之间的局部相似性。
- 使用这个马尔可夫矩阵，然后我们将非负最小二乘回归 (NNLS) 应用到 GCS。这使我们能够将 GCS 表示为马尔可夫矩阵中捕获的不同转录邻域的函数。
- 应用扩散过程，根据马尔可夫过程的概率结构迭代调整 GCS。注意：这不是 GCS，而是经过 NNLS 调整的 GCS。
- 结果值在 0 和 1 之间进行排序和缩放，代表细胞的相对分化状态的预测顺序（0，分化程度更高；1，分化程度较低）。
CytoTRACE 的来龙去脉（摘抄）
下面的长篇大论，可以不看：
2020年1月24日，来自斯坦福大学的Newman团队的研究人员，在Science杂志上发表了基于scRNA-seq多样性来对细胞发育潜能进行鉴定和预测的方法，其题目为“Single-cell transcriptional diversity is a hallmark of developmental potential”。该文章证明了每个细胞表达的基因数可以作为一种更加简单但更强健的发育潜能的关键标志，并且研究者通过权衡转录组多样性，开放了一种计算方法--CytoTRACE，用于从scRNA-seq数据中预测分化状态。他们用CytoTRACE在对不同组织类型和物种的细胞进行了分化状态分析，解析了52个实验中定向发育轨迹的结果，CytoTRACE的表现比之前的分化轨迹鉴定方法和大约19,000注释基因集的预测结果更佳。此外，CytoTRACE有助于实现对静态干细胞的鉴定和揭示乳腺癌形成的诱导基因。该研究还建立了关键基于RNA的发育潜能特征来描绘细胞层次结构的分析平台。
该研究系统地评估了包括19,000个注释基因集的基于RNA的特征，以鉴定出不依赖于组织类型、物种和平台，并且更加准确预测细胞分化状态的主要因素。在此基础上，综合实验结果开发了一种非监督式框架，通过单细胞转录组预测相关分化状态。最后，将该预测方法与其他最先进的算法相比，探讨其在健康和癌症组织研究中，鉴定与干细胞和分化相关的关键基因，从而验证该预测方法的功能和优势。
该研究的目的是摆脱对发育方向或中间态等先验知识的依赖，独立地、更稳健地鉴定RNA发育潜能相关的关键要素。研究者先利用scRNA-seq数据，评估了19,000多个细胞间的潜在关联，这些数据几乎囊括了Molecular Signature Database(n=17,810)所有可用的基因集，和来自ENCODE、ChEA和mRNAsi的896个覆盖转录因子结合位点的基因集。参考了StemID、SCENT和SLICE三个衡量转录熵值来推断干性的计算技术，再进一步探讨了基因数（或每个细胞能检测到的表达基因），反映分化状态的实际效能。为此，研究者创建了一套训练集，由9个验证实验分化轨迹的金标准scRNA-seq数据库组成，将这些数据集与先前研究常用的标志数据集进行优先级排序，同时确保从哺乳动物受精卵到终末分化细胞的各个层次的广泛采样。
写在文末
我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com
如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：
We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.
十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你
腾讯云开发者
