视频链接:https://www.youtube.com/watch?v=XmHDexCtjyw&list=PLjiXAZO27elC_xnk7gVNM85I2IQl5BEJN&index=11
练习地址:https://github.com/NBISweden/excelerate-scRNAseq/blob/master/session-trajectories/session-trajectories.md
一、为什么要做轨迹分析以及轨迹分析的定义
在进行标准的单细胞分析流程得到聚类结果后,可以进行轨迹分析,并进一步进行基因表达分析。
但并不是只有聚类结果才可进行轨迹分析,标准分析中的很多步骤都可以follow with轨迹分析。
在整个生命生长发育过程中,细胞都在不断从一种功能“状态”过渡到另一种功能“状态”(如下图)。处于不同状态的细胞表达不同的基因,产生蛋白质和代谢物的动态重复,从而完成它们的工作。当细胞在不同的状态间转化时,会经历转录重组的过程,一些基因被沉默,而另一些则被激活。这些瞬态通常很难描述,因为在更稳定的端点状态之间净化细胞可能很困难或不可能。但由于这个过程是连续发生的,我们可以使用轨迹推断(TI,trajectory inference)的方法可以根据测序的细胞(瞬时状态)之间表达模式的相似性对单细胞沿着轨迹进行排序,以此来模拟细胞动态变化的过程。也就是重建分化轨迹或者拟时间轴。
但值得注意的是,并不是所有样本都适合进行轨迹分析。
比如我们明确知道骨髓中存在分化中间态的细胞(存在从干细胞到成熟细胞的分化过程),因此骨髓的样本可以进行轨迹分析。但一些样本如PBMC中几乎都是分化成熟的细胞,虽然也可以做轨迹分析,但是意义不大。此外,存在分化过程的细胞也有差别。比如B细胞到浆细胞的分化更像是一个线形过程,不存在分支(branch),所以B细胞不适合做branching。但Th细胞向其亚型的分化可能就存在多个分支。
自2014年以来,TI的各种算法得到飞速发展,至2018年已有接近60种方法。
Trajectory Inference主要方法的pipeline总结:主要包括两个step:降维和轨迹建模
二、降维方法
降维的方法包括线性降维PCA,ICA等,和非线性降维TSNE,UMAP,DF等。在学习轨迹分析之前,先来了解两种之前接触的比较少的降维方法:ICA和DF。
1. ICA (Independent Component Analysis)独立成分分析
ICA是数据结构的一种方法(A method for decomposing the data)。monocle1使用的就是ICA方法。
ICA与PCA比较类似,PCA(对高斯分布的数据效果较好)是将高变基因分配到主要的主成分中,用主成分来进行后续分析。而ICA是将数据解构,从混杂的信号中分离原始的多个生物信号。
PCA和ICA的区别:
主成分分析假设源信号间彼此非相关,独立成分分析假设源信号间彼此独立。
主成分分析认为主元之间彼此正交,样本呈高斯分布;独立成分分析则不要求样本呈高斯分布。
ICA的缺点:
- ICA假设它找出来的生物信号都是相互独立的。
- 每个信号的来源都是非高斯分布。举例来说,在教室中放一些麦克风,在很多人都同时讲话的时候,我们可以使用ICA来对混杂信号进行解构以判断是谁在讲话。这些混杂的声音信号就是非高斯分布的。但是很多的生物学信号都是高斯分布的。对单细胞数据来说,也很难区分是高斯分布还是非高斯分布。
总结:ICA和PCA一样,是一种线性降维方法。常被用于评估数据的原始组成。在ICA中,这些原始信号被认为是互相独立的,而且,ICA会先假定单细胞数据是非高斯分布的,实际上往往不是这样。不同的信号在ICA分析中同等重要,但ICA不能确定实际有多少个信号源。
2. DF (Diffusion Maps)扩散映射
Diffusion maps是一种非线性降维方法。
Diffusion maps原理讲解视频:https://www.bilibili.com/video/av38891467/
Diffusion Map用的是Diffusion Process的方法。如果两个点距离较近,则从一个点随机行走到一个点的概率就大。反之,如果两个点距离较远,则从一个点随机行走到一个点的概率就小。Deffusion Map就是这样将两个点之间的距离转换成它们之间能够产生随机行走的概率 ,并用这个随机行走的过程去捕捉数据的neighborhood结构,从而将一个高维的扭曲的数据展开,变成一个低维的visualization。
简单来说,为了把可能性转化成距离,DM可以计算B到C的可能性,再计算A到C的可能性。根据公式,如果两种可能性差不多大,那么他们的差值就趋于0。说明A到B的过程可以通过C来很好的连接起来。
DM是一种非线性降维(UMAP和tSNE也是非线形降维)。点和点之间(也就是细胞和细胞之间的距离)是通过probability来计算的。
三、定义轨迹
在学习了ICA和DM两种降维方法后,现在我们想要建立细胞之间的关系,定义轨迹应该从哪里开始,在哪里结束。
1. MST (minimum spamming tree)
举个例子:下图中有很多点,每个点之间的距离都可以计算(比如使用DM来计算点和点之间的possibility)。将点连线,寻找一个所有的点之间距离加和最小的连接方式,得到的结果如黑色的粗线所示,这就是最小生成树。细胞数目越多,MST的轨迹构建越准确。
monocle1中使用的就是这种方式。如下图a:每个细胞都代表了高维空间中的一个点,将高维空间降维(使用PCA/ICA或UMAP/TSNE),随后使用MST定义细胞轨迹,并将细胞按照MST构建的生成树排序,标注上细胞类型,就可以得到细胞轨迹。
但值得注意的是,MST只能构建细胞轨迹,但不能告诉你端点是转录起始点还是终止点,也就是不知道是从哪个方向向哪个方向分化。所以如果有先验知识(比如干细胞向别的细胞分化),就会容易很多。
此外,由于MST没有循环,所以不适用于增殖细胞(细胞周期)样本。
2. RGE(Reverse graph embedding, i.e. DDRTree and others)
第二种方法叫做反向图嵌入。
如下图A和B:在使用MST进行轨迹推断时,由于最小生成树高度依赖于每个点的位置和点与点之间的距离,仅仅只是某个点的位置有些微变化就会得到完全不同的细胞轨迹。而REG的方法(图C)则是先对细胞进行聚类,再对细胞群的平均值进行轨迹构建。
Monocle2中使用的就是RGE方法(DDRTree)。
上图显示的是RGE的工作原理。每个细胞都代表了高维空间中的一个点,使用PCA或其他方法来对细胞进行降维后,根据假定的细胞cluster的中心点来对轨迹进行构建。随后计算细胞到假设轨迹的距离,并将细胞分配到距离细胞最近的轨迹cluster上,分配完成后对中心点进行更新,重建轨迹,再将二维轨迹投射到多维空间里,比较与原始数据的契合度,如果match的不好,就重新降维和构建轨迹,循环这个过程,直到细胞轨迹能充分反映原始data(类似降维中的TSNE和UMAP循环)。这时就可以选择轨迹的root(需要先验知识),并对拟时间轴或者发育轨迹进行定义。根据轨迹图中的分叉,还可以定义cell fate。
REG还衍生出了许多方法,比如PAGA、Slingshot、TSCAN、CellRouter等。
Monocle3进行聚类的原理与Monocle2类似
Monocle3的工作流程:scRNAseq数据 --> 预处理(标准化+PCA)--> 降维 --> 聚类 --> 拟时间轴的建立(DDRTree、SimplePPT、L1-graph)--> 差异分析
和Monocle2相比,Monocle3的主要update:
3. RNA velocity (gene expression trajectory)
RNA velocity是基于真实的转录动力学,可用于细胞基因表达的动态分化的研究。
如上左图,刚转录出的mRNA包含外显子和内含子,经过splicing切除内含子后,得到用于编码蛋白的spliced mRNA。spliced mRNA的丰度由未成熟mRNA的splicing速度和降解速率共同决定。如上中图:每个点代表一个细胞,在拟时间轴上,未经过剪切的mRNA的出现始终早于经过剪切的mRNA。如上右图:红色代表未经过剪切的mRNA,蓝色代表经过剪切的mRNA,可以看出,这些细胞的应该是从左往右分化的,因此Velocity可以用于定义轨迹的起点分支和终点。也就是说,Velocity可以在不知发育过程的前提下,预测谱系的方向(如下图)。
Velocity可以用于周期的轨迹
总结:
Which method should I use?