单细胞分析基础流程理论 - 洪瑜的博客

主要参考学习生信星球公众号，仅个人学习使用。

整体流程介绍

首先是数据导入
紧接着进行质量控制和归一化，这个步骤包含了矫正实验因素以及其他影响因素，目的是从 raw count 得到 clean count
有了 clean count，就要挑取导致表达量差异的基因们（用于后面的降维）。因为降维的过程其实就是去繁存简，每个基因的变化都对整体有影响，对细胞来讲都是一个变化维度。但这么多维度我们看不了也处理不了，需要尽可能保证真实差异的前提下减少维度的数量，因此需要挑出那些更能代表整体差异的基因
分群（clustering）：意在探索如何把 scRNA 数据集给拆分掉，变成一小块一小块的数据，这每一小块基因变化都是相似的
差异分析（differential expression）：不同组的细胞之间表达量差异是如何产生的
数据集的整合(Integrating datasets)：scRNA 数据集越来越多，数据集之间的比较也日益显著
处理大型数据（large scale data）：这部分仅仅依靠内存存储的数据是不够的，还有更好的办法
还有：轨迹推断、细胞周期推断等等特定需求

# 准备好表达矩阵后使用list构建出SingleCellExperiment对象
sce <- SingleCellExperiment(assays = list(counts = counts_matrix))
# 获得singlecellexperiment中的表达矩阵
assay(sce, "counts")
counts(sce)

# 扩展assays
sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)
logcounts(sce) #提取对应矩阵

# Normalization归一化：去除样本文库差异(测序深度不一致)

# 自定义扩展assays
counts_100 <- assay(sce, "counts") + 100
assay(sce, "counts_100") <- counts_100

# 关于注释信息metadata
# 细胞注释信息
# 直接构建
sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
                         colData = cell_metadata)
# 后续添加
colData(sce) <- DataFrame(cell_metadata) # 数据框格式
colData(sce) # 获取metadata
sce <- scater::addPerCellQC(sce) # 自动计算一些metadata
# 基因注释信息
# rowData：是一个数据框的结构，它就存储了核心assays矩阵的基因相关信息

# 归一化--文库大小因子
sce <- scran::computeSumFactors(sce)
sizeFactors(sce)
sizeFactors(sce) <- scater::librarySizeFactors(sce)
sizeFactors(sce)

# 降维结果--reducedDims
# 数值型矩阵的list，行名是原来矩阵的列名（就是细胞、样本），列就是各种维度信息
reducedDim(sce, "PCA")
reducedDim(sce, "TSNE")
reducedDim(sce)

# metadata接口
# 存储任意类型的数据，只要给它一个名字。
metadata(sce) <- list(favorite_genes = my_genes)
metadata(sce)

# 对列添加label
colLabels(sce) <- LETTERS[1:3]
colLabels(sce)

想要发现罕见细胞类群，就要多获得细胞；想要探索潜在的微小差异，就要加大测序深度。

10X 的数据：使用 CellRanger 软件，基于 STAR 比对到参考基因组，然后统计每个基因的 UMIs 数量
Pseudo-alignment 方法（如 alevin）：就像之前用的 salmon、kallisto 意思一样，不需要比对参考基因组，节省时间、内存
对于一些高度 multiplexed 的方法：可以使用 scPipe 包：提供了一套综合的分析流程，利用 Rsubread 比对，然后统计每个基因的 UMIs 数量
CEL-seq、CEL-seq2 数据：scruff 包可以专门分析
read-based 方法：可以使用常规 bulk 转录组定量的流程（比如smartseq2 就可以用 hisat2+featureCounts）
任何包含 spike-in 转录本的数据：spike-in 序列都要在比对、定量之前加到参考基因组中

细胞质控

文库大小（library size）：指的是每个细胞中所有基因的表达量之和。细胞的文库如果很小，说明在文库制备过程中存在 RNA 的损失，要么是由于细胞裂解，要么是 cDNA 捕获不足导致后续的扩增量少
每个细胞中表达基因的数目（number of expressed features in each cell）：指的是细胞中非 0 表达量的基因数量。如果细胞中基本没有基因表达，很可能是转录本捕获失败
比对到 spike-in 转录本的 reads 比例（proportion of reads mapped to spike-in transcripts）：计算方法是：spike-in counts / all features (including spike-ins) for each cell。每个细胞都应该加入等量的外源转录本（spike-in），如果哪个细胞的 spike-in 比例提高了，说明它的内源 RNA 含量减少（比如在细胞分选阶段出现的细胞裂解或者 RNA 降解）
比对到线粒体基因组的 reads 比例（proportion of reads mapped to genes in the mitochondrial genome）：如果没有 spike-in，那么使用线粒体指标也是能说明问题的。比对到线粒体的 reads 增多，说明细胞质中的 RNA 减少，可能存在细胞穿孔的情况，而这个孔的大小，可能只是将细胞质中存在的 mRNA 流出去，但线粒体的体积超过了孔的大小，所以还留在了细胞中，造成一定程度的富集，导致线粒体 RNA 占比升高。

scater 包的 perCellQCMetrics()函数，其中 sum 这一列表示每个细胞的文库大小；detected 这一列表示检测到的基因数量；subsets_Mito_percent 这一列表示比对到线粒体基因组的 reads 占比；altexps_ERCC_percent 表示比对到 ERCC spike-in 的 reads 占比

scRNA 的归一化 normalization

scRNA 数据中文库的差异还是很常见的，来源于：不同细胞的起始底物浓度不同，导致 cDNA 捕获或 PCR 扩增效率差异。归一化的目的就是去除细胞间与真实表达量无关的技术因素，方便后续比较。

归一化与批次处理还是不同的。**归一化不管实验的批次因素，只考虑细胞中存在的技术误差（比如测序深度），而批次处理既要考虑实验批次，又要考虑技术误差（比如不同实验时间、不同细胞系、不同文库制备方法、不同测序方法、不同测序深度）。技术误差对不同基因的影响是相似的，比如有的技术误差会影响基因长度、GC 含量，那么基本上所有基因都会受影响。但不同的批次产生的影响是不同的，比如不同时间做的实验效果就是不一样，那么基因的表达量可能也会受到这个影响。尽管有的 R 包可以同时处理技术误差和批次效应（例如 zinbwave），但大部分的分析还是各顾各的。

归一化，归的是技术误差，而不是批次效应

normalization 为归一化，scale 为标准化,一般来说，应该先进行 normalization，再进行 scale，而 scale 的结果会用来后续的降维和聚类

Normalization “normalizes” within the cell for the difference in sequenicng depth / mRNA throughput. 主要着眼于样本的文库大小差异
Scaling “normalizes” across the sample for differences in range of variation of expression of genes . 主要着眼于基因的表达分布差异

文库大小就是指每个细胞中所有基因表达量之和

文库相关的 size factor 基于的假设是：任意两对细胞间的差异表达基因都是平衡的。也就是说，一组基因的上调表达程度，势必会被另一组基因的下调表达程度相抵消。但是在 scRNA 中差异基因表达并不是平衡的（也称为：存在 composition biases），因此 library size normalization 可能不会得到准确的归一化结果，仅仅是能用，但不是最精确。精确的归一化结果也不是 scRNA 数据探索的主要任务。归一化的目的是更好地为下游展示数据做铺垫。因此使用 library size normalization 对细胞分群和鉴定每群的 marker 基因也是足够的。

composition biases

现在想象一个例子：两个细胞 A、B，一个基因 X 在 A 细胞相比于 B 细胞表达量更高。这种高表达意味着：

当细胞文库大小一定（例如实验采用 library quantification/文库定量的方法），更多的测序资源偏向 X 基因，势必会降低其他本来不是差异表达基因的覆盖度【简言之，一个盘子就这么多水果，一个异类胃口大开，吃了很多，其他人本来也能吃饱，但现在也要挨饿】
当不限制细胞文库大小时，X 基因的 reads 或 UMIs 增加，也会增加总体的文库大小，导致计算的 library size factor 增加，原本非差异基因的表达量也会由于 size factor 的增加而相对之前降低【简言之，一个异类拉高了整体平均成绩，导致其他本来成绩平平的学生看上去成绩发生了下降】

这两种情况都会导致一个结果：细胞 A 的其他非差异基因都会由于 X 基因的”过度表现“而不经意地”下调“

那么如何去除这种 composition biases 呢？之前在 bulk RNA-Seq 中研究非常透彻了。例如 DESeq2 使用的 estimateSizeFactorsFromMatrix()，edgeR 使用的 calcNormFactors() 都在归一化过程中考虑到了这点。它的假设是：大部分基因都不是差异基因。

去卷积方法计算 size factor,得到更准确的结果，虽然对分群改善不大，但对于每个基因的表达量数据估计与解释还是很重要的，因为它考虑了 composition biases 的影响。

使用 spike-in 计算 size factor,基于的假设是：外源 RNA（spike-in）添加到每个细胞时的量都是一样的。如果出现同一 spike-in 转录本表达量在各个细胞之间的差异，问题只能是细胞相关，例如不同细胞捕获效率、不同细胞测序深度等等。

添加到每个细胞的 spike-in 的量都是一定的
对偏差的反应，和内源基因是一样的

对于归一化方法的选择，取决于生物学假设。大多数情况下，总体 RNA 含量的变化我们并不感兴趣，可以直接用 library size 或 deconvolution factors 当做系统性偏差去除。但如果总体 RNA 含量变化与某个感兴趣的生物学过程相关，例如细胞周期活性或 T 细胞激活，spike-in 的归一化方法就可以把这种变化保留下来，而去除其他系统性偏差。

得到 Size Factors 后进行数据转换

logNormCounts: 用某个细胞中每个基因或 spike-in 转录本的表达量除以这个细胞计算的 size factor，最后还进行一个 log 转换，得到一个新的矩阵：logcounts
至于为什么取 log？就是为了不受真实值的影响，而关注变化倍数。这里举个例子，X 基因在细胞 A 的表达量是 50，在 B 细胞的表达量是 10；而 Y 基因在细胞 A 的表达量是 1000，在 B 细胞的表达量是 1100。哪个基因更能吸引你的注意呢？我们是不是会关注变化的倍数？
另外，既然要取 log，就应该要注意存在很多 0 表达量的数值，因此 log 前还需要给表达量加上 1，log1p = log(x + 1)，这个 1 就称作：pseudo-count。

挑选高变化基因

挑选哪些基因进行聚类分群和降维分析，这是非常关键的，挑选的基因一定要有代表性，尽可能多的包含生物信息而剔除其他技术噪音。
一般来说，如果一个基因在细胞群体中变化幅度很大，就是受关注对象，我们也会认为是生物因素导致了这么大的差异。这样的基因叫做：高度变化基因（highly variable genes ，HVGs）

衡量每个基因的变化程度/方差(Variance)

使用 log-counts 衡量变化程度基因 log 后的表达量在细胞间变化幅度越大，就说明细胞间距离越远（这个距离指欧氏距离 Euclidean distances）
计算每个基因的变化程度很简单，但要基于这么多基因的结果去挑选最有价值的基因（feature），还需要构建一个模型。这个模型叫：mean-variance relationship
使用变异系数(CV)衡量变化程度

如果两组数据的观测值在一个数量级，那么可以用标准差来比较它们的离散度
如果两个数据差别太大（比如一群老鼠和一群大象的重量），然后就要用 CV 或 CV^2，这样可以去除不同量纲均值差异的影响如果 CV^2 越大，偏离拟合曲线越远，就可能包含更多的生物偏差，更具有生物学意义

技术误差估计

基于 spike-in 有 spike-in 的话，它和内源基因不同，不会受到生物因素的影响。利用 spike-in 画的拟合线，更能代表技术偏差
没有 spike-in 可以考虑利用数据分布来表示技术噪音，例如只考虑技术噪音的话，UMI counts 通常会呈现近似泊松分布

批次效应的影响

绘制拟合曲线时加入批次信息（Fitting block-specific trends）在一个批次中的 HVGs 信息。一般可以一个批次一个批次地去处理，看看不同批次的差异
利用实验设计矩阵（design matrix）类似于 DEseq2 中的实验设计矩阵

挑选高变化基因 HVGs

基于方差最大化提取在生物学因素中方差最大的（也就是对波动贡献最大）基因。优点就是可以自己控制基因的数量，方便对后面的分析复杂度有一个大体的估计。比较主观。
基于假设检验原假设是：基因的方差与之前拟合的曲线相符（也就是对生物学角度的波动没什么贡献）。由此可以设置 adjusted p-values 的阈值，
适用于 HVGs 的确大部分都是与研究点相关，使用 FDR 只会让结果更准确；但缺点是：比第一种指定数量的方法更难以预测。既然是假设检验，就可能存在假阳性，可能得到更多不相干的基因。慎用。
保留之前拟合曲线上面全部的基因拟合曲线以下的基因都是不感兴趣的，因为它们的生物学偏差很小，于是把它们全部剔除，留下剩下的。为了偏差最小化，选择了噪音最大化（保留了很多与研究问题无关的基因）。

降维

每个基因现在都作为数据的一个维度存在。基因的表达决定了细胞的分布
PCA 的作用就是：对超过 4 维的数据降维到一个 2D/3D 平面图中，并且这个图中”相似相聚”
为什么可以通过 PCA 可以看批次效应？因为 PCA 图中的每个点都是一个样本，这个点中包含了大量的基因表达量信息；如果说本来属于生物学重复的几个样本在 PCA 图上离得很远，那么就意味着它们包含的基因表达量差异很大，这是不符合实际的，因此可能存在批次效应

SVD 是一种矩阵分解方法，相当于因式分解，目的纯粹就是将一个矩阵拆分成多个矩阵相乘的形式。

选择合适的 PC 数

如果选的 PC 多，可以避免丢弃一些生物信号，但同时又增加了噪音的风险
很多有经验的人会设置 PC 数在一个“合理”区间，例如 10-50 之间

利用 elbow plot 确定 PC 数基于的假设是：每个 PCs 都能捕获一些生物差异，而且前面的 PC 比后面的 PC 包含的差异信息更多，更有价值。就像“二八准则”，仅有 20%的变因操纵着 80%的局面，PCA 中也一样。于是当我们从头到尾观察每个 PC 的差异贡献时，会有一个明显的“落差”，“落差”之前的那些 PC 就是“二八准则”中 20%的变因。
利用技术噪音设置一个差异贡献阈值，将每个 PC 的贡献率加起来，直到达到这个阈值，保留其中的所有 PC。例如设置 80%，也就是保留能够解释 80%差异的 PCs。通过计算数据有多少比例的差异与生物因素相关，从而得到这个阈值。
利用细胞分群的推断
自定义 PC 数量

非负矩阵分解

非负矩阵分解是一个用来替代主成分分析的方法。它利用两个低维的矩阵（W 和 H，一个代表基因，一个代表细胞）来估计整个表达矩阵，并且其中所有的数据都要求是非负。它的想法和 PCA 类似，都是“缩小”数据，NMF 是利用小的矩阵来归纳大矩阵的主要特性，最终降噪、压缩数据。相比于 PCA 的坐标，NMF 坐标更容易解释，因为更大的值表示相应因子中基因的表达量更高。
相比于 PCA 的优点是：可以通过找到基因相关的小矩阵（W）每个因子中表达量高的基因，然后通过这些基因推测相应的细胞，再看细胞相关的小矩阵（H）中相应的高表达细胞，从而将细胞与细胞类型对应起来

聚类分群

分类：类别是已知的，通过对已知分类的数据进行训练和学习，找到这些不同类的特征，再对未分类的数据进行分类。属于有监督学习。
聚类：事先不知道数据会分为几类，通过聚类分析将数据聚合成几个群体。聚类不需要对数据进行训练和学习。属于无监督学习。

基于图形的聚类(Graph-based clustering)

最基础的想法是：首先构建一个图，其中每个节点都是一个细胞，它与高维空间中最邻近的细胞相连。连线基于细胞之间的相似性计算权重，权重越高，表示细胞间关系更密切。如果一群细胞之间的权重高于另一群细胞，那么这一群细胞就被当做一个群体 “communiity”
这种方法使得每个细胞都被强制连接到一定数量的相邻细胞，这减少了仅由一两个离群细胞组成的无意义群体的风险。

Graph-based 聚类使用

需要关注的问题：

构建这个图需要设置几个邻居？
用什么方法确定边的权重？
用什么检测 community 的方法来定义 cluster？

重要的参数是 k ，含义是：the number of nearest neighbors used to construct the graph。如果 k 设置越大，得到的图之间联通程度越高，cluster 也越大。
权重计算：type=”number” ：根据近邻数量计算；type=”jaccard” ：根据 Jaccard index of the two sets of neighbors 计算

当使用 graph-based 方法时，模块化（modularity）是一个评价 cluster 的指标。modularity 值越大，表示 cluster 内部的细胞与细胞之间连线（edge）数最多，内部分散效果越好，从而避免了不同 cluster 之间邻近的细胞形成连线。

k-means 聚类

通过不断迭代找出 k 个中心(centroid)，然后将各个细胞分配至最近的中心点。

k 值需要提前定义。假如真实有 10 种细胞类型，但我现在只定义 k=9，那么有两个细胞类型就会硬生生地”合二为一“。但其他方法，例如上面的基于图形的，即使分辨率设置的很低（意思就是看的很模糊），也是能看出这两种不是同一种类型
根据它的计算方法，需要不断重复多运行几次，保证得到的 cluster 结果是稳定的
既然 k-means 是找 k 个中心，那么有中心就会有区域半径，就会限定数据的范围。但实际上，很多分组并没有常规的数据大小和结构，也就是说，一组数据不会这么”乖乖地“落在我们画好的范围中

可以用 within-cluster sum of squares (WCSS) 方法对每个亚群进行评价。
基于 WCSS 又可以计算 root-mean-squared deviation (RMSD)，反映的是偏离平均位置的程度。就是用每个亚群的 WCSS 值除以该亚群的细胞数量然后再开方。
RMSD 值低的好。如果 RMSD 值低，就表示和其他亚群的混杂程度更小，更像是一个纯粹的亚群。

t-SNE 就是一个可视化的工具，不要试图量化 t-SNE 结果上的 cluster。上面的点不能说明什么问题，t -SNE 的算法会让稠密的点变膨胀，让原本稀疏的点变紧凑，不可以用 cluster 的大小来衡量亚群的异质性。另外它并没有义务帮我们保留 cluster 的相对位置，因此我们不能根据点的位置远近来确定 cluster 之间的相似程度，每运行一次点的位置都是变动的。

一般就是通过层次聚类来判断 cluster 之间的关系。

层次聚类

思路：贪婪地将样本聚集成簇，然后将这些簇聚集成更大的簇，依此类推，直到所有样本都属于一个簇为止。
层次聚类包含的算法也有很多种，差别在于如何层层聚集。
因为需要计算细胞间的距离矩阵所以运算量较大，一般只使用于小型的 scRNA 数据。

可以借助”轮廓图“（silhouette width）检查分群的质量。会对每个细胞都计算一个 silhouette width 值，如果一个细胞的 width 值为正并且越大，表示相对于其他亚群的细胞，这个细胞和它所在亚群中的细胞更接近，分群效果越好；如果 width 为负，就表示这个亚群的这个细胞和其他亚群的细胞更接近，即分群效果不太理想。

评价分群的稳定性

分群的稳定性：分群之前操作的小小波动，也不会影响到最后的分群结果；高稳定性的分群结果其实也方便了别人进行重复。scran 利用自助法（bootstrapping）进行评价。
自助法（Bootstrap Method，Bootstrapping，或自助抽樣法）是一种从给定训练集中有放回的均匀抽样，也就是说，每当选中一个样本，它等可能地被再次选中并被再次添加到训练集中。
细胞有放回的抽样，得到一个”自助重复“数据集，然后用他进行聚类，看看是否能得到相同的亚群。
自助法（Bootstrapping）是一个通用的评价 cluster 稳定性的方法，适用于各种聚类算法。
高稳定性的亚群不一定是高质量的，因为如果一个亚群质量一直很差，它也很稳定。

Marker 基因检测

导致 cluster 出现差异的那部分基因中，尤其是那些差异最显著的基因中，最可能包含 marker 基因。因此一个首要任务就是去进行 cluster 之间的差异分析，最后把每个 cluster 中最显著的前多少基因拿出来，放在一起。

findMarkers() 会对两两比较结果做一个排名，然后选择 p 值比较显著的一些作为 Top 基因，返回的结果包括了所有 cluster 的 logFC 情况，比如这里有 18 个 cluster，那么返回的结果也包含 18 个 cluster 的 logFC。以上调差异表达为例：某个基因在 cluster9 与 cluster1 相比下上调，但这个基因在 cluster9 和 cluster2 中不上调，依然会把这个基因列出来。

寻找 cluster 特异的 marker 基因,只选在某个 cluster 与其他各个 clusters 相比都差异表达的基因。

细胞类型注释

最常见的就是通过去了解参与特定生物学过程的经过验证的基因，利用 GO、KEGG 数据库；另外还可以将表达谱与之前发布的参考数据集做对比，因为参考数据集中样本或细胞的生物学状态都经过证实。

使用参考数据集

SingleR会计算数据中的细胞与参考数据集中的细胞之间 Spearman 相关性，找到最相关的属性赋予细胞。为了减少数据的噪音，SingleR 包先鉴定 marker 基因，然后根据 marker 基因的表达量去计算相关性。

使用 marker 基因集

高表达的 marker 基因与每个细胞的生物学功能相关，因此可以用 marker 基因作为桥梁，连接测试数据集与参考数据集

基于 marker 基因的富集分析

先得到每个 cluster 的标志性 marker 基因，然后对这些代表该 cluster 的基因们进行富集分析，看看它们集中在哪些通路

数据集整合后的批次矫正

一般不同批次的数据前期处理分开进行，可以根据各自的特点先对数据进行过滤。
常规的前期处理步骤还是：质控=》归一化=》找 HVGs=》降维=》聚类

检查批次效应的严重性

最简单的方法是把两个数据合并，然后 PCA+t-SNE。

矫正批次效应之线性回归

对每个基因表达量拟合一个线性模型。例如 limma 的 removeBatchEffect()、sva 的 comBat() 。如果要使用这类方法，就需要假设：批次间的细胞组成相同。另外的一个假设是：批次效应的累积的，对于任何给定的基因，在不同亚群中经过任何因素诱导的表达变化倍数是相同的。

MNN 矫正

Mutual nearest neighbors，是对一个批次的细胞找到它们在另一个批次中对应的最相邻细胞集合。
简单解释就是：batchA 中有一个细胞 a，然后想根据挑选的特征基因（如 HVGs）去在 batchB 中鉴定与 a 细胞空间相近的细胞们。同样的，对 batchB 中的细胞 b 重复这个过程，也鉴定它在 batchA 中的近邻。
原理：在进行批次矫正前，MNN 找到的一对对细胞都是生物学状态相似的，因此如果再有不同，那么就姑且认为是外部的批次效应导致的，也就方便了去除。
与线性回归方法相比，MNN 方法不会假设细胞群组成相同或者事先已知。MNN 会自己学习细胞群的结构并进行估计。

如果要探索基因表达相关，最好还是要在未矫正数据基础上进行，并且还要把批次信息封锁掉。

多样本间差异分析

多样本的 scRNA 差异分析，主要分成两大类：差异表达（differential expression）和差异丰度（differential abundance），前者检测相同类型细胞在不同条件下表达的变化，后者检测不同条件下细胞类型(或状态等)组成的变化。

不同处理间差异基因表达分析 DE analysis

差异分析基于 edgeR 的 quasi-likelihood (QL)方法，使用负二项广义线性模型(NB GLM) 来处理过度分散的表达量数据。

差异细胞丰度分析 DA analysis

检测不同条件下细胞类型(或状态等)组成的变化；如果说 DE analysis 重点在基因的表达量变化，那么 DA analysis 就着眼于细胞数量的变化。

细胞周期推断

许多细胞周期相关的关键节点（比如通过检查点）属于转录后调控，在转录组数据中不可见，但通过表达量的变化也可以进行初步的推断。
Cell cycle

先找细胞周期蛋白基因，再对照数据集

细胞周期蛋白（cyclins）控制着整个细胞周期的进展，并在细胞周期各阶段的表达模式具有一定特点。

cyclin A 在 S 和 G2 期表达；
cyclin B 在 G2 后期和 M（mitosis 有丝分裂期）表达量最高
cyclin D 在整个过程都表达，但在 G1 期表达最高；
cyclin E 在 G1/S 过渡期表达最高

先找 marker 基因，再对应细胞周期蛋白

可以利用常规的寻找差异基因的方法，找找有没有基因在细胞周期蛋白中上调，来帮助判断某个 cluster 中是否有更多细胞属于某个细胞周期。

辅以参考数据

使用别人已经做好的细胞周期注释数据，来辅助我们自己的数据进行注释。

用已发表的分类器

scran 包的 cyclone 分类器。如果 G1 分数高于 0.5 并且大于 G2/M 的分数，那么就判断为 G1 期；同样的方法可以判断 G2/M 期；如果二者分数都不大于 0.5，那就是 S 期。

移除细胞周期导致的差异

如果要去掉这部分影响，可以先按照上面的方法推断出细胞周期，然后将每个周期当成一个独立的批次，再进行批次矫正。最简单的方法是使用线性模型（如 regressBatches()）来移除这个影响。
这个操作并非常规的分析流程。因为大部分数据中，细胞周期导致的差异是次要的，不会超过细胞真实生物差异的影响。

细胞轨迹推断

做这个分析之前，最好先问几个问题：
确定数据会体现发育轨迹吗？也就是研究的样本是不是和发育相关的？
数据中的细胞会体现出中间态吗？
是否认为轨迹会出现分支？
并且要注意：
任何数据都可以强行画出轨迹，但不一定都有生物学意义！
先要保证目前找到的 HVGs 和降维结果符合我们的预期，才能继续向下分析