单细胞轨迹推断流程，记录谱系构建、伪时间分析和概念。

Source: Slingshot: Trajectory Inference for Single-Cell Data

本篇为slingshot摘要笔记。简介全单细胞谱系分析流程，特别是谱系重建过程和伪时序分析。 主要使用slingshot包。

1 Introduction

  1. how to understand unsupervised or semi-supervised manner?
  2. what is minimum spanning tree?

最小的输入是一个在降维空间内代表细胞的matrix，以及一个cluster labels的向量。通过这两个输入，我们可以：

  • 使用 getLineages 函数，在聚类上构建最小生成树（MST），从而确定全局谱系结构。

  • 利用 getCurves 函数拟合同步主曲线，构建平滑的线序并推断伪时间变量。

  • 使用内置可视化工具评估每个步骤的输出结果。

  • 第一个数据集（称为 "单一轨迹 "数据集）是在下文中生成的，旨在代表一个单一谱系，其中三分之一的基因与转变相关。

  • 第二个数据集（"分叉 "数据集）由坐标矩阵（如通过 PCA、ICA、扩散图等得到的坐标矩阵）和 k 均值聚类生成的聚类标签组成。这个数据集代表了一个分叉轨迹，我们可以利用它来演示slingshot提供的一些附加功能。

2 Upstream Analysis

2.1 Gene Filtering

第一步通常是数据降维和过滤，这可以提高下游分析的速度，同时保证信息的损失最小。

1
2
3
# Filer出表达量≥3，表达的cluster数量≥10的Gene
geneFilter <- apply(assays(sce)$counts,1,function(x){ sum(x >= 3) >= 10 }) 
sce <- sce[geneFilter, ]

2.2 Normalization

Another important early step in most RNA-Seq analysis pipelines is the choice of normalization method. This allows us to remove unwanted technical or biological artifacts from the data, such as batch, sequencing depth, cell cycle effects, etc.

Recommend normalization technique: scone.

The order of these steps may change depending on the choice of method:

  • ZINB-WaVE: 考虑技术变量的同时进行了降维处理
  • MNN: 在降维处理后校正批次效应。

因为我们使用的是模拟数据，我们不必担心批次效应或其他混杂因子，因此，我们的模拟数据采用完全量化归一化方法(full quantile normalization)，这是一种行之有效的方法，可强制每个细胞具有相同的表达值分布。

1
2
3
4
5
6
7
8
9
10
11
12
FQnorm <- function(counts){ 
	rk <- apply(counts,2,rank,ties.method='min') 
	# 对每列向量进行排名操作，返回排名的结果。结果是一个与`counts`相同大小的矩阵，其中的元素是对应列的排名值。
	counts.sort <- apply(counts,2,sort) 
	# 对每列的元素进行升序排序
	refdist <- apply(counts.sort,1,median) 
	# 计算每行的中位数
	norm <- apply(rk,2,function(r){ refdist[r] }) 
	rownames(norm) <- rownames(counts) 
	return(norm) 
} 
assays(sce)$norm <- FQnorm(assays(sce)$counts)

2.3 Dimensionality Reduction

slingshot 的基本假设是，转录相似的细胞将在一些降维的空间中彼此靠近。由于我们在构建谱系和测量伪时间时使用 Euclidean distances，因此数据的低维表示很重要。

有许多方法可用于此任务，我们将有意避免确定哪种方法是“最佳”方法的问题，因为这可能取决于数据类型、收集方法、上游计算选择和许多其他因素。我们将演示两种降维方法：主成分分析（PCA）和一致流形近似与投影（UMAP，通过uwot包）。

The "best" method depends on the type of data, method of collection, upstream computational choices, and many other factors.

在进行 PCA 时，我们不按基因的方差对其进行调整(scale)scale=FALSE，因为我们认为所有基因的信息量并不相同。我们希望在强表达、高变异的基因中发现信号，而不是通过强制基因间方差相等来抑制这种信号。在绘图时，我们确保设置长宽比，以免扭曲感知距离。

1
2
3
4
# PCA :通过PCA将数据降维到二维空间，并通过散点图展示数据在这两个主成分上的分布情况。帮助观察数据中的模式、群集或者异常值。
pca <- prcomp(t(log1p(assays(sce)$norm)), scale. = FALSE) 
rd1 <- pca$x[,1:2] 
plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)

以下是详细解释：

1
pca <- prcomp(t(log1p(assays(sce)$norm)), scale. = FALSE)
  • assays(sce)$norm：提取sce中的norm字段，即上一个代码段中通过FQnorm函数处理过的归一化数据。
  • log1p：对数据取对数，并加上1。这是一个常见的处理步骤，特别是对于计数数据，可以减小数据中的变异性。
  • t：对矩阵进行转置。
  • prcomp：进行主成分分析，得到数据的主成分。
  • scale. = FALSE：表示在进行PCA时不进行缩放。
1
rd1 <- pca$x[,1:2]
  • 提取PCA结果中的前两个主成分，存储在rd1中。这样，rd1包含了数据在第一和第二主成分上的投影。
1
plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)
  • 使用 plot 函数绘制散点图。
  • rd1 是包含了数据在前两个主成分上的投影的矩阵，col = rgb(0,0,0,.5)设置点的颜色为黑色，透明度为0.5，pch=16表示使用实心圆作为点的标记，asp = 1表示保持横纵坐标的比例一致。
1
2
3
4
5
# UMAP from uwot
library(uwot)
rd2 <- uwot::umap(t(log1p(assays(sce)$norm))) 
colnames(rd2) <- c('UMAP1', 'UMAP2') 
plot(rd2, col = rgb(0,0,0,.5), pch=16, asp = 1)

我们将向SingleCellExperiment对象添加两个降维（即PCA和UMAP），但是后续我们的分析则重点关注PCA结果。

2.4 CLustering Cells

slingshot的输入最好是一个有cluster labels的向量，因此我们最好聚类（如果没有，会被识别为1个整类）

这一步骤中识别的聚类将会被用来决定潜在谱系的全局结构（即，他们的数量number，他们在何时从哪里分支when they branch off from one another, 这些分支事件的大概位置approximate locations of those branching events），这是与经典的单细胞聚类的目的（识别在数据集中的所有生物相关细胞类型）不同的。例如，当决定Global lineage structure时，没必要去区分不成熟和成熟的神经元细胞，因为这俩可能落在同一个谱系上面。

我们有两种聚类方法，其假设都是相似的，假设低维空间的欧氏距离反映细胞间的生物学差异，这两种方法分别是Gaussian mixture modeling和k-means

Gaussian mixture modeling是在mclust包提出的，是基于Bayesian information criterion (BIC)方法的自动确定聚类数量的方法。

1
2
3
4
5
6
7
8
9
10
library(mclust, quietly = TRUE)

cl1 <- Mclust(rd1)$classification
colData(sce)$GMM <- cl1

library(RColorBrewer)
plot(rd1, col = brewer.pal(9,"Set1")[cl1], pch=16, asp = 1)
# - `col = brewer.pal(9,"Set1")[cl1]`：使用`brewer.pal`函数生成颜色调色板，共有9种颜色，然后根据`cl1`中的聚类结果为每个点选择对应的颜色。
# - `pch=16`：使用实心圆作为点的标记。
# - `asp = 1`：保持横纵坐标的比例一致。

k-means的聚类数量则取决于k值的选取，并且principal curves are quite robust to the choice of k。

1
2
3
4
cl2 <- kmeans(rd1, centers = 4)$cluster
colData(sce)$kmeans <- cl2

plot(rd1, col = brewer.pal(9,"Set1")[cl2], pch=16, asp = 1)

3 Using Slingshot

至此，我们已经具备了在模拟数据集上运行 slingshot 所需的一切。这是一个分两步走的过程，包括用基于聚类的最小生成树（MST）识别全局谱系结构，以及拟合描述每个谱系的同步主曲线。

slingshot的功能主要是通过2步实现的：

  • Identifing the global lineage structure with a cluster-based minumum spanning tree(MST)
  • fitting simultaneous principal curves to describe each lineage.

这两个步骤可以使用 getLineages 和 getCurves 函数单独运行，也可以使用包装函数 slingshot（推荐）一起运行。我们将使用包装函数来分析单轨迹数据集，但稍后会在分叉数据集上演示单个函数的用法。

即，可以单独地调用这两步，分别通过getLineages和getCurves，但是我们更推荐一次性使用slingshot。

我们输入的需要包括一个降维的坐标矩阵，一个聚类标签的向量集。在我们的单轨迹数据中，这些元素包含在了一个SingleCellExperiment对象中。

对于slingshot ，采用PCA降维方法和Gaussian mixture modeling识别聚类标签的运行示例如下：

1
2
3
4
5
6
# sce <- slingshot(sce, clusterLabels = 'GMM', reducedDim = 'PCA') 
# 也就是前面提到的：进行slingshot需要一个降维数据rd1，一个分类向量cl1
# 因此，这里提供的是这俩量，而非得到这两的方法PCA/GMM，这需要在前面的步骤做，俩结果做输入
sce <- slingshot(sce, clusterLabels = cl1, reducedDim = rd1 )

BiocManager::install("DelayedMatrixStats") # 根据报错补包

输出也是一个SingleCellExperiment对象，并且附有slingshot结果，即所有的结果都会存储在一个PseudotimeOrdering对象中，这个结果会被添加在原始 SingleCellExperiment 对象的colData里，可以通过如下面这样的语句访问：

1
2
# 查看slingshot result
colData(sce)$slingshot

此外，所有推断的拟时序变量（每个谱系1个），本例为slingPseudotime_1，也会被加到colData

1
2
3
4
5
6
colData(sce)
# DataFrame with 300 rows and 2 columns
#                 slingshot slingPseudotime_1
#      <PseudotimeOrdering>         <numeric>
# c1                   ####          0.685367
# c2                   ####          0.000000

要提取所有的 slingshot 结果到一个单独的对象中，我们可以使用以下两种函数，使用哪一种取决于我们想要哪一种形式：

1
2
3
4
5
# PseudotimeOrdering object: a SummarizedExperiment extension,
# which are flexible containers that will be useful for most purposes.
as.PseudotimeOrdering
# SlingshotDataSet object: primarily used for visualization, as several plotting methods are included with the package
as.SlingshotDataSet

下面是用拟时序标注颜色点的单轨迹数据的推断谱系的可视化：

1
2
3
4
5
6
7
8
9
10
summary(sce$slingPseudotime_1)

library(grDevices)
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)

plotcol <- colors[cut(sce$slingPseudotime_1, breaks=100)]

# plot(reducedDims(sce)$PCA, col = plotcol, pch=16, asp = 1)
plot(rd1, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2, col='black')

我们可以看到谱系结构是首先通过基于聚类的最小生成树估计的，通过使用type参数

1
2
3
# plot(reducedDims(sce)$PCA, col = brewer.pal(9,'Set1')[sce$GMM], pch=16, asp = 1)
plot(rd1, col = brewer.pal(9,'Set1')[cl1], pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2, type = 'lineages', col = 'black')

补充待解决：

  1. 原例中把pca的结果保存在了reducedDims(sce)$PCA，我们得到的slingshot结果中只有一个reducedDimNames(1): slingReducedDim，其value = NULL，所以pca值应如何恰当的保存？
  2. 如何定制想要的点/线的颜色？

4 Downstream Analysis

4.1 Identifying temporally dynamic genes

在运行slingshot以后，我们通常会对在整个发育过程中表达有变化的基因产生兴趣。我们将使用tradeSeq包演示这种类型的分析。

对于每个基因，我们将用一个负二项噪点分布拟合一个 general additive model (GAM)来对基因表达与拟时序的关系（这可能是非线性的）；然后我们检查其关系的显著性。

1
2
3
4
5
6
7
library(tradeSeq)

# fit negative binomial GAM
sce <- fitGAM(sce)

# test for dynamic expression
ATres <- associationTest(sce)

how to understand general additive model (GAM)?

我们可以基因p值选择top gene，然后用热图可视化其基因表达与伪时间的关系。 以下绘制了p值最小的前250个基因的表达热图

1
2
3
4
5
6
7
topgenes <- rownames(ATres[order(ATres$pvalue), ])[1:250]
pst.ord <- order(sce$slingPseudotime_1, na.last = NA)

heatdata <- assays(sce)$counts[topgenes, pst.ord]
heatclus <- sce$GMM[pst.ord]

heatmap(log1p(heatdata), Colv = NA, ColSideColors = brewer.pal(9,"Set1")[heatclus])

5 Detailed Slingshot Functionality

在此，我们将提供更多细节，并重点介绍 slingshot 软件包的一些附加功能。 我们将使用附带的 slingshotExample (即我们前面准备的第二个数据集）进行说明。

该数据集旨在表示低维空间中的细胞，并附带一组由 k-均值聚类生成的聚类标签集。我们将直接使用低维坐标矩阵，并将聚类标签作为附加参数提供，而不是构建需要基因水平数据的完整 SingleCellExperiment 对象。

5.1 Identifying global lineage structure

getLineages 函数的输入：一个 n × p 矩阵，一个长度为 n 的聚类结果向量。它使用最小生成树（MST）映射相邻聚类之间的连接，并识别通过这些连接表示谱系的路径。该函数的输出是一个 PseudotimeOrdering（伪时间排序），其中包含输入以及推断出的 MST（用 igraph 对象表示）和 lineages（聚类名称的有序向量）。

这种分析可以完全在无监督的情况下进行，也可以通过指定已知的起点和终点聚类在半监督的情况下进行。如果我们没有指定起点，slingshot 会根据解析性选择一个起点，最大限度地增加分裂前各系共享的聚类数量。如果没有分裂或多个聚类产生相同的解析得分，则任意选择起始聚类。在我们的模拟数据中，slingshot 选择簇 1 作为起始簇。

  • 我们一般建议根据先验知识（样本采集时间或已确定的基因标记）指定初始聚类。
  • 这种指定不会影响 MST 的构建，但会影响分支曲线的构建。
1
2
3
4
5
lin1 <- getLineages(rd, cl, start.clus = '1') 
lin1

plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16) 
lines(SlingshotDataSet(lin1), lwd = 3, col = 'black')

在这一步，slingshot还允许指定已知端点。在构建 MST 时，被指定为终端单元状态的群集将被限制为只有一个连接（即它们必须是叶节点）。这种限制可能会影响树的其他部分的绘制，如下一示例所示，我们将第 3 个簇指定为端点。

1
2
3
4
lin2 <- getLineages(rd, cl, start.clus= '1', end.clus = '3') 

plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16) 
lines(SlingshotDataSet(lin2), lwd = 3, col = 'black', show.constraints = TRUE)

这种类型的监督有助于确保结果与先前的生物学知识相一致。具体来说，它可以防止已知的末期细胞命运被归类为过渡状态。

为了更好地控制，我们还可以向 getLineages 传递一些额外的参数：

  • dist.method: Character，指定如何计算聚类之间的距离。
    • 默认是 "slingshot"，它使用的是聚类中心之间的距离，并通过它们的全联合协方差矩阵进行归一化(normalized by their full, joint covariance matrix)。
    • 如果聚类较小（细胞数少于数据维数），则会自动切换到使用对角联合协方差矩阵(diagonal joint covariance matrix)。其他选项包括 simple（欧氏）、scaled.full、scaled.diag 和 mnn（mutual nearest neighbor-based distance, 基于互近邻的距离）。
  • omega: numeric/logical，是一个粒度参数，允许用户设置连接距离的上限。表示每个真实聚类与人工 .OMEGA 聚类(在拟合 MST 后会删除)之间的距离。这对于识别不属于任何谱系的离群集群或分离不同的轨迹非常有用。
    • 如果值为 "true"（"真"），则表示应使用 1.5（median edge length of unsupervised MST, 无监督 MST 的中位边长）的启发式，这意味着最大允许距离为该距离的 3 倍。

构建 MST 后，getLineages 会识别穿过树的路径，并将其指定为谱系。在这一阶段，谱系将由一组有序的簇名组成，从根簇开始，以叶簇结束。getLineages 的输出是一个伪时间排序PseudotimeOrdering，它包含了从根簇到叶子簇的所有伪时间排序。

5.2 Constructing smooth curves and ordering cells

为了模拟这些不同支系的发展过程，我们将使用函数 getCurves 构建平滑曲线。使用基于所有细胞的平滑曲线可以消除细胞投影到片断线性轨迹顶点上的问题，并使slingshot对聚类结果中的噪声更加稳健(robust)。

为了构建平滑的线段，getCurves 采用了类似于（Hastie 和 Stuetzle，1989 年）中主曲线的迭代过程。当只有一条直线时，得到的曲线就是通过数据中心的主曲线，但有一个调整：初始曲线是根据聚类中心之间的线性连接而不是数据的PC1构建的。这一调整增加了稳定性，通常会加快算法的收敛速度。

Hastie, Trevor, and Werner Stuetzle. 1989. “Principal Curves.” Journal of the American Statistical Association 84 (406): 502–16.

当有两条或更多的谱系时，我们会在算法中增加一个额外的步骤：平均共享细胞附近的曲线(averaging curves near shared cells)。在尚未分化的细胞上，两条系谱应该能达成一致，因此在每次迭代时，我们都会对这些细胞附近的曲线进行平均。这样可以提高算法的稳定性，并产生平滑的分支线。

1
2
3
4
5
crv1 <- getCurves(lin1) 
crv1

plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16) 
lines(SlingshotDataSet(crv1), lwd = 3, col = 'black')

getCurves 的输出结果是更新后的PseudotimeOrdering，其中包含了同步主曲线以及关于如何拟合这些曲线的附加信息。slingPseudotime 函数按细胞提取每个细胞的伪时间值矩阵，其中的 NA 值表示未分配到特定谱系的细胞。slingCurveWeights 函数提取类似的权重矩阵，将每个细胞分配到一个或多个谱系。

可以使用 slingCurves 函数访问曲线对象，得到 principal_curve 对象列表。 这些对象由以下slot组成：

  • s：构成曲线的点矩阵。它们对应于数据点的正交投影。
  • ord：指数，用于将曲线上的细胞按投影顺序排列。
  • lambda：从起点到每个细胞的投影之间沿曲线的长度。值矩阵由 slingPseudotime 函数返回。
  • dist_ind：数据点与它们在曲线上的投影之间的距离平方。
  • dist：投影距离的平方和。
  • w：沿该谱系的权重向量。如果单元格位于分支事件之前，则有可能多个分支的权重都是 11（或非常接近 11）。返回值的矩阵。

5.3 Running Slingshot on large datasets

对于大型数据集，建议使用 slingshot/getCurves时使用 approx_points 参数。

  • 这样可以指定曲线的分辨率（即唯一点的数量, the number of unique points）。
  • 虽然 MST 构造是在聚类上运行的，但随着数据集规模的扩大，将所有点迭代投影到一条或多条曲线上的过程可能会增加计算负担。 因此，将 approx_points 的默认值设为 150150 或数据集中的细胞数（较小者为准）将大大降低探索性分析的计算成本，同时对结果轨迹的影响也最小。

对于最大化 "密集(dense) "曲线，设置 approx_points = FALSE，曲线的点数将与数据集中的细胞数相同。不过，需要注意的是，迭代曲线拟合过程中的每个投影步骤的计算复杂度将与 n2n^{2} 成正比（其中 nn 为细胞数）。在存在分支系的情况下，这些密集的曲线也会影响曲线平均和收缩的复杂性(curve averaging and shrinkage)。

我们建议的值为 100−200100 - 200，因此默认值为 150150。请注意，限制曲线上唯一点的数量并不会对唯一伪时间值的数量造成类似限制，如下所示。即使 approx_points 的值低至不切实际的 55，我们仍然可以看到伪时间值的全色渐变：

1
2
3
4
5
6
sce5 <- slingshot(sce, clusterLabels = 'GMM', reducedDim = 'PCA', approx_points = 5) 

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100) 
plotcol <- colors[cut(sce5$slingPseudotime_1, breaks=100)] 
plot(reducedDims(sce5)$PCA, col = plotcol, pch=16, asp = 1) 
lines(SlingshotDataSet(sce5), lwd=2, col='black')

5.4 Multiple Trajectories

在某些情况下，我们需要识别多个互不关联的轨迹。

Slingshot 在初始 MST 构建中通过引入一个名为omega的人工簇(artificial cluster)来处理这个问题。这个人造簇与每个真实簇之间有一个固定长度的间隔，使得任何两个真实簇之间的最大距离是这个长度的2倍。

实际上，这对 MST 中允许的最大边长设置了限制。设置 omega = TRUE 将实现一条经验法则，即允许的最大边长等于不含人工簇的 MST 边长中值的 3 倍（注：这相当于说 omega_scale 的默认值是 1.5）。

1
2
3
4
5
6
rd2 <- rbind(rd, cbind(rd[,2]-12, rd[,1]-6)) 
cl2 <- c(cl, cl + 10) 
pto2 <- slingshot(rd2, cl2, omega = TRUE, start.clus = c(1,11)) 

plot(rd2, pch=16, asp = 1, col = c(brewer.pal(9,"Set1"), brewer.pal(8,"Set2"))[cl2]) 
lines(SlingshotDataSet(pto2), type = 'l', lwd=2, col='black')

拟合 MST 后，slingshot会像往常一样拟合同步主曲线，并分别处理每条轨迹。

1
2
plot(rd2, pch=16, asp = 1, col = c(brewer.pal(9,"Set1"), brewer.pal(8,"Set2"))[cl2])
lines(SlingshotDataSet(pto2), lwd=2, col='black')

5.5 Projecting Cells onto Existing Trajectories

有时，我们可能只想使用细胞的一个子集来确定轨迹，或者我们可能会获得新数据，并将其投射到现有轨迹上。无论哪种情况，我们都需要一种方法来确定新细胞格沿着之前构建的轨迹的位置。为此，我们可以使用预测函数（由于该函数并非slingshot的原生函数，请参阅 "?predict,PseudotimeOrdering-method "以获取相关文档）。

1
2
3
4
5
6
7
8
# our original PseudotimeOrdering 
pto <- sce$slingshot 

# simulate new cells in PCA space 
newPCA <- reducedDim(sce, 'PCA') + rnorm(2*ncol(sce), sd = 2) 

# project onto trajectory 
newPTO <- slingshot::predict(pto, newPCA)

这将产生一个新的混合对象，其中包含原始数据中的轨迹（曲线），以及细胞的伪时间值和权重。作为参考，原始细胞在下面以灰色显示，但它们并不包含在predict的输出中。

1
2
3
4
newplotcol <- colors[cut(slingPseudotime(newPTO)[,1], breaks=100)] 
plot(reducedDims(sce)$PCA, col = 'grey', bg = 'grey', pch=21, asp = 1, xlim = range(newPCA[,1]), ylim = range(newPCA[,2])) 
lines(SlingshotDataSet(sce), lwd=2, col = 'black') 
points(slingReducedDim(newPTO), col = newplotcol, pch = 16)