URL：https://www.bioconductor.org/packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html Slingshot：单细胞数据的轨迹推断 3.1 介绍：该部分将展示完整的单细胞谱系分析工作流程，特别强调谱系重建和拟时序推断的过程。利用(Street al. 2017 )slingshot中提出的软件包。
slingshot 的目标是利用细胞簇来揭示全局结构，并将这种结构转换为由一维变量（称为 “ 伪时间 ”）表示。该部分内容提供了一些工具，用于以无监督或半监督方式学习细胞簇关系，并构建代表每条脉络的平滑曲线，以及每一步的可视化方法。3.1.1 概括 Slingshot 的输入是一个表示降维空间中单元格的矩阵和一个群组标签向量。
有了这两个输入，就可以：使用 getLineages 函数在集群上构建最小生成树（MST），从而确定全局线序结构。使用 getCurves 函数拟合同步主曲线，构建平滑的线状结构并推断伪时间变量。使用内置可视化工具评估每个步骤的输出结果。3.1.2 数据载入 接下来，我们将教学中的两个模拟数据跑一遍流程：data1：第一个数据集（称为“单轨迹”数据集）在下面生成，旨在表示单个谱系，其中三分之一的基因与转变有关。
该数据集将包含在一个SingleCellExperiment 对象（Lun and Risso 2017）中，并将用于演示完整的工作流程。
library (Seurat) library (SingleCellExperiment) library (tidyverse) library (RColorBrewer) # generate synthetic count data representing a single lineage means <- rbind ( # non-DE genes matrix ( rep ( rep ( c ( 0.1 , 0.5 , 1 , 2 , 3 ), each = 300 ), 100 ), ncol = 300 , byrow = TRUE ), # early deactivation matrix ( rep ( exp ( atan ( (( 300 : 1 ) - 200 ) / 50 )), 50 ), ncol = 300 , byrow = TRUE ), # late deactivation matrix ( rep ( exp ( atan ( (( 300 : 1 ) - 100 ) / 50 )), 50 ), ncol = 300 , byrow = TRUE ), # early activation matrix ( rep ( exp ( atan ( (( 1 : 300 ) - 100 ) / 50 )), 50 ), ncol = 300 , byrow = TRUE ), # late activation matrix ( rep ( exp ( atan ( (( 1 : 300 ) - 200 ) / 50 )), 50 ), ncol = 300 , byrow = TRUE ), # transient matrix ( rep ( exp ( atan ( c (( 1 : 100 ) / 33 , rep ( 3 , 100 ), ( 100 : 1 ) / 33 ) )), 50 ), ncol = 300 , byrow = TRUE ) ) counts <- apply (means, 2 , function (cell_means){ total <- rnbinom ( 1 , mu = 7500 , size = 4 ) rmultinom ( 1 , total, cell_means) }) rownames (counts) <- paste0 ( 'G' , 1 : 750 ) colnames (counts) <- paste0 ( 'c' , 1 : 300 ) sce <- SingleCellExperiment ( assays = List ( counts = counts)) 结果展示：> head (counts)[ 1 : 5 , 1 : 10 ] c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 G1 1 1 0 0 0 0 1 0 2 1 G2 7 1 4 0 6 4 5 4 9 8 G3 11 2 4 5 8 8 5 3 12 13 G4 16 5 8 8 22 9 14 20 26 19 G5 21 6 20 17 32 26 20 27 30 > sce class : SingleCellExperiment dim : 750 300 metadata ( 0 ) : assays ( 1 ) : counts rownames ( 750 ) : G1 G2 ... G749 G750 rowData names ( 0 ) : colnames ( 300 ) : c1 c2 ... c299 c300 colData names ( 0 ) : reducedDimNames ( 0 ) : mainExpName : NULL altExpNames ( 0 ) : > data2：第二个数据集（“分叉”数据集）由一个坐标矩阵（通过 PCA、ICA、扩散图等获得）以及由以下函数生成的聚类标签组成。
该数据集代表了一条分叉轨迹，它将使我们能够展示所提供的一些附加功能。
library (slingshot, quietly = FALSE ) data ( "slingshotExample" ) rd <- slingshotExample $ rd cl <- slingshotExample $ cl dim (rd) # data representing cells in a reduced dimensional space 结果展示：dim (rd) # [1] 140 2 length (cl) # vector of cluster labels # [1] 140 3.2 上游分析：3.2.1 基因过滤 在单细胞系谱数据集分析时，我们首先需要降低数据的维度，而筛选掉不提供有用信息的基因是常见的第一步。
这将显著提升后续分析的速度，同时将信息损失降至最低。在基因筛选阶段，我们保留了在足够多的细胞中稳健表达的基因，这些基因至少足以构成一个细胞簇，使它们成为潜在的有趣的细胞类型标记基因。我们将最小簇大小设定为10个细胞，并定义一个基因至少有3个模拟读数为“稳健表达”。
geneFilter <- apply ( assays (sce) $ counts, 1 , function (x){ sum (x >= 3 ) >= 10 }) table (geneFilter) sce <- sce[geneFilter, ] > table (geneFilter) geneFilter FALSE TRUE 7 743 3.2.2 标准化 在实践中，数据标准化可以帮助我们从数据中去除不需要的技术或生物学误差，例如批次、测序深度、细胞周期效应等。
此处，推荐使用软件包scone( Cole and Risso 2018 )，ZINB-WaVE (Risso et al. 2018 ) 在考虑技术变量的同时进行降维，而 MNN (Haghverdi et al. 2018 )在降维后校正批次效应。本教学使用的是模拟数据，因此我们无需担心批次效应或其他潜在混杂因素。因此，将继续进行全分位数标准化，可强制每个细胞具有相同的表达值分布。
FQnorm <- function (counts){ rk <- apply (counts, 2 ,rank, ties.method= 'min' ) 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) 3.2.3 降维 Slingshot的核心假设：转录相似的细胞在降维空间中位置接近。
为了构建谱系和测量伪时间，数据的低维表示非常关键。在实际操作中，通常采用两种降维方法：主成分分析（PCA）和均匀流形逼近与投影（UMAP）。在进行PCA时，选择不按基因的方差进行缩放，因为认为并非所有基因的信息量都是相同的，更关注那些表达强烈且变异度高的基因，而不是通过强制基因间方差相等来削弱这些信号。在绘制图表时，会设置合适的长宽比，以避免感知上的距离扭曲。
UMAP是一种较新的降维技术，它在保持全局数据结构的同时，还能揭示局部结构特征。通过uwot软件包，我们可以应用UMAP来进行降维处理。
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 ) 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 结果。
reducedDims (sce) <- SimpleList ( PCA = rd1, UMAP = rd2) 3.2.4 细胞聚类 Slingshot分析的最终输入是细胞的聚类标签向量。如果未提供此信息，Slingshot将把数据视为单一聚类，并拟合一个标准的主曲线。然而，即使在预期只有单一谱系的数据集中，程序也推荐对细胞进行聚类，有助于发现新的分支事件。
这一步骤中识别的聚类将用于确定潜在谱系的全局结构（即它们的数量、它们彼此分支的时间点，以及这些分支事件的大致位置）。这与单细胞数据聚类的典型目标不同，后者旨在识别数据集中所有生物学上相关的细胞类型。在本分析中，将采用了两种聚类方法，它们都假设低维空间中的欧几里得距离反映了细胞之间的生物学差异：高斯混合模型和k均值聚类。
前者在mclust包中实现（Scrucca et al. 2016），并具有基于贝叶斯信息准则（BIC）自动确定聚类数量的方法。
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 ) cl2 <- kmeans (rd1, centers = 4 ) $ cluster colData (sce) $ kmeans <- cl2 plot (rd1, col = brewer.pal ( 9 , "Set1" )[cl2], pch= 16 , asp = 1 ) 至此，所有的准备数据已经完成，接下来对模拟数据集运行Slingshot分析。
这个过程包括两个步骤：首先是使用基于聚类的最小生成树（MST）来识别全局谱系结构，其次是拟合同时主曲线来描述每个谱系。这两个步骤可以分别使用getLineages和getCurves函数运行，也可以通过包装函数slingshot一起运行（推荐），接下来先使用包装函数展示轨迹推断，随后展示在分叉数据集（data2）上展示各个独立函数的使用。Slingshot包装函数在一次调用中执行轨迹推断的两个步骤。
所需的输入包括一个降维后的坐标矩阵和一组聚类标签。这些可以是独立的对象，或者在单细胞测序数据的情况下，是包含在SingleCellExperiment对象中的元素。
要使用PCA降维和高斯混合模型识别的聚类标签运行Slingshot，我们将执行以下操作：sce <- slingshot (sce, clusterLabels = 'GMM' , reducedDim = 'PCA' ) 如上所述，如果没有提供聚类结果，将假定所有细胞都属于同一个聚类，并将构建一个单一曲线。如果没有提供降维结果，Slingshot将使用reducedDims返回列表的第一个元素。
输出是一个包含Slingshot结果的SingleCellExperiment对象。所有结果都存储在一个PseudotimeOrdering对象中，该对象被添加到原始对象的colData中，可以通过colData(sce)$slingshot访问。此外，所有推断出的伪时间变量（每个谱系一个）也被单独添加到colData中。
要将所有Slingshot结果提取为一个对象，我们可以使用as.PseudotimeOrdering或as.SlingshotDataSet函数，具体取决于我们想要的形式。PseudotimeOrdering对象是SummarizedExperiment对象的扩展，这是一个灵活的存储对象。SlingshotDataSet对象主要用于可视化。如下，我们用伪时间着色的点来可视化单细胞数据的推断谱系。
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 ) lines ( SlingshotDataSet (sce), lwd= 2 , col= 'black' ) 还可以看到如何使用type参数通过基于聚类的最小生成树初步估计谱系结构。
plot ( reducedDims (sce) $ PCA, col = brewer.pal ( 9 , 'Set1' )[sce $ GMM], pch= 16 , asp = 1 ) lines ( SlingshotDataSet (sce), lwd= 2 , type = 'lineages' , col = 'black' ) 3.3 下游分析 3.3.1 识别时间动态基因 运行后slingshot，通常使用tradeSeq包对在发育过程中改变其表达的基因进行探索。
对于每个基因，使用负二项噪声分布拟合一般加性模型 (GAM)，以模拟基因表达和伪时间之间的（潜在非线性）关系。然后，计算基因表达和伪时间之间的显著关联 associationTest。进一步可以根据 p 值选出最显著的基因，并用热图可视化它们在发育过程中的表达情况。这里我们使用了表达最活跃的前 250 个基因。
library (tradeSeq) sce <- fitGAM (sce) ATres <- associationTest (sce) 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]) 3.4 功能详细解说 以下是对Slingshot包的详细功能说明和一些附加功能的强调：Slingshot将使用包含的slingshotExample数据集作为示例。
该数据集旨在表示降维的细胞基因表达数据，并附带了由k均值聚类生成的一组聚类标签。3.4.1 确定全局谱系结构 getLineages函数接受一个n×p的矩阵和一个长度为n的聚类结果向量作为输入。它使用最小生成树（MST）映射相邻聚类之间的连接，并识别通过这些连接的路径，这些路径代表谱系。该函数的输出是一个PseudotimeOrdering对象，其中包含输入数据以及推断出的MST（由igraph对象表示）和谱系（聚类名称的有序向量）。
这种分析通过完全无监督地进行，也可以通过指定已知的起始和终止点聚类以半监督的方式进行。如果不指定起始点，Slingshot会基于简洁性选择一个，最大化谱系分裂前共享的聚类数量。如果没有分裂或多个聚类产生相同的简洁性得分，起始聚类将被任意选择。在模拟数据的情况下，Slingshot选择聚类1作为起始聚类。然而，通常建议根据先验知识（无论是样本收集时间还是已建立的基因标记）指定一个初始聚类。
lin1 <- getLineages (rd, cl, start.clus = '1' ) lin1 # class: PseudotimeOrdering # dim: 140 2 # metadata(3): lineages mst slingParams # pathStats(2): pseudotime weights # cellnames(140): cell-1 cell-2 ... cell-139 cell-140 # cellData names(2): reducedDim clusterLabels # pathnames(2): Lineage1 Lineage2 # pathData names(0): plot (rd, col = brewer.pal ( 9 , "Set1" )[cl], asp = 1 , pch = 16 ) lines ( SlingshotDataSet (lin1), lwd = 3 , col = 'black' ) 在此步骤中，slingshot还允许指定已知端点。
在构建 MST 时，指定为终端单元状态的簇将被限制为只有一个连接。如下一个示例所示，我们将簇 3 指定为端点。
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是一个字符参数，用于指定如何计算聚类之间的距离。默认值是“slingshot”，它使用聚类中心之间的距离，并由它们的完整联合协方差矩阵进行归一化。在小聚类（细胞数量少于数据维度）存在的情况下，这将自动切换到使用对角联合协方差矩阵。
其他选项包括simple（欧几里得）距离、scaled.full（完整协方差矩阵的缩放版本）、scaled.diag（对角协方差矩阵的缩放版本）和mnn（基于相互最近邻的距离）。omega是一个粒度参数，允许用户设置连接距离的上限。它代表每个真实聚类与一个人工的.OMEGA聚类之间的距离，该聚类在拟合最小生成树（MST）后会被移除。这对于识别不属于任何谱系的离群聚类或分离不同的轨迹非常有用。这个参数可以是一个数值参数或逻辑值。
值为TRUE表示应使用启发式方法，即1.5倍的无监督MST中位数边长，这意味着允许的最大距离将是该距离的3倍。在构建MST之后，getLineages会识别通过树的路径，将其指定为谱系。在这个阶段，一个谱系将由一个有序的聚类名称集合组成，从根聚类开始，以叶聚类结束。getLineages的输出是一个PseudotimeOrdering对象，它保存了所有输入以及一系列谱系和一些关于它们如何构建的额外信息。
然后，这个对象被用作getCurves函数的输入。3.4.2 构建平滑曲线并排序细胞 为了更好地模拟细胞沿着不同谱系的发展路径，我们将使用一个名为getCurves的函数来绘制平滑的曲线 getCurves函数通过一个迭代过程来构建平滑的谱系曲线，这个过程与Hastie和Stuetzle在1989年提出的主曲线方法类似。
如果只有一个谱系，那么构建的曲线将穿过数据的中心点，类似于主曲线，但有所改进：初始曲线是基于聚类中心之间的直线连接来构建的，而不是基于数据的第一主成分。这种改进使得算法更稳定，并且通常可以加快算法的收敛速度。当存在多个谱系时，算法会增加一个步骤：在共同细胞附近对曲线进行平均处理。由于这些共同细胞尚未分化，不同的谱系在这些细胞上应该有很好的一致性。
因此，在每次迭代中，我们会对这些细胞附近的曲线进行平均，这样可以提高算法的稳定性，并生成平滑的分支谱系曲线。
crv1 <- getCurves (lin1) crv1 # class: PseudotimeOrdering # dim: 140 2 # metadata(4): lineages mst slingParams curves # pathStats(2): pseudotime weights # cellnames(140): cell-1 cell-2 ... cell-139 cell-140 # cellData names(2): reducedDim clusterLabels # pathnames(2): Lineage1 Lineage2 # pathData names(0): 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函数来得到曲线对象，该函数将返回一个主曲线对象列表。这些对象包含以下几个部分：s：构成曲线的点的矩阵。这些点对应于数据点的正交投影。
ord：可以用来根据细胞在曲线上的投影顺序排列细胞的索引。lambda：从曲线的起点到每个细胞投影的弧长。slingPseudotime函数返回这些值的矩阵。dist_ind：数据点与其在曲线上的投影之间的平方距离。dist：平方投影距离的总和。w：沿着这条谱系的权重向量。被明确分配到这条谱系的细胞将有一个权重为1，而被分配到其他谱系的细胞将有一个权重为0。如果细胞位于分支事件之前，它们可能对多个谱系有权重为1（或非常接近1）。
slingCurveWeights函数返回这些值的矩阵。3.4.3 在大型数据集上运行 Slingshot 对于大型的数据集，我们推荐在使用slingshot或getCurves时加入approx_points参数，以便控制曲线的精细程度，即曲线上有多少个独特的点。尽管最小生成树是基于聚类构建的，但是随着数据集变大，将所有数据点投影到曲线上的过程可能会变得计算成本很高。
因此，我们将approx_points的默认值设为150或者数据集中细胞数量的较小值，这样可以在对结果影响不大的情况下，显著降低计算成本。如果你想要得到密集的曲线，可以将approx_points设置为FALSE，这样曲线上的点数将与数据集中的细胞数相同。但请注意，这样做会使得每次迭代拟合过程中的计算复杂度增加到与细胞数量的平方成正比。此外，如果数据中存在分支谱系，这种密集的曲线还会增加曲线平均和收缩过程的复杂度。
我们建议将approx_points的值设定在100到200之间，这也是我们选择150作为默认值的原因。限制曲线上独特点的数量并不会影响到伪时间值的独特数量。即使将approx_points的值设得很低，比如5，我们仍然可以从伪时间值中观察到完整的颜色变化。
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' ) 3.4.4 多条轨迹 在某些情况下，可能需要识别多个感兴趣的轨迹。
Slingshot通过在最初的最小生成树（MST）构建中引入一个人造聚类，称为omega——来处理这个问题。这个人造聚类与每个真实聚类之间相隔一个固定的长度，这意味着任何两个真实聚类之间的最大距离是这个长度的两倍。实际上，这为MST中允许的最大边长设定了一个限制。设置omega = TRUE将采用一种经验规则，即允许的最大边长等于没有人造聚类时构建的MST的中位数边长的3倍（即：omega_scale的默认值为1.5）。
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对每条轨迹分别处理。
plot (rd2, pch= 16 , asp = 1 , col = c ( brewer.pal ( 9 , "Set1" ), brewer.pal ( 8 , "Set2" ))[cl2]) lines ( SlingshotDataSet (pto2), lwd= 2 , col= 'black' ) 3.4.5 将细胞投射到现有轨迹上 有时，可能只想使用细胞的子集来确定轨迹，或者可能会获得想要投射到现有轨迹上的新数据。
无论哪种情况，都需要一种方法来确定新细胞沿先前构建的轨迹的位置。为此，我们可以使用函数predict。
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) 这将产生一个新的组合对象，其中包含来自原始数据的轨迹（曲线），但包含新细胞的伪时间值和权重。
作为参考，原始细胞在下方显示为灰色。
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 )