1、简介 测试文件：请在Slingshot拟时序分析学习手册推送中获取。2、官方教学解读 URL：https://www.bioconductor.org/packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html Slingshot：单细胞数据的轨迹推断 2.1 介绍：该部分将展示完整的单细胞谱系分析工作流程，特别强调谱系重建和拟时序推断的过程。
利用(Street al. 2017 )slingshot中提出的软件包。slingshot 的目标是利用细胞簇来揭示全局结构，并将这种结构转换为由一维变量（称为 “ 伪时间 ”）表示。该部分内容提供了一些工具，用于以无监督或半监督方式学习细胞簇关系，并构建代表每条脉络的平滑曲线，以及每一步的可视化方法。2.1.1 概括 Slingshot 的输入是一个表示降维空间中单元格的矩阵和一个群组标签向量。
有了这两个输入，就可以：使用 getLineages 函数在集群上构建最小生成树（MST），从而确定全局线序结构。使用 getCurves 函数拟合同步主曲线，构建平滑的线状结构并推断伪时间变量。使用内置可视化工具评估每个步骤的输出结果。2.1.2 数据载入 接下来，我们将教学中的两个模拟数据跑一遍流程：data1：第一个数据集（称为“单轨迹”数据集）在下面生成，旨在表示单个谱系，其中三分之一的基因与转变有关。
该数据集将包含在一个SingleCellExperiment 对象（Lun and Risso 2017）中，并将用于演示完整的工作流程。
Sys.time # [1] "2025-04-28 14:12:56 CST" library (Seurat) library (SingleCellExperiment) library (tidyverse) library (RColorBrewer) # 设置工作目录 setwd ( "/data/06_scRNA-seq_psedutime/" ) # 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 0 1 1 0 1 0 1 0 0 1 # G2 3 2 0 2 6 3 1 0 4 6 # G3 3 8 1 4 10 3 1 7 6 2 # G4 7 15 1 16 14 5 3 5 17 10 # G5 12 34 5 12 21 11 10 20 25 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 # [1] 140 2 length (cl) # [1] 140 2.2 上游分析：2.2.1 基因过滤 在单细胞系谱数据集分析时，我们首先需要降低数据的维度，而筛选掉不提供有用信息的基因是常见的第一步。
这将显著提升后续分析的速度，同时将信息损失降至最低。在基因筛选阶段，我们保留了在足够多的细胞中稳健表达的基因，这些基因至少足以构成一个细胞簇，使它们成为潜在的有趣的细胞类型标记基因。我们将最小簇大小设定为10个细胞，并定义一个基因至少有3个模拟读数为“稳健表达”。
geneFilter <- apply ( assays (sce) $ counts, 1 , function (x){ sum (x >= 3 ) >= 10 }) table (geneFilter) # geneFilter # FALSE TRUE # 21 729 sce <- sce[geneFilter, ] 2.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) 2.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) 2.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) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 0.000 9.083 21.616 21.766 34.805 44.369 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' ) 2.3 下游分析 2.3.1 识别时间动态基因 运行后slingshot，通常使用tradeSeq包对在发育过程中改变其表达的基因进行探索。
对于每个基因，使用负二项噪声分布拟合一般加性模型 (GAM)，以模拟基因表达和伪时间之间的（潜在非线性）关系。然后，计算基因表达和伪时间之间的显著关联 associationTest。进一步可以根据 p 值选出最显著的基因，并用热图可视化它们在发育过程中的表达情况。这里我们使用了表达最活跃的前 250 个基因。
# BiocManager::install("tradeSeq") 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], col = colorRampPalette ( c ( " , " , " , " , " ))( 100 )) 2.4 功能详细解说 以下是对Slingshot包的详细功能说明和一些附加功能的强调：Slingshot将使用包含的slingshotExample数据集作为示例。
该数据集旨在表示降维的细胞基因表达数据，并附带了由k均值聚类生成的一组聚类标签。2.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函数的输入。2.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函数返回这些值的矩阵。2.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' ) 2.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' ) 2.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 ) 3、外部数据实操 接下来，我们从GEO数据库下载真实单细胞转录组测序数据，并使用R语言Slingshot包进行单细胞转录组测序 拟时序分析；
全流程包括：原始数据下载———测序数据预处理———降维聚类———细胞注释———利用关键细胞进行拟时序分析。3.1 单细胞数据下载：GSE163973 该数据集主要探讨了皮肤纤维化疾病中的成纤维细胞异质性，特别是以瘢痕疙瘩（KD）为研究模型。该研究使用真皮进行scRNA-seq分析，因为瘢痕疙瘩代表皮肤真皮纤维化疾病。
经过严格的质量控制，我们获得了 40655 个细胞的转录组（瘢痕疙瘩（nsample=3）：21488;正常疤痕（nsample=3）：19167）。我们下载该数据集中每个样本（其中包含3个KD样本，3个NC样本）的原始文件，即barcodes.tsv、genes.tsv和matrix.mtx文件。
3.2 测序数据预处理：接下来，我们使用R语言加载数据并构建疾病组与正常组样本的Seurat对象：3.2.1 载入相关的包 if ( ! require (Seurat)) install.packages ( "Seurat" ) if ( ! require (dplyr)) install.packages ( "dplyr" ) if ( ! require (ggplot2)) install.packages ( "ggplot2" ) if ( ! require (scran)) install.packages ( "scran" ) if ( ! require (patchwork)) install.packages ( "patchwork" ) if ( ! require (tidyr)) install.packages ( "tidyr" ) if ( ! require (scater)) install.packages ( "scater" ) if ( ! require (stringr)) install.packages ( "stringr" ) if ( ! require (tibble)) install.packages ( "tibble" ) if ( ! require (SingleCellExperiment)) install.packages ( "SingleCellExperiment" ) library (hdf5r) library (Matrix) library (openxlsx) library (clustree) library (paletteer) 3.2.2 处理各组样本 # 正常样本 rm ( list= ls ) normal_scar <- Read10X ( data.dir = "data/NS1/" ) class (normal_scar) normal_scar_obj <- CreateSeuratObject ( counts = normal_scar, project = "NS1" ) normal_scar_obj[[ "percent.mt" ]] <- PercentageFeatureSet (normal_scar_obj, pattern = "^MT-" ) normal_scar_obj <- subset (normal_scar_obj, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10 ) VlnPlot1 <- VlnPlot (normal_scar_obj, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) for (sample in c ( "NS2" , "NS3" )) { normal_scar1 <- Read10X ( data.dir = paste0 ( "data/" ,sample, "/" )) class (normal_scar1) normal_scar_obj1 <- CreateSeuratObject ( counts = normal_scar1, project = sample) normal_scar_obj1[[ "percent.mt" ]] <- PercentageFeatureSet (normal_scar_obj1, pattern = "^MT-" ) normal_scar_obj1 <- subset (normal_scar_obj1, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10 ) normal_scar_obj <- merge (normal_scar_obj,normal_scar_obj1) } # saveRDS(normal_scar_obj,file = "result_slingshot/normal_Object.rds") # 疾病样本 rm ( list= ls ) keloid_data <- Read10X ( data.dir = "data/KL1/" ) keloid_obj <- CreateSeuratObject ( counts = keloid_data, project = "KL1" ) keloid_obj[[ "percent.mt" ]] <- PercentageFeatureSet (keloid_obj, pattern = "^MT-" ) keloid_obj <- subset (keloid_obj, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10 ) VlnPlot1 <- VlnPlot (keloid_obj, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) for (sample in c ( "KL2" , "KL3" )) { keloid_data1 <- Read10X ( data.dir = paste0 ( "data/" ,sample, "/" )) class (keloid_data1) keloid_data_obj1 <- CreateSeuratObject ( counts = keloid_data1, project = sample) keloid_data_obj1[[ "percent.mt" ]] <- PercentageFeatureSet (keloid_data_obj1, pattern = "^MT-" ) keloid_data_obj1 <- subset (keloid_data_obj1, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10 ) keloid_obj <- merge (keloid_obj,keloid_data_obj1) } # saveRDS(keloid_obj,file = "result_slingshot/keloid_Object.rds") 3.2.3 合并两组数据，构建seurat_object rm ( list= ls ) normal_object <- readRDS ( "result_slingshot/normal_Object.rds" ) keloid_object <- readRDS ( "result_slingshot/keloid_Object.rds" ) seurat_object <- merge (keloid_object, y = normal_object) seurat_object[[ "percent.mt" ]] <- PercentageFeatureSet (seurat_object, pattern = "^MT-" ) VlnPlot1 <- VlnPlot (seurat_object, features = c ( "nFeature_RNA" , "nCount_RNA" , "percent.mt" ), ncol = 3 ) 3.2.4 数据标准化 # 数据标准化 seurat_object <- NormalizeData (seurat_object, normalization.method = "LogNormalize" , scale.factor = 10000 ) seurat_object <- FindVariableFeatures (seurat_object, selection.method = "vst" , nfeatures = 2000 ) top10 <- head ( VariableFeatures (seurat_object), 10 ) plot1 <- VariableFeaturePlot (seurat_object) plot2 <- LabelPoints ( plot = plot1, points = top10, repel = TRUE ) plot1 + plot2 # 数据缩放，对高变基因进行标准化 seurat_object <- ScaleData (seurat_object, features = VariableFeatures (seurat_object)) # 降维 seurat_object <- RunPCA (seurat_object, features = VariableFeatures ( object = seurat_object)) dims = 1:15, cells = 1000, balanced = TRUE) DimPlot (seurat_object, reduction = 'pca' , group.by = 'orig.ident' ) seurat_object <- JackStraw (seurat_object, num.replicate = 100 ) seurat_object <- ScoreJackStraw (seurat_object, dims = 1 : 20 ) PCA_Score <- JackStrawPlot (seurat_object, dims = 1 : 20 ) 3.3 细胞聚类 seurat_object <- FindNeighbors (seurat_object, reduction = "pca" , dims = 1 : 20 ) seurat_object <- FindClusters (seurat_object, resolution = 0.3 ) # Modularity Optimizer version 1.3.0 by Ludo Waltmand Nees Jan van Eck # # Number of nodes: 41167 # Number of edges: 1367889 # # Running Louvain algorithm... # Maximum modularity in 10 random starts: 0.9481 # Number of communities: 18 # Elapsed time: 25 seconds seurat_object <- RunUMAP (seurat_object, reduction = "pca" , dims = 1 : 20 ) head ( Idents (seurat_object), 5 ) # AAACCCAAGAACAGGA-1_1_1 AAACCCAAGATCCTAC-1_1_1 AAACCCAAGCCTGACC-1_1_1 # 0 4 0 # AAACCCAAGCGCAATG-1_1_1 AAACCCACAAACTCGT-1_1_1 # 0 0 # Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 DimPlot (seurat_object, reduction = "umap" , label = TRUE , group.by = "RNA_snn_res.0.3" ) 3.3.1 细胞簇注释 seurat_object <- JoinLayers (seurat_object) seurat_object.markers <- FindAllMarkers (seurat_object, group.by = "RNA_snn_res.0.3" , min.pct = 0.5 , logfc.threshold = 0.5 , only.pos = T) top10 <- seurat_object.markers %>% group_by (cluster) %>% dplyr :: filter (avg_log2FC > 1 ) %>% slice_head ( n = 5 ) %>% ungroup mesenchymal_stem_cell <- c ( "THY1" , "TAGLN" , "COL1A2" , "DCN" , "PECAM1" , "CDH5" , "VWF" , "KRT5" , "KRT14" , "KRT1" , "CD68" , "CSF1R" , "CYBB" , "CPA3" , "HPGDS" , "TPSB2" , "MLANA" , "CD3D" , "CD3E" , "KRT19" , "KRT7" , "KRT8" ) 3.3.2 定义细胞类型 new.cluster.ids <- c ( "Endothelial cells" , "Fibroblast cells" , "Fibroblast cells" , "Epidermal stem cells" , "Fibroblast cells" , "Fibroblast cells" , "Endothelial cells" , "Epidermal stem cells" , "Macrophages" , "Fibroblast cells" , "Endothelial cells" , "Mast cells" , "T cells" , "Melanocytes" , "Endothelial cells" , "Epidermal stem cells" , "Endothelial cells" , "Keratinocyte" ) names (new.cluster.ids) <- levels (seurat_object) new.cluster.ids # 0 1 2 # "Endothelial cells" "Fibroblast cells" "Fibroblast cells" # 3 4 5 # "Epidermal stem cells" "Fibroblast cells" "Fibroblast cells" # 6 7 8 # "Endothelial cells" "Epidermal stem cells" "Macrophages" # 9 10 11 # "Fibroblast cells" "Endothelial cells" "Mast cells" # 12 13 14 # "T cells" "Melanocytes" "Endothelial cells" # 15 16 17 # "Epidermal stem cells" "Endothelial cells" "Keratinocyte" seurat_object <- RenameIdents (seurat_object,new.cluster.ids) seurat_obj_Annotated <- seurat_object # saveRDS(seurat_obj_Annotated,file = "result_slingshot/seurat_obj_Annotated.rds") 3.3.3 可视化注释结果 rm ( list= ls ) seurat_object <- readRDS ( "result_slingshot/seurat_obj_Annotated.rds" ) seurat_object <- JoinLayers (seurat_object) seurat_object.markers <- FindAllMarkers (seurat_object, min.pct = 0.5 , logfc.threshold = 1 , only.pos = T) DimPlot (seurat_object, reduction = "umap" , label = TRUE ) 3.3.4 显著表达基因可视化 mesenchymal_stem_cell <- c ( "SELE" , "DCN" , "LGALS7" , "CD86" , "TPSAB1" , "PTPRC" , "SOX10" , "DCD" ) VlnPlot (seurat_object, features = mesenchymal_stem_cell, ncol = 2 ) 3.4 关键细胞 在下游分析中，我们可以选择一些关键细胞群进行拟时序分析，关键细胞可以通过生物标志物鉴定，也可以根据现有研究来定义，在这里我们选择”Fibroblast cells”,“Endothelial cells”,“Epidermal stem cells”做进一步的研究分析：seurat_object @ meta.data $ cell_type <- seurat_object @ active.ident seurat_object $ group <- ifelse ( grepl ( "^KL" , seurat_object $ orig.ident), "KD" , ifelse ( grepl ( "^NS" , seurat_object $ orig.ident), "NC" , "unknown" )) metadata <- data.frame (seurat_object @ meta.data) unique (seurat_object $ cell_type) subset_clusters <- subset (seurat_object, idents = c ( "Fibroblast cells" , "Endothelial cells" , "Epidermal stem cells" )) # saveRDS(subset_clusters,file = "result_slingshot/subset_clusters.rds") 3.5 拟时序分析 在这里，我们基于上面筛选的关键细胞簇进行后续拟时序分析：3.5.1 载入相关包 rm ( list= ls ) library (Seurat) library (SingleCellExperiment) library (tidyverse) library (RColorBrewer) 3.5.2 可视化细胞簇 subset_clusters <- readRDS ( "result_slingshot/subset_clusters.rds" ) DimPlot (subset_clusters, reduction = "umap" , label = TRUE ) mesenchymal_stem_cell <- c ( "PECAM1" , "CDH5" , "VWF" , "THY1" , "TAGLN" , "COL1A2" , "DCN" , "KRT5" , "KRT14" , "KRT1" ) DotPlot (subset_clusters, features = mesenchymal_stem_cell) + RotatedAxis + scale_color_gradientn ( colours = c ( ' , ' , ' , ' )) seurat_object <- subset_clusters 3.5.3 安装slingshot # if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("slingshot") 接下来，我们从现有的Seurat对象中，根据cluster_list指定的簇标识符，提取对应的细胞，并将这部分数据转换成SingleCellExperiment对象，以便进行进一步的分析或操作。
library (slingshot) cluster_list <- unique (seurat_object $ cell_type) select_cluster = as.SingleCellExperiment ( subset (seurat_object, idents = cluster_list), assay = 'RNA' ) select_cluster # class: SingleCellExperiment # dim: 33538382 # metadata(0): # assays(2): counts logcounts # rownames(33538): MIR1302-2HG FAM138A ... AC213203.1 FAM231C # rowData names(0): # colnames(38382): AAACCCAAGAACAGGA-1_1_1 AAACCCAAGATCCTAC-1_1_1 ... # TTTGTTGTCGAACCAT-1_2 TTTGTTGTCTTCTAAC-1_2 # colData names(9): orig.ident nCount_RNA ... group ident # reducedDimNames(2): PCA UMAP # mainExpName: RNA # altExpNames(0): colnames ( colData (select_cluster)) # [1] "orig.ident" "nCount_RNA" "nFeature_RNA" "percent.mt" # [5] "RNA_snn_res.0.3" "seurat_clusters" "cell_type" "group" # [9] "ident" 3.5.4 轨迹推断 sim <- slingshot (select_cluster, clusterLabels = 'cell_type' , reducedDim = 'UMAP' , start.clus= "Epidermal stem cells" , # 选择起点 end.clus = NULL # 这里我们不指定终点 ) 这里我们定义”Epidermal stem cells”为细胞发育起点，并且不规定终点，让算法进行模拟计算。
slingshot函数，用于单细胞数据分析中的谱系推断，帮助从高维数据中推测细胞发展轨迹。下面是各个参数的含义和作用：1. data 类型：数据对象，包含用于谱系推断的坐标矩阵。说明：这可以是 matrix、SingleCellExperiment、SlingshotDataSet 或 PseudotimeOrdering 类型的数据对象，通常是高维数据的降维结果（如PCA、t-SNE、UMAP等）。
2. clusterLabels 类型：向量或矩阵 说明：每个细胞的簇分配标签。如果是向量，则每个细胞对应一个标签；如果是矩阵，则为每个细胞与簇的加权关系，其中“-1”表示未分配的细胞。3. reducedDim 类型：字符或矩阵（可选） 说明：指定要使用的降维数据。如果提供了该参数，则可以指定哪一个降维结果（例如PCA、UMAP或t-SNE）用于谱系推断。
4. start.clus 类型：字符（可选） 说明：指定谱系推断的起始簇，通常是谱系树的根簇。5. end.clus 类型：字符（可选） 说明：指定强制作为叶节点的簇，通常是谱系树的终点簇。6. dist.method 类型：字符（可选） 说明：指定计算簇之间距离的方法，默认值为“slingshot”。可以使用不同的距离度量方法。
3.5.5 可视化计算结果 colnames ( colData (sim)) # [1] "orig.ident" "nCount_RNA" "nFeature_RNA" # [4] "percent.mt" "RNA_snn_res.0.3" "seurat_clusters" # [7] "cell_type" "group" "ident" # [10] "slingshot" "slingPseudotime_1" SlingshotDataSet (sim) # class: SlingshotDataSet # # Samples Dimensions # 38382 2 # # lineages: 1 # Lineage1: Epidermal stem cells Endothelial cells Fibroblast cells # # curves: 1 # Curve1: Length: 47.865 Samples: 38382 plot ( reducedDims (sim) $ UMAP, pch= 16 , asp = 1 ) lines ( SlingshotDataSet (sim), lwd= 2 , col= brewer.pal ( 9 , "Set1" )) legend ( "right" , legend = paste0 ( "lineage" , 1 ), col = unique ( brewer.pal ( 6 , "Set1" )), inset= 0.8 , pch = 16 ) 美化图形，添加渐变色。
colors <- colorRampPalette ( brewer.pal ( 11 , 'Spectral' )[ - 6 ])( 100 ) plotcol <- colors[ cut (sim $ slingPseudotime_1, breaks= 100 )] plotcol[ is.na (plotcol)] <- "lightgrey" plot ( reducedDims (sim) $ UMAP, col = plotcol, pch= 16 , asp = 1 ) lines ( SlingshotDataSet (sim), lwd= 2 , col= brewer.pal ( 9 , "Set1" )) legend ( "right" , legend = paste0 ( "lineage" , 1 ), col = unique ( brewer.pal ( 6 , "Set1" )), inset= 0.8 , pch = 16 ) 从模拟计算结果我们可以看到，3种关键细胞，由”Epidermal stem cells”向”Endothelial cells”，再向”Fibroblast cells”分化。
3.6 探究分化过程中的基因表达 3.6.1 tradeSeq 提供了一种灵活的回归模型拟合方法，可用于查找沿轨迹中的一个或多个谱系差异表达的基因。基于拟合模型，发现表达与伪时间相关的基因，或沿轨迹差异表达（在特定区域）的基因。它为每个基因拟合一个负二项式广义加性模型 （GAM），并对 GAM 的参数进行推理。
# if (!require("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("tradeSeq") 3.6.2 数据准备 counts <- sim @ assays @ data $ counts crv <- SlingshotDataSet (sim) counts：从模拟数据对象sim中提取基因表达矩阵（计数数据）。
sim@assays @data$counts提取计数矩阵，每行代表一个基因，每列代表一个细胞。crv：将sim转换为SlingshotDataSet对象。Slingshot是一种分析单细胞发育轨迹的方法，SlingshotDataSet对象描述了细胞间的发育轨迹及其分支结构。
3.6.3 确定最佳模型的参数 k set.seed ( 1234 ) icMat <- evaluateK ( counts = counts, sds = crv, k = 3 : 5 , nGenes = 20 , # 这里只做演示，实际计算请设置200-500 verbose = T) evaluateK：评估不同模型中参数k的最佳值，k表示基因表达模型中的平滑参数数量。参数说明：counts：输入的基因表达矩阵，用于拟合模型。
sds：发育轨迹数据（SlingshotDataSet对象），定义了细胞的轨迹和分支信息。k = 3:10：尝试从3到10个不同的平滑参数数量。nGenes = 200：随机选择200个基因进行评估，加快分析速度。verbose = T：显示详细的日志信息，方便跟踪运行状态。输出：icMat：存储不同k值下模型的评估结果，通常基于信息准则（如AIC、BIC）选择最优k值。
这一计算时间较长，结果展示如下：下图中，我们需要找到快速下降曲线的转折点；如果这个点小于6，则选择对应的数字；如果这个点大于等于6，则选择6。
这里我们选择 nknots = 6，进行下一步计算：内容解读可参考：https://statomics.github.io/tradeSeq/articles/tradeSeq.html 3.6.4 计算伪时间 pseudotime <- slingPseudotime (crv, na = FALSE ) slingPseudotime(crv, na = FALSE)：计算每个细胞的伪时间值。
crv 是一个 SlingshotDataSet 对象，包含了细胞的发育轨迹和分支结构。伪时间（Pseudotime） 是一种用于量化细胞在发育过程中所处“时间”点的概念。尽管我们没有真实的时间轴，伪时间可以帮助揭示细胞的发育或分化状态。该函数计算每个细胞沿轨迹的相对位置（伪时间），并将这些值存储在 pseudotime 中。na = FALSE：表示如果在计算过程中有缺失值，是否应跳过这些细胞。
设为FALSE时，缺失值将被处理（默认行为）。3.6.5 计算细胞权重 cellWeights <- slingCurveWeights (crv) slingCurveWeights(crv)：计算每个细胞在各个轨迹曲线上的权重。在发育轨迹分析中，细胞的表达模式可以影响它们在特定轨迹上的位置。该函数计算每个细胞在各个发育轨迹（或曲线）上的“权重”，即它们在不同轨迹分支的概率分配。
3.6.6 拟合广义加性模型（GAM） top_genes <- head ( order ( rowVars (counts), decreasing = TRUE ), 200 ) counts_subset <- counts[top_genes, ] # 约7min system.time ({ sce <- fitGAM ( counts = counts_subset, pseudotime = pseudotime, cellWeights = cellWeights, nknots = 6 , verbose = FALSE ) }) # user system elapsed # 393.259 66.290 473.043 fitGAM：拟合一个广义加性模型（GAM）来描述基因在伪时间中的表达变化。
counts：基因表达矩阵，其中行是基因，列是细胞。pseudotime：之前计算的每个细胞的伪时间值。它提供了每个细胞在发育轨迹中的相对时间点，作为GAM拟合的一个自变量。cellWeights：每个细胞在不同轨迹曲线上的权重，作为GAM拟合中的额外信息，用于提升模型的准确性。nknots = 6：指定在拟合GAM时使用的“结点”数量。结点数量影响模型的灵活性，通常越多结点，模型越灵活，但也可能导致过拟合。这里指定了6个结点。
3.6.7 基因表达与拟时序的相关性 assoRes <- associationTest (sce) head (assoRes) # waldStat df pvalue meanLogFC # COL1A1 99151.88 5 0 2.1863965 # MALAT1 10867.17 5 0 0.8055341 # COL1A2 107704.27 5 0 2.1889041 # COL3A1 84147.77 5 0 2.1775600 # RPLP1 50220.89 5 0 0.7993864 # RPS18 27798.16 5 0 0.6785392 3.6.8 起止点相关性最高的基因 startRes <- startVsEndTest (sce) head (startRes) # waldStat df pvalue logFClineage1 # COL1A1 29795.940 1 0 8.5624909 # MALAT1 4889.731 1 0 1.5013682 # COL1A2 30976.256 1 0 8.1200847 # COL3A1 24290.809 1 0 8.0383031 # RPLP1 5866.138 1 0 -0.6587992 # RPS18 8062.045 1 0 -0.7984236 startVsEndTest 是一个用于评估单细胞转录组数据中发育轨迹（lineage）起始点（start）和终止点（end）之间差异表达的函数。
该函数的主要目的是通过比较发育轨迹起始点与终止点的基因表达差异，来识别在发育过程中可能有重要作用的基因。startVsEndTest 的结果是一个包含差异表达基因（DEGs）分析的输出，通常包括统计量（如Wald统计量）、自由度、p值、对数变化值（logFC）等。你提供的输出是一个数据框，其中包含以下列：waldStat：Wald统计量用于检验基因在起始点和终止点之间是否存在显著的表达变化。
它是一个用于假设检验的统计量，通常用于广义线性模型中。数值越大，表示该基因在轨迹的起始和终止点之间的差异越显著。df：这是自由度（degrees of freedom），通常用于检验统计显著性的分布。由于每个基因的表达模型是基于伪时间拟合的，因此它的自由度是1，表示每个基因有一个独立的估计参数（即线性模型中的斜率）。pvalue：p值用于检验基因表达的差异是否显著。p值越小，表示该基因在起始点和终止点之间的表达差异越显著。
logFClineage1：logFClineage1 是基因在发育轨迹的第1条线路（lineage 1）中起始和终止点之间的对数变化（log2 fold change）。这是基因在起始点和终止点之间表达变化的对数变动。
3.6.9 按照相关性排序 oStart <- order (startRes $ waldStat, decreasing = TRUE ) # 挑选相关性最强的基因，并可视化 sigGeneStart <- names (sce)[oStart[ 1 ]] # > sigGeneStart # [1] "DCN" plotSmoothers (sce, counts, gene = sigGeneStart) 在UMAP图中展示基因的差异表达：plotGeneCount (crv, counts, gene = sigGeneStart) 3.6.10 使用ggplots可视化基因表达趋势 coldata <- data.frame ( celltype = sim @ colData $ cell_type, plotcol = sim @ colData $ group) rownames (coldata) = colnames (sim) # sce中对应信息取出 filter_coldata <- coldata[ colnames (sce),] # 添加拟时序信息 filter_coldata $ Pseudotime = sce $ crv $ pseudotime.Lineage1 # top6 genes top6 <- names (sce)[oStart[ 1 : 6 ]] top6_exp = sce @ assays @ data $ counts[top6,] top6_exp = log2 (top6_exp + 1 ) %>% t # 获得最终的清洁数据 plt_data = cbind (filter_coldata, top6_exp) colnames (plt_data) # [1] "celltype" "plotcol" "Pseudotime" "CD63" "SPARC" # [6] "S100A6" "VIM" "COL1A2" "LGALS1" 3.6.11 绘图 library (ggplot2) library (ggsci) library (ggpubr) library (RColorBrewer) getPalette = colorRampPalette ( brewer.pal ( 12 , "Paired" )) mycolors = getPalette ( length ( unique (plt_data $ celltype))) plt_list = list for (gene in top6) { print (gene) p = ggscatter ( data = plt_data, x = 'Pseudotime' , y = gene, color = 'celltype' , size = 0.6 ) + geom_smooth ( se = F, color = 'orange' ) + theme_bw + scale_color_manual ( values = mycolors) + theme ( legend.position = 'none' ) plt_list[[gene]] = p } # [1] "CD63" # [1] "SPARC" # [1] "S100A6" # [1] "VIM" # [1] "COL1A2" # [1] "LGALS1" library (patchwork) wrap_plots (plt_list) ggsave ( filename = 'result_slingshot/top6_genes.pdf' , width = 9 , height = 5 ) 4、下游探索 4.1 DimPlot和FeaturePlot # 提取拟时间 sl.info <- as.data.frame ( slingPseudotime (sim)) head (sl.info) # Lineage1 # AAACCCAAGAACAGGA-1_1_1 17.469736 # AAACCCAAGATCCTAC-1_1_1 43.886954 # AAACCCAAGCCTGACC-1_1_1 17.944789 # AAACCCAAGCGCAATG-1_1_1 21.178838 # AAACCCACAAACTCGT-1_1_1 20.885255 # AAACCCACAAGACGGT-1_1_1 1.142355 # 为了避免与Seurat原有的注释重合，我们给这些变量名加一个后缀 colnames (sl.info) <- "slingshot_pseudotime" colnames (sl.info) # [1] "slingshot_pseudotime" seurat_object <- AddMetaData (seurat_object, metadata = sl.info) # 可以利用Seurat原有的降维图观察细胞类型、拟时序中state、拟时间分布：library (ggsci) color <- c ( pal_npg ( 10 ), pal_jco ( 10 ), pal_lancet ( 9 ), pal_aaas ( 10 )) DimPlot (seurat_object, group.by = 'cell_type' , cols = color) | FeaturePlot (seurat_object, features = 'slingshot_pseudotime' , cols = c ( "lightgrey" , 'darkred' )) 4.2 VlnPlot library (ggplot2) # VlnPlot可以用于绘制Pseudotime这样的数值型变量 VlnPlot (seurat_object, features = 'slingshot_pseudotime' , group.by = 'cell_type' , pt.size = 0 ,, cols= color) + theme ( axis.text.x = element_text ( angle = 40 , hjust = 1 )) | VlnPlot (seurat_object, features = 'slingshot_pseudotime' , group.by = 'group' , pt.size = 0 , cols= color) + theme ( axis.text.x = element_text ( angle = 40 , hjust = 1 )) 4.3 ggplot2绘图 cellname <- rownames (sl.info) data4plot <- cbind (seurat_object @ meta.data, as.data.frame (seurat_object @ reductions $ umap @ cell.embeddings)[cellname,]) head (data4plot) # orig.ident nCount_RNA nFeature_RNA percent.mt # AAACCCAAGAACAGGA-1_1_1 KL1 8791 2335 3.810716 # AAACCCAAGATCCTAC-1_1_1 KL1 15880 3589 3.425693 # AAACCCAAGCCTGACC-1_1_1 KL1 12648 2989 3.913662 # AAACCCAAGCGCAATG-1_1_1 KL1 7076 1932 7.405314 # AAACCCACAAACTCGT-1_1_1 KL1 5235 1994 5.386819 # AAACCCACAAGACGGT-1_1_1 KL1 15944 2853 3.712995 # RNA_snn_res.0.3 seurat_clusters cell_type # AAACCCAAGAACAGGA-1_1_1 0 0 Endothelial cells # AAACCCAAGATCCTAC-1_1_1 4 4 Fibroblast cells # AAACCCAAGCCTGACC-1_1_1 0 0 Endothelial cells # AAACCCAAGCGCAATG-1_1_1 0 0 Endothelial cells # AAACCCACAAACTCGT-1_1_1 0 0 Endothelial cells # AAACCCACAAGACGGT-1_1_1 3 3 Epidermal stem cells # group slingshot_pseudotime umap_1 umap_2 # AAACCCAAGAACAGGA-1_1_1 KD 17.469736 2.496052 -9.570583 # AAACCCAAGATCCTAC-1_1_1 KD 43.886954 -2.775920 9.123234 # AAACCCAAGCCTGACC-1_1_1 KD 17.944789 2.970328 -8.379124 # AAACCCAAGCGCAATG-1_1_1 KD 21.178838 7.411248 -8.299144 # AAACCCACAAACTCGT-1_1_1 KD 20.885255 6.304723 -7.642037 # AAACCCACAAGACGGT-1_1_1 KD 1.142355 -13.288835 -3.727686 library ( "ggplot2" ) ggplot (data4plot, aes ( x= umap_1, y= umap_2)) + geom_point ( aes ( x= umap_1, y= umap_2, size= slingshot_pseudotime, colour = cell_type), data = data4plot) + scale_color_manual ( values = color) + scale_size_continuous ( range = c ( 0.3 , 1.5 )) + theme_classic + theme ( axis.line = element_line ( size= 0.6 , colour = "black" , arrow = ggplot2 :: arrow ( length = ggplot2 :: unit ( 0.1 , "inches" ), ends = "last" , type = 'open' ))) 4.4 细胞比例探索 这部分内容其实我们以前出过专门的教程，详情见沉浸式统计细胞比例。
在这里也一样，我们可以去统计拟时序的state和celltype之间的关系。
4.4.1 饼图 # 统计细胞类型 sl.info.sorted <- sl.info[ order (sl.info $ slingshot_pseudotime), , drop = FALSE ] # 获取行名（细胞名） cell_names <- rownames (sl.info.sorted) # 按长度划分成三等份 n <- length (cell_names) group1 <- cell_names[ 1 : floor (n / 3 )] group2 <- cell_names[( floor (n / 3 ) + 1 ) : ( floor ( 2 * n / 3 ))] group3 <- cell_names[( floor ( 2 * n / 3 ) + 1 ) : n] sl.info.sorted[group1, "new" ] <- "group1" sl.info.sorted[group2, "new" ] <- "group2" sl.info.sorted[group3, "new" ] <- "group3" sl.info.sorted[, c ( "cell_type" , "group" )] <- seurat_object @ meta.data[ rownames (sl.info.sorted), c ( "cell_type" , "group" )] sl.info.sorted $ cell_type <- factor (sl.info.sorted $ cell_type, levels = c ( "Epidermal stem cells" , "Endothelial cells" , "Fibroblast cells" )) mynames <- table (sl.info.sorted $ cell_type) %>% names mynames # [1] "Epidermal stem cells" "Endothelial cells" "Fibroblast cells" # 统计细胞数目 myratio <- table (sl.info.sorted $ cell_type) %>% as.numeric myratio # [1] 6017 13975 18390 # 整理细胞标签 pielabel <- paste0 (mynames, " (" , round (myratio / sum (myratio) * 100 , 2 ), "%)" ) pielabel # [1] "Epidermal stem cells (15.68%)" "Endothelial cells (36.41%)" # [3] "Fibroblast cells (47.91%)" cols <- pal_npg ( "nrc" )( 10 ) cols # [1] " " " " " " # [7] " " " " pie (myratio, labels= pielabel, radius = 1.0 , clockwise= F, main = "Celltype" , col = cols) mynames <- table (sl.info.sorted $ new) %>% names myratio <- table (sl.info.sorted $ new) %>% as.numeric pielabel <- paste0 (mynames, " (" , round (myratio / sum (myratio) * 100 , 2 ), "%)" ) pielabel # [1] "group1 (33.33%)" "group2 (33.33%)" "group3 (33.33%)" pie (myratio, labels= pielabel, radius = 1.0 , clockwise= F, main = "State" , col = cols) 4.4.2 堆积柱形图 cellnum <- table (sl.info.sorted $ cell_type,sl.info.sorted $ new) cellnum <- cellnum[ ! apply (cellnum, 1 , sum) == 0 ,] # 按列统计细胞比例 cell.prop <- as.data.frame ( prop.table (cellnum, margin= 2 )) colnames (cell.prop) <- c ( "Celltype" , "group" , "Proportion" ) head (cell.prop) # Celltype group Proportion # 1 Epidermal stem cells group1 0.4702204158 # 2 Endothelial cells group1 0.5292324527 # 3 Fibroblast cells group1 0.0005471315 # 4 Epidermal stem cells group2 0.0000000000 # 5 Endothelial cells group2 0.5399405972 # 6 Fibroblast cells group2 0.4600594028 p.bar <- ggplot (cell.prop, aes (group,Proportion, fill= Celltype)) + geom_bar ( stat= "identity" , position= "fill" ) + scale_fill_manual ( values= cols) + ggtitle ( "cell proportion" ) + theme_bw + theme ( axis.ticks.length= unit ( 0.5 , 'cm' )) + guides ( fill= guide_legend ( title= NULL )) p.bar 4.4.3 累计密度分布图 densyplot <- ggplot (sl.info.sorted, aes ( x= slingshot_pseudotime, group= cell_type, fill= cell_type)) + geom_density ( alpha= 0.3 , adjust= 1.5 ) + scale_fill_manual ( values = c ( pal_npg ( "nrc" )( 10 )) ) + theme_bw densyplot 4.5 基因表达量的探索 4.5.1 在轨迹图上绘制表达量 # 用所有的基因Scale,便于后续提取矩阵基因都在 seurat_object <- ScaleData (seurat_object, features = rownames (seurat_object)) plotGeneCount (sim, counts, gene = "CD63" ) library (circlize) mycolor <- colorRampPalette ( colors = c ( ' , 'white' , ' ))( 100 ) plotGeneCount (sim, counts, gene = "CD63" ) + scale_colour_gradientn ( colours = mycolor, guide = "colorbar" ) 4.5.2 热图 其实热图的教程我们做过很多：热图绘制天花板：ComplexHeatmap|
一、初探 热图绘制天花板| 二、两万字长文探索注释功能 ComplexHeatmap|
三、用富集分析结果做成词云注释热图 在理解了monocle对象的数据结构和绘制热图的语法后，大家完全可以自行拼凑出下面的图，流程有些长，有点耐心哦：首先整理一下注释信息和表达矩阵 library (ComplexHeatmap) # 注意这里绘制热图如果出现报错 "Image must have at least 1 frame to write a bitmap",因为完整的矩阵是有3w多个细胞，当绘制矩阵过大时候就会出现这样的情况 sce.obj_sub <- subset (seurat_object, downsample= 500 ) expr4plot <- GetAssayData (sce.obj_sub, slot = "scale.data" ) sl.info <- sl.info[ rownames (sl.info) %in% colnames (expr4plot),,drop = F] sl.info <- sl.info[ order (sl.info $ slingshot_pseudotime),,drop = F] expr4plot <- expr4plot[, rownames (sl.info)] dim (sl.info) # [1] 1500 1 dim (expr4plot) # [1] 33538 1500 # 设置顶部的psedotime、celltype的注释 celltype_colors <- c ( "Endothelial cells" = " , "Fibroblast cells" = " , "Epidermal stem cells" = " ) pseudotime_colors <- colorRamp2 ( range (sl.info $ slingshot_pseudotime, na.rm = TRUE ), c ( " , "green" ) # 从浅黄到红 ) top_anno <- HeatmapAnnotation ( Pseudotime= sl.info $ slingshot_pseudotime, Celltype= sce.obj_sub @ meta.data[ rownames (sl.info), "cell_type" ], col = list ( Pseudotime = pseudotime_colors, Celltype = celltype_colors)) # 可以先看一下注释的效果 draw (top_anno) # 找一些感兴趣的基因集，比如各种细胞的marker，这里以每种细胞的Top20为例，实际过程中大家可以选择更多的基因，便于富集到有意义的通路：# 先展示基因来源 sub.marker <- FindAllMarkers (sce.obj_sub, only.pos = TRUE , min.pct = 0.25 , logfc.threshold = 0.25 ) top20gene <- sub.marker %>% group_by (cluster) %>% top_n ( n = 20 , wt = avg_log2FC) left_ann <- rowAnnotation ( cellmarker= top20gene $ cluster) draw (left_ann) 来做个富集分析用于词云注释，富集分析的教程可以看这里：手把手教你做单细胞测序数据分析
（七）——基因集富集分析 library (clusterProfiler) library (org.Hs.eg.db) pathway4wc <- data.frame for (runcell in unique (top20gene $ cluster)) { try ({pathway4wc}) gene.df <- filter (top20gene,cluster == runcell) %>% pull (gene) gores <- enrichGO (gene.df, org.Hs.eg.db, keyType = "SYMBOL" , ont = "BP" , pvalueCutoff = 0.05 , pAdjustMethod = "BH" , qvalueCutoff = 0.2 , minGSSize = 10 , maxGSSize = 500 , readable = FALSE , pool = FALSE ) go.all.res <- as.data.frame (gores @ result) pathway.dt <- top_n (go.all.res, - 10 , p.adjust) pathway.dt <- cbind (pathway.dt, data.frame ( cluster= rep (runcell, nrow (pathway.dt)))) pathway4wc <- rbind (pathway4wc,pathway.dt) } head (pathway4wc) # ID Description GeneRatio BgRatio # GO:0045446 GO:0045446 endothelial cell differentiation 5/15 121/18903 # GO:0003158 GO:0003158 endothelium development 5/15 139/18903 # GO:0061028 GO:0061028 establishment of endothelial barrier 3/15 48/18903 # GO:0043114 GO:0043114 regulation of vascular permeability 3/15 49/18903 # GO:0001885 GO:0001885 endothelial cell development 3/15 65/18903 # GO:0002064 GO:0002064 epithelial cell development 4/15 213/18903 # pvalue p.adjust qvalue geneID # GO:0045446 2.821312e-08 1.704073e-05 9.503368e-06 PECAM1/CLDN5/SOX17/KDR/ROBO4 # GO:0003158 5.661023e-08 1.709629e-05 9.534354e-06 PECAM1/CLDN5/SOX17/KDR/ROBO4 # GO:0061028 6.843439e-06 1.100228e-03 6.135812e-04 PECAM1/CLDN5/ROBO4 # GO:0043114 7.286277e-06 1.100228e-03 6.135812e-04 CLDN5/PLVAP/TACR1 # GO:0001885 1.714324e-05 1.953748e-03 1.089577e-03 PECAM1/CLDN5/ROBO4 # GO:0002064 1.940809e-05 1.953748e-03 1.089577e-03 PECAM1/CLDN5/KDR/ROBO4 # Count cluster # GO:0045446 5 Endothelial cells # GO:0003158 5 Endothelial cells # GO:0061028 3 Endothelial cells # GO:0043114 3 Endothelial cells # GO:0001885 3 Endothelial cells # GO:0002064 4 Endothelial cells # 构建词云的text # 遍历每个唯一的聚类 mytext <- list for (i in 1 : length ( unique (pathway4wc $ cluster))) { # 选择p.adjust最小的前5个通路 temp.text <- filter (pathway4wc,cluster == unique (pathway4wc $ cluster)[i]) %>% top_n ( - 2 ,p.adjust) %>% pull (Description) # 将当前聚类的通路描述信息添加到列表中 mytext[[i]] <- temp.text } # 为列表元素命名，使用聚类的唯一值 names (mytext) <- unique (pathway4wc $ cluster) # 依据显著性设置词云大小，越显著的通路字体越大 # 创建一个空向量用于存储每个聚类的p.adjust值 myfontsize <- c # 遍历每个唯一的聚类 for (i in 1 : length ( unique (pathway4wc $ cluster))) { # 筛选出当前聚类且p.adjust小于0.05的数据 temfs <- filter (pathway4wc,cluster == unique (pathway4wc $ cluster)[i]) %>% filter (p.adjust < 0.05 ) %>% # 选择p.adjust最小的前5个通路 top_n ( - 2 ,p.adjust) %>% pull (p.adjust) # 将当前聚类的p.adjust值添加到向量中 myfontsize <- c (myfontsize,temfs) } # 加载scales包，用于数据缩放 library (scales) # 对p.adjust取负对数，然后进行缩放，使字体大小范围在0到10之间 myfontsize <- - log2 (myfontsize) %>% scales :: rescale (., to= c ( 0 : 1 )) * 10 # 创建行注释对象 # 使用anno_textbox函数绘制文本框注释 max_line_length <- 30 # 定义换行函数，保留 list 结构 wrap_list_text <- function (lst, width = max_line_length) { lapply (lst, function (vec) { sapply (vec, function (txt) { paste ( strwrap (txt, width = width), collapse = " \n " ) }) }) } # 应用函数 mytext_wrapped <- wrap_list_text (mytext) right_ann <- rowAnnotation ( textbox = anno_textbox (top20gene $ cluster, mytext_wrapped, # 允许文本换行 add_new_line = T, by = "anno_block" , # 设置文本框的样式 gp = gpar ( col = c ( pal_npg ( 10 ), pal_nejm ( 8 ), pal_aaas ( 4 )) , fontsize= myfontsize), # 设置文本框的圆角 round_corners = TRUE , # 设置文本框的最大宽度 max_width= unit ( 3 , 'mm' ))) col_fun <- colorRamp2 ( quantile (expr4plot) %>% as.numeric , c ( ' , ' , 'white' , ' , ' )) finalhmp <- Heatmap (expr4plot[top20gene $ gene,], show_column_names = F, show_row_names = F, col = col_fun, top_annotation = top_anno, cluster_columns = F, row_split = top20gene $ cluster, row_title = NULL , right_annotation = right_ann, border = T, left_annotation = left_ann ) draw (finalhmp) save.image ( "result_slingshot/slingshot.Rdata" ) Sys.time # [1] "2025-04-28 15:15:31 CST"