单细胞转录组学轨迹分析解析2-Slingshot代码解析

通过上节对Slingshot文献的基本讲解，对这个拟时序的分析方法有了基本的了解，作者也公布了流程的代码，并分享在https://bioconductor.org/packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html上。

作者的文章发地很早（2018年），但是最近10月份（2023年10月）对代码进行了更新。

图片

1 Introduction

这个流程将展示一个完整的单细胞谱系分析工作流程，特别强调谱系重建和拟时序推理的过程。我们将利用Street等中提出的Slingshot(Street, Risso et al. 2018)，并对他们的应用内容进行展示。

Slingshot的目标是使用细胞亚群的结果对全部发育结果进行可视化，并将这种结构转换为由一维变量表示的平滑谱系，称为“拟时序”。我们提供的工具用于以无监督或半监督的方式学习聚类关系，并构建代表每个谱系的平滑曲线，以及每个步骤的可视化方法。

1.1 Overview

Slingshot分析需要两个输入文件：细胞的坐标矩阵和聚类标签向量。随后，

l使用 getLineages 函数在集群上构造最小生成树 （MST） 来识别全局谱系结构。

l使用 getCurves 函数拟合同步主曲线来构造平滑的谱系并推断伪时间变量。

l使用内置的可视化工具对结果进行可视化。

1.2 Datasets

我们将在此流程分析中使用两个test数据集。第一个（称为“单轨迹”数据集）在下面生成，旨在表示单个谱系，其中三分之一的基因与过渡相关。该数据集将包含在 SingleCellExperiment 对象 （Lun and Risso 2017） 中，并将用于演示完整的“从头到尾”工作流程。

Slingshot主要的分析矩阵是基于SingleCellExperiment的，如果大家一般都是用seurat矩阵的，可以先去学习一下SingleCellExperiment矩阵，便于后面的理解。

2 Upstream Analysis

2.1 Gene Filtering

上一节提到的Slingshot也可以用于前面的数据过滤等流程分析，作者也是提供了相应的代码。

这里为了进行数据过滤的目的是为了为了开始对单个谱系数据集进行分析，我们需要降低数据的维度，过滤没有信息的基因，从而提高下游分析的速度，将后期的分析错误率降到最低。

对于基因过滤步骤，我们保留了在至少足够多的细胞中稳定表达的任何基因以构成一个簇，使它们成为潜在的细胞类型标记基因。我们将最小簇大小设置为10个细胞，如果一个基因的模拟计数至少为 3 个reads，则将其定义为“robustly expressed”。

library(slingshot, quietly = FALSE)    
data("slingshotExample")# filter genes down to potential cell-type markers
# at least M (15) reads in at least N (15) cells
geneFilter <- apply(assays(sce)$counts,1,function(x){
    sum(x >= 3) >= 10
})
sce <- sce[geneFilter, 

2.2 Normalization

在大多数RNA-Seq分析流程中，选择归一化的方法是早期数据过滤中比较重要的步骤，主要是去除技术噪音，如批次效应、测序深度、细胞周期效应等。这使我们能够从数据中去除不需要的技术或生物伪影，例如批次、测序深度、细胞周期效应等。基于此我们推荐使用scone package 。

由于现在是用test数据进行后续的分析，因此我们不考虑批次效应的问题，直接对数据进行quantile normalization处理。

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.3 Dimensionality Reduction

Slingshot的基本假设是：转录水平类似的细胞在某个降维水平中距离较近。由于我们在构建谱系和测量伪时间时使用欧几里得距离，因此对数据进行低维表示非常重要。

我们将演示两种降维方法：主成分分析（PCA）和均匀流形近似和投影（UMAP，通过uwot包）。

在进行PCA时，我们不会根据基因的方差来缩放基因，主要是因为不是所有基因的表达量都相同。我们希望在强烈表达、高度可变的基因中找到信号。在绘图时，我们确保设置纵横比，以免扭曲计算距离。

pca <- prcomp(t(log1p(assays(sce)$norm)), scale. = FALSE)rd1 <- pca$x[,1:2]plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)

图片

library(uwot)## Loading required package: Matrix#### Attaching package: 'Matrix'## The following object is masked from 'package:S4Vectors':    ####     expandrd2 <- 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.4 Clustering Cells

Slingshot的最终输入是细胞的簇标签向量。如果未提供，则 slingshot 会将数据视为单个聚类并拟合标准主曲线。但是，我们建议即使在仅存在单个谱系的数据集中也对细胞进行聚类，因为它允许潜在地发现新的分支事件。

在此流程中确定的聚类将用于确定基础谱系的全局结构（即它们的数量、它们何时分支，以及这些分支事件的大致位置）。这与聚类单细胞数据的典型目标不同，后者是识别数据集中存在的所有生物学相关细胞类型。例如，在确定全局谱系结构时，无需区分未成熟神经元和成熟神经元，因为这两种细胞类型可能都属于谱系的同一部分。

在我们的分析中，我们实现了两种聚类方法，它们同样假设低维空间中的欧几里得距离反映了细胞之间的生物学差异：高斯混合建模和 k-means方法。前者在 mclust 包中实现，并具有一种基于贝叶斯信息准则（BIC）确定聚类数量的自动化方法。

library(mclust, quietly = TRUE)## Package 'mclust' version 6.0.0## Type 'citation("mclust")' for citing this R package in publications.#### Attaching package: 'mclust'## The following object is masked from 'package:mgcv':####     mvncl1 <- Mclust(rd1)$classificationcolData(sce)$GMM <- cl1library(RColorBrewer)plot(rd1, col = brewer.pal(9,"Set1")[cl1], pch=16, asp = 1)

图片

虽然 k-means 没有类似的功能，同步主曲线对 k 的选择非常robust，所以我们选择一个 k的 4 比较随意。如果这个数字太低，我们可能会错过一个真正谱系结果，如果它太高或有大量的小簇，我们可能会看到很多假的分支事件。

图片

总结

目前是对前面的数据过滤的内容进行解析，但是一般我们在做轨迹分析的时候基本上前面的数据过滤和降维基本上已经分析完了，因此主要是后面的轨迹分析，所以主要还是明天下一节的内容，我们主要是对这个软件前面的内容进行一下了解。

参考文献

Street, K., et al. (2018). "Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics." BMC Genomics 19(1).