全球SARS-COV-2基因组中与菌匹司病毒相关的突变签名

　　我们在突变类别的计算中使用了三个现有的基因组数据来源（以及后来的上下文突变光谱），与已知的使用molnupiravir以及在没有molnupiravir的情况下与典型的SARS-COV-2进化相关。我们分析了Alteri等人17的数据集，该数据集包含两个用Molnupiravir和未经处理的个体治疗的个体的纵向数据。为此，我们从BioProject No下载了FastQ文件。ERP142142使用FASTQ-DL46。我们使用miniMAP2将这些读取映射到HU-1参考基因组，然后在每个位置提取每个底座的呼叫数量。我们确定了与第0天序列相比的突变，计算该位点具有100或更多读数的变体，其中5％或以上是第0天共识的变体。作为次要数据集，我们使用了敏捷试验中的数据（Donovan-Banfield等人18，Bioprojectno。PrJNA854613）。关于Alteri等人17与敏捷试验之间的Molnupiravir突变的性质有普遍的共识，除了仅在敏捷试验中看到的高G-TO-T突变率。先前关于Molnupiravir突变类别的证据，以及在敏捷数据中未经处理的个体中也可以看到高g-T-T率的事实，这使我们得出结论，敏捷数据中的G-TO-T信号代表了技术人工制品。 　　我们使用了先前由Ruis等人计算的BA.1突变频谱作为在相关时间段内典型的SARS-COV-2演化下突变类别和光谱的示例。要比较用molnupiravir和未经治疗的个体治疗的个体之间的突变负担比较突变负担，我们缩放了Ruis等人。19典型BA.1的数据集具有与Alteri等人中未经处理的个体相同数量的总突变。 　　为了确定哪些突变类是使用molnupiravir使用的诊断，我们首先计算了Alteri等人的每个突变类别来自每个突变的比例。17Molnupiravir数据集和Ruis等人。19BA.1 BA.1数据集。然后，我们计算了molnupiravir和（天真的）BA.1数据集之间这些比例的比率。为了在这些比率中建立不确定性，我们从每组突变（带有1,000个bootstrap重复序列）进行了引导重新采样。这些数据如图2b所示。 　　为了识别具有类似洋吡韦的突变模式的全球数据库中的序列，我们分析了由Ultrafast样品放置在现有树上（USHER）Team47上建造的定期更新的突变通知树，该树使用INSDC和Gisaid的几乎所有全局数据，这是McBroome等人的版本。树22。我们使用最初从tackoniumtools48改编的脚本提取数据，然后修改以使用大树资源管理器（BTE）49。该脚本将元数据从测序数据库中添加到每个节点，然后使用简单的启发式方法将这些元数据传递给父节点：（1）如果将所有后代都注释为当年，则将父节点注释一年；（2）如果将其所有后代都注释到该国，则对父母节点进行注释；（3）父母的平均年龄（年龄通知）后代的平均年龄注释。跨越多年或国家的后代的节点很少见。我们还使用Chronumental50计算了节点的更细微的时间估计。我们在调查了感兴趣的各个分支中，使用了Talkonium48，Usher Web Interface47，NextStrain51和NextClade52。 　　我们将“高g-to A分支”定义为至少十个突变的分支机构，其中90％以上是过渡，超过25％的是G-TO-A突变，其中C-TO-A突变超过20％。这样的阈值产生的特异性非常高，这是根据随着时间的推移检测罕见事件率（Molnupiravir治疗）的明显变化的能力来判断的。我们还使用Ruis等人的突变类型的分布创建了敏感性和特异性的模拟度量。19和Alteri等人17我们对不同的分支长度（N）进行了这些计算，从10到20。在每种情况下，我们从每个天真和molnupiraviraviraviravir raviravir相关的突变类别分布中进行了10,000个N突变。然后，我们评估了这些绘制的哪些比例满足上面定义的标准。在摩尔班司迪维尔相关的类别分布的情况下，该比例代表了灵敏度。在典型的BA.1分布的情况下，该比例代表1-特定性。我们获得了分支长度10的灵敏度为46％，特异性为98.9％，分支长度13的灵敏度为63％，特异性为98.6％，分支长度15的灵敏度为71％，特异性为98.6％，敏感性为99.8％，对分支长度的敏感性为64％，特异性为99.8％。 　　为了衡量高g-to A分支是否显示突变率的统计学显着增加，我们在时间上使用了Chronumental的分支长度估计，并进行了统计测试，并在2022年开始对节点上的两侧t检验进行了统计测试，仅查看至少十个突变的节点。 　　为了根据高g-to-a签名的存在来测试年龄元数据是否有所不同，我们从2022年开始将所有我们的节点占据了最小分支长度（10或更多），并根据高g-to-a签名的存在或不存在。我们进行了两侧t检验，以测试所见效果的重要性。为了验证效果并非基本上是由采用后代节点平均值的启发式驱动的，我们重复了仅考虑具有单个后代的分支的分析，并发现了高度相似的结果。 　　为了鉴定基于molnupiravir的诱变的首选核苷酸环境，我们计算了单基替代光谱。对于高G-TO-A分支，我们从Usher突变的树上提取突变路径。使用Wuhan-Hu-1基因组（登录号NC_045512.2）鉴定了每个突变的上下文，并结合了路径早些时候获得的突变。根据基因组含量，通过将突变的数量除以武汉-HU-1基因组中的起始三胞胎的数量来重塑突变计数。muttui（https://github.com/chrisruis/muttui）用于重建和绘制突变光谱。 　　为了计算Alteri等人的单基替代频谱。17数据集，我们使用了BioProject no中的映射读数。PRJNA854613，再次采用具有100个或更多读数的站点，其中5％或更多的位置与第0天共识不同。我们通过将每个突变计数除以Wuhan-Hu-1基因组中的起始三重态数（登录号NC_045512.2），将突变计数重新计算为突变负担。 　　我们对Donovan-Banfield等人进行了类似的分析。数据18。我们使用了从第一天（治疗前）和第五天（治疗后）收集的样品的深层测序数据，该数据接受了65例接受安慰剂治疗的患者和58例接受Molnupiravir治疗的患者。对于每个患者，我们将第一天样本的共有序列用作参考序列，并将突变鉴定为第五天的变体样本，从至少100倍覆盖率的基因组位点中至少5％的读取中的患者参考序列远离患者参考序列。从患者参考序列中鉴定出每个突变的周围核苷酸上下文。 　　为了确保安慰剂和molnupiravir处理之间的任何光谱差异并不是由于先前观察到的SARS-COV-2变体之间的频谱差异19,20，我们比较了处理之间的变体分布（扩展数据图3）。分布非常相似。 　　我们比较了每种过渡突变类型中的上下文模式，评估了从高G到A系统发育分支的16个可能的三核苷酸上下文的相似性与Alteri等人的17数据集的相似性，并分别来自Donovan-Banfield等人。17DataSet和Ruis et al.19 Control.19 Control.19对照组。对于每个数据集组合，计算每个过渡突变类别的余弦相似性。我们在长分支数据中执行了相同的相关分析，将G-TO-A子集与C-T子集进行了比较，将每个G-TO-A上下文与C-TO-T数据集中的反向补充相匹配。 　　鉴于在构造usher突变的树木高度发散序列的过程中，我们决定执行二级分析以识别具有莫纳皮拉维尔签名的发散序列。我们将NextClade52用于此任务。我们提供了完整的全长FastA序列数据集；包括每个可以与NextClade对齐的序列。NextClade将每个序列放在稀疏的参考系统发育树上。它的输出包括一个“未标记的私人突变”列，该列包含相对于树的节点处的私人突变，不包括在其他进化枝中非常常见的恢复突变和突变。我们分析了这组突变，以实现造菌病的突变类别分布。 　　我们选择了具有20个或更多突变的序列，其中20％或更多的是G-TO-A，C-TO-T是C-TO-T或90％或更多的序列是转变。同样，在推出Molnupiravir后，这些日期大量丰富了日期。我们使用usher.bio将确定的序列放在倒下的全局树上，并使用NextStrain51可视化了这棵树。 　　为了测试图4C中所示的序列中的100个或更多突变的上下文与我们在本研究中确定的Molnupiravir Spectrum或典型的BA.1频谱更兼容，我们使用多项式模型进行了分析。我们的目的是忽略突变类别本身的信号（因为这些信号被用来选择序列有趣），而仅考虑通过发生过渡突变的上下文中添加的额外信息。对于每个过渡类别（G-TO-A，C-T，A-TO-G，T-TO-C），我们创建了两个多项式模型的三核苷酸上下文，一种使用我们在这项研究中定义的长分支molnupiravir Spectrum，一个使用Ruis等人的BA.1 Spectrum。在每种情况下，我们都乘以基因组中三核苷酸上下文的次数，以消除对此参数的先前归一化。我们评估了在两个模型下观察感兴趣顺序的上下文计数的可能性，并计算了每个贝叶斯因子（g-to-a：35,017，c-to-t：6,068：6,068：6,068，a-to-g-g：53，t-to-to-c-c：1.22）。这些结合得出的贝叶斯系数为1.4×1010。 　　我们检查了组成这些分支的突变类型。我们使用BTE来确定观察到的每个突变是否是同义词。根据该分支是否短（少于十个突变）或长（十个或更多突变），以及它是否具有高g-to-a签名，对突变进行了分组。我们计算了每种情况下非同义词的突变的比例，计算了这些比例的二项式置信区间，并使用使用R. 　　我们绘制了高g-TO-A分支的突变分布，根据突变是否是同义词，同时还针对四个或更多高的G-TO-A分支中最常见的非同步突变绘制了分布。用高斯核和500 bp的带宽进行核密度估计。 　　使用GGTREE53从Usher树中绘制了补充剂中的散装树木。 　　对于图4A中所示的20个个体的群集，我们观察到了乌瑟尔（Usher）在某些位置缺少覆盖范围而导致的群集中突变的树中的小瑕疵。因此，我们重新计算了我们在这里显示的树。我们在群集中取了20个序列，使用NextClade52对齐三个最接近的外群序列，使用IQ-Tree54计算了一棵树，并使用Treetime55重建了突变的树。我们使用Figtree56可视化了树。 　　英国卫生安全局（UKHSA）分析了来自英格兰的高g-to-a签名的49个序列，分为35个集群。基于NHS编号，将序列与Blueteq处理记录34链接。为所有序列建立了链接。分析发现，在35个不同的簇中，有11个涉及单个处方的莫纳皮拉维尔（Molnupiravir），簇命中率为31％。仅计算治疗日期后采样的序列，没有上限。 　　我们的工作有一些局限性。将特定的分支识别为拥有类似洋匹司病毒的特征是一种概率，而不是绝对判断：莫纳皮拉维尔只会产生少数突变（通常是试验数据表明的情况），分支长度太小，无法置信地分配突变的原因。因此，我们将分析局限于长分支。这种方法也可能无法检测到具有大量摩尔班司米拉维尔诱导的突变的分支，以及来自其他原因的大量突变（可能发生在慢性感染中）。我们鉴定与摩尔班司蒂拉维尔相关序列的方法使用了简单的阈值，该阈值是在具有不同类别突变类别的分支上的突变比例。这种方法的简单性不会使检测概率成为分支长度的函数，使我们能够执行分析，例如在不同条件下查看分支长度的分布；但是，未来通过更细微的统计方法提高灵敏度的未来分析（我们对此进行了实验，发现在第一种情况下以其灵活性为优点的简单方法），并考虑上下文突变谱本身作为检测的信号，在未来的工作中都将是有价值的。 　　我们根据国家发现了与摩尔努维韦相关的序列的截然不同的速率，这部分反映了在不同地理区域中是否使用molnupiravir以及如何使用。但是，在使用molnupiravir的设置中测序基因组的速率也将有贡献。例如，如果主要用于这些设施的老年护理设施中的Molnupiravir比一般社区中的病毒更有可能被测序，则这将提高此类序列的确定率。此外，很可能将某些序列专门分析为特定研究的一部分，因为样品表现出在莫诺皮迪拉维尔治疗后持续的测试阳性。根据不同位置的测序优先级，此类效果可能会有所不同。我们确定了从高g到A节点降低的序列簇。在许多情况下，详细且独特的元数据表明，一个特定的簇是由不同患者的序列组成的，表明莫纳皮迪拉维尔诱导的突变传播。但是，在没有此类数据的情况下，群集也与表示从单个个体获取的多个样本兼容。 　　我们的分析研究了共识序列，这意味着要检测到突变，必须达到宿主中人群的很高比例。在共识序列中查看深度测序数据和混合基础呼叫的分析将是有价值的。 　　有关研究设计的更多信息可在与本文有关的自然投资组合报告摘要中获得。