如何优化一个传统分析方法还发了14分

Benchmarking principal component analysis for large-scale single-cell RNA-sequencing大规模单细胞RNA测序的基准主成分分析

一. 研究背景

单细胞RNA测序(scRNA-seq)技术的出现使许多类型的细胞异质性检查成为可能。例如，由各种组织组成的细胞亚群，稀有细胞和干细胞壁，与细胞周期进程有关的连续基因表达变化，空间坐标， 许多scRNA-seq研究已经捕获了分化成熟度的差异。由于细胞异质性的测量高度依赖于同时测量的细胞数量，因此已经开发了各种各样的大规模scRNA测序技术，包括那些使用细胞分选设备的技术，Fludigm C1 ，基于液滴的技术和ingle-cell组合索引 RNA测序(Sci-RNA-seq)。主成分分析(PCA)是一种适用于许多情况的常用UML算法。尽管其用途广泛，但仍有很多原因导致不清楚如何对大规模scRNA-seq进行PCA。

而主成分分析(PCA)是分析单细胞RNA序列的基本方法(scRNA-seq)数据集，但是对于大规模的scRNA-seq数据集，计算时间长且消耗大量内存，而一些快速PCA算法可能会忽略一些重要的差异基因表达模式。因此，需要一些快速且高效存储的PCA算法，实际的计算时间和存储效率在很大程度上取决于特定的实现，包括编程语言，加载输入文件的方法以及数据格式。在本研究中，作者检查了用于大型scRNA-seq数据集的快速且内存高效的PCA算法的实用性。

二. 分析流程

三. 结果解读

作者先回顾了一下现有的PCA算法及其实现(图1)，图1是对这项工作中的基准测试概述。

图1a. 这项工作的示意图

图1b. 使用实际数据集进行基准测试的评估指标

图1c. 使用综合数据集进行基准测试的评估指标

1. 审查PCA算法和实现

PCA被广泛用于数据可视化，数据质量控制，特征选择，去噪，插补，批量确认和删除效应，细胞周期效应的确认和估计，稀有细胞类型检测，细胞类型和细胞状态相似性搜索，伪时间推断， 空间重建和其他非线性降维，聚类方法的输入等，并节省内存空间。各种各样的scRNA-seq数据分析工具都包括PCA作为内部功能，或者利用PC分数作为下游分析的输入。

作者回顾了现有的PCA算法和实现，并将算法分为六类，即基于相似度转换(SimT)，基于下采样(DS)，基于奇异值分解(SVD)基于更新(SU) ，基于Krylov子空间(Krylov)，基于梯度下降(GD)和基于随机投影(Rand)。作者列出了21种PCA实施方案(包含10种算法)，表1总结了已审查的PCA实施与scRNA-seq研究的对应关系。

表1. scRNA-seq研究中PCA实施的用例

2. 真实数据集

考虑到在有限的时间，计算资源和人力下评估的大量方法之间的权衡，作者选择了真实数据集进行基准测试。最新的scRNA-seq方法分为两个类别，即全长scRNA-seq方法和具有特定细胞解离和细胞/分子条形码技术的高通量scRNA-seq方法，例如基于液滴的和分离的池实验。因为通过后一种技术增加了用scRNA-seq测量的细胞数量，所以作者选择了以下四种通过这种技术生成的数据集：人外周血单核细胞(PBMC)，人胰腺细胞(Pancreas)，小鼠脑和脊髓 (BrainSpinalCord)，以及来自皮层，海马和心室区的小鼠细胞(Brain)(表2)。

表2. 用于基准测试的真实数据集

3. PCA算法的准确性

在这里，作者通过使用四个实际数据集来评估各种PCA算法的准确性。为了分析PBMC和Pancreas数据集，作者将prcomp的结果设置为黄金标准， 将其他实现与该结果进行比较(1b和图2)。

需要说明的是，对于Brain-SpinalCord和Brain数据集分析，LAPACK的完整SVD在计算上很困难。当无法定义基本事实时，作者将方法相互比较。因此，作者只是使用几种不同的标准(例如特征值的大小和聚类精度)相互比较了方法的结果。

其中，图2为其他实现与prcomp结果进行比较结果汇总：图2a为作者的文献综述总结了理论性质；图2b为prcomp结果与每个实现相关的属性；图2c通过对真实数据集和综合数据集进行基准测试评估展示性能；图2d为通过一些指标评估的用户友好性。

图2. 其他实现与prcomp结果进行比较

首先，作者对每种PCA算法的结果进行了t随机邻居嵌入(t-SNE)和统一流形逼近和投影(UMAP)，并对每种PCA算法的结果进行比较，并比较了原始研究检测到的簇结构的清晰度(图1b和3)。

图3比较t随机邻居嵌入(t-SNE)图。使用经验数据集执行的多个主成分分析(PCA)比较：PBMC(102个细胞)，Pancreas (103个细胞)，BrainSpinalCord(105个细胞)和Brain数据集(106个细胞)。使用每个PCA实施的结果执行t-SNE。由原始研究确定的一些不同的聚类被错误地合并为单个聚类(红色圆圈聚类)。其中：IncrementalPCA，主要是为了解决单机内存限制的。

图3. 比较t随机邻居嵌入(t-SNE)图

作者还对PCA实施的所有结果执行了Louvain聚类(图1b和图4)，ARI值表明，与黄金标准或其他实现相比，下采样和sgd(OnlinePCA.jl)的结果更差。

其中，图4通过Louvain聚类结果的调整后的Rand指数(ARI)评估聚类准确性。针对PBMC(102个细胞)，Pancreas (103个细胞)，BrainSpinalCord(105个细胞)和Brain数据集(106个细胞)执行了多个主成分分析(PCA)，并对PCA结果进行Louvain聚类。对于每个PCA结果，Louvin聚类进行十次，并计算平均值，分类标签与相应的原始论文相同。

图4. 聚类精度比较

接下来，作者对金标准的PC与其他PCA实现方案进行了全面比较(图1b和5a)。因为PC是单位向量，所以当两个PC指向相同或相反方向时，它们的叉积分别为1或− 1。作者进行了PCA优化，为所有比较计算了叉积的绝对值，范围为0到1，并评估了是否准确计算出了对应于较低特征值的较高PC。得出结论：特征向量，加载矢量中，下采样，orthiter / gd / sgd(OnlinePCA.jl)和PCA(dask-ml)都会随着PC尺寸的增加而变得不准确，PCA实施不正确可能会降低加载矢量的准确性。

其中，图5a计算了金标准方法的特征向量与其他主成分分析(PCA)实现的特征向量之间所有组合的叉积的绝对值。该值越接近1(黄色)，则两个对应的特征向量彼此越接近。如果两个PCA结果相等但不考虑符号差异，则该图中的矩阵将变为单位矩阵。图5b计算了金标准方法的前500个最大绝对值元素与其他PCA实现的相交元素数量。该值越接近500(黄色)，则两个相应的加载向量彼此越接近。如果两个PCA结果相等但不考虑其符号差异，则该图中矩阵的所有对角元素变为500。

图5. 特征向量所有组合，所有加载向量组合的比较

下采样，IncrementalPCA(sklearn)和sgd(OnlinePCA.jl)的特征值的分布也与其他实现不同(图6)。其中，图6a表示PBMCs，图6b表示Pancreas，图6c表示BrainSpincalCord，图6d表示 Brain；x轴表示每个主成分分析(PCA)实现的主要成分，y轴表示对应的特征值。

图6. 四个真实数据集特征值分布比较

4. 计算时间，内存使用率和可伸缩性

作者比较了所有PCA实现的计算时间和内存使用情况(图7)。对于Brain-SpinalCord数据集，降采样本身比大多数PCA实现要快，但是其他预处理步骤速度较慢；对于Brain数据集，下采样变得比大多数PCA实现要慢，并且随着数据矩阵大小的增加，这种趋势很明显，因为下采样基于LAPACK中的完整SVD。

其中，图7a展示了诸如二进制化和规范化之类的预处理步骤的经过时间(橙色条)以及每次PCA计算本身的经过时间(绿色条)。仅当对Brain数据集执行PCA实现时，作者才使用内部的Julia脚本进行预处理。这是因为无法将此数据集作为R语言的data.frame加载到内存空间。图7b展示了每个经验数据集计算的所有主成分分析(PCA)实现的内存使用情况(蓝色条形)。

图7. 比较经验数据集的经过时间和最大内存使用量

5. 文件格式与性能之间的关系

在辅助文件中，作者计算了oocPCA_CSV(R，oocRPCA)，IncrementalPCA(Python，sklearn)和orthiter / gd / sgd / halko / algorithm971(Julia，OnlinePCA.jl)等核心实现中大脑矩阵的通过次数。在oocPCA_CSV(R，oocRPCA)和IncrementalPCA(Python，sklearn)中，数据矩阵以CSV格式传递给这些函数，在其他核心实现中，首先将数据矩阵二进制化并压缩得到Zstd文件格式。结果显示计算时间与实施的通过次数相关，二进制化和数据压缩大大加快了计算时间。并且，稀疏矩阵的使用可能加快scRNA-seq数据集的计算时间。

6. 用户指南

基于所有基准测试结果以及作者在这项工作中的实施，作者提出了一些用户指南(图8)。考虑到生物信息学研究结合了多种工具来构建用户的特定工作流程，因此编程语言是选择正确的PCA实施方案的重要因素。因此，作者根据语言对PCA实施进行了分类(即R [111]，Python [112]和Julia [113]；图8，按列)。除了数据矩阵的大小，作者还根据实现加载数据的方式(内存中或内核外)以及输入矩阵格式(密集或稀疏，按行，图8)对实现方式进行了分类。

其中，图8推荐的PCA实现基于书面语言和矩阵大小进行分类，还描述了每个PCA实施的推荐参数(红色)。在这里，作者将数据矩阵的GC值定义为基因数×细胞数。

图8. 用户指南

小结

在这项工作中，作者回顾了现有的快速且高效存储的PCA算法和实现，并评估了它们在大规模scRNA-seq数据集中的实际应用。这项工作提供了四个关键贡献：

• 首先，作者回顾一下现有的PCA算法及其实现；

• 其次，作者介绍了使用所选PCA算法和实现的基准测试；

• 第三，作者为质量控制(QC)，过滤和特征选择提供了一些PCA算法和实用程序功能的一些原始实现。所有命令都在快速且高效内存的Julia包中实现；

• 最后，作者为最终用户和软件开发人员提出准则作者的基准测试表明，以选择合适的PCA实施方案，一些基于Krylov子空间和随机奇异值分解的PCA算法比其他算法更快，内存效率更高且更准确。