Slingshot 拟时序分析学习手册
⚠️ 本页用于教学学习，不构成临床诊疗建议。请结合实际数据和论文原文审慎使用。

拟时序分析（Pseudotime Analysis） 是一种在单细胞转录组学中常用的方法，它通过对细胞在特定生物过程中的转录组数据进行排序，推断其发展或变化的时间顺序。尽管单细胞数据通常是静态的（在一个时间点采集），拟时序分析通过基因表达的变化模式来重建细胞的发育轨迹或动态过程，从而模拟细胞随时间推移的变化，揭示细胞在不同状态下的演变路径和关键转折点。

拟时序分析作用 推断细胞的发育阶段：通过分析基因表达变化，将处于不同分化阶段的细胞排序，模拟出一个拟时间轴。识别关键的基因调控事件：通过比较不同拟时间阶段的细胞，识别在特定转折点上调或下调的基因，找出与细胞分化或功能转换相关的关键分子。推断细胞状态转换：分析细胞从一种状态向另一种状态过渡的动态变化，尤其在干细胞分化或疾病模型中应用广泛。常用工具 在单细胞转录组分析中，使用R语言进行拟时序分析时，有多种工具包可供选择。

以下是一些常用的R语言拟时序分析工具包，以及它们的特点和使用方法：1. Monocle是单细胞RNA-Seq数据分析的流行工具，能够进行细胞聚类、分类、计数以及构建单细胞轨迹。它使用一种称为“反向图嵌入”的方法来重建复杂的单细胞轨迹。Monocle 3是该工具的最新版本，提供了改进的算法和更多的功能。

https://bioconductor.org/packages/release/bioc/html/monocle.html 2. Slingshot是一个用于单细胞数据的拟时序分析工具，它通过聚类和轨迹推断来揭示细胞群体之间的关系和潜在的发育路径。它使用最小生成树（MST）来识别全局谱系结构，并构建平滑的谱系轨迹。Slingshot能够处理多个轨迹，并允许用户指定起点和终点，以确保分析结果与已知的生物学信息一致。

https://www.bioconductor.org/packages/release/bioc/html/slingshot.html 3. TSCAN（Time-series Clustering Analysis for Single-cell and bulk RNA-seq data）是一个专门用于单细胞和批量RNA-seq数据的时间序列聚类分析工具包。

它能够准确识别基因表达在时间上的变化模式，并将样本或细胞划分为不同的聚类。TSCAN支持伪时间分析，有助于理解细胞类型的演变和发展轨迹。https://github.com/zji90/TSCAN

二、Slingshot 拟时序分析方法学 【 参考文献 】：https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-018-4772-0 2.1 研究背景与目标：单细胞转录组技术能够解析复杂细胞群体的异质性，尤其在干细胞分化、发育轨迹等研究中具有重要价值。

然而，准确推断多分支细胞谱系（如发育过程中的分化路径）和伪时间（pseudotime）（反映细胞分化进程的连续变量）仍是计算生物学中的挑战。传统方法如Monocle、TSCAN等在处理多分支或噪声数据时存在局限性。Slingshot是一种从单细胞基因表达数据推断细胞谱系和伪时间的新方法。Slingshot被证明可以正确识别一到三个分支轨迹的生物信号。该模拟研究表明，Slingshot推断出的伪时间比其他领先方法更准确。

2.2 Slingshot方法的核心创新 Slingshot是一种两阶段算法，旨在从单细胞转录组数据中推断细胞谱系和伪时间：谱系结构推断：基于细胞聚类构建最小生成树（MST），识别全局谱系分支结构。用户可半监督地指定起始或终止聚类，提升生物学解释性。伪时间估计：通过 同步主曲线（simultaneous principal curves）拟合每个谱系的平滑轨迹，将高维数据映射为一维伪时间变量。该方法在噪声数据中表现稳健，支持多分支轨迹。

2.2.1 谱系结构推断（Lineage Identification） 1. 输入与预处理 输入数据：细胞表达矩阵：n个细胞 × J个基因的标准化表达矩阵（如log归一化后的读数值）。降维矩阵（推荐）：n×J’ 的降维坐标（如PCA、t-SNE、UMAP），J’ << J，避免高维空间距离失效。细胞聚类标签：通过k-means、GMM等算法获得的K个细胞簇（每个簇对应一种细胞状态）。

2. 最小生成树（MST）构建 目标：将K个细胞簇视为节点，构建全局分化路径。关键步骤：距离度量：采用考虑簇形状的马氏距离（Mahalanobis-like distance）：软聚类支持：若聚类结果为概率分布（如GMM），加权计算均值与协方差矩阵。MST生成：基于上述距离构建最小生成树，默认以最大直径簇为根节点，或允许用户指定根节点（半监督模式）。

3. 生物学监督机制 初始状态指定：用户指定根节点（如干细胞簇），MST从根节点出发生成所有可能路径。终末状态约束（可选）：用户指定终末簇（如成熟细胞簇），Slingshot将终末簇连接到最近的非终末簇，优化路径生物学合理性。分支检测：MST中每个分支节点对应一个分化事件（如二叉、三叉分支），形成多谱系结构。

2.2.2 伪时间估计（Pseudotime Inference） 1. 同步主曲线（Simultaneous Principal Curves） 目标：为每个谱系拟合平滑曲线，将细胞投影到曲线上获得伪时间。技术流程：初始化：基于谱系的MST路径（簇中心连线）定义初始曲线，而非传统主成分（加速收敛）。迭代优化。投影与伪时间计算：将细胞投影到当前曲线上，以弧长作为伪时间（起点为0）。

曲线更新：单谱系：对每个降维轴，用平滑样条函数（smoothing spline）拟合伪时间与坐标的关系，更新曲线形状。多谱系同步拟合：引入收缩（shrinkage）机制，使共享路径的曲线向平均曲线收缩，确保伪时间一致性。收敛判断：当细胞投影距离的平方和变化小于阈值时停止迭代。2. 多谱系伪时间分配 细胞权重：根据细胞到各谱系曲线的投影距离计算权重，权重高的谱系主导该细胞的伪时间。

伪时间一致性：共享路径细胞的伪时间在不同谱系中高度一致（收缩机制保证），分支后逐渐分化。2.2.3 实验验证与参数设置 Slingshot工具被应用于三个单细胞RNA-Seq数据集（HSMM、qNSC和OE），以验证其在谱系推断和伪时间推断方面的性能。这些数据集涵盖了不同的生物学背景和细胞类型，为评估Slingshot的通用性和准确性提供了丰富的测试场景。

1. HSMM数据集包含212个人类骨骼肌成肌细胞，用于研究其发育成成熟肌管的过程。该数据集是一个单一谱系的示例，分析中去除了一个污染的间质细胞群。研究中使用了Monocle生成的聚类标签以及通过k均值聚类获得的标签，并将数据表示为通过ICA获得的二维数据。Slingshot在该数据集上的应用展示了其对单一谱系数据的处理能力，能够准确地推断出细胞的发育轨迹。

2. qNSC数据集包含101个成年小鼠海马静止神经干细胞及其直系后代，这些细胞参与神经发生。该数据集的分析重点是细胞异质性和连续时间发育动态。Slingshot在该数据集中不仅关注了主要的IPC谱系，还描述了两种可能的细胞命运的发育轨迹。这表明Slingshot能够处理复杂的谱系结构，并在存在多个分支时提供全面的分析。3. OE数据集包含616个成年小鼠嗅上皮细胞，追踪静止干细胞发育成三种独特的终末细胞命运。

该数据集的分析展示了Slingshot在处理多谱系数据时的能力，能够准确地识别出不同的发育路径。2.2.4 模拟研究 为了进一步检验Slingshot的性能，研究者使用Bioconductor R包splatter生成了人工单细胞RNA测序数据集。模拟研究分为两部分：第一部分包含两个分支谱系，细胞数在120到1500之间变化；第二部分包含五个分支谱系，细胞数在220到1320之间变化。

通过调整基因在路径上的差异表达概率，研究者模拟了不同的信噪比场景，以评估Slingshot在各种条件下的表现。在聚类方法的选择上，Slingshot展示了对层次聚类、k均值聚类和高斯混合模型的稳健性。此外，研究者还通过Kendall等级相关系数评估了Slingshot与其他谱系推断方法的性能，结果表明Slingshot在推断伪时间与真实伪时间之间的一致性方面表现出色。

2.3 文章结果介绍 2.3.1 在真实数据集上的噪声鲁棒性：首先使用人类骨骼肌成肌细胞（HSMM）数据集的子集对一些知名方法的稳定性进行了检验。其中，Monocle 方法在单个细胞上构建 MST，并根据 PQ 树沿着 MST 的最长路径对细胞进行排序，是所比较的方法中最不稳定的。其绘制的路径高度可变，甚至对少量噪声也很敏感。相比之下，基于簇的 MST 方法和主曲线方法在类似自举的样本上表现出稳定性。

基于簇的 MST 由于分段线性路径的顶点，多个细胞通常会被分配相同的伪时间，对应于顶点的值。主曲线方法是最稳定的方法，但在更复杂的数据集上，它有明显的局限性，即只能表征单个谱系。这也是将主曲线扩展以适应多个分支谱系的原因。注：图a展示了三种不同的谱系推断方法在原始数据集和50个数据子样本上的细胞排序情况：Monocle方法（红色）：通过在所有细胞上构建最小生成树（MST），识别最长路径来推断细胞排序。

Waterfall和TSCAN方法（紫色）：通过聚类细胞并用MST连接聚类中心来推断细胞排序。这里聚类是通过k-means方法进行的，k=5。Embeddr和Slingshot方法（绿色）：使用主曲线方法，即通过非线性拟合数据来推断细胞排序。所有方法都使用了独立成分分析（ICA）进行降维处理。图b展示了基于50个数据子样本的伪时间与原始数据集伪时间的散点图。从图中可以看出：左图：红色表示Monocle方法的伪时间比较。

可以看到，伪时间在原始数据和子样本之间的相关性较差，说明Monocle方法在子样本上的稳定性较差。中图：紫色表示Waterfall和TSCAN方法的伪时间比较。这些方法的伪时间在原始数据和子样本之间的相关性较好，说明这些方法在子样本上的稳定性较好。右图：绿色表示Embeddr和Slingshot方法的伪时间比较。这些方法的伪时间在原始数据和子样本之间的相关性最好，说明这些方法在子样本上的稳定性最好。

2.3.2 多谱系推理 在分析嗅觉上皮（OE）数据时，Slingshot利用局部监督标记成熟的支持细胞（mSus）、微绒毛细胞（MV）和成熟的嗅觉感觉神经元（mOSN）作为终端状态，尽管只有第一个对最终的基于簇的MST有影响。Slingshot推断出的谱系结构建立了两个分支的顺序，并随后得到了验证。具体来说，证明了支持细胞是通过水平基底细胞（HBC）的直接转换产生的，而微绒毛细胞和神经元细胞则需要一个中间的、增殖状态。

Slingshot还允许使用一种形状敏感的距离度量，这种度量受到马氏距离的启发，它根据两个簇的协方差结构来缩放簇中心之间的距离。这与Waterfall和TSCAN中使用的基于簇中心的标准欧几里得距离不同，后者未能利用簇的形状，导致在HBC分化的早期错误地识别了一个虚假的分支事件。而Slingshot能够正确识别与先前生物学知识一致的谱系，而其他谱系检测方法则没有。

例如，Monocle 2只识别出两个谱系，其中一个以球状基底细胞（GBC）为终点，这是一个已知的过渡状态，而两个谱系都包含支持细胞和微绒毛细胞，这些是已知的独立谱系的终点。注：嗅觉上皮（OE）数据集。在三谱系 OE 数据集上，Slingshot 和 Monocle 2 推断出的每个谱系的伪时间变量。图（a）：细胞类型之间已知的生物学关系。

图（b）：对于 Monocle 2，我们使用 DDRTree 算法获得数据的二维表示，并根据 HBC 细胞簇中细胞占比最高的情况选择起始状态。图（c）：对于 Slingshot，我们使用前五个主成分（PCs）并通过 RSEC 对细胞进行聚类。将 HBC 细胞簇指定为起点，mSus 细胞簇指定为一个终点；其他终点在无监督的情况下确定。

2.3.3 上游计算选择的稳健性 数据集模拟：为了对不同的谱系推断方法进行更定量的比较，并检验Slingshot对上游计算选择的鲁棒性，研究者们使用Bioconductor R包splatter生成的合成数据集进行了一项模拟研究。在研究的第一部分，所有模拟的数据集都包含一个初始路径，该路径分化为两个不同的谱系（见图4a）。在研究的第二部分，每个数据集都是从一个更复杂的分支结构模拟出来的，包含五个不同的谱系（见图4b）。

对于两谱系部分的模拟研究，生成了1200个合成数据集，而对于五谱系部分，模拟了300个数据集。推断的伪时间准确性的测量方法如下：对于每个真实的谱系，根据真实和推断的伪时间变量之间的Kendall秩相关系数，识别所有推断的谱系中的最佳匹配。将这些值平均到所有真实谱系上，就得到了特定方法在特定数据集上的准确性得分。与标准的皮尔逊相关系数一样，Kendall秩相关系数的值在-1到1之间，值越接近1表示推断的伪时间和真实伪时间之间的一致性越好。

准确性测量：依据真实和推断的伪时间变量之间的肯德尔等级相关系数，为每个真实谱系在所有推断谱系中找到最佳匹配。对所有真实谱系的这些值求平均，得到特定方法在特定数据集上的准确度分数，肯德尔等级相关系数值在 - 1 到 1 之间，越接近 1 表明推断和真实伪时间之间一致性越好。注：图a展示了两谱系情况下的分支结构。图b展示了五谱系情况下的更复杂分支结构。图c和图d：这两个图展示了不同方法在两谱系和五谱系数据集上的准确性得分分布。

每个密度图旁边的条形图表示在该方法返回错误数据集的百分比。以下是对这些结果的解释：两谱系情况 Monocle：大多数策略表现良好，通常产生双峰分布的准确性得分，一个峰值接近0，另一个较大峰值在0.5或以上。然而，Monocle比其他任何方法更频繁地返回错误。Monocle 2：比Monocle更一致，但总体准确性较低。它很少返回接近0的得分，并且显示出较少的双峰性，特别是在使用四维或五维RGE时。

总体准确性得分较低可能部分是由于它识别出的大量虚假分支事件。Diffusion Pseudotime (DPT)：由于大量细胞缺少分支分配，导致准确性得分人为降低。在两谱系和五谱系情况下，不使用分支信息的DPT-1策略达到了最高的中位数准确性得分。TSCAN：使用全分位数标准化的TSCAN产生的准确性得分与Monocle 2相当。当运行推荐的预处理步骤时，TSCAN的表现稍差，特别是在没有全分位数标准化的情况下。

混合方法（使用TSCAN进行降维和聚类，然后使用Slingshot进行伪时间推断）产生了最高的中位数准确性得分。五谱系情况 Monocle和Monocle 2：继续识别大量虚假谱系，并且样本大小与其推断的谱系数量之间仍然存在强烈的相关性。Slingshot和TSCAN：在所有谱系类型中通常表现优于其他方法。Slingshot在两谱系和五谱系数据集上的表现相对稳定，准确性得分分布较窄，表明其在不同数据集上的表现较为一致。

这些结果表明，Slingshot和TSCAN可能是处理复杂谱系结构时更可靠的谱系推断工具。

2.3.4 对聚类（方法）的鲁棒性 该研究通过模拟数据集验证了Slingshot算法在单细胞轨迹推断中的核心优势与局限性，其核心结论可总结为以下三点：1. 聚类方法鲁棒性 Slingshot通过”同步主曲线（simultaneous principal curves）“对基于聚类的MST（最小生成树）进行平滑处理，使其结果对具体聚类方法（层次聚类/k-means/高斯混合模型）的选择具有较强鲁棒性。

相比之下，TSCAN等直接将细胞投影到聚类中心构建MST的方法，因依赖聚类中心位置而结果波动较大。

2. 聚类数量K的关键影响 实验表明，Slingshot的表现更敏感于 聚类数量K 而非具体聚类方法：K过小：可能导致分支事件漏检 K过大：可能引入虚假分支 适度K值：即使降维质量不同（如4D-PCA与3D-PCA），仍能保持结果稳定性 3. 方法局限性 虽然主曲线平滑增强了稳定性，但本质上仍属于”基于聚类的MST方法”，无法完全规避此类方法的固有缺陷：正确识别全局谱系结构的前提是初始聚类能近似反映真实发育轨迹 极端K值导致的过分割/欠分割问题与其他同类方法存在共性 研究通过两系谱拓扑模拟数据验证，当初始聚类能大致捕捉发育轨迹时，Slingshot通过主曲线实现了对局部聚类波动的强容错能力，这使其在单细胞轨迹推断中具有更优的实用价值。

注：在模拟的双谱系数据集上，Slingshot的伪时间推断对不同的聚类方法表现出良好的稳健性。研究使用了分层聚类、k-means和高斯混合建模（GMM）等不同的聚类方法，并在不同的聚类数量（K）下进行了测试。结果显示，Slingshot在不同的聚类方法和K值下均能生成相似的准确性得分分布。当K=3时，Slingshot通常无法检测到分支事件，导致生成的伪时间与真实谱系不完全匹配。

随着聚类数量的增加，Slingshot能够更准确地推断出谱系结构和伪时间。当K值较高时，准确性得分开始缓慢下降，这可能是由于Slingshot开始过度拟合并识别出更多的虚假分支事件。此外，研究还发现，使用四维PCA进行降维时，Slingshot能够产生更一致的高准确性得分，而三维PCA则导致得分高度可变。这表明在使用Slingshot时，选择合适的降维方法和聚类数量对于获得准确的伪时间推断至关重要。

2.4 文章讨论 Slingshot方法总结 1. 方法概述：Slingshot通过将谱系推断问题分为两个步骤来解决：首先是全局谱系结构的识别，其次是每个谱系中细胞的伪时间推断。这种方法结合了高度稳定的最小生成树（MST）和灵活的多分支主曲线拟合技术，能够在噪声数据中稳定地识别复杂的谱系结构。全局谱系结构识别：通过在细胞簇上构建最小生成树（MST），Slingshot能够识别谱系的数量和分支位置。

伪时间推断：使用同时主曲线方法拟合每个谱系的平滑曲线，从而为每个细胞推断出伪时间。2. 性能表现：准确性：在模拟数据集和真实数据集上，Slingshot能够准确推断出单个谱系和多分支谱系的结构。鲁棒性：Slingshot对噪声和不同的聚类方法表现出较高的鲁棒性，即使在不同的降维技术和聚类数量下，也能保持稳定的伪时间推断结果。

灵活性：Slingshot允许用户在分析中加入局部监督，例如指定起始簇和终端簇，从而在不牺牲灵活性的情况下提高结果的生物学一致性。3. 应用场景：不同数据类型：Slingshot适用于多种单细胞数据类型，包括RNA测序数据和质谱细胞术数据。现有分析管道集成：该方法设计灵活，能够与现有的单细胞数据分析管道无缝集成。4. 优势总结：Slingshot在处理单细胞基因组学数据的谱系和伪时间推断中表现出色，特别是在处理复杂分支结构和噪声数据时。

它提供了一种灵活且鲁棒的方法，能够准确推断细胞的发育轨迹，并且可以与多种降维和聚类技术结合使用。5. 实验结果：在模拟数据集和真实数据集上的实验结果表明，Slingshot能够准确识别复杂的多分支谱系结构，并且在不同噪声水平和细胞数量下表现出较高的鲁棒性。例如，在嗅觉上皮数据集中，Slingshot能够正确识别复杂的三谱系结构，而其他方法则无法做到。Slingshot是一种强大且灵活的工具，适用于单细胞基因组学数据的谱系和伪时间推断。

它结合了处理噪声数据所需的稳定性以及识别复杂结构的灵活性，能够为动态基因表达分析提供关键的谱系推断步骤。

三、官方教学解读 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 )

四、外部数据实操 接下来，我们从GEO数据库下载真实单细胞转录组测序数据，并使用R语言Slingshot包进行单细胞转录组测序 拟时序分析；全流程包括：原始数据下载———测序数据预处理———降维聚类———细胞注释———利用关键细胞进行拟时序分析。4.1 单细胞数据下载：GSE163973 该数据集主要探讨了皮肤纤维化疾病中的成纤维细胞异质性，特别是以瘢痕疙瘩（KD）为研究模型。

该研究使用真皮进行 scRNA-seq 分析，因为瘢痕疙瘩代表皮肤真皮纤维化疾病。经过严格的质量控制，我们获得了 40655 个细胞的转录组（瘢痕疙瘩（nsample=3）：21488;正常疤痕（nsample=3）：19167）。我们下载该数据集中每个样本（其中包含3个KD样本，3个NC样本）的原始文件，即barcodes.tsv、genes.tsv和matrix.mtx文件。

4.2 测序数据预处理：接下来，我们使用R语言加载数据并构建疾病组与正常组样本的Seurat对象：4.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) 4.2.2 处理各组样本 rm ( list= ls ) normal_scar <- Read10X ( data.dir = "./00_rawdata/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 ( "./00_rawdata/" ,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 = "normal_Object.rds" ) rm ( list= ls ) keloid_data <- Read10X ( data.dir = "./00_rawdata/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 ( "./00_rawdata/" ,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 = "keloid_Object.rds" ) 4.2.3 合并两组数据，构建seurat_object rm ( list= ls ) normal_object <- readRDS ( "../00_rawdata/normal_Object.rds" ) keloid_object <- readRDS ( "../00_rawdata/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 ) 4.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 ) variable_genes <- 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) seurat_object <- JackStraw (seurat_object, num.replicate = 100 ) seurat_object <- ScoreJackStraw (seurat_object, dims = 1 : 20 ) PCA_Score <- JackStrawPlot (seurat_object, dims = 1 : 20 ) p <- DimPlot (seurat_object, reduction = 'pca' , group.by = 'orig.ident' ) 4.3 细胞聚类 seurat_object <- RunUMAP (seurat_object, reduction = "pca" , dims = 1 : 20 ) seurat_object <- FindNeighbors (seurat_object, reduction = "pca" , dims = 1 : 20 ) seurat_object <- FindClusters (seurat_object, resolution = 0.3 ) head ( Idents (seurat_object), 5 ) plot <- DimPlot (seurat_object, reduction = "umap" , label = TRUE , group.by = "RNA_snn_res.0.3" ) 4.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" ) 4.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 seurat_object <- RenameIdents (seurat_object,new.cluster.ids) seurat_obj_Annotated <- seurat_object saveRDS (seurat_obj_Annotated, file = "seurat_obj_Annotated.rds" ) 4.3.3 可视化注释结果 rm ( list= ls ) seurat_object <- readRDS ( "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) plot <- DimPlot (seurat_object, reduction = "umap" , label = TRUE ) 4.3.4 显著表达基因可视化 mesenchymal_stem_cell <- c ( "SELE" , "DCN" , "LGALS7" , "CD86" , "TPSAB1" , "PTPRC" , "SOX10" , "DCD" ) plot <- VlnPlot (seurat_object, features = mesenchymal_stem_cell, ncol = 2 ) image 4.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 = "subset_clusters.rds" ) 4.5 拟时序分析 在这里，我们基于上面筛选的关键细胞簇进行后续拟时序分析：4.5.1 载入相关包 rm ( list= ls ) library (Seurat) library (SingleCellExperiment) library (tidyverse) library (RColorBrewer) 4.5.2 可视化细胞簇 subset_clusters <- readRDS ( "subset_clusters.rds" ) plot <- DimPlot (subset_clusters, reduction = "umap" , label = TRUE ) mesenchymal_stem_cell <- c ( "PECAM1" , "CDH5" , "VWF" , "THY1" , "TAGLN" , "COL1A2" , "DCN" , "KRT5" , "KRT14" , "KRT1" ) plot <- DotPlot (subset_clusters, features = mesenchymal_stem_cell) + RotatedAxis + scale_color_gradientn ( colours = c ( ' , ' , ' , ' )) seurat_object <- subset_clusters 4.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 colnames ( colData (select_cluster)) 结果展示：> 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" "RNA_snn_res.0.3" [6] "seurat_clusters" "cell_type" "group" "ident" > 4.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”。可以使用不同的距离度量方法。

4.5.5 可视化计算结果 colnames ( colData (sim)) SlingshotDataSet (sim) 结果展示：> colnames(colData(sim)) [1] "orig.ident" "nCount_RNA" "nFeature_RNA" "percent.mt" "RNA_snn_res.0.3" [6] "seurat_clusters" "cell_type" "group" "ident" "slingshot" [11] "slingPseudotime_1" > SlingshotDataSet(sim) class: SlingshotDataSet Samples Dimensions 38382 2 lineages: 1 Lineage1: Epidermal stem cells Endothelial cells Fibroblast cells curves: 1 Curve1: Length: 48.537 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”分化。

4.6 探究分化过程中的基因表达 4.6.1 tradeSeq 提供了一种灵活的回归模型拟合方法，可用于查找沿轨迹中的一个或多个谱系差异表达的基因。基于拟合模型，发现表达与伪时间相关的基因，或沿轨迹差异表达（在特定区域）的基因。它为每个基因拟合一个负二项式广义加性模型 （GAM），并对 GAM 的参数进行推理。

if ( ! require ( "BiocManager" , quietly = TRUE )) install.packages ( "BiocManager" ) BiocManager :: install ( "tradeSeq" ) 4.6.2 数据准备 counts <- sim @ assays @ data $ counts crv <- SlingshotDataSet (sim) counts：从模拟数据对象sim中提取基因表达矩阵（计数数据）。

sim@assays @data$counts提取计数矩阵，每行代表一个基因，每列代表一个细胞。crv：将sim转换为SlingshotDataSet对象。Slingshot是一种分析单细胞发育轨迹的方法，SlingshotDataSet对象描述了细胞间的发育轨迹及其分支结构。

4.6.3 确定最佳模型的参数 k set.seed ( 1234 ) icMat <- evaluateK ( counts = counts, sds = crv, k = 3 : 5 , nGenes = 20 , 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 4.6.4 计算伪时间 pseudotime <- slingPseudotime (crv, na = FALSE ) slingPseudotime(crv, na = FALSE)：计算每个细胞的伪时间值。

crv 是一个 SlingshotDataSet 对象，包含了细胞的发育轨迹和分支结构。伪时间（Pseudotime） 是一种用于量化细胞在发育过程中所处“时间”点的概念。尽管我们没有真实的时间轴，伪时间可以帮助揭示细胞的发育或分化状态。该函数计算每个细胞沿轨迹的相对位置（伪时间），并将这些值存储在 pseudotime 中。na = FALSE：表示如果在计算过程中有缺失值，是否应跳过这些细胞。

设为FALSE时，缺失值将被处理（默认行为）。4.6.5 计算细胞权重 cellWeights <- slingCurveWeights (crv) slingCurveWeights(crv)：计算每个细胞在各个轨迹曲线上的权重。在发育轨迹分析中，细胞的表达模式可以影响它们在特定轨迹上的位置。该函数计算每个细胞在各个发育轨迹（或曲线）上的“权重”，即它们在不同轨迹分支的概率分配。

4.6.6 拟合广义加性模型（GAM） top_genes <- head ( order ( rowVars (counts), decreasing = TRUE ), 200 ) counts_subset <- counts[top_genes, ] system.time ({ sce <- fitGAM ( counts = counts_subset, pseudotime = pseudotime, cellWeights = cellWeights, nknots = 6 , verbose = FALSE ) }) fitGAM：拟合一个广义加性模型（GAM）来描述基因在伪时间中的表达变化。

counts：基因表达矩阵，其中行是基因，列是细胞。pseudotime：之前计算的每个细胞的伪时间值。它提供了每个细胞在发育轨迹中的相对时间点，作为GAM拟合的一个自变量。cellWeights：每个细胞在不同轨迹曲线上的权重，作为GAM拟合中的额外信息，用于提升模型的准确性。nknots = 6：指定在拟合GAM时使用的“结点”数量。结点数量影响模型的灵活性，通常越多结点，模型越灵活，但也可能导致过拟合。这里指定了6个结点。

4.6.7 基因表达与拟时序的相关性 assoRes <- associationTest (sce) head (assoRes) > head(assoRes) waldStat df pvalue meanLogFC COL1A1 76617.452 2 0 2.5987026 MALAT1 6995.086 2 0 0.8395357 COL1A2 83287.236 2 0 2.6722302 COL3A1 65739.148 2 0 2.6731177 RPLP1 41623.380 2 0 0.8002227 RPS18 25055.523 2 0 0.5648929 > 4.6.8 起止点相关性最高的基因 startRes <- startVsEndTest (sce) head (startRes) > head (startRes) waldStat df pvalue logFClineage1 COL1A1 58068.684 1 0 6.8558583 MALAT1 4258.177 1 0 0.7492381 COL1A2 63639.221 1 0 6.8190368 COL3A1 50532.819 1 0 6.7238449 RPLP1 22381.552 1 0 - 0.6871365 RPS18 22745.007 1 0 - 0.6963100 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）。这是基因在起始点和终止点之间表达变化的对数变动。

4.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) 4.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) 结果展示：> head (plt_data) celltype plotcol Pseudotime DCN CD63 COL1A2 CCDC80 COL6A2 AAACCCAAGAACAGGA - 1_1_1 Endothelial cells KD 17.6562951 0.000000 2.807355 1.000000 0.000000 0.000000 AAACCCAAGATCCTAC - 1_1_1 Fibroblast cells KD 44.1862229 8.266787 5.754888 8.471675 6.672425 4.954196 AAACCCAAGCCTGACC - 1_1_1 Endothelial cells KD 18.1090814 0.000000 3.459432 0.000000 0.000000 0.000000 AAACCCAAGCGCAATG - 1_1_1 Endothelial cells KD 21.4856377 0.000000 1.584963 0.000000 0.000000 0.000000 AAACCCACAAACTCGT - 1_1_1 Endothelial cells KD 21.1379977 0.000000 1.584963 1.000000 0.000000 0.000000 AAACCCACAAGACGGT - 1_1_1 Epidermal stem cells KD 0.9835969 0.000000 2.000000 0.000000 0.000000 0.000000 CST3 AAACCCAAGAACAGGA - 1_1_1 3.321928 AAACCCAAGATCCTAC - 1_1_1 6.714246 AAACCCAAGCCTGACC - 1_1_1 3.906891 AAACCCAAGCGCAATG - 1_1_1 4.169925 AAACCCACAAACTCGT - 1_1_1 2.584963 AAACCCACAAGACGGT - 1_1_1 2.321928 > colnames (plt_data) [ 1 ] "celltype" "plotcol" "Pseudotime" "DCN" "CD63" "COL1A2" "CCDC80" "COL6A2" [ 9 ] "CST3" > 4.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 } library (patchwork) wrap_plots (plt_list) ggsave ( filename = 'top6_genes.pdf' , width = 9 , height = 5 )

五、结论 本期教学，我们讲解并带领大家完成了单细胞转录组测序的全流程，从数据预处理、质量控制，到基因表达定量分析、聚类分析，再到使用 Slingshot 进行拟时序分析，了解了如何运用 Slingshot 工具进行细胞发育轨迹的推断，探讨了基因在不同发育阶段中的表达动态。拟时序分析是单细胞转录组学中的一项重要技术，它允许研究人员在没有时间序列数据的情况下推断细胞的动态变化。

Slingshot 作为一种强大的工具，通过结合聚类和轨迹推断，为理解复杂的细胞发育过程提供了新的视角。通过 Slingshot，研究人员可以揭示细胞状态的连续变化，识别关键的基因调控事件，并推断细胞状态转换的动态过程。在未来的研究中，拟时序分析将会在更多的生物医学领域中得到应用，包括癌症研究、免疫反应、组织再生等方面，推动精准医疗的发展。

希望通过本期教学，大家不仅掌握了单细胞转录组学分析的基本方法，也能够运用这些技术探索细胞发育和疾病发生的奥秘。我们鼓励学员们继续深入学习，并将这些技能应用到自己的科研工作中，推动学术进步与创新。

本教程用到的环境 大家如果遇到报错，请先检查一下各包的版本是否与我一致：sessionInfo R version 4.3.2 (2023-10-31 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045) Matrix products: default locale: [1] LC_COLLATE=Chinese (Simplified)_China.utf8 LC_CTYPE=Chinese (Simplified)_China.utf8 [3] LC_MONETARY=Chinese (Simplified)_China.utf8 LC_NUMERIC=C [5] LC_TIME=Chinese (Simplified)_China.utf8 time zone: Asia/Shanghai tzcode source: internal attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets methods base other attached packages: [1] paletteer_1.6.0 clustree_0.5.1 ggraph_2.2.1 [4] openxlsx_4.2.7.1 hdf5r_1.3.11 scater_1.30.1 [7] scran_1.30.2 scuttle_1.12.0 patchwork_1.3.0 [10] ggpubr_0.6.0 ggsci_3.2.0 future_1.34.0 [13] Matrix_1.6-5 tradeSeq_1.16.0 slingshot_2.10.0 [16] TrajectoryUtils_1.10.1 princurve_2.1.6 RColorBrewer_1.1-3 [19] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1 [22] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 [25] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 [28] tidyverse_2.0.0 SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0 [31] Biobase_2.62.0 GenomicRanges_1.54.1 GenomeInfoDb_1.38.8 [34] IRanges_2.36.0 S4Vectors_0.40.2 BiocGenerics_0.48.1 [37] MatrixGenerics_1.14.0 matrixStats_1.4.1 Seurat_5.1.0 [40] SeuratObject_5.0.2 sp_2.1-4 loaded via a namespace (and not attached): [1] spatstat.sparse_3.1-0 bitops_1.0-9 httr_1.4.7 [4] tools_4.3.2 sctransform_0.4.1 backports_1.5.0 [7] R6_2.5.1 lazyeval_0.2.2 uwot_0.2.2 [10] mgcv_1.9-0 withr_3.0.2 gridExtra_2.3 [13] progressr_0.15.1 cli_3.6.3 textshaping_0.4.1 [16] spatstat.explore_3.3-3 fastDummies_1.7.4 labeling_0.4.3 [19] spatstat.data_3.1-4 ggridges_0.5.6 pbapply_1.7-2 [22] systemfonts_1.1.0 parallelly_1.40.1 limma_3.58.1 [25] rstudioapi_0.17.1 generics_0.1.3 ica_1.0-3 [28] spatstat.random_3.3-2 car_3.1-3 zip_2.3.1 [31] ggbeeswarm_0.7.2 abind_1.4-8 lifecycle_1.0.4 [34] yaml_2.3.10 edgeR_4.0.16 carData_3.0-5 [37] SparseArray_1.2.4 Rtsne_0.17 grid_4.3.2 [40] promises_1.3.2 dqrng_0.4.1 crayon_1.5.3 [43] miniUI_0.1.1.1 lattice_0.21-9 beachmat_2.18.1 [46] cowplot_1.1.3 pillar_1.10.0 knitr_1.49 [49] metapod_1.10.1 future.apply_1.11.3 codetools_0.2-19 [52] leiden_0.4.3.1 glue_1.8.0 spatstat.univar_3.1-1 [55] data.table_1.16.4 vctrs_0.6.5 png_0.1-8 [58] spam_2.11-0 gtable_0.3.6 rematch2_2.1.2 [61] cachem_1.1.0 xfun_0.49 S4Arrays_1.2.1 [64] mime_0.12 tidygraph_1.3.1 survival_3.5-7 [67] statmod_1.5.0 bluster_1.12.0 fitdistrplus_1.2-1 [70] ROCR_1.0-11 nlme_3.1-163 bit64_4.5.2 [73] RcppAnnoy_0.0.22 irlba_2.3.5.1 vipor_0.4.7 [76] KernSmooth_2.23-22 colorspace_2.1-1 tidyselect_1.2.1 [79] bit_4.5.0.1 compiler_4.3.2 BiocNeighbors_1.20.2 [82] DelayedArray_0.28.0 plotly_4.10.4 scales_1.3.0 [85] lmtest_0.9-40 digest_0.6.37 goftest_1.2-3 [88] spatstat.utils_3.1-1 rmarkdown_2.29 XVector_0.42.0 [91] htmltools_0.5.8.1 pkgconfig_2.0.3 sparseMatrixStats_1.14.0 [94] fastmap_1.2.0 rlang_1.1.4 htmlwidgets_1.6.4 [97] shiny_1.10.0 DelayedMatrixStats_1.24.0 farver_2.1.2 [100] zoo_1.8-12 jsonlite_1.8.9 BiocParallel_1.36.0 [103] BiocSingular_1.18.0 RCurl_1.98-1.16 magrittr_2.0.3 [106] Formula_1.2-5 GenomeInfoDbData_1.2.11 dotCall64_1.2 [109] munsell_0.5.1 Rcpp_1.0.13-1 viridis_0.6.5 [112] reticulate_1.40.0 stringi_1.8.4 zlibbioc_1.48.2 [115] MASS_7.3-60 plyr_1.8.9 listenv_0.9.1 [118] ggrepel_0.9.6 deldir_2.0-4 graphlayouts_1.2.1 [121] splines_4.3.2 tensor_1.5 hms_1.1.3 [124] locfit_1.5-9.10 igraph_2.1.2 spatstat.geom_3.3-4 [127] ggsignif_0.6.4 RcppHNSW_0.6.0 reshape2_1.4.4 [130] ScaledMatrix_1.10.0 evaluate_1.0.1 BiocManager_1.30.25 [133] tzdb_0.4.0 tweenr_2.0.3 httpuv_1.6.15 [136] RANN_2.6.2 polyclip_1.10-7 scattermore_1.2 [139] ggforce_0.4.2 rsvd_1.0.5 broom_1.0.7 [142] xtable_1.8-4 RSpectra_0.16-2 rstatix_0.7.2 [145] later_1.4.1 viridisLite_0.4.2 ragg_1.3.3 [148] memoise_2.0.1 beeswarm_0.4.0 cluster_2.1.4 [151] timechange_0.3.0 globals_0.16.3