3′与5′单细胞RNA测序中协议依赖性基因表达偏差的表征与缓解策略

《Frontiers in Bioinformatics》：Characterizing and mitigating protocol-dependent gene expression bias in 3′ and 5′ single-cell RNA sequencing

： 单细胞RNA测序（scRNA-seq）已实现细胞异质性的大规模表征；然而，整合通过不同文库制备方案生成的数据集仍具挑战性。例如，10X Genomics 3′与5′化学法之间的比较因转录本末端捕获和扩增差异所引入的协议依赖性技术偏差而变得复杂。尽管归一化及通常的批次校正（batch correction）是scRNA-seq数据集预处理中的必要步骤，但对于可靠的跨协议比较而言，何种校正方法最为合适甚至是否必要，目前仍不明确。本研究利用来自36名受试者、涵盖六个数据集及多种组织的匹配3′与5′ scRNA-seq数据，系统表征了与协议相关的表达差异。研究发现，基因表达差异并非遍及整个转录组，而是由一组相对较小且可重复的协议偏差基因所驱动。这些协议特异性差异并非主要由GC含量或基因长度解释。排除受影响基因可改善跨协议一致性，表明大多数基因无需激进校正即可直接比较。随后，研究人员对十一种常用归一化方法进行了基准测试，结果显示，尽管fastMNN等方法在细胞群体匹配良好时可改善统计对齐，但在细胞类型重叠不完全的生物学现实场景中，这些方法可能扭曲基因水平信号并夸大差异表达。综合来看，本研究结果表明3′与5′ scRNA-seq之间的协议偏差范围有限，针对偏差基因的靶向处理为跨协议数据整合与比较提供了一种替代归一化或批次校正策略的新途径。

一、研究背景与问题

单细胞RNA测序（single-cell RNA sequencing, scRNA-seq）是近年来发展迅速的基因表达分析技术，能够在单细胞分辨率下解析组织与疾病状态中的细胞异质性，广泛应用于细胞分化、免疫应答及肿瘤异质性等关键生物学问题的研究。当前scRNA-seq技术需要单细胞分离、mRNA捕获与条形码标记、cDNA合成、片段化及扩增等步骤，最终生成约300碱基对长度的条形码化转录本片段。由于在逆转录过程中需添加唯一分子标识符（Unique Molecular Identifiers, UMIs）以降低扩增噪声，条形码仅附着于每条转录本的一端，即3′端或5′端，这意味着只有包含相应末端的片段才会被扩增和测序。

最初，scRNA-seq方案将条形码置于用于捕获mRNA分子的poly(dT)寡核苷酸上，UMIs因此附着于转录本的3′端。然而，这种方法在解析T细胞受体（T-cell receptors, TCRs）和B细胞受体（B-cell receptors, BCRs）重排及可变区方面存在局限，因为只能测序到恒定区。为实现在单细胞分辨率下表征TCR和BCR的重排及可变区，方案被修改为将UMIs和细胞特异性条形码置于转换寡核苷酸上，从而附着于转录本的5′端。因此，许多较新的scRNA-seq研究，尤其是聚焦免疫景观与免疫应答的研究，采用5′ scRNA-seq方案。然而，这种方法因3′与5′技术在基因表达测量上的技术偏差，使全面图谱构建和荟萃分析的开展变得复杂。

已有基准研究评估了不同计算方法在共享低维空间整合多种方案的效果，但这些方法仅局限于改善聚类、轨迹推断等仅利用scRNA-seq数据低维表征的算法，无法实现不同细胞群体或生物学条件间基因表达水平的比较。截至目前，尚无系统性研究评估可用于校正跨scRNA-seq方案技术差异的基因表达空间归一化方法。尽管有研究表明某些技术效应可通过回归掉基因长度等简单调整部分校正，但这些方法因细胞质转录后调控和解离相关应激反应等混杂生物学效应，无法完全调和scRNA-seq与单核RNA测序（snRNA-seq）数据。相比之下，10X Genomics 3′与5′ scRNA-seq方案间的差异被认为主要源于逆转录和下游扩增过程中转录本末端捕获的技术差异，这些协议相关效应不反映细胞间真实的生物学差异，因此可能适合使用现有归一化和批次校正方法进行校正。然而，校正的必要程度以及何种方法能在不扭曲生物学意义信号的前提下最有效地调和3′与5′ scRNA-seq数据，仍不明确。

二、研究内容与结论

本研究系统表征了3′与5′ scRNA-seq之间的协议依赖性基因表达差异，并评估了现有归一化和批次校正方法的性能。利用来自六个公开数据集、36名受试者的匹配3′与5′ scRNA-seq数据，研究人员首先鉴定了一组在这些化学法之间表达持续偏差的基因，并考察了其复现性、表达分布、生物学功能、细胞类型特异性及转录本特征。随后，研究人员评估了十一种选定校正方法在缓解上述偏差方面的能力，通过比较校正后3′与5′数据中细胞类型特异性表达谱的相似性、评估校正后3′与5′样本的混合程度，以及检验细胞类型特异性标记基因的一致性来评价性能。最后，研究人员在旨在整合部分重叠3′与5′数据集的合成用例中测试了表现最佳的方法。结果表明，排除受协议差异强烈影响的基因是所测试方法中最有效的策略，可减少3′与5′ scRNA-seq数据集间比较表达水平时的假阳性发现。

三、关键技术方法

研究人员从CELLxGENE数据库获取了六个包含36名受试者匹配3′与5′ Chromium scRNA-seq数据的人类数据集，涵盖胸腺、骨髓、肝脏、脾脏、肾脏、皮肤等多种组织。数据预处理使用Seurat（v5.1.0）进行文库大小归一化至每细胞10,000计数并自然对数转换，通过vst模型鉴定前3,000个高变基因（Highly Variable Genes, HVGs）。协议偏差基因通过Wilcoxon秩和检验对每个受试者进行差异表达分析，最终在至少五个受试者中以1e^?50阈值显著的867个基因被选为共识偏差基因集。研究还使用了七个额外的匹配3′/5′验证样本评估偏差基因集的泛化性。归一化方法包括SCTransform回归与独立模式、M3Drop、mnnCorrect、fastMNN、Scanorama、ComBat、limma、基于管家基因的Z变换、scVI及scArches。性能评估采用Pearson相关系数、余弦相似度、均方误差（MSE）、欧氏距离和Jensen-Shannon散度（JSD）等统计指标，以及UMAP可视化、Seurat MixingMetric和标记基因重叠评分等生物学指标。此外，研究人员通过构建部分重叠（25%重叠）和完全重叠（100%重叠）的合成场景，使用马修斯相关系数（Matthews correlation coefficient, MCC）评估各方法恢复细胞类型标记基因的能力。

四、研究结果

3.1 36个匹配样本中10X Genomics 3′与5′协议偏差基因的探索

研究人员对36名具有匹配10X 3′和5′ scRNA-seq数据的受试者进行差异基因表达分析，在不同显著性阈值下鉴定出609至2,598个在平台间持续偏差的基因。移除各偏差基因集均改善了匹配3′与5′数据在低维空间的对齐，UMAP显示原始数据中明显的检测法相关分离在排除偏差基因后大幅减少，细胞类型一致性得到改善。MixingMetric在移除仅609个基因后显著增加。在表达水平上，609和867基因偏差集在匹配3′与5′细胞类型谱间的余弦相似度显著低于随机基因集，表明偏差基因捕获了协议相关表达差异的绝大部分。867基因集因在较小基因集导致对齐不足和较大基因集过度移除之间取得平衡而被选为最终偏差基因集。在七个未参与发现过程的验证比较中，移除预定义的867基因在所有验证比较中均改善了3′/5′混合度，且偏差基因集显示出比表达匹配随机基因集更低的3′/5′余弦相似度，表明这些基因在保留数据中捕获了复发性协议相关表达差异。偏差基因倾向于在多个数据集中高表达，75%位于大多数数据集的最高表达四分位，并富集于管家功能，包括细胞翻译和氧化磷酸化。偏差基因在GC含量模式上无差异，表明协议差异并非主要由针对GC富集序列的PCR扩增偏差驱动。偏差基因略短于预期，但这可能由其管家功能富集所解释。内含子读数纳入将867基因集中显著基因数从525个减少至390个，主要减少了5′文库中高表达显著基因的数量，但大多数偏差基因仍保持显著，表明内含子读数纳入减弱了部分偏差基因调用，但不足以解释所有复发性3′/5′协议偏差。

3.2 十一种常见scRNA-seq批次校正方法缓解10X 3′/5′协议偏差的基准测试

研究人员比较了十一种常用归一化方法在校正3′和5′ scRNA-seq基因表达值同时保留生物学变异方面的性能。所有归一化方法均独立应用于36名受试者各自的匹配3′和5′样本，以避免受试者间遗传异质性的混杂效应。

3.3 校正后3′与5′相似性的定量测量

通过比较校正后3′和5′ scRNA-seq数据中先前鉴定的偏差基因和3,000个高变基因的细胞类型特异性表达来评估性能。ComBat和mnnCorrect在提高3′和5′表达间的Pearson相关性及降低JSD和欧氏距离方面持续产生最大改善。基于管家基因的Z变换表现最差，无证据表明其改善3′和5′一致性。fastMNN在欧氏距离和JSD方面表现相似但Pearson相关性表现可变。limma表现与ComBat相似但相关性较弱。Scanorama改善了Pearson相关性、余弦相似度和JSD，但未持续优于mnnCorrect或fastMNN。scVI和scArches在相关性和JSD方面与其他方法相当，但显著增加了3′和5′间的欧氏距离。M3Drop和SCTransform在偏差基因集上的表现优于高变基因，尤其在欧氏距离和MSE方面。SCTransform的“split”选项（独立拟合各样本回归）在所有指标和基因集上均显著优于用于校正3′与5′差异的回归模型方法。

3.4 校正性能的生物学测量

研究人员评估了每种归一化方法在低维空间对齐3′和5′数据集的能力。表现最佳的方法为ComBat、fastMNN和M3Drop，它们实现了3′和5′ scRNA-seq数据集的完全重叠且不干扰细胞类型聚类。Scanorama在从校正后细胞×基因表达矩阵生成UMAP时表现可变，但从其嵌入生成UMAP时在多个数据集中显示出更强的3′/5′对齐，表明Scanorama作为潜空间整合方法优于作为下游基因水平分析的校正表达值来源。limma、Z变换和SCTransform回归未能整合3′和5′数据集。所有方法均能保留细胞类型区分，在十个数据集中的九个中细胞类型标签的MixingMetric为0。在细胞类型特异性标记基因一致性方面，fastMNN、scVI和scArches达到最高标记重叠，但这可能由显著标记基因数量膨胀6至10倍所驱动。Scanorama产生的3′/5′相交标记基因少于任何其他方法。有趣的是，移除867个偏差基因对细胞类型特异性标记基因鉴定的一致性影响可忽略不计，同时显著减少了低维空间中的批次效应。

3.5 校正3′和5′基因表达的计算成本

大多数归一化方法扩展性良好，所需内存低于20 GB，在两个AMD EPYC 9654 CPU上运行时间低于5分钟。mnnCorrect例外，需要30.55 GB内存和51分钟运行时间。深度学习方法scVI和scArches需要更长运行时间，但在不同硬件上运行。

3.6 应用批次校正方法整合部分重叠的3′和5′数据集

研究人员在三个选定数据集（胸腺、肝脏和骨髓）中引入人工稀疏性，使3′和5′样本间仅约25%的细胞共享，部分细胞类型仅存在于单一检测法中。在稀疏（25%重叠）和完全重叠设置中，移除867个偏差基因后的对数归一化数据产生最高MCC分数，其次是未校正的对数归一化表达。这表明在此背景下，批次校正可能无法改善甚至可能损害标记基因检测的准确性。在校正方法中，fastMNN在25%重叠场景中表现突出，达到次高MCC分数，但仍低于简单对数归一化。fastMNN仅略微增加DE基因总数，表明其具有良好的假阳性控制。Scanorama和limma在25%重叠设置中表现相似并膨胀了DEG数量，scVI的膨胀最为强烈。校正方法在此设置中的不佳表现并非归因于所选地面真值标记基因在3′和5′检测法间缺乏差异，因为在所有三个案例中至少10%的地面真值标记与偏差基因集重叠。校正方法确实显著改变了这些基因的表达，但未能准确移除3′与5′技术引入的特异性偏差。fastMNN和scVI在校正之外还平滑和插补了表达值，成功保留了差异明显的基因的DE模式，但有时导致分布重叠基因的过度平滑。fastMNN和Scanorama还独特地压缩了基因表达值范围。limma和M3Drop倾向于以与预期DE模式相矛盾的方式改变零值，导致FOXP3和KLRD1等基因失去显著性。

五、讨论与结论

3′和5′ scRNA-seq化学法均被广泛使用，5′方案因能捕获重排TCR和BCR序列常用于免疫组库分析，而3′方案因其成本效益和历史使用在大规模单细胞图谱中占主导。公开可用的scRNA-seq数据集日益包含3′和5′化学法的混合，为跨协议比较和荟萃分析创造了机遇与挑战。大多数现有批次校正方法仅在降维空间中对齐数据集，而非直接校正基因水平表达值。本研究通过表征协议特异性表达差异并在生物学现实场景中评估归一化和校正策略，填补了这一空白。

研究发现大多数协议诱导变异由相对较小的基因子集驱动，这被认为源于转录本特异性捕获和扩增效率差异。偏差基因主要相对于3′数据在5′中上调，且其5′和3′末端间GC含量差异不大于无偏基因，不支持GC依赖性扩增偏差作为主要解释。内含子与外显子定量也无法完全解释观察到的协议偏差基因。由于偏差基因的身份和偏差方向在不同数据集中高度保守，研究人员预期该偏差源于独立的技術效应，且由于10X 3′与5′化学法除转录本末端捕获和条形码/UMI掺入设计外使用相同试剂，这不太可能是特定方案的结果，而是基因本身的特征。预期该偏差在其他单细胞RNA测序平台中保守，但物种间保守性较低。

在统计基准中，fastMNN、mnnCorrect和ComBat持续改善匹配3′和5′表达谱间的一致性，与先前显示互近邻方法和线性批次校正方法在中等批次效应和共享细胞群体设置中表现稳健的基准研究一致。mnnCorrect虽常实现略强于fastMNN的对齐，但计算成本显著更高。ComBat在统计指标上表现良好，但其全局线性调整的依赖性限制了其在细胞类型不平衡或批次间重叠不完全场景中的适用性。SCTransform变体行为存在一致差异，“split”方法优于回归模型。在生物学相关场景中，性能发生转变：fastMNN保持稳健的标记恢复，而scVI或M3Drop等统计指标相对较好的方法显示基因检测膨胀或表达动力学扭曲。重要的是，在分析前排除867个协议偏差基因的小集合可显著改善下游对齐，在某些场景中可能提供比复杂校正流程更可靠的替代方案。