LolliPop：利用时间序列正则化解码废水测序数据，追踪SARS-CoV-2基因组变体动态

《PLOS Computational Biology》：Tracking SARS-CoV-2 genomic variants in wastewater sequencing data with LolliPop

这是一篇关于开发新方法（LolliPop）以解决废水测序（WBE）中关键挑战——病毒变体反卷积问题的研究论文。作者提出了一种针对时间序列数据进行正则化处理的统计方法，能够从存在大量缺失数据（如覆盖率极低的样本）的混合样本中，稳健地估算SARS-CoV-2等病原体基因组变体的相对丰度，并提供了高效计算置信区间（CI）的方法。该研究在瑞士国家废水监测项目和模拟数据中验证了其方法的准确性和实用性，为大规模废水流行病学监测提供了高效可靠的计算工具。

　　

方法

变体反卷积

研究将废水测序数据中观察到的突变频率建模为变体谱的线性组合。具体而言，假设在时间点t，携带特定突变的读数比例期望值等于携带该突变的变体相对丰度之和。由此定义了变体反卷积问题：在给定变体定义矩阵B和突变频率时间序列y_t的条件下，求解每个时间点所有变体的相对丰度π_t。

时间正则化

针对废水数据噪声高、缺失值多（常导致设计矩阵秩亏）的问题，研究引入了时间正则化。该方法基于变体丰度具有时间连续性的假设，对相邻时间步的相对丰度向量之间的差异施加融合岭（fused ridge）惩罚。研究发现，这种惩罚在数学上等价于核平滑（kernel smoothing），使用高斯核（Gaussian kernel）及其带宽超参数。时间正则化使得丰度估计对极高水平的缺失数据具有鲁棒性。

求解反卷积问题

研究使用Python科学计算库Scipy中的优化例程来求解带正则化的损失最小化问题。损失函数选用软L₁损失（SL1），它在最小二乘（LS）损失和最小绝对偏差损失之间进行插值，以平衡鲁棒性和效率。

置信区间

为了提供决策所需的估计不确定性，研究提出了两种计算置信区间（CI）的策略。一是基于对突变位点进行重采样的自助法（bootstrap），该方法概念直接但计算量大。二是基于准二项分布（quasibinomial）假设推导的渐近标准误差，并在对数比（logit）尺度上构建Wald型置信区间，以避免区间超出[0,1]范围，同时通过过离散（overdispersion）因子进行调整以更准确地捕捉读数间的依赖性。与自助法相比，该方法能提供相似的结果但计算速度快约30倍。

实施与应用

上述方法被集成到Python软件包LolliPop中，可从GitHub和Bioconda获取。研究使用的瑞士废水监测数据来源于2021年1月至9月期间从六个主要污水处理厂（WWTP）收集的1295个NGS数据集，并使用V-pipe 3.0进行处理。变体定义基于Cov-Spectrum数据库中临床序列的突变谱。

结果

与临床数据比较

研究将LolliPop从废水数据推断的变体相对丰度时间序列与基于临床测序的估计值进行了比较。即使在病毒发病率极低、测序覆盖率极差（例如许多月份仅有个位数比例的扩增子被覆盖）的时期，该方法仍能准确地进行反卷积。废水推断的动态与临床测序结果高度一致（）。通过加权交叉相关分析发现，废水信号平均领先临床信号约3.3天，相关性高达0.91。

置信区间

研究展示了通过自助法（）和Wald法（）计算的置信区间。两种方法得到的区间有很大重叠（Wald区间平均覆盖了约89%的自助法区间），但Wald区间通常更保守。在SARS-CoV-2 RNA浓度较低的时期，不确定性一致性地更高，反映了数据质量的影响。

模拟研究

研究通过模拟数据评估了方法对缺失数据和变体相似性的鲁棒性。模拟了多种场景，包括Delta变体取代Alpha变体、多个高度相似的Omicron亚变体共循环、模型过指定（包含不存在的变体）和模型欠指定（忽略存在的变体），以及由人工设计的、具有最大重叠突变的5个高度相关变体的混合。

结果表明，在没有正则化的情况下，缺失值水平的增加会严重影响反卷积的准确性。然而，引入时间正则化后，即使在缺失数据比例极高（>90%）的情况下，仍能以高精度（R²> 0.9）恢复真实的变体丰度信号（）。

在模型过指定的情况下，当覆盖率中等或较高时影响不大，但在低覆盖率下可能导致零星假阳性检测。模型欠指定则会系统性地导致被忽略变体的丰度被错误地归因于相关变体。与时间正则化相比，从LS损失切换到SL1损失对鲁棒性的提升较小，仅在无正则化时略有改善。时间正则化虽然会在时间序列的起点和终点引入少量平滑偏差，但显著降低了估计方差，总体上提升了准确性。

超参数

研究评估了平滑带宽（σ）和损失尺度参数（α）的影响。增加带宽最初能降低均方根误差（RMSE），但过高的值会因过度平滑而增加偏差。在存在高比例缺失数据时，使用非零带宽的益处尤为明显。损失尺度参数对RMSE的影响微乎其微。在真实废水数据上分析发现，增加带宽可提高与临床数据拟合的R²，而损失尺度参数在0.1到0.25之间时拟合效果最佳。

运行时间

在一台标准笔记本电脑上，使用LS损失函数对全部1295个废水样本（225个时间点，6个地点，5个变体，94个突变）进行反卷积耗时约2秒。使用SL1损失函数则耗时约1分钟。计算重新参数化的Wald置信区间耗时约1分钟，而基于1000次自助样本构建置信区间则需约30分钟。

讨论

研究开发的LolliPop方法有效解决了废水测序中的变体反卷积问题。其核心优势在于通过时间正则化处理，使方法对废水数据中常见的高噪声和大量缺失值具有极强的鲁棒性。关键超参数是控制平滑程度的带宽σ，即使设置很小的非零值也能显著提升估计的稳健性和准确性。另一个超参数损失尺度α影响较小，为兼顾计算效率与鲁棒性，推荐使用LS损失或α值在0.1至0.25之间的SL1损失。

该方法虽然基于SARS-CoV-2数据开发，但其原理具有普适性，已成功应用于流感A和呼吸道合胞病毒（RSV）的废水监测中。其计算效率高，适合大规模国家监测项目。

LolliPop的局限性在于依赖预定义的变体谱进行反卷积，新兴变体需要外部临床数据确定其特征性突变谱后才能被有效追踪。因此，将其嵌入集成的基因组监测系统中效果最佳。变体选择是关键步骤，COJAC工具可用于在反卷积前确认变体的存在。此外，当前实现将基因组位置视为独立，可能降低对极低丰度变体的检测灵敏度，未来可考虑纳入局部单体型信息。

方法本身具有扩展性：可探索其他损失函数（如引入L₁正则化以在候选谱系众多时增加稀疏性）；可根据病毒在废水中的脱落动态知识定制非对称核函数；未来还可纳入空间连续性假设，对多个监测点的信息进行部分池化以提高估计稳健性。

总之，LolliPop通过结合时间序列平滑与反卷积，有效应对了废水测序数据的高噪声和大量缺失挑战，并能通过快速计算的解析置信区间量化不确定性，为大规模废水流行病学项目中的基因组变体追踪提供了强大工具。