作者：Diya Das ，Kelly Street ，Davide Risso 最后修改：6月28日。
8.1概述
8.1.1说明
本研讨会将作为实验室会议（简要介绍，然后是动手编码）提供，指导参与者在Bioconductor工作流程中分析单细胞RNA测序数据，分为三部分：
维数减少，可以解释零膨胀，过度分散和批量效应
采用基于重采样的方法的小区聚类，从而产生稳健且稳定的聚类
谱系轨迹分析揭示了连续的分支发育过程
我们将提供实验室会话的工作示例，以及此存储库中的一组独立注释。
注意：本次研讨会的前一版本在BioC 上得到了很好的参与，但所提供的工具已经针对互操作性进行了大量更新（最值得注意的是，通过使用SingleCellExperiment该类），我们收到了许多提供更新工作流程的请求。我们计划将此研讨会的反馈意见纳入我们F1000工作流程的修订版。
8.1.2先决条件
我们期望R语法的基本知识。熟悉S4对象可能会有所帮助，但不是必需的。更重要的是，参与者应熟悉RNA测序实验的概念和设计。不需要直接使用单细胞RNA-seq，并且将说明单细胞RNA-seq与大量RNA-seq相比的主要挑战。
8.1.3参与
这将是一个动手实践的研讨会，每个学生使用他们的笔记本电脑将分析提供的示例数据集。研讨会将是教师将向学生展示的示例代码（通过此存储库提供）和简短练习的混合。
8.1.4使用 R / Bioconductor封装
zinbwave：https：///packages/zinbwave
clusterExperiment：https：///packages/clusterExperiment
弹弓：https：///packages/slingshot
8.1.5时间大纲
活动时间
单细胞RNA-seq分析简介15米
zinbwave（降维）30米
clusterExperiment（群集）30米
弹弓（谱系轨迹分析）30米
问题/扩展15米
8.1.6研讨会的目标和目标
学习目标
描述了单细胞RNA-seq分析的目标
确定典型的单细胞RNA-seq分析的主要步骤
根据模型拟合评估每个步骤的结果
合成连续步骤的结果以解释生物学意义并开发生物模型
应用此工作流程对其他单细胞RNA-seq数据集进行完整分析
学习目标
计算和解释单细胞数据的低维表示
从数据中识别并删除技术变化的来源
识别细胞亚群（簇）并评估其稳健性
推断对应于区分细胞的谱系轨迹
通过发育“伪时间”命令细胞
识别在细胞分化中起重要作用的基因
8.2入门
本次研讨会介绍的工作流程包括四个主要步骤：
使用zinbwaveBioconductor软件包，考虑零膨胀和过度分散以及调整基因和细胞水平协变量的维数减少;
使用基于重采样的顺序集成聚类的强大且稳定的细胞聚类，如clusterExperimentBioconductor包中所实现的;
使用slingshotR包，通过沿着谱系的发育进展推断细胞谱系和细胞的排序; 和
沿着谱系的DE分析。
图8.1：分析scRNA-seq数据集的工作流程。在右侧，工作流生成的主要图表。
在整个工作流程中，我们使用单个SingleCellExperiment对象来存储scRNA-seq数据以及实验中可用的任何基因或细胞级元数据。
8.2.1数据
图8.2：小鼠嗅上皮细胞的干细胞分化。这个数字是在Fletcher等人的许可下复制的。（）。
该研讨会使用来自小鼠嗅上皮（OE）中干细胞分化的scRNA-seq研究的数据（Fletcher等人，）。嗅上皮包含成熟的嗅觉感觉神经元（mOSN），其通过神经发生在上皮细胞中通过球状基底细胞（GBC）的分化而不断更新，所述球状基底细胞是上皮中的活跃增殖细胞。当整个组织发生严重损伤时，嗅上皮细胞可以从称为水平基底细胞（HBC）的正常静止干细胞中再生，其被激活以分化并重建上皮细胞中的所有主要细胞类型。
我们使用scRNA-seq数据集作为案例研究来研究HBC干细胞在嗅上皮中存在的不同细胞类型的差异。为了绘制由HBC产生的多种细胞谱系的发育轨迹，使用Fluidigm C1微流体细胞捕获平台对FACS纯化的细胞进行scRNA-seq，然后进行Illumina测序。通过计数映射到它的读数的总数来量化给定细胞中每个基因的表达水平。然后使用类似于本工作流程中的统计分析管道将细胞分配到不同谱系。最后，使用体内谱系示踪通过实验验证结果。有关数据生成和统计方法的详细信息，请参阅（Fletcher等人; Risso等人b; Street等人; Risso等人a）。
在本次研讨会中，我们描述了一系列步骤，以恢复原始研究中发现的谱系，从基因x细胞矩阵的原始计数公开可用https://www.ncbi.nlm.nih.gov/geo/ query / acc.cgi？acc = GSE95601。
需要以下包。
suppressPackageStartupMessages({ # Bioconductor library(BiocParallel) library(SingleCellExperiment) library(clusterExperiment) library(scone) library(zinbwave) library(slingshot) # CRAN library(gam) library(RColorBrewer)})#> Warning: replacing previous import 'SingleCellExperiment::weights' by#> 'stats::weights' when loading 'slingshot'#> Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display#> Warning: 'rgl_init' failed, running with rgl.useNULL = TRUEset.seed(20)
8.2.2并行计算
该BiocParallel包可用于允许并行计算zinbwave。在这里，我们使用单个CPU来运行该功能，register进行串行模式BiocParallel。鼓励在系统中访问多个核心的用户使用多个核心来提高速度。
register(SerialParam())
8.3本SingleCellExperiment类
每个细胞中所有基因的计数可作为GitHub R包的一部分获得drisso/fletcherdata。在过滤之前，数据集具有849个细胞和28,361个检测到的基因（即具有非零读数的基因）。
library(fletcherdata)data(fletcher)fletcher#> class: SingleCellExperiment #> dim: 28284 849 #> metadata(0):#> assays(1): counts#> rownames(28284): Xkr4 LOC102640625 ... Ggcx.1 eGFP#> rowData names(0):#> colnames(849): OEP01_N706_S501 OEP01_N701_S501 ... OEL23_N704_S503#> OEL23_N703_S502#> colData names(19): Experiment Batch ... CreER ERCC_reads#> reducedDimNames(0):#> spikeNames(0):
在整个研讨会期间，我们使用该类SingleCellExperiment来跟踪单个对象中的计数及其相关元数据。
图8.3：SingleCellExperiment类的示意图。
细胞级元数据包含质量控制措施，序列批次ID以及原始出版物中的簇和谱系标签（Fletcher等人，）。具有群集标签的单元-2未分配给原始出版物中的任何群集。
colData(fletcher)#> DataFrame with 849 rows and 19 columns#> Experiment Batch publishedClusters NREADS#> <factor> <factor> <numeric> <numeric>#> OEP01_N706_S501 K5ERRY_UI_96HPT Y01 1 3313260#> OEP01_N701_S501 K5ERRY_UI_96HPT Y01 1 2902430#> OEP01_N707_S507 K5ERRY_UI_96HPT Y01 1 2307940#> OEP01_N705_S501 K5ERRY_UI_96HPT Y01 1 3337400#> OEP01_N704_S507 K5ERRY_UI_96HPT Y01 -2 117892#> ... ... ... ... ...#> OEL23_N704_S510 K5ERP63CKO_UI_14DPT P14 -2 2407440#> OEL23_N705_S502 K5ERP63CKO_UI_14DPT P14 -2 2308940#> OEL23_N706_S502 K5ERP63CKO_UI_14DPT P14 12 2215640#> OEL23_N704_S503 K5ERP63CKO_UI_14DPT P14 12 2673790#> OEL23_N703_S502 K5ERP63CKO_UI_14DPT P14 7 2450320#> NALIGNED RALIGN TOTAL_DUP PRIMER#> <numeric> <numeric> <numeric> <numeric>#> OEP01_N706_S501 3167600 95.6035 47.9943 0.0154566#> OEP01_N701_S501 2757790 95.0167 45.015 0.0182066#> OEP01_N707_S507 2178350 94.3852 43.7832 0.0219196#> OEP01_N705_S501 3183720 95.3952 43.2688 0.0183041#> OEP01_N704_S507 98628 83.6596 18.0576 0.0623744#> ... ... ... ... ...#> OEL23_N704_S510 2305060 95.7472 47.1489 0.0159111#> OEL23_N705_S502 2203300 95.4244 62.5638 0.0195812#> OEL23_N706_S502 2108490 95.1637 50.6643 0.0182207#> OEL23_N704_S503 2568300 96.0546 60.5481 0.0155611#> OEL23_N703_S502 2363500 96.4567 48.4164 0.0122704#> PCT_RIBOSOMAL_BASES PCT_CODING_BASES PCT_UTR_BASES#> <numeric> <numeric> <numeric>#> OEP01_N706_S501 2e-06 0.20013 0.230654#> OEP01_N701_S501 0 0.182461 0.1#> OEP01_N707_S507 0 0.152627 0.207897#> OEP01_N705_S501 2e-06 0.169514 0.207342#> OEP01_N704_S507 1.4e-05 0.110724 0.199174#> ... ... ... ...#> OEL23_N704_S510 0 0.287346 0.314104#> OEL23_N705_S502 0 0.337264 0.297077#> OEL23_N706_S502 7e-06 0.244333 0.262663#> OEL23_N704_S503 0 0.343203 0.338217#> OEL23_N703_S502 8e-06 0.259367 0.238239#> PCT_INTRONIC_BASES PCT_INTERGENIC_BASES PCT_MRNA_BASES#> <numeric> <numeric> <numeric>#> OEP01_N706_S501 0.404205 0.165009 0.430784#> OEP01_N701_S501 0.465702 0.150027 0.384271#> OEP01_N707_S507 0.511416 0.12806 0.360524#> OEP01_N705_S501 0.457556 0.165586 0.376856#> OEP01_N704_S507 0.489514 0.73 0.309898#> ... ... ... ...#> OEL23_N704_S510 0.250658 0.147892 0.60145#> OEL23_N705_S502 0.230214 0.135445 0.634341#> OEL23_N706_S502 0.355899 0.137097 0.506997#> OEL23_N704_S503 0.174696 0.143885 0.68142#> OEL23_N703_S502 0.376091 0.126294 0.497606#> MEDIAN_CV_COVERAGE MEDIAN_5PRIME_BIAS MEDIAN_3PRIME_BIAS#> <numeric> <numeric> <numeric>#> OEP01_N706_S501 0.843857 0.061028 0.521079#> OEP01_N701_S501 0.91437 0.03335 0.373993#> OEP01_N707_S507 0.955405 0.014606 0.49123#> OEP01_N705_S501 0.81663 0.101798 0.525238#> OEP01_N704_S507 1.19978 0 0.706512#> ... ... ... ...#> OEL23_N704_S510 0.698455 0.198224 0.419745#> OEL23_N705_S502 0.830816 0.105091 0.398755#> OEL23_N706_S502 0.805627 0.103363 0.431862#> OEL23_N704_S503 0.745201 0.118615 0.38422#> OEL23_N703_S502 0.711685 0.196725 0.377926#> CreER ERCC_reads#> <numeric> <numeric>#> OEP01_N706_S501 1 10516#> OEP01_N701_S501 3022 9331#> OEP01_N707_S507 2329 7386#> OEP01_N705_S501 717 6387#> OEP01_N704_S507 60 992#> ... ... ...#> OEL23_N704_S510 659 0#> OEL23_N705_S502 1552 0#> OEL23_N706_S502 0 0#> OEL23_N704_S503 0 0#> OEL23_N703_S502 2222 0
8.4预处理
使用Bioconductor软件包scone，我们根据功能中实施的质量控制过滤器metric_sample_filter并根据以下标准去除低质量的细胞（图8.4）:( 1）过滤掉总读数低或比对百分比低的样品（ 2）过滤掉管家基因检测率低的样本。有关过滤程序的详细信息，请参阅scone vignette。
8.4.1样本过滤
# QC-metric-based sample-filteringdata("housekeeping")hk = rownames(fletcher)[toupper(rownames(fletcher)) %in% housekeeping$V1]mfilt <- metric_sample_filter(counts(fletcher), nreads = colData(fletcher)$NREADS, ralign = colData(fletcher)$RALIGN, pos_controls = rownames(fletcher) %in% hk, zcut = 3, mixture = FALSE, plot = TRUE)
图8.4：SCONE：过滤低质量的细胞。
# Simplify to a single logicalmfilt <- !apply(simplify2array(mfilt[!is.na(mfilt)]), 1, any)filtered <- fletcher[, mfilt]dim(filtered)#> [1] 28284 747
样品过滤后，我们留下了747个优质细胞。
最后，为了提高计算效率，我们只保留了1000个变量最大的基因。这似乎是该工作流程的说明性目的的合理选择，因为我们能够恢复已发表的分析中发现的生物信号（（Fletcher等人，））。但是，一般而言，我们建议谨慎选择基因过滤方案，因为适当的选择是依赖于数据集。
我们可以使用clusterExperiment包中的函数来计算基于方差（makeFilterStats）的过滤器统计信息，并将过滤器应用于数据（filterData）。
filtered <- makeFilterStats(filtered, filterStats="var", transFun = log1p)filtered <- filterData(filtered, percentile=1000, filterStats="var")filtered#> class: SingleCellExperiment #> dim: 1000 747 #> metadata(0):#> assays(1): counts#> rownames(1000): Cbr2 Cyp2f2 ... Rnf13 Atp7b#> rowData names(1): var#> colnames(747): OEP01_N706_S501 OEP01_N701_S501 ... OEL23_N704_S503#> OEL23_N703_S502#> colData names(19): Experiment Batch ... CreER ERCC_reads#> reducedDimNames(0):#> spikeNames(0):
在最初的工作中（Fletcher等人，），细胞聚集成14个不同的聚类，151个细胞未分配到任何聚类（即聚类标记-2）。
publishedClusters <- colData(filtered)[, "publishedClusters"]col_clus <- c("transparent", "#1B9E77", "antiquewhite2", "cyan", "#E7298A", "#A6CEE3", "#666666", "#E6AB02", "#FFED6F", "darkorchid2", "#B3DE69", "#FF7F00", "#A6761D", "#1F78B4")names(col_clus) <- sort(unique(publishedClusters))table(publishedClusters)#> publishedClusters#> -2 1 2 3 4 5 7 8 9 10 11 12 14 15 #> 151 90 25 54 35 93 58 27 74 26 21 35 26 32
8.5归一化和降维：ZINB-WaVE
在scRNA-seq分析中，维度降低通常用作下游分析之前的初步步骤，例如聚类，细胞谱系和假定位排序，以及DE基因的鉴定。这使得数据变得更容易处理，无论是从统计（参见维数的诅咒）还是从计算的角度来看。另外，可以减少技术噪声，同时保留通常固有的低维信息（Dijk等人; Pierson和Yau ; Risso等人b）。
在这里，我们使用在Bioconductor R包中实现的零膨胀的基于负二项式的期望变化提取（ZINB-WaVE）方法来执行降维zinbwave。该方法适合ZINB模型，该模型考虑零膨胀（丢失），过度分散和数据的计数性质。该模型可以包括细胞级截距，其用作全局缩放归一化因子。用户还可以指定基因水平和细胞水平协变量。包含观察到的和未观察到的细胞水平协变量使得能够归一化复杂的非线性效应（通常称为批次效应），而基因水平协变量可用于调整序列组成效应（例如，基因长度和GC-内容效果）。图中提供了ZINB-WaVE模型的示意图8.5。有关ZINB-WaVE模型和估算程序的更多详细信息，请参阅原始手稿（Risso等人，b）。
图8.5：ZINB-WaVE：ZINB-WaVE模型的示意图。这个数字是在Risso等人的许可下复制的。（）。
与大多数降维方法一样，用户需要指定新的低维空间的维数。在这里，我们使用K = 50尺寸并通过矩阵调整批次效果X。
clustered <- zinbwave(filtered, K = 50, X = "~ Batch", residuals = TRUE, normalizedValues = TRUE)))
请注意，fletcherdata程序包包含clustered已包含ZINB-WaVE因子的对象。我们可以加载这些对象以避免等待计算。
data(clustered)
8.5.1规范化
该函数zinbwave返回一个SingleCellExperiment对象，该对象包括规范化表达式度量，定义为ZINB-WaVE模型与用户指定的基因和细胞级协变量拟合的偏差残差。这些残差可用于可视化目的（例如，在热图，箱形图中）。注意，在这种情况下，低维矩阵W不包括在残差的计算中以避免去除感兴趣的生物信号。
assayNames(clustered)#> [1] "normalizedValues" "residuals" "counts"norm <- assay(clustered, "normalizedValues")norm[1:3,1:3]#> OEP01_N706_S501 OEP01_N701_S501 OEP01_N707_S507#> Cbr2 4.531898 4.369185 -4.142982#> Cyp2f2 4.359680 4.324476 4.124527#> Gstm1 4.724216 4.621898 4.403587
8.5.2降维
该zinbwave函数的主要用途是执行降维。得到的低维矩阵W存储在reducedDim名为的槽中zinbwave。
reducedDimNames(clustered)#> [1] "zinbwave"W <- reducedDim(clustered, "zinbwave")dim(W)#> [1] 747 50W[1:3, 1:3]#> W1 W2 W3#> OEP01_N706_S501 0.5494761 1.1745361 -0.93175747#> OEP01_N701_S501 0.4116797 0.3015379 -0.46922527#> OEP01_N707_S507 0.7394759 0.3365600 -0.07959226
W通过使用欧几里德距离执行多维缩放（MDS），可以在二维中可视化低秩矩阵。为了验证W确实捕获了感兴趣的生物信号，我们在散点图中显示MDS结果，其颜色对应于原始发布的簇（图8.6）。
W <- reducedDim(clustered)d <- dist(W)fit <- cmdscale(d, eig = TRUE, k = 2)plot(fit$points, col = col_clus[as.character(publishedClusters)], main = "", pch = 20, xlab = "Component 1", ylab = "Component 2")legend(x = "topleft", legend = unique(names(col_clus)), cex = .5, fill = unique(col_clus), title = "Sample")
图8.6：ZINB-WaVE：低维矩阵W的MDS，其中每个点代表一个细胞，细胞通过原始发表的聚类进行颜色编码。
8.6细胞聚类：RSEC
下一步是根据W在前一步骤中计算的低维矩阵对单元进行聚类。我们使用RSECBioconductor R包中的函数中实现的基于重采样的顺序集合聚类（RSEC）框架clusterExperiment。具体而言，给定一组用户提供的基本聚类算法和相关的调整参数（例如，k -means，具有k的值范围），RSEC生成候选聚类的集合，可选择重新采样单元格并使用顺序紧密集群程序（Tseng and Wong ）。基于候选聚类上的样本的共聚类水平获得共识聚类。通过合并类似的集群进一步压缩共识聚类，这通过创建集群层次结构，处理树以及测试姐妹节点之间的差异表达来完成，其中DE不足的节点崩溃。与监督学习一样，重采样极大地提高了聚类的稳定性，并且考虑到方法和调整参数的集合使我们能够利用基本算法的不同优势并避免调整参数的主观选择。
clustered <- RSEC(clustered, k0s = 4:15, alphas = c(0.1), betas = 0.8, reduceMethod="zinbwave", clusterFunction = "hierarchical01", minSizes=1, ncores = NCORES, isCount=FALSE, dendroReduce="zinbwave", subsampleArgs = list(resamp.num=100, clusterFunction="kmeans", clusterArgs=list(nstart=10)), verbose=TRUE, consensusProportion = 0.7, mergeMethod = "none", random.seed=424242, consensusMinSize = 10)
同样，先前加载的clustered对象已包含RSEC上面运行的结果，因此我们不在此处评估上面的块。
clustered#> class: ClusterExperiment #> dim: 1000 747 #> reducedDimNames: zinbwave #> filterStats: var #> -----------#> Primary cluster type: makeConsensus #> Primary cluster label: makeConsensus #> Table of clusters (of primary clustering):#> -1 c1 c2 c3 c4 c5 c6 c7 #> 184 145 119 105 100 48 33 13 #> Total number of clusterings: 13 #> Dendrogram run on 'makeConsensus' (cluster index: 1)#> -----------#> Workflow progress:#> clusterMany run? Yes #> makeConsensus run? Yes #> makeDendrogram run? Yes #> mergeClusters run? No
请注意，RSEC函数的结果是ClusterExperiment类的对象，它SingleCellExperiment通过添加有关聚类结果的其他信息来扩展类。
is(clustered, "SingleCellExperiment")#> [1] TRUEslotNames(clustered)#> [1] "transformation" "clusterMatrix" #> [3] "primaryIndex" "clusterInfo" #> [5] "clusterTypes" "dendro_samples" #> [7] "dendro_clusters" "dendro_index" #> [9] "dendro_outbranch" "coClustering" #> [11] "clusterLegend" "orderSamples" #> [13] "merge_index" "merge_dendrocluster_index"#> [15] "merge_method" "merge_demethod" #> [17] "merge_cutoff" "merge_nodeProp" #> [19] "merge_nodeMerge" "int_elementMetadata" #> [21] "int_colData" "int_metadata" #> [23] "reducedDims" "rowRanges" #> [25] "colData" "assays" #> [27] "NAMES" "elementMetadata" #> [29] "metadata"
可以使用plotClusters函数可视化生成的候选聚类（图8.7），其中列对应于不同聚类的单元格和行。每个样本都根据其对该行的聚类进行颜色编码，其中已选择颜色以尝试匹配跨行显示大的重叠的聚类。第一行对应于跨所有候选聚类的一致聚类。
plotClusters(clustered)
图8.7：RSEC：使用clusterExperiment包中的RSEC函数找到的候选聚类。
该plotCoClustering函数生成共聚类矩阵的热图，该矩阵为每对细胞记录它们在候选聚类中聚集在一起的时间比例（图8.8）。
plotCoClustering(clustered)
图8.8：RSEC：共聚类矩阵的热图。
整个共识簇的细胞分布可以在图8.9中看到，如下所示：
table(primaryClusterNamed(clustered))#> #> -1 c1 c2 c3 c4 c5 c6 c7 #> 184 145 119 105 100 48 33 13
plotBarplot(clustered, legend = FALSE)
图8.9：RSEC：我们工作流的RSEC群集中每个群集的小区数量的条形图。
我们的聚类中的单元格分布总体上与原始发布的聚类中的细胞分布一致（图8.10），主要区别在于几个已发布的聚类在此合并为单个聚类。这种差异可能是由于我们从前1000个基因开始，这可能不足以区分密切相关的聚类。
clustered <- addClusterings(clustered, colData(clustered)$publishedClusters, clusterLabel = "publishedClusters")## change default color to match with Figure 7clusterLegend(clustered)$publishedClusters[, "color"] <- col_clus[clusterLegend(clustered)$publishedClusters[, "name"]]plotBarplot(clustered, whichClusters=c("makeConsensus", "publishedClusters"), xlab = "", legend = FALSE,missingColor="white")
图8.10：RSEC：每个群集的小区数量的条形图，用于我们的工作流程的RSEC群集，由原始发布的群集进行颜色编码。
plotClustersTable(clustered, whichClusters=c("makeConsensus","publishedClusters"))
图8.11：RSEC：每个集群的单元数量的混淆矩阵，用于我们的工作流程的RSEC集群和原始发布的集群。
图8.12显示了1000个变异基因的标准化表达量度的热图，其中细胞根据RSEC共识进行聚类。
# Set colors for additional sample dataexperimentColors <- bigPalette[1:nlevels(colData(clustered)$Experiment)]batchColors <- bigPalette[1:nlevels(colData(clustered)$Batch)]metaColors <- list("Experiment" = experimentColors, "Batch" = batchColors)plotHeatmap(clustered, whichClusters = c("makeConsensus","publishedClusters"), clusterFeaturesData = "all", clusterSamplesData = "dendrogramValue", breaks = 0.99, colData = c("Batch", "Experiment"), clusterLegend = metaColors, annLegend = FALSE, main = "")
图8.12：RSEC：标准化表达的热图测量1,000个最可变的基因，其中行对应于RSEC簇排序的细胞的基因和列。
最后，我们可以使用低维矩阵的MDS在二维空间中可视化细胞，W并根据新发现的RSEC簇对细胞着色（图8.13）; 对于原始发布的簇，这与图8.6是同源的。
plotReducedDims(clustered,whichCluster="primary",reducedDim="zinbwave",pch=20, xlab = "Component1", ylab = "Component2",legendTitle="Sample",main="", plotUnassigned=FALSE)
图8.13：RSEC：低维矩阵W的MDS，其中每个点代表一个小区，并且通过RSEC聚类对小区进行颜色编码。
8.7细胞谱系和伪时间推断：弹弓
我们现在演示如何使用Bioconductor包slingshot来推断分支细胞谱系并通过沿着每个谱系的发育进程来命令细胞。（Street et al.）中提出的方法包括两个主要步骤：（1）使用最小生成树（MST）推断全局谱系结构（即，谱系数和它们分支的位置）以上确定的集群RSEC（2）使用同时主曲线的新方法推断每个谱系的细胞伪时间变量。（1）中的方法允许识别任何数量的新谱系，同时还适应使用领域特定知识来监督树的部分（例如，已知的终端状态）; （2）中的方法为平滑的分支谱系产生了强大的伪时间。
该分析由slingshot函数执行，结果存储在SlingshotDataSet对象中。该函数的最小输入是单元的低维表示和一组簇标签; 这些可以是单独的对象（即矩阵和向量），或者如下所示，是SingleCellExperiment对象的组件。当提供SingleCellExperiment对象作为输入时，输出将是包含SlingshotDataSet作为int_metadata列表元素的更新对象，可以通过该SlingshotDataSet函数访问该对象。对于谱系推断过程的更低级控制，这两个步骤可以通过函数getLineages和单独运行getCurves。
从最初发表的工作中，我们知道起始簇应该对应于HBC，并且末端簇对应于MV，mOSN和mSUS单元。此外，我们知道在MV和mOSN细胞之间的差异之前，GBC应该处于结点（图8.2）。我们在这里找到的聚类与原始聚类之间的对应关系如下。
table(data.frame(original = publishedClusters, ours = primaryClusterNamed(clustered)))#> ours#> original -1 c1 c2 c3 c4 c5 c6 c7#> -2 55 45 30 5 6 1 3 6#> 1 44 46 0 0 0 0 0 0#> 2 1 0 0 0 24 0 0 0#> 3 2 2 1 0 49 0 0 0#> 4 2 2 31 0 0 0 0 0#> 5 44 47 2 0 0 0 0 0#> 7 4 0 54 0 0 0 0 0#> 8 26 1 0 0 0 0 0 0#> 9 0 0 1 65 1 0 0 7#> 10 1 0 0 0 0 25 0 0#> 11 2 2 0 0 17 0 0 0#> 12 0 0 0 35 0 0 0 0#> 14 2 0 0 0 2 22 0 0#> 15 1 0 0 0 1 0 30 0
群集名称描述颜色对应
C1HBC红色原来的1,5
C2管理支持蓝色原创4,7
C3mOSN绿色原9,12
C4GBC橙子原来的2,3,11
C5不成熟的神经元紫色原10,14
C6MV棕色原来的15
C7mOSN浅蓝原来9
为了推断谱系和伪时间，我们将Slingshot应用于低维矩阵的4维MDS W。我们发现Slingshot结果对于MDS 的维数k是稳健的（我们尝试k从2到5）。在这里，我们使用Slingehot的半监督版本，我们只提供起始集群的身份，但不提供终端集群的身份。
pseudoCe <- clustered[,!primaryClusterNamed(clustered) %in% c("-1")]X <- reducedDim(pseudoCe,type="zinbwave")mds <- cmdscale(dist(X), eig = TRUE, k = 4)lineages <- slingshot(mds$points, clusterLabels = primaryClusterNamed(pseudoCe), start.clus = "c1")#> Using full covariance matrix
在讨论同时主要曲线之前，我们通过绘制群集上的MST来检查谱系的全局结构。这表明我们的实现已经恢复了已发表作品中的谱系（图8.14）。该slingshot软件包还包括使用软件包中的功能进行三维可视化的功能，如图8.2所示。plot3drgl
colorCl<-convertClusterLegend(pseudoCe,whichCluster="primary",output="matrixColors")[,1]pairs(lineages, type="lineages", col = colorCl)
图8.14：弹弓：在4维MDS空间中通过聚类进行颜色编码的单元格，聚类中心之间的连接线表示推断的全局谱系结构。
找到全局谱系结构后，slingshot构建一组平滑的分支曲线，以推断出相应的伪时间变量。从每个谱系的单个细胞构建同时主曲线，而不是细胞簇。在此迭代过程中，如果单元格明显更接近于相应的曲线，则甚至可以将单元格重新分配给不同的谱系。这slingshot减少了对原始聚类的依赖，并且通常更稳定。最终曲线如图8.15所示。
pairs(lineages, type="curves", col = colorCl)
图8.15：弹弓：细胞在4维MDS空间中通过簇进行颜色编码，平滑曲线代表每个推断的谱系。
lineages#> class: SlingshotDataSet #> #> Samples Dimensions#> 563 4#> #> lineages: 3 #> Lineage1: c1 c4 c5 c7 c3 #> Lineage2: c1 c4 c6 #> Lineage3: c1 c2 #> #> curves: 3 #> Curve1: Length: 9.5238 Samples: 361.21#> Curve2: Length: 7.8221 Samples: 234.14#> Curve3: Length: 4.2828 Samples: 254.85
作为替代方案，我们可以将MDS结果合并到clustered对象中并slingshot直接应用于它。在这里，我们需要指定我们想要使用MDS结果，因为slingshot否则将使用reducedDims列表的第一个元素（在这种情况下，来自10维W矩阵zinbwave）。
reducedDim(pseudoCe, "MDS") <- mds$pointspseudoCe <- slingshot(pseudoCe, reducedDim = "MDS", start.clus = "c1")#> Using full covariance matrixpseudoCe#> class: ClusterExperiment #> dim: 1000 563 #> reducedDimNames: zinbwave MDS #> filterStats: var #> -----------#> Primary cluster type: makeConsensus #> Primary cluster label: makeConsensus #> Table of clusters (of primary clustering):#> c1 c2 c3 c4 c5 c6 c7 #> 145 119 105 100 48 33 13 #> Total number of clusterings: 14 #> No dendrogram present#> -----------#> Workflow progress:#> clusterMany run? Yes #> makeConsensus run? Yes #> makeDendrogram run? No #> mergeClusters run? NocolData(pseudoCe)#> DataFrame with 563 rows and 23 columns#> Experiment Batch publishedClusters NREADS#> <factor> <factor> <numeric> <numeric>#> OEP01_N706_S501 K5ERRY_UI_96HPT Y01 1 3313260#> OEP01_N701_S501 K5ERRY_UI_96HPT Y01 1 2902430#> OEP01_N707_S507 K5ERRY_UI_96HPT Y01 1 2307940#> OEP01_N705_S501 K5ERRY_UI_96HPT Y01 1 3337400#> OEP01_N702_S508 K5ERRY_UI_96HPT Y01 -2 525096#> ... ... ... ... ...#> OEL23_N704_S510 K5ERP63CKO_UI_14DPT P14 -2 2407440#> OEL23_N705_S502 K5ERP63CKO_UI_14DPT P14 -2 2308940#> OEL23_N706_S502 K5ERP63CKO_UI_14DPT P14 12 2215640#> OEL23_N704_S503 K5ERP63CKO_UI_14DPT P14 12 2673790#> OEL23_N703_S502 K5ERP63CKO_UI_14DPT P14 7 2450320#> NALIGNED RALIGN TOTAL_DUP PRIMER#> <numeric> <numeric> <numeric> <numeric>#> OEP01_N706_S501 3167600 95.6035 47.9943 0.0154566#> OEP01_N701_S501 2757790 95.0167 45.015 0.0182066#> OEP01_N707_S507 2178350 94.3852 43.7832 0.0219196#> OEP01_N705_S501 3183720 95.3952 43.2688 0.0183041#> OEP01_N702_S508 484847 92.3349 18.806 0.0248804#> ... ... ... ... ...#> OEL23_N704_S510 2305060 95.7472 47.1489 0.0159111#> OEL23_N705_S502 2203300 95.4244 62.5638 0.0195812#> OEL23_N706_S502 2108490 95.1637 50.6643 0.0182207#> OEL23_N704_S503 2568300 96.0546 60.5481 0.0155611#> OEL23_N703_S502 2363500 96.4567 48.4164 0.0122704#> PCT_RIBOSOMAL_BASES PCT_CODING_BASES PCT_UTR_BASES#> <numeric> <numeric> <numeric>#> OEP01_N706_S501 2e-06 0.20013 0.230654#> OEP01_N701_S501 0 0.182461 0.1#> OEP01_N707_S507 0 0.152627 0.207897#> OEP01_N705_S501 2e-06 0.169514 0.207342#> OEP01_N702_S508 0 0.130247 0.230848#> ... ... ... ...#> OEL23_N704_S510 0 0.287346 0.314104#> OEL23_N705_S502 0 0.337264 0.297077#> OEL23_N706_S502 7e-06 0.244333 0.262663#> OEL23_N704_S503 0 0.343203 0.338217#> OEL23_N703_S502 8e-06 0.259367 0.238239#> PCT_INTRONIC_BASES PCT_INTERGENIC_BASES PCT_MRNA_BASES#> <numeric> <numeric> <numeric>#> OEP01_N706_S501 0.404205 0.165009 0.430784#> OEP01_N701_S501 0.465702 0.150027 0.384271#> OEP01_N707_S507 0.511416 0.12806 0.360524#> OEP01_N705_S501 0.457556 0.165586 0.376856#> OEP01_N702_S508 0.477167 0.161738 0.361095#> ... ... ... ...#> OEL23_N704_S510 0.250658 0.147892 0.60145#> OEL23_N705_S502 0.230214 0.135445 0.634341#> OEL23_N706_S502 0.355899 0.137097 0.506997#> OEL23_N704_S503 0.174696 0.143885 0.68142#> OEL23_N703_S502 0.376091 0.126294 0.497606#> MEDIAN_CV_COVERAGE MEDIAN_5PRIME_BIAS MEDIAN_3PRIME_BIAS#> <numeric> <numeric> <numeric>#> OEP01_N706_S501 0.843857 0.061028 0.521079#> OEP01_N701_S501 0.91437 0.03335 0.373993#> OEP01_N707_S507 0.955405 0.014606 0.49123#> OEP01_N705_S501 0.81663 0.101798 0.525238#> OEP01_N702_S508 1.13937 0 0.671565#> ... ... ... ...#> OEL23_N704_S510 0.698455 0.198224 0.419745#> OEL23_N705_S502 0.830816 0.105091 0.398755#> OEL23_N706_S502 0.805627 0.103363 0.431862#> OEL23_N704_S503 0.745201 0.118615 0.38422#> OEL23_N703_S502 0.711685 0.196725 0.377926#> CreER ERCC_reads slingClusters slingPseudotime_1#> <numeric> <numeric> <character> <numeric>#> OEP01_N706_S501 1 10516 c1 NA#> OEP01_N701_S501 3022 9331 c1 1.17234765019331#> OEP01_N707_S507 2329 7386 c1 1.05857338852636#> OEP01_N705_S501 717 6387 c1 1.60463630065258#> OEP01_N702_S508 6 1218 c1 1.15934111346639#> ... ... ... ... ...#> OEL23_N704_S510 659 0 c2 NA#> OEL23_N705_S502 1552 0 c2 NA#> OEL23_N706_S502 0 0 c3 8.14552572709918#> OEL23_N704_S503 0 0 c3 8.53595122837615#> OEL23_N703_S502 2222 0 c2 NA#> slingPseudotime_2 slingPseudotime_3#> <numeric> <numeric>#> OEP01_N706_S501 NA 0.692789182245903#> OEP01_N701_S501 1.16364108574816 1.14696263207989#> OEP01_N707_S507 1.0611741023521 1.03755897807688#> OEP01_N705_S501 1.61033070603652 1.4465916227885#> OEP01_N702_S508 1.1674954947514 1.4206584204892#> ... ... ...#> OEL23_N704_S510 NA 2.01842841303396#> OEL23_N705_S502 NA 3.75230631265741#> OEL23_N706_S502 NA NA#> OEL23_N704_S503 NA NA#> OEL23_N703_S502 NA 2.74576551163381
slingshot应用于ClusterExperiment对象的结果仍然是类ClusterExperiment。请注意，我们没有指定一组群集标签，这意味着slingshot应该使用默认primaryClusterNamed向量。
在工作流程中，我们使用很大程度上无监督的版本来恢复群集的合理排序slingshot。但是，在其他一些情况下，我们注意到我们需要为算法提供更多指导以找到正确的顺序。getLineages用户可以选择提供已知的端集群，它表示对MST的约束，要求这些集群为叶节点。以下是slingshot在监督设置中使用的代码，我们知道集群c3，c6并c2代表终端单元命运。
pseudoCeSup <- slingshot(pseudoCe, reducedDim = "MDS", start.clus = "c1", end.clus = c("c3", "c6", "c2"))
8.8沿谱系的差异表达分析
在将细胞分配到谱系并在谱系中对它们进行排序之后，我们感兴趣的是找到在伪时间内具有非恒定表达模式的基因。
更正式地说，对于每个谱系，我们使用稳健的局部回归方法黄土以灵活，非线性的方式模拟基因的标准化表达量度和假性时间之间的关系。然后，我们可以使用gam包测试每个基因随时间没有变化的零假设。我们对神经元谱系实施这种方法，并在图8.16的热图中通过p值显示前100个基因的表达量度。
t <- colData(pseudoCe)$slingPseudotime_1y <- transformData(pseudoCe)gam.pval <- apply(y,1,function(z){ d <- data.frame(z=z, t=t) tmp <- gam(z ~ lo(t), data=d) p <- summary(tmp)[4][[1]][1,5] p})
topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:100]pseudoCe1 <- pseudoCe[,!is.na(t)]orderSamples(pseudoCe1)<-order(t[!is.na(t)])plotHeatmap(pseudoCe1[topgenes,], clusterSamplesData = "orderSamplesValue", breaks = .99)
图8.16：DE：标准化表达的热图测量神经元谱系的100个最显着的DE基因，其中行对应于假定的有序细胞的基因和列。
