用于使用神经网络进行种系和体细胞变体调用的系统和方法

摘要

本公开提供利用神经网络诸如卷积神经网络来分析由测序仪生成的基因组序列数据并且生成鉴别和描述所述序列数据内的种系和/或体细胞变体的准确预测数据的系统和方法。

说明书

相关申请的交叉引用

本申请要求 2018 年 8 月 13 日提交的美国临时申请号 62/718,338 和 2019年 5 月 31 日提交的美国临时申请号 62/855,541 的优先权，这些临时申请各自借此以全文引用的方式并入。

通过引用并入

本说明书中提到的所有出版物和专利申请都以相同的程度通过引用并入本文，就好像每个单独的出版物或专利申请都被具体且单独地指出通过引用并入一样。

技术领域

本发明的实施例总体上涉及用于变体调用的系统和方法，并且更具体地涉及使用机器学习算法进行变体调用。

背景技术

使用测序数据鉴别遗传突变是生物信息学中的一项至关重要的任务，它对包括癌症在内的多种疾病的诊断、预后和治疗具有影响。这是一项不平凡的任务，尤其是当用于生成测序数据的测序技术具有较高的错误率（这是单分子测序技术所常见的）时或当突变以较低的频率（癌症突变尤其如此）或在复杂的基因组区域内发生时。测序样本中癌症突变的频率较低可能是由于样本污染（因为肿瘤样本中可能含有一些来自正常细胞的 DNA）或肿瘤异质性所致。

传统方法通过限制特定技术和待测突变类型来缩窄问题的范围。因此，存在用于鉴别种系和癌症突变的不同成套工具和框架，其根据突变是大还是小进一步分层。另外，用于低错误测序技术（例如，短读取测序）的方法不同于用于高错误测序技术（例如，单分子测序）的方法。差异主要体现在数据标志 (data signature) 以及用于将错误与真实突变区分开的统计模型中。显式数据标志和统计模型的使用限制了传统方法对多种技术的通用性。此外，即使是传统方法，也需要调整参数来捕获样本的特定特征，这可能取决于多种因素，包括所使用的技术、测序方案、测序样本、样本纯度等。

深度学习是一种机器学习技术，只要有正确的训练数据，就可以使用多个数据层隐式地捕获数据标志和统计模型。就突变检测问题的通用性而言，这使得深度学习成为潜在的有吸引力的解决方案。已知几种用于对种系突变进行分类的深度学习解决方案。

Poplin, R. 等人 (“Creating a Universal SNP and Small Indel VariantCaller With Deep Neural Networks”, bioRxiv, 092890 (2016)) 提出了一种方法，其中首先通过基因组扫描鉴别候选变体定位，然后将在每个所鉴别的候选周围比对的读取的堆叠图像作为输入提供到卷积神经网络 (CNN) 中，然后该卷积神经网络预测候选变体的基因型。这种方法不是很实用，因为它需要复杂的神经网络架构和大量内存。

Luo, R. 等人 (“Clairvoyante: a Multi-task Convolutional Deep NeuralNetwork For Variant Calling In Single Molecule Sequencing. bioRxiv, 310458(2018)) 使用数据汇总方法，该方法包括扫描基因组以鉴别候选变体位点，并为每个候选位置准备多个输入矩阵，每种变体类型一个矩阵。Chin, J.（“Simple ConvolutionalNeural Network for Genomic Variant Calling with TensorFlow”，https://towardsdatascience.com/simple- convolution-neural-network-for-genomic-variant-calling-with-tensorflow-c085dbc2026f，2017 年 7 月 16 日出版）中讨论了类似的技术。这些方法也不是最优的，因为它们不能清楚地表示整个窗口中的所有插入序列。例如，Clairvoyante 系统为插入事件分配了一个独立的输入矩阵，该矩阵仅捕获变体定位处的插入，而忽略窗口内其他定位处的真实/错误插入。有关整个窗口中插入序列精确定位的全图信息对于将真实插入变体与背景噪声区分开来至关重要。

Torracinta, R. 等人 (“Adaptive somatic mutations calls with deeplearning and semi- simulated data. bioRxiv, 79087 (2016)) 描述了将六层全连接神经网络应用于一组手动提取的特征。但是仅具有全连接层的系统不能利用卷积神经网络提供的功能，即，使用见于候选变体的局部基因组情境中的模式直接从原始序列数据中学习特征表示。此外，由于全连接的网络更加复杂，因此与使用 CNN 所能实现的相比，该方法的通用性和可扩展性较差。

发明内容

本发明的总体上涉及用于变体调用的系统和方法，并且更具体地涉及使用机器学习算法进行变体调用。

在一些实施例中，提供了用于种系变体调用的方法。该方法可以包括获得参考序列、多条序列读取和候选变体在序列读取内的位置；通过将一个或多个空格插入一条或多条序列读取中而获得经扩充序列读取；通过将一个或多个空格插入参考序列中而获得经扩充参考序列；将候选变体周围的经扩充序列读取的片段转换为样本矩阵；将候选变体周围的经扩充参考序列的片段转换为参考矩阵；将样本矩阵和参考矩阵提供给受训神经网络；以及，在受训神经网络的输出端，获得关于多条序列读取内的变体的预测数据。

在一些实施例中，该方法进一步包括检测多条序列读取内的一个或多个插入碱基，其中扩充序列读取和参考序列包括：对于在任何序列读取中的任一条中检测到的每个插入碱基，在该插入碱基的位置处将空格插入参考样本中。

在一些实施例中，该方法进一步包括：对于在所述序列读取中的任一条中检测到的每个插入碱基，在该插入碱基的位置处将空格插入每一条在该插入碱基的位置处没有检测到插入的序列读取中。

在一些实施例中，样本矩阵包括至少四个表示四种类型的核苷酸碱基的行，每一行表示在经扩充序列读取的片段内不同位置处的相应核苷酸碱基类型的碱基数；和至少一个表示在经扩充序列读取的片段内不同位置处的插入的空格数的行。

在一些实施例中，参考矩阵具有与样本矩阵相同的维数，并且其中参考矩阵提供经扩充参考序列内不同核苷酸碱基和空格的具体定位的完全表示。

在一些实施例中，受训神经网络包括受训卷积神经网络。

在一些实施例中，该方法进一步包括将以下至少一项提供给受训神经网络：变体位置矩阵，其表示经扩充序列读取的片段内候选变体的位置；覆盖率矩阵，其表示经扩充序列读取的片段的覆盖率或深度；比对特征矩阵，其表示经扩充序列读取的比对特征；知识性碱基矩阵，其表示关于一个或多个变体的公开已知信息的信息。

在一些实施例中，关于变体的预测数据包括以下至少一项：变体的预测类型；变体的预测位置；变体的预测长度；以及变体的预测基因型。

在一些实施例中，关于变体的预测数据包括变体的预测类型，并且其中神经网络配置为产生用于变体的预测类型的多个值中的一个，该多个值包括：第一值，其指示变体为假阳性的概率；第二值，其指示变体为单核苷酸多态性变体的概率；第三值，其指示变体为缺失变体的概率；以及第四值，其指示变体为插入变体的概率。

在一些实施例中，提供了用于体细胞变体调用的方法。该方法可以包括获得多条正常序列读取和多条肿瘤序列读取；将正常序列读取的片段和肿瘤序列读取的片段分别转换为正常样本矩阵和肿瘤样本矩阵；将正常样本矩阵和肿瘤样本矩阵馈入受训卷积神经网络中；以及在受训卷积神经网络的输出端，获得多条肿瘤序列读取内的体细胞变体的预测类型。

在一些实施例中，多条肿瘤序列读取表示患者的肿瘤样本的遗传信息，并且多条正常序列读取表示患者的正常样本的遗传信息。

在一些实施例中，将正常序列读取转换为正常样本矩阵包括通过将一个或多个空格插入一条或多条正常序列读取中而扩充正常序列读取的片段；并且将肿瘤序列读取的片段转换为肿瘤样本矩阵包括通过将一个或多个空格插入一条或多条肿瘤序列读取中而扩充肿瘤序列读取的片段。

在一些实施例中，肿瘤样本矩阵包括至少一个用于每种核苷酸碱基类型的行，每一行表示在肿瘤序列读取的片段内每个位置处的相应核苷酸碱基类型的出现次数；和至少一个表示在肿瘤序列读取的片段内每个位置处的插入空格数的行。

在一些实施例中，该方法进一步包括向受训卷积神经网络提供一个或多个矩阵，该一个或多个矩阵表示从一个或多个其他变体调用者获得的特征，该调用者已经分析了多个肿瘤序列读取和/或多个正常序列读取。

在一些实施例中，该方法进一步包括获得参考序列；将参考序列转换为参考矩阵；以及将参考矩阵与正常样本矩阵和肿瘤样本矩阵一起馈入受训卷积矩阵中。

在一些实施例中，提供一种非暂时性计算机可读介质，其包含指令，该指令当由计算机系统的一个或多个处理器执行时，使计算机系统实施包括以下各项的操作：获得多条正常序列读取和多条肿瘤序列读取；将正常序列读取的片段和肿瘤序列读取的片段分别转换为正常样本矩阵和肿瘤样本矩阵；将正常样本矩阵和肿瘤样本矩阵馈入受训卷积神经网络中；以及在受训卷积神经网络的输出端，获得多条正常序列读取内的体细胞变体的预测类型。

在一些实施例中，提供一种计算机系统，其包含一个或多个处理器并且耦接至一个或多个存储有指令的非暂时性计算机可读存储器，该指令当由所述计算机系统执行时，使该计算机系统实施包括以下各项的操作：获得多条肿瘤序列读取；通过将一个或多个空格插入一条或多条肿瘤序列读取中获得经扩充肿瘤序列读取；将肿瘤序列读取的片段转换为肿瘤样本矩阵；将正常样本矩阵和肿瘤样本矩阵馈入受训神经网络中；以及在受训神经网络的输出端，获得多条肿瘤序列读取内的体细胞变体的预测类型。

在一些实施例中，提供了用于变体调用的方法。该方法可以包括：获得参考序列和多条序列读取；任选地对多条序列读取与参考序列实施第一比对，除非所获得的多条序列读取和参考序列是以已经比对的配置获得的；从经比对序列读取和参考序列中鉴别候选变体位置；在候选变体位置周围扩充序列读取和/或参考序列，以实现多条序列读取与参考序列的第二比对；从经扩充参考序列生成针对候选变体位置的参考矩阵，并从多条经扩充序列读取生成针对候选变体位置的样本矩阵；将参考矩阵和样本矩阵输入神经网络中；以及用该神经网络确定候选变体位置处是否存在变体类型。

在一些实施例中，扩充序列读取和/或参考序列的步骤包括将一个或多个空格引入序列读取和/或参考序列，以解释序列读取中的插入和/或缺失。

在一些实施例中，该方法进一步包括从训练数据集生成多个训练矩阵，其中训练矩阵具有与样本矩阵和参考矩阵相对应的结构，其中训练数据集包含序列数据，该序列数据包含多个突变，该多个突变包含单核苷酸变体、插入和缺失；以及使用多个训练矩阵训练神经网络。

在一些实施例中，训练数据集包含多个子集，其中每个子集包含从 0% 到 100%范围内的肿瘤纯度水平，其中该子集中的至少两个子集各自具有不同的肿瘤纯度水平。

在一些实施例中，该子集中的至少三个子集各自具有不同的肿瘤纯度水平。

在一些实施例中，多个子集包括具有小于约 30% 的肿瘤纯度水平的第一子集、具有约 30% 和 70% 之间的肿瘤纯度水平的第二子集和具有至少约 70% 的第三肿瘤纯度水平的第三子集。

在一些实施例中，多个子集包括具有小于约 40% 的肿瘤纯度水平的第一子集、具有约 40% 和 60% 之间的肿瘤纯度水平的第二子集和具有至少约 60% 的肿瘤纯度水平的第三子集。

在一些实施例中，多个子集包括具有小于约 10% 的肿瘤纯度水平的子集。

在一些实施例中，多个子集包括具有小于约 5% 的肿瘤纯度水平的子集。

在一些实施例中，训练数据集包括合成数据。

在一些实施例中，合成数据包括人工生成的突变，其中人工生成的突变包含单核苷酸变体、插入和缺失。

在一些实施例中，训练数据集包括真实数据，其中真实数据包括真实突变，其中真实突变包括单核苷酸变体、插入和缺失。

在一些实施例中，训练数据集包含多个子集，其中每个子集包含从 0% 到 100%范围内的变体等位基因频率，其中该子集中的至少两个子集各自具有不同的变体等位基因频率。

在一些实施例中，该子集中的至少三个子集各自具有不同的变体等位基因频率水平。

在一些实施例中，该子集中的至少一个子集具有至少 2.5% 的变体等位基因频率。

在一些实施例中，该子集中的至少一个子集具有至少 5% 的变体等位基因频率。

在一些实施例中，该子集中的至少一个子集具有至少 10% 的变体等位基因频率。

在一些实施例中，该方法进一步包括将来自至少一个突变调用算法的至少一个预测输入神经网络中。

在一些实施例中，至少一个预测包括来自至少三个独立的突变调用算法的至少三个预测。

在一些实施例中，至少一个预测包括来自至少五个独立的突变调用算法的至少五个预测。

在一些实施例中，训练数据集包括合成数据和真实数据的集合体。

在一些实施例中，训练数据集包括至少 5% 的合成数据。

在一些实施例中，训练数据集包括至少 10% 的合成数据。

在一些实施例中，训练数据集包括全基因组测序数据。

在一些实施例中，训练数据集包括全外显子组测序数据。

在一些实施例中，训练数据集包括靶向测序数据。

在一些实施例中，训练数据集包括从福尔马林固定石蜡包埋样本获得的数据。

在一些实施例中，训练数据集包括全基因组测序数据、全外显子组测序数据、靶向测序数据和从福尔马林固定石蜡包埋样本获得的数据中的至少两者。

在一些实施例中，训练数据集包括全基因组测序数据、全外显子组测序数据、靶向测序数据和从福尔马林固定石蜡包埋样本获得的数据中的至少三者。

在一些实施例中，训练数据集包括全基因组测序数据、全外显子组测序数据、靶向测序数据和从福尔马林固定石蜡包埋样本获得的数据。

在一些实施例中，提供了用于变体调用的方法。该方法可以包括获得参考序列、多条肿瘤序列读取和多条正常序列读取；任选地实施多条肿瘤序列读取和多条正常序列读取与参考序列的第一比对，除非所获得的多条肿瘤序列读取和多条正常序列读取及参考序列是以已经比对的配置获得的；从经比对肿瘤序列读取、正常序列读取和参考序列中鉴别候选变体位置；在候选变体位置周围扩充肿瘤序列读取和/或正常序列读取和/或参考序列以实现多条肿瘤序列读取和多条正常序列读取与参考序列的第二比对；从经扩充参考序列生成针对候选变体位置的参考矩阵，从多条经扩充肿瘤序列读取生成针对候选变体位置的肿瘤矩阵，并且从多条经扩充正常序列读取生成针对候选变体位置的正常矩阵；将参考矩阵、肿瘤矩阵和正常矩阵输入神经网络中；以及用该神经网络确定候选变体位置处是否存在变体类型。

在一些实施例中，该方法进一步包括从训练数据集生成多个训练矩阵，其中训练矩阵具有与肿瘤矩阵、正常矩阵和参考矩阵相对应的结构，其中训练数据集包括肿瘤序列数据和正常序列数据；以及使用多个训练矩阵训练神经网络。

在一些实施例中，肿瘤序列数据和正常序列数据两者均包括多个突变，该突变包括单核苷酸变体、插入和缺失。

在一些实施例中，正常序列数据包括最多 5% 的肿瘤序列数据。

在一些实施例中，正常序列数据包括最多 10% 的肿瘤序列数据。

在一些实施例中，肿瘤序列数据包括约 10% 至 100% 之间的肿瘤纯度水平。

在一些实施例中，训练数据集包括多个肿瘤序列数据子集，其中每个肿瘤序列数据子集包含从 10% 到 100% 范围内的肿瘤纯度水平，其中肿瘤序列数据子集中的至少两个子集各自具有不同的肿瘤纯度水平。

在一些实施例中，肿瘤序列数据子集中的至少三个子集各自具有不同的肿瘤纯度水平。

在一些实施例中，多个肿瘤序列数据子集包括具有小于约 30% 的肿瘤纯度水平的第一肿瘤序列数据子集、具有约 30% 和 70% 之间的肿瘤纯度水平的第二肿瘤序列数据子集和具有至少约 70% 的肿瘤纯度水平的第三肿瘤序列数据子集。

在一些实施例中，多个肿瘤序列数据子集包括具有小于约 40% 的肿瘤纯度水平的第一肿瘤序列数据子集、具有约 40% 和 60% 之间的肿瘤纯度水平的第二肿瘤序列数据子集和具有至少约 60% 的肿瘤纯度水平的第三肿瘤序列数据子集。

在一些实施例中，训练数据集包括合成数据。

在一些实施例中，合成数据包括人工生成的突变，其中人工生成的突变包含单核苷酸变体、插入和缺失。

在一些实施例中，训练数据集包括真实数据，其中真实数据包括真实突变，其中真实突变包括单核苷酸变体、插入和缺失。

在一些实施例中，训练数据集包括全基因组测序数据。

在一些实施例中，训练数据集包括全外显子组测序数据。

在一些实施例中，训练数据集包括靶向测序数据。

在一些实施例中，训练数据集包括从福尔马林固定石蜡包埋样本获得的数据。

在一些实施例中，提供了一种系统。该系统可以包括处理器，其配置成执行实施权利要求 18 至 64 中任一项中所述的步骤。

附图简要说明

本发明的新颖特征在所附的权利要求书中具体阐述。通过参考以下具体实施方式并结合附图，可以更好地理解本发明的特征和优点，具体实施方式阐述了其中利用本发明原理的说明性实施例，在附图中：

图 1 示出了根据一些实施例的包括测序装置和计算机系统的示例性系统；

图 2a 和图 2b 示出了根据一些实施例的用于种系变体调用的示例性方法；

图 3 示出了根据一些实施例的图 2 的示例性方法的扩充步骤；

图 4 示出了根据一些实施例的可以提供给神经网络和从神经网络获得的示例性数据；

图 5 示出了根据一些实施例的神经网络的示例性架构；

图 6 示出了根据一些实施例的用于体细胞变体调用的示例性方法；

图 7 示出了根据一些实施例的可作为图 6 的示例性方法的一部分实施的扩充步骤；

图 8 示出了根据一些实施例的可以提供给神经网络和从神经网络获得的数据的另一示例；

图 9 示出了根据一些实施例的神经网络的示例性架构；

图 10a 示出了针对给定的候选体细胞 SNV 的输入矩阵的生成；图 10b 示出了输入矩阵和网络架构；图 10c 示出了针对给定的候选体细胞缺失的输入矩阵的生成；图10d 示出了针对给定的候选体细胞插入的输入矩阵的生成；

图 11a 至 图 11c 示出了不同检测方法在铂双样本混合物数据集上的性能；

图 12a 至 图 12c 示出了不同检测方法在 Dream 阶段 3 数据集上的性能；

图 13a 至 图 13c 示出了不同检测方法在 Dream 阶段 4 数据集上的性能；

图 14a 至 图 14c 示出了两种不同检测方法在 PacBio 数据集上的性能；

图 15a 至图 15b 示出了序列覆盖率对于不同检测方法在全外显子组样本混合物数据集上的性能的影响；

图 16a 至图 16d 示出了不同检测方法在铂肿瘤掺入数据集上的性能；

图 17a 和图 17c 示出了不同的体细胞突变检测方法在全外显子组数据集上的性能；图 17b 和 图 17d 示出了体细胞突变检测方法在靶向组数据集上的性能；

图 18a 示出了使用基于全基因组和全外显子组训练的模型对外显子组混合物数据集进行测试的性能分析；图 18b 示出了使用基于全基因组和全外显子组训练的模型对靶向组混合物数据集进行测试的性能分析；

图 19 示出了在 Dream 阶段 3、Dream 阶段 4、铂双样本混合物、铂肿瘤掺入、PacBio 和外显子组数据集中的基准真值 (ground truth) INDEL 的尺寸分布；

图 20a 和图 20b。图 20a 示出了针对各种数据集，各种检测方法基于预测的体细胞突变位置和类型对 INDEL 进行测试的性能分析；

图 21 示出了针对 Dream 挑战阶段 3 数据集的交叉样本训练的性能分析；

图 22a 示出了不同体细胞突变检测方法在真实数据集 COLO-829 上的性能；图22b 示出了不同体细胞突变检测方法在真实数据集 CLL1 上的性能；图 22c 示出了不同体细胞突变检测方法在真实数据集 TCGA-AZ-6601 上的性能；

图 23 是用于 Microsoft Azure 实验的 261 个 TCGA 癌症样本的列表；

图 24a 至图 24m 示出了所测试的不同网络架构；

图 25 示出了图 24a 至图 24m 所示的不同网络架构的性能分析；

图 26a 和图 26b 示出了不同体细胞突变检测算法的运行时比较；

图 27 示出了在由 SEQC-II 联合体分类为四个置信度水平（高、中、低和未分类）的针对 HCC1395 的调用结果超集合中，基准真值 SNV 和 INDEL 体细胞突变的 VAF 分布；

图 28a 至图 28e 显示了 NeuSomatic 在 SEQC-II 数据集中的 123 个复本上的总体性能；

图 29a 至图 29d 显示了各种模型在 21 个 WGS 复本数据集上的性能比较。图29a 示出了 F1 得分 (%)，图 29b 示出了跨不同体细胞突变调用者的精确率-召回率比较。在这里，SEQC-WGS-GT50-掺入WGS10 受训模型用于 NeuSomatic 和 NeuSomatic-S。图29c 显示了各种 NeuSomatic 受训模型的 F1 得分 (%) 比较。图 29d 显示了不同体细胞突变检测方法和基于 WGS 数据集的 NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 30a 至图 30d 显示了肿瘤纯度数据集的性能比较。图 30a 显示了跨不同覆盖率（10× 至 300×）和肿瘤纯度（5% 至 100%）的不同体细胞突变调用者的 SNV 和INDEL F1 得分 (%) 比较。图 30b 显示了匹配的正常样本对肿瘤污染的鲁棒性：在 80×覆盖率和 5% 至 100% 肿瘤纯度下，当匹配的正常样本被 5% 的肿瘤污染时，F1 得分 (%)与纯正常样本的 F1 得分 (%) 相比有所改变。对于图 30a 和图 30b，SEQC-WGS-GT50-掺入WGS10 受训模型用于 NeuSomatic 和 NeuSomatic-S。图 30c 显示了各种 NeuSomatic受训模型的 F1 得分 (%) 比较。图 30d 显示了针对不同 INDEL 尺寸的性能分析。负INDEL 尺寸反映缺失；

图 30e 和图 30f 显示了基于肿瘤纯度数据集的精确率-召回率分析。图 30e 显示了 SNV，并且图 30f 显示了跨不同覆盖率（10× 至 300×）和肿瘤纯度（5% 至 100%）的不同体细胞突变调用者的 INDEL 准确性比较；

图 30g 显示了基于肿瘤正常滴定数据集的各种 NeuSomatic-S 受训模型的精确率-召回率比较；

图 30h 显示了不同体细胞突变检测方法和基于肿瘤-正常滴定数据集的NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 31a 至图 31f 显示了文库制备、WES 和 AmpliSeq 数据集的性能比较。图31a 显示了跨不同文库-试剂盒和 DNA 量的不同体细胞突变调用者的 SNV 和 INDEL F1得分 (%) 比较。图 31b 显示了针对不同文库-试剂盒和 DNA 量的各种 NeuSomatic 受训模型的精确率-召回率比较。图 31c 和图 31e 显示了针对不同体细胞突变调用者的 SNVF1 得分 (%) 比较。图 31d 和图 31f 显示了各种 NeuSomatic 受训模型的 F1 得分 (%)比较。对于图 31a、图 31c 和图 31e，SEQC-WGS-GT50-掺入WGS10 受训模型用于NeuSomatic 和 NeuSomatic-S；

图 31g 显示了基于文库制备数据集的各种 NeuSomatic-S 受训模型的精确率-召回率比较；

图 31h 显示了不同体细胞突变检测方法和基于文库制备数据集的 NeuSomatic网络受训模型的 F1 得分 (%) 性能；

图 31i 显示了不同体细胞突变检测方法和基于 WES 数据集的 NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 31j 显示了不同体细胞突变检测方法和基于 AmpliSeq 数据集的NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 32a 至图 32d 显示了 FFPE 数据集的性能比较。图 32a 和图 32b 显示了基于 16 个 FFPE WGS 复本的性能。图 32c 和图 32d 显示了基于 14 个 FFPE WGS 复本的性能。图 32a 和图 32c 显示了跨 FFPE 和新鲜正常样本的不同体细胞突变调用者的 F1得分 (%) 比较。图 32b 和图 32d 显示了各种 NeuSomatic 受训模型的 F1 得分 (%) 比较。图 32b 显示了使用 FFPE 和新鲜正常样本进行的不同体细胞突变调用者的 16 个FFPE WGS 复本的精确率-召回率图。对于图 32a 和图 32c，SEQC-WGS-GT50-掺入WGS10 受训模型用于 NeuSomatic 和 NeuSomatic-S；

图 32e 显示了不同体细胞突变检测方法和基于 FFPE WGS 数据集的NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 32f 显示了不同体细胞突变检测方法和基于 FFPE WES 数据集的NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 33a 和图 33b 是分析使用 NeuSomatic（图 33a）和 NeuSomatic-S（图 33b）针对来自不同数据集的九个复本进行训练对基于靶标样本训练的影响的图表。对于每个样本，使用两个样本特异性模型：一个仅基于靶标样本训练，另一个则使用来自 SEQC-WGS 掺入模型的另外 10% 的训练数据进行训练；

图 33c 显示了 SEQC-WGS-掺入模型和两个样本特异性 NeuSomatic 网络受训模型的 F1 得分 (%) 性能；

图 34a 至图 34d 显示了不同 INDEL 尺寸、不同 VAF 分布以及在困难区域上的性能分析。图 34a 和图 34b 显示了使用不同调用者和不同训练方法对跨 SEQC-II 数据集的不同 INDEL 尺寸的性能分析。负 INDEL 尺寸反映缺失。图 34c 显示了使用不同调用者对跨 SEQC 数据集的具有不同 VAF 范围的突变的性能分析。图 34d 显示了不同体细胞调用者在困难区域上的性能分析，该困难区域包括不同尺寸的串联重复序列 (TR) 和片段复制；

图 35 示出了使用不同调用者对跨 SEQC 数据集的具有不同 VAF 范围的突变的性能分析。

图 36a 和图 36b 显示了使用不同的针对 NeuSomatic（图 36a）和 NeuSomatic-S（图 36b）的训练方法对跨 SEQC 数据集的具有不同 VAF 范围的突变的性能分析；

图 37a 至图 37c 显示了使用不同的调用者（图 37a）和针对 NeuSomatic（图37b）和 NeuSomatic-S（图 37c）的训练方法（图 37b 和图 37c）对肿瘤/正常滴定数据集的具有不同 VAF 范围的突变的性能分析；

图 38 显示了使用不同调用者在全基因组与困难基因组区域上的性能比较，该困难基因组区域包括不同尺寸的串联重复序列 (TR) 和片段复制；

图 39a 和图 39b 显示了不同的针对 NeuSomatic（图 39a）和 NeuSomatic-S（图39b）的训练方法在困难区域上的性能分析，该困难区域包括不同尺寸的串联重复序列(TR) 和片段复制；

图 40a 和图 40b 显示了使用不同的针对 NeuSomatic（图 40b）和 NeuSomatic-S（图 40c）的训练方法在全基因组和困难基因组区域上的性能比较，该困难区域包括不同尺寸的串联重复序列 (TR) 和片段复制；

图 41a 至图 41c 显示了使用不同的调用者（图 41a）和针对 NeuSomatic（图41b）和 NeuSomatic-S（图 41c）的训练方法（图 41b 和图 41c）对肿瘤/正常滴定数据集上困难区域的性能分析，该困难区域包括不同尺寸的串联重复序列 (TR) 和片段复制；以及

图 42 显示了基于 WGS 数据集的私有 FN 调用结果的 VAF 分布的小提琴图比较。

具体实施方式

此外，本公开（除了描述其它以外）还描述一种用于种系变体调用的方法，该方法可以包括获得参考序列、多条序列读取和候选变体在序列读取内的位置；通过将一个或多个空格插入一条或多条序列读取中而获得经扩充序列读取；通过将一个或多个空格插入参考序列中而获得经扩充参考序列；将候选变体周围的经扩充序列读取的片段转换成样本矩阵；将候选变体周围的经扩充参考序列的片段转换成参考矩阵；将样本矩阵和参考矩阵提供给受训神经网络；以及，在受训神经网络的输出端，获得关于多条序列读取内的变体的预测数据。所公开的系统和方法允许直接从原始数据捕获重要的变体信号，并且对于不同的测序技术、样本纯度和测序策略（诸如全基因组与靶标富集）一致地实现高精确率。

图 1 阐述了说明性系统 100，其包括通信耦接至计算机系统 102 的测序装置110。测序装置 110 可以直接（例如，通过一条或多条通信电缆）或通过网络 130 耦接至计算机系统 102，网络 130 可以是互联网或广域网络、局域网络、有线网络和/或无线网络的任何其他组合。在一些实施例中，计算机系统 102 可以包括在或整合在测序装置 110 中。在一些实施例中，测序装置 110 可以对包含遗传物质的样本进行测序，并产生所得的测序数据。测序数据可以被发送到计算机系统 102（例如，通过网络 130），或存储在存储装置上并在以后的阶段被发送到计算机系统 102（例如，通过网络 130）。在一些实施例中，计算机系统 102 可以包括或可以不包括显示器 108 和一个或多个输入装置（未示出），用于从使用者或操作者（例如，技术人员或遗传学家）接收命令。在一些实施例中，计算机系统 102和/或测序装置 110 可以由使用者或其他装置通过网络 130 远程访问。因此，在一些实施例中，本文讨论的各种方法可以在计算机系统 102 上远程运行。

计算机系统 102 可以包括一个计算机装置或任何类型的多个计算机装置的组合，诸如个人计算机、膝上型计算机、网络服务器（例如，本地服务器或公共/私有/混合云中包括的服务器）、移动装置等等，其中一些或所有装置可以互连。计算机系统 102 可以包括一个或多个处理器（未示出），每个处理器可以具有一个或多个核心。在一些实施例中，计算机系统 102 可以包括一个或多个通用处理器（例如，CPU）、专用处理器诸如图形处理器(GPU)、数字信号处理器、或这些和其他类型的处理器的任意组合。在一些实施例中，可以使用定制或可定制的电路诸如专用集成电路 (ASIC) 或现场可编程门阵列 (FPGA) 来实现计算机系统中的一些或所有处理器。在一些实施例中，计算机系统 102 还可以检索并执行存储在整合到或以其他方式通信耦接至计算机系统 102 的一个或多个存储器或存储装置（未示出）中的非暂时性计算机可读指令。存储器/存储装置可以包括非暂时性计算机可读存储介质的任意组合，包括各种类型的半导体存储芯片（DRAM、SRAM、SDRAM、闪存、可编程只读存储器）等等。也可以使用磁盘和/或光盘。存储器/存储装置也可以包括可读取和/或可写入的可移动存储介质；此类介质的示例包括光盘 (CD)、只读数字通用光盘（例如，DVD-ROM、双层 DVD-ROM）、只读和可读 Blu-ray® 磁盘、超密度光盘、闪存存储卡（例如，SD 卡、小型 SD 卡、微型 SD 卡等）等等。在一些实施例中，数据和其他信息（例如测序数据）可以存储在一个或多个远程定位（例如，云存储）中，并与系统 100 的其他部件同步。

图 2a 示出了根据一些实施例的用于种系变体调用的示例性方法 200。方法 200可以例如以软件的形式（即，存储在一个或多个非暂时性介质中的可由例如计算机系统102 的一个或多个处理器访问和执行的一组指令）或以固件、硬件或其任意组合的形式来实现。

方法 200 可以在步骤 210 处开始，在步骤 210 处可以获得参考序列。参考序列可以例如从一个或多个私有或公共资源库获得，诸如参考序列（RefSeq，由美国国家生物技术信息中心 (NCBI) 建立的核苷酸碱基序列的开放获取、注释和管理集合）或 NCBIGenomes FTP 网站（其存储了一组不同生物的完整基因组）。应当理解，在一些实施例中，参考序列的特定副本可以本地存储（例如，存储在计算机系统 102 的存储器中），而在其他实施例中，参考序列可以从远程服务器获得，例如，通过网络 130 获得。此外，在一些实施例中，可以获得完整的参考序列，而在其他实施例中，可以获得参考序列的一个或多个区段，例如，仅与具体测定法相关联的一个或多个区段。因此，本文所用的“参考序列”通常是指参考序列的一个或多个区段，其可以包括或可以不包括完整的参考序列。

在步骤 220，可以获得与由测序装置 110 测序的遗传样本（例如，包含患者的DNA 或 RNA 物质的样本）相对应的多条序列读取。如上所述，可以直接从测序装置 110 或从通信耦接至计算机系统 102 的一个或多个本地或远程易失性或非易失性存储器、存储装置或数据库中获得序列读取。所获得的序列读取可以已经预处理的（例如，预先比对的），或者它们可以是“原始的”，在这种情况下，方法 200 还可以包括预处理（例如，预先比对）步骤（未示出）。而且，尽管在一些实施例中，可以获得完整的序列读取（由测序装置 110 生成），但在其他实施例中，可以仅获得序列读取的区段。因此，如本文所用，“获得序列读取”通常是指获得一条或多条（例如，相邻）序列读取的一个或多个区段。

在步骤 230，可以获得序列读取内的多个候选变体位置。在一些实施例中，这包括扫描获得的多条序列读取，将它们与获得的参考序列进行比较，并确定该序列读取内的似乎包括一些类型的变体（例如，插入或缺失变体 (INDEL)、单核苷酸变体 (SNV) 或结构变体 (SV)）的一组一个或多个位置。查找候选变体位置可以包括例如检查参考中的所有位置，并确定至少一个序列读取区别于参考的位置。在确定候选位置时，也可以使用一些过滤器，例如，对于待视为候选变体位置的特定位置而言，针对必须不同于参考的读取数或百分比的过滤器。在其他实施例中，序列读取可能已经被扫描和分析（例如，通过独立的软件和/或硬件模块），在这种情况下，可以从耦接至计算机系统 102 的存储器或数据库中获得多个候选位置。在步骤 240，该方法可以从多个鉴别出的候选变体位置前进到下一个候选变体位置。

在步骤 250，在当前候选变体位置周围扩充序列读取和参考序列，以实现精确的多序列比对 (MSA)。结合图 3 可以更好地说明这一步骤，该图显示了多个示例性序列读取311、示例性参考序列 310 和示例性候选变体位置 350。如图 3 所示，扩充参考序列 310和序列读取 311 可包括将一个或多个空位或空格（例如，空格 325）插入参考序列 310 内和/或一个或多个序列读取 311 中，从而产生经扩充参考序列 320 和与经扩充参考序列320 精确比对的多条序列读取 321。

具体地，如图 3 所示，在一些实施例中，扩充可包括：确定或检测具体序列读取中的插入碱基；在该插入碱基位置处将空格插入参考样本中；在相同位置处将空格插入每一条在该相同位置处没有检测到插入的序列读取中（在一些情况下，这对应于除检测到插入的具体序列读取以外的所有序列读取）；以及对每个在每条序列读取中检测到的插入碱基重复该过程。应当理解，在检测到两个碱基或更长的潜在插入的位置中，可以分别插入两个或更多个空格。应当理解，可以使用其他合适的技术将空格插入序列读取和参考序列中，从而可以最终实现精确的多序列比对 (MSA)，这意味着经扩充序列读取内的所有碱基（或至少在没有出现缺失或插入的那些定位）都可以彼此精确地比对，并与经扩充参考序列精确比对。

在一些实施例中，扩充还可以包括检测具体序列读取中的缺失碱基，并在潜在缺失变体位置处将空格插入该序列读取中。应当理解，在检测到缺失碱基的情况下，不需要扩充参考序列。

还应理解，“插入空格”可以包括或可以不包括对存储了序列读取和参考序列的存储器的实际修改，并且可以代之以其他操作的固有或隐含部分（例如，方法 200 的其他步骤）以实现相同的最终结果。例如，代替修改序列读取和参考序列的原始定位，可以以将会造成新副本将空格存储在所有正确位置处的方式将序列读取和参考序列选择性地从原始定位处复制到独立的存储器定位处。如本文所用，“存储空格”可以意味着存储除对应于四种类型的核苷酸碱基（A、C、T 和 G）的四种值以外的任何随机或预定值（例如，0）。作为另一个示例，如下所述，当将序列读取和参考序列转换成计数矩阵时，可能不是在任何步骤中物理地插入空格，而是可以在算法上考虑空格。

在一些实施例中，在候选变体位置周围扩充序列读取和参考序列可以意味着仅扩充包括候选变体位置的可变长度或预定长度的一部分。然而，在其他实施例中，可以一次性地扩充可用序列读取和参考序列的完整长度，在这种情况下，步骤 250 可以仅实施一次（例如，在步骤 240 之前或甚至在步骤 230 之前），而不是针对每一个候选变体位置重复步骤 250。

在步骤 260，可以将经扩充序列读取的片段（即，窗口）转换（或“总结”）成样本矩阵。该片段可以具有固定的长度，并且可以将候选变体位置放置在中间或可以不放置在中间。例如，在图 3 所示的示例中，扩充片段长 16 列（位置），并且包括候选变体位置 350的左侧 7 列和右边 8 列，从而将候选变体位置 350 基本放置在扩充片段的中间。

在一些实施例中，样本矩阵可以是具有与正在被转换的片段相同的列数的计数矩阵。对于每种类型的核苷酸碱基（A、C、G 和 T），矩阵可以包括一个横行，每一横行表示相应核苷酸碱基在经扩充序列读取的每个位置（列）上出现的次数。参考图 3 的示例，样本矩阵331 中的横行“A”包括值“6000000306520000”，指示碱基“A”在经扩充读取的片段内的第一位置出现 6 次；在第八位置出现 3 次；在第十位置出现 6 次；在第十一位置出现 5 次；在第十二位置出现两次；并且在片段内的任何其他位置出现 0 次（未出现）。

在一些实施例中，样本矩阵还可以包括一个横行，该横行表示空格（在图 3 中表示为“-”）在经扩充序列读取内的每个位置出现的次数。例如，在图 3 中，样本矩阵 331 的第一横行（标记为“-”）包括值“0000014000140000”，指示经扩充序列读取在片段内的位置6、7、11 和 12 处分别包括 1、4、1 和 4 个空格。应当理解，尽管例如在图 3 的示例中将空格横行“-”显示为矩阵的第一横行，但是在其他实施例中，空格横行可以放置在矩阵内的其他地方（例如，最后）。还应当理解，横/列表示法是任意的，并且在其他实施例中，该表示法可以颠倒，这意味着样本矩阵可以具有 5 列（A、C、T、G、-），并且横行数对应于片段的长度。据此，通常可以将矩阵描述为具有至少五个“行”（横行或列），一个行表示每种核苷酸类型的计数，并且至少一个行表示经扩充序列读取内的空格计数。

在步骤 270，可以类似地将经扩充参考序列的与序列读取的片段相对应的片段转换成具有与样本矩阵相同维数的参考矩阵。例如，在图 3 中，扩充参考 320 转换为参考矩阵 330。类似于样本矩阵 331，参考矩阵 330 可以具有四个横行，指示每种类型的核苷酸在经扩充参考序列内的每个位置出现的次数，以及空格横行，指示在扩充参考序列内的每个位置处的空格数目。然而，应当理解，在其他实施例中，参考矩阵 330 可以具有与图 3所示不同的参考序列不同碱基的总结/表示法。

在图 3 所示的实施例中，参考矩阵 320 也被归一化以具有与样本矩阵 331 相同的值范围。如图 3 所示，归一化可以包括将每个计数乘以归一化因子，诸如序列读取的总数（在该示例中乘以 6）。可以理解，在其他实施例中，代替对参考矩阵进行归一化，可以对样本矩阵或对两个矩阵实施归一化。在其他实施例中，矩阵可能根本不需要被归一化。还应当理解，步骤 270 不需要在步骤 260 之后实施，并且在一些实施例中可以在步骤 260之前或与之并行实施。

在步骤 280，可以将参考矩阵和样本矩阵作为输入提供给受训深度神经网络，并且在步骤 290，可以获得受训深度神经网络的输出，其中，输出可以包括各种关于由样本矩阵表示的序列读取内包括的变体的预测，如下文更详细讨论。

本领域技术人员还应该理解，除非明确相反地指示，否则方法 200 中的步骤的数目和顺序不必限于如图 2a 所示的步骤的数目和顺序。因此，在不同的实施例中，方法 200的一些步骤可以被重新排序或并行实施；某些步骤可以设为可选或省略；并且可以添加一些附加步骤而不脱离本公开的精神和范围。

例如，图 2b 示出了另一个实施例，其中在步骤 210' 处，获得参考序列和多条相对应的序列读取。可以以任何顺序（即顺序地或同时地）获得参考序列和多条相对应的序列读取。接下来，在步骤 220' 中，将参考序列和多条相对应的序列读取进行比对。在一些实施例中，可以以比对的格式获得参考序列和多条相对应的序列读取，这可以消除实施第一比对步骤的需要。在步骤 230' 处，从经比对序列读取和参考序列中鉴别候选变体位置。在步骤 240' 处，以与上述相同的方式（即，通过插入空位以解释插入和/或缺失）在候选变体位置周围扩充多条序列读取和/或参考序列，以实现多条序列读取与参考序列的第二比对。

在一些实施例中，可对任一候选变体位置上的约 2 至 5、10、15、20、25、30、35、40、45、50、55、60、65、70、75、80、85、90、95 或 100 个碱基的窗口进行扩充和比对。在一些实施例中，可对候选变体位置任一侧的至少约 5、10、15、20、25、30、35、40、45、50、55、60、65、70、75、80、85、90、95 或 100 个碱基的窗口进行扩充和比对。在一些实施例中，可对候选变体位置任一侧的不超过约 5、10、15、20、25、30、35、40、45、50、55、60、65、70、75、80、85、90、95 或 100 个碱基的窗口进行扩充和比对。

在步骤 250' 处，从经扩充参考序列生成针对候选变体位置的参考矩阵，并从多条经扩充序列读取生成针对候选变体位置的样本矩阵。如上所述，矩阵包括在每个序列位置的每个碱基或空位（即，A、C、G、T 或空位）的频率。在步骤 260' 处，将参考矩阵和样本矩阵输入受训神经网络中。

在步骤 270' 处，神经网络确定在候选变体位置处是否存在变体类型。例如，神经网络可以包括多个卷积层，这些卷积层处理参考矩阵和样本矩阵，并生成由一个或多个分类器和回归器（即，突变类型的分类器和长度分类器，以及位置回归器）处理的输出，确定变体类型、尺寸和位置。在一些实施例中，神经网络可以包括多达约 5、10、15、20、25、30、35、40、45、50、55、60、65、70、75、80、85、90、95 或 100 层。在一些实施例中，卷积层可以被构造成块。在一些实施例中，块的数目是层的数目的大约一半。

如果参考序列和相对应的序列读取包括超过一个变体候选位置，则对于每个变体候选位置可以重复步骤 230' 至 260'。另外，可以通过获得额外的参考序列和额外的相对应序列读取来重复步骤 210' 至 260'。

尽管在图 2a 和图 2b 中示出的示例在针对种系变体调用的章节中进行了描述，但这些示例可以适用于体细胞变体调用。相反，在针对体细胞变体调用的章节中描述的实施例也可以适用于种系变体调用。例如，图 2b 中描述的多条序列读取可以细分为两个组：正常序列读取和肿瘤序列读取。正常序列读取可以从受试者的正常组织样本获得，而肿瘤序列读取可以从受试者的肿瘤组织样本获得。正常矩阵和肿瘤矩阵将会与参考矩阵一起生成，而不是生成单个样本矩阵。该方法的其余部分将以类似的方式进行。

图 4 示出了示例性实施例，其中参考矩阵 330 和样本矩阵 331（连同下面讨论的其他可选输入）被馈入受训深度神经网络 350 中，然后，深度神经网络 350 输出关于样本矩阵 331 所表示的序列读取的一个或多个预测。

在一些实施例中，网络 350 可以是卷积神经网络 (CNN)，其说明性示例在图 5中显示。如图 5 所示，网络 350 可以包括多个串联连接的卷积层 (510A-510I)，其最终馈入第一全连接层 520A，其输出馈入四个独立的输出各种预测的全连接层 520A-520C。在图5 的示例中，不同的卷积层 510 可以具有不同的过滤器尺寸（例如 1x3、3x3 或 5x5）并可以通过额外的处理层（未明确显示）互连，诸如整流线性单元 (ReLU)、池化层、批归一化(BN) 层等。此外，如图 5 所示，在一些实施例中，有时可以经由标识捷径连接 515 将输入到特定层中的输入连接到后续层的输出，以便帮助将信号保持在网络的更深层中。

在一些实施例中，网络 350 可以输出与经扩充序列读取的片段内的候选变体相关联的一个或多个预测。如图 5 的示例中所示，这样的预测之一可以在全连接层 520B 的输出端获得并且可以指示候选变体位置是否对应于变体，如果是，则指示该变体的类型。更具体地，在一些实施例中，预测可以包括与四种可能类型的变体相关联的至少四个概率值：NONE（无变体，即假阳性调用结果）；SNP/SNV（单核苷酸多态性/变体）；INS（插入变体）；和DEL（缺失变体）。基于这些概率值，可以确定哪种类型的变异是最可能的，以及具体预测的置信度。

在一些实施例中，神经网络 350 还可以输出（例如，在全连接层 520C 的输出端）由样本矩阵 331 表示的经扩充序列读取的片段内的变体的预测位置。在其他实施例中，可以省略该输出，并且可以基于针对给定的候选变体而选择序列读取的片段的方式，假设变体的位置是已知的（例如，在片段的中心），如上所述。

在一些实施例中，神经网络 350 还可以输出（例如，在全连接层 520D 的输出端）变体的预测长度。例如，如果预测的变体是一个碱基长（例如，SNP 变体或一碱基 DEL 或INS 变体），则可以输出“1”的长度；如果预测的变体是两碱基长的 DEL 或 INS 变体，则可以输出“2”的长度；以此类推。在一些实施例中，如果输出指示变体是 2 个碱基长或更长，则后处理步骤可以解析正被插入或缺失的确切序列。

在一些实施例中，神经网络 350 还可以输出（例如，在全连接层 520E 的输出端）与该变体相关联的预测基因型，例如，表示该变体的概率为 (1) 纯合参考（非变体）；(2)杂合变体（只有母系或父系副本具有变体）；(3) 纯合变体（其中两个副本具有相同变体）；或 (4) 其他（每个副本都有不同的变体）。应当理解，在其他实施例中，可以省略这些输出中的一些和/或可以添加额外的输出。

现在参考网络的输入，除了上面讨论的参考矩阵 330 和样本矩阵 331 之外，在一些实施例中，可以提供其他输入以进一步改善预测的准确性和/或扩展可以在网络的输出端获得的信息量。如图 4 的示例中所示，额外的输入可以包括，例如，变体位置矩阵340、覆盖率矩阵 343 和一个或多个比对特征矩阵 344。

在一些实施例中，所有输入矩阵可以作为一个大的三维矩阵提供到网络 350 中。矩阵的维数可以是例如 5 x s x k，其中 5 对应于横行/列的数目（-、A、C、T、G）；s 对应于片段的长度（例如，16）；并且 k 对应于上述的不同二维（例如，5 x s）矩阵的数目。在一些实施例中，k 可以高达 30，或者甚至更高。

变体位置矩阵 340 可以包括参考矩阵和样本矩阵内的候选变体位置的二维表示。例如，矩阵可以具有五个横行（-、A、C、T、G），并且每一列可以表示片段内的位置。例如，变体位置矩阵 340 可以在对应于候选变体位置的列的所有横行中均包括一个值（例如，1），并且在所有其他列的所有横行中均包括另一个值（例如，0）。

覆盖率矩阵 343 可以表示序列读取的具体片段的覆盖率或深度。例如，在一些实施例中，覆盖率矩阵 343 可以在其所有元素中包括相同的值，该值表示片段内不同读取的覆盖率/深度（例如，平均覆盖率/深度）。在其他实施例中，覆盖率矩阵 343 可在不同列处包括不同的值，每个值表示每一列的相应覆盖率/深度。比对特征矩阵 344 可以表示各种与序列读取的质量及其比对相关联的度量。此类度量可以包括例如碱基质量、映射质量(mapping quality)、链偏向性、剪切信息等。这些比对特征可通过为序列读取提供更多情境来帮助改善预测准确性，例如，通过将序列碱基和/或其映射 (mapping) 的质量告知网络。网络 350 的额外输入（为简洁起见未示出）可以包括关于已知变体的各种数据。可以从例如 dbSNP、COSMIC、ExAC 等公共和/或私有知识库中获取此类数据。

应当理解，图 5 仅显示了网络 350 的一种示例性配置，并且具有其他合适的配置/架构的神经网络可以用于分析经扩充序列读取（由样本矩阵表示）并预测变体的类型及其其他特征，而不脱离本公开的范围和精神。这样的替代配置可以包括更少的层、额外的层、具有不同参数的层、额外的输入或输出、更少的输入或输出等等。此外，在一些实施例中，网络 350 可以不是卷积神经网络 (CNN)，但可以是另一种类型的深度神经网络(DNN)，即，在输入层与输出层之间具有多个层的另一种类型的人工神经网络 (ANN)，而不脱离本公开的范围和精神。

在一些实施例中，神经网络 350 也可以被训练以实施体细胞变体调用。图 6 示出了使用受训神经网络 350 进行体细胞变体调用的示例性方法 600。

在步骤 610 处，方法 600 包括例如以上文结合方法 200 的步骤 210 所述者类似的方式获得参考序列。在步骤 620 处，方法 600 包括获得多条肿瘤序列读取和多条相对应的正常（非肿瘤）序列读取。例如，在一些实施例中，肿瘤序列读取可以包括患者的肿瘤（癌性）组织的测序结果，而正常序列读取可以包括取自相同患者的正常（非癌性）组织的测序结果。在其他实施例中，正常序列读取可以包括取自不同患者的正常组织的测序结果。可以例如以上文结合方法 200 的步骤 220 所述类似的方式获得两种类型的序列读取。

在步骤 630 处，可以通过从存储器或公共/私有数据库获得预定的变体位置即一组白名单候选位置，或通过实施在肿瘤序列读取与参考序列之间的比较（例如，以上文结合方法 200 的步骤 230 所述的方式），来获得多个候选体细胞变体位置。

在步骤 640 处，该方法进行到多个候选体细胞变体位置之一。在步骤 660 处，在候选体细胞变体位置周围将正常序列读取的片段转换为正常样本矩阵。类似地，在步骤665 处，在候选体细胞变体位置周围将肿瘤序列读取的片段（在长度和定位上基本上对应于正常序列读取的片段）转化为肿瘤样本矩阵。在步骤 670 处，在候选体细胞变体位置周围将参考序列的片段（在长度和定位上基本上对应于正常序列读取的片段）转化为参考矩阵。

在一些实施例中，在将以上三个片段转换成它们的相应参考矩阵之前，可以使用上文结合方法 200 的步骤 250 描述的扩充技术来扩充它们中的每一个。在其他实施例中，可以在将片段转换成矩阵之前对它们实施不同类型的扩充技术；例如，可以使用 k-mer塌缩方法，其中将长度 L > k 的均聚物塌缩成长度 k，其中 k 为大于或等于 1 的任何数字。在其他实施例中，可以完全省略扩充步骤，并且可以将原始的、未处理的序列读取和参考序列直接转换为它们的相应矩阵。

在图 7 所示的示例中，实施类似于在方法 200 的步骤 250 实施的扩充。在该示例中，首先扩充参考序列 310 的片段、肿瘤序列读取 311-A 的片段和正常序列读取 311-B 的片段，以分别产生经扩充参考序列 320、经扩充肿瘤序列读取 321-A 和经扩充正常序列读取 321-B，然后分别转换为参考矩阵 330、肿瘤样本矩阵 331-A 和正常样本矩阵331-B。

现在回到方法 600，在步骤 680 处，将三个矩阵馈入（作为输入提供）受训神经网络中，并且在步骤 690 处，获得受训神经网络的输出，其中，输出至少包括对肿瘤序列读取（例如，在候选体细胞变体位置或其附近）中含有的体细胞变体类型的预测（概率估计）。

图 8 显示了用于体细胞变体调用的受训神经网络 350 的示例性图示。在该示例中，网络 350 接收参考矩阵 330、肿瘤样本矩阵 331-A 和正常样本矩阵 331-B 作为其输入。

应当理解，尽管图 8 显示了输入至网络 350 中的所有三种类型的矩阵，但在一些实施例中，参考矩阵 330 和/或正常样本矩阵 331-B 可以被省略并且不提供给网络350，而且网络 350 仍然可以实施体细胞变体调用。在省略正常样本矩阵的实施例中，可以从方法 600 的步骤 620 中省略获得正常序列读取，并且可以完全省略将正常序列读取转换成矩阵的步骤 660。

在一些实施例中，为了进一步改善预测的准确性和/或扩展可以在网络的输出端获得的信息量，网络还可以接收以下的额外矩阵作为输入：诸如表示候选体细胞变体位置的体细胞变体位置矩阵 340；表示肿瘤序列读取的片段和正常序列读取的片段的相应覆盖率的肿瘤覆盖率矩阵 343-A 和正常覆盖率矩阵 343-B；以及一个或多个比对特征矩阵344，该矩阵以类似于上文结合图 4 所述者类似的方式表示序列读取的特征。另外，尽管在图 8 中为了简洁起见未示出，但在一些实施例中，可以从诸如 dbSNP、COSMIC、ExAC 等公共和/或私有知识库获得关于已知变体的数据，并将其作为输入提供给网络 350。

在一些实施例中，为了进一步改善体细胞变体调用的准确性，还可以在其输入端向神经网络 350 提供一个或多个其他调用者的特征矩阵 345。这些矩阵可以表示（例如，总结）从已经处理序列读取的相同片段的一个或多个其他（例如，第三方）变体调用应用程序或算法获得的一个或多个特征。此类特征可以包括由其他应用程序/算法预测的变体类型、长度和位置，以及质量得分、映射得分、变体显著性得分和其他比对特征。如上文结合图4 所讨论的，在一些实施例中，馈入网络 350 的所有输入矩阵可以被“组合”并且提供为一个大的三维矩阵。

图 9 显示了示例性的受训神经网络 350，其接收上述各种输入并输出预测，该预测至少指示在其输入端由体细胞样本矩阵 331-A 表示的体细胞序列读取的片段中包含的体细胞变体的变体类型、位置和长度。在一些实施例中，网络 350 还可以包括输出变体的基因型的全连接层 520E（如图 5 所示），但是在体细胞变体调用的情况下，该基因型并不重要，因此可以从网络中省略。

基于图 9 和图 5 的示例，应当理解，在一些实施例中，具有相同或基本相同的架构的网络 350 可以用于预测种系变体和体细胞变体，取决于已经对其进行训练的数据类型（如下文所讨论的），并且取决于提供给它的输入类型。然而，在一些实施例中，可以针对种系变体调用或体细胞变体调用来修改和优化网络 350 的架构。

本领域技术人员应当理解，在神经网络 350 可以开始实施准确的种系变体调用或体细胞变体调用之前，首先需要分别基于种系或体细胞训练数据（即，训练序列）对其进行训练。训练可以包括，例如，在大量训练序列上实施方法 200 或 600 的所有步骤，但随后还可以向网络提供“基准真值”数据（例如，实际已知的变体类型及其位置长度、基因型等），以便在每次运行之后能够通过被称为“反向传播”的过程来调整其可训练的过滤器和其他可训练的参数，从而使网络逐渐将其输出误差最小化。

在一些实施例中，网络 350 可以基于具有充分表征的基准真值变体的基因组诸如 NA12878 基因组进行训练。在一些实施例中，代替真实基因组或除了真实基因组之外，可以使用各种基于模拟的策略来训练网络。例如，为了训练用于种系变体调用的网络，可以使用 Mu, J. C. 等人 (“VarSim: a high-fidelity simulation and validationframework for high-throughput genome sequencing with cancer applications.”Bioinformatics 31, 1469–1471 (2015)) 中描述的方法模拟具有给定的一组变体的合成样本。作为另一个示例，为了训练用于体细胞变体调用的网络，可以例如使用 Eving, A.等人 (“Combining tumor genome simulation with crowdsourcing to benchmarksomatic single- nucleotide-variant detection,” Nature methods (2015)) 中描述的方法，将其中已经掺入随机变体的正常样本馈入该网络中。替代地或额外地，可以将两个不同的正常种系样本以不同比例混合以生成合成的肿瘤/正常训练样本。作为另一个示例，可以通过在两个种系样本的变体定位之间切换读取来生成具有期望的等位基因频率分布的虚拟肿瘤/正常训练样本。

体细胞突变是癌症发生、进展和治疗中的关键标志。由于肿瘤正常样本的交叉污染、肿瘤异质性、测序伪影和覆盖率，因此准确检测体细胞突变具有挑战性。通常，有效过滤由上述问题引入的假阳性调用结果并精确保留难以捕捉的真阳性调用结果（这些调用结果可能在等位基因频率 (AF) 较低的情况下发生或在复杂性较低的区域发生），对于精确的体细胞突变检测算法至关重要。

迄今为止，已经开发了一系列工具来解决体细胞突变检测问题，包括 MuTect2

这里，我们通过利用深度卷积神经网络 (CNN) 解决了肿瘤测序数据统计模型的通用性和复杂性的局限性。最近，CNN 在不同结构域的分类问题中表现出了显著的性能，包括种系变体调用

我们引入了 NeuSomatic，这是第一个基于 CNN 的体细胞突变检测方法，它可以有效地利用源自序列比对以及其他方法的信号来准确地鉴别体细胞突变。与其他专注于种系变体的基于深度学习的方法不同，NeuSomatic 正在解决由于肿瘤样本的复杂性而在准确性方面尚未得到满足的更大需求。它可以直接从原始数据有效地捕获重要的突变信号，并且对于不同的测序技术、样本纯度和测序策略（诸如全基因组与靶标富集）一致地实现高精确率。

NeuSomatic 网络的输入是通过扫描肿瘤样本以及匹配的正常样本的序列比对而鉴别出的候选体细胞突变（图 10a 至图 10d）。通过其他方法报告的体细胞突变也可以包括在此候选列表中。对于每个候选基因座，我们构造一个 3 维特征矩阵 M（尺寸为 k × 5× 32），由 k 个通道组成，每个通道的尺寸为 5 × 32，以捕获来自以该基因座为中心的区域的信号。矩阵 M 的每个通道都有五个表示四种核苷酸碱基和空位碱基（“-”）的横行，以及 32 个表示候选定位周围的比对列的列。

图 10a 示出了用于给定的候选体细胞 SNV 的输入矩阵制备的虚拟示例。提取候选体细胞突变周围七个碱基的窗口中的序列比对信息。然后通过添加空位以解释读取中的插入而扩充参考序列。然后将经扩充比对总结为参考矩阵、肿瘤计数矩阵和正常计数矩阵。计数矩阵在比对的每一列中记录 A/C/G/T 和空位（“-”）字符的数目，而参考矩阵在每一列中记录参考碱基。然后通过覆盖率对计数矩阵进行归一化，以反映每一列中的碱基频率。保留了独立的通道来记录肿瘤和正常覆盖率。

图 10b 示出了输入 3 维矩阵和所提出的 NeuSomatic 网络架构。输入矩阵由参考通道、肿瘤和正常频率通道、覆盖率和位置通道以及其后的几个总结比对特征的通道组成。在集成模式下使用时，NeuSomatic 还包括用于其他各个方法特征的额外通道。NeuSomatic 网络架构由九个卷积层组成，这些卷积层由四个具有捷径标识连接的块构成。我们在最后一层使用两个 softmax 分类器和一个回归器来预测突变类型、尺寸和位置。

图 10c 示出了用于给定的候选体细胞缺失的输入矩阵制备的虚拟示例。提取候选体细胞突变周围 7 个碱基的窗口中的序列比对信息。然后通过添加空位以解释读取中的插入而扩充参考序列。然后将经扩充比对总结为参考矩阵、肿瘤计数矩阵和正常计数矩阵。计数矩阵在比对的每一列中记录 A/C/G/T 和空位（“-”）字符的数目，而参考矩阵在每一列中记录参考碱基。然后通过覆盖率对计数矩阵进行归一化，以反映每一列中的碱基频率。保留了独立的通道来记录肿瘤和正常覆盖率。

图 10d 示出了用于给定的候选体细胞插入的输入矩阵制备的虚拟示例。提取候选体细胞突变周围 7 个碱基的窗口中的序列比对信息。然后通过添加空位以解释读取中的插入而扩充参考序列。然后将经扩充比对总结为参考矩阵、肿瘤计数矩阵和正常计数矩阵。计数矩阵在比对的每一列中记录 A/C/G/T 和空位（“-”）字符的数目，而参考矩阵在每一列中记录参考碱基。然后通过覆盖率对计数矩阵进行归一化，以反映每一列中的碱基频率。保留了独立的通道来记录肿瘤和正常覆盖率。

前三个通道分别是参考通道、肿瘤频率通道和正常频率通道，它们总结了候选基因座周围的参考碱基以及该区域中不同碱基的频率。为了捕获插入，我们在候选基因座周围使用与读取比对中的插入相对应的空位扩充了参考序列（图 10a、图 10c 和图 10d）。因此，肿瘤频率矩阵和正常频率矩阵的每一列分别代表肿瘤样本和正常样本的相对应多序列比对 (MSA) 列中 A/C/G/T/空位碱基的频率。其余的通道总结了针对支持不同碱基的读取的其他特征，诸如覆盖率、碱基质量、映射质量、链偏向性和剪切信息。如果在集成模式下使用 NeuSomatic，则我们还将使用了用于由单独的体细胞突变检测方法报告的特征的额外通道。通过这种简洁而又全面的结构化表示，NeuSomatic 可以利用肿瘤、正常和参考中的必要信息将具有低 AF 的难以捕获的体细胞突变与种系变体以及测序错误区分开来。这种设计还可以在 CNN 中使用卷积过滤器来捕获矩阵子块中的情境模式。

为了与基因组问题中使用的其他 CNN 方法进行比较，DeepVariant

NeuSomatic 采用新颖的 CNN 结构，该结构可在给定的特征矩阵 M 的情况下预测候选体细胞突变的类型和长度（图 10b）。所提出的 CNN 由九个卷积层组成，该卷积层由四个具有捷径标识连接的块构成，其灵感来自 ResNet

由于 NeuSomatic 可以在独立模式和集成模式下使用，因此我们使用NeuSomatic-S 表示独立模式，而保留 NeuSomatic 来表示集成模式。我们将 NeuSomatic和 NeuSomatic-S 与最新的体细胞突变检测方法（包括 MuTect2

对于第一个合成数据集，如先前的研究

图 11a 和图 11b 是用于铂双样本混合物数据集的性能分析的图表。在此数据集中，使用了四种肿瘤和正常样本纯度的情况（50% T:100% N、70% T:95% N、50% T:95% N 和25% T:95% N）。图 11a 是用于精确率-召回率分析的图：置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在曲线上，并用实心圆标记。图11b 是用于针对不同 INDEL 尺寸的 INDEL 准确性（F1 得分）的性能分析的曲线图。图11c 是显示不同体细胞突变检测方法在铂双样本混合数据集上的性能的表。对于每种方法，我们报告精确率-召回率曲线中达到最高 F1 得分的质量得分阈值的精确率、召回率和F1 得分。（RC：召回率，PR：精确率，F1：F1 得分）。

对于第二个合成数据集，我们使用 ICGC-TCGA DREAM 挑战阶段 3 和阶段 4 数据集

图 12a 和图 12b 是用于 DREAM 阶段 3 数据集的性能分析的曲线图。在此数据集中，使用了五种肿瘤和正常样本纯度的情况（100% T:100% N、50% T:100% N、70% T:95%N、50% T:95% N 和 25% T:95% N）。图 12a 是用于精确率-召回率分析的曲线图：置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在曲线上，并用实心圆标记。图 12b 是用于针对不同 INDEL 尺寸的 INDEL 准确性（F1 得分）的性能分析的曲线图。图 12c 是显示不同体细胞突变检测方法在 Dream 挑战阶段 3 数据集上的性能的表。对于每种方法，我们报告精确率-召回率曲线中达到最高 F1 得分的质量得分阈值的精确率、召回率和 F1 得分。（RC：召回率，PR：精确率，F1：F1 得分）。

图 13a 和图 13b 是用于 DREAM 阶段 4 数据集的性能分析的曲线图。在此数据集中，使用了五种肿瘤和正常样本纯度的情况（100% T:100% N、50% T:100% N、70% T:95%N、50% T:95% N 和 25% T:95% N）。图 13a 是用于精确率-召回率分析的曲线图：置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在曲线上，并用实心圆标记。图 13b 是用于针对不同 INDEL 尺寸的 INDEL 准确性（F1 得分）的性能分析的曲线图。图 13c 是显示不同体细胞突变检测方法在 Dream 挑战阶段 4 数据集上的性能的表。对于每种方法，我们报告精确率-召回率曲线中达到最高 F1 得分的质量得分阈值的精确率、召回率和 F1 得分。（RC：召回率，PR：精确率，F1：F1 得分）。

对于第三个合成数据集，如以前的研究

为了评定 NeuSomatic 在不同靶标富集上的性能，我们使用了来自 AshkenaziJewish trio

图 17a 示出了各种方法在外显子组样本混合物上的性能分析。在此数据集中，使用了四种肿瘤和正常样本纯度的情况（50% T:100% N、70% T:95% N、50% T:95% N 和 25%T:95% N）。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高F1 得分打印在曲线上，并用实心圆标记。这里的训练是 NeuSomatic、NeuSomatic-S 和SomaticSeq 基于外显子组数据的训练。

图 17b 示出了各种方法在靶标组样本混合物上的性能分析。在此数据集中，使用了四种肿瘤和正常样本纯度的情况（50% T:100% N、70% T:95% N、50% T:95% N 和 25% T:95% N）。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1得分打印在曲线上，并用实心圆标记。这里的训练是 NeuSomatic、NeuSomatic-S 和SomaticSeq 基于外显子组数据的训练。

图 17c 示出了不同体细胞突变检测方法在全外显子组样本混合数据集上的性能。对于每种方法，我们报告精确率-召回率曲线中达到最高 F1 得分的质量得分阈值的精确率、召回率和 F1 得分。（RC：召回率，PR：精确率，F1：F1 得分）。

图 17d 示出了不同体细胞突变检测方法在靶向组数据集上的性能。对于每种方法，我们报告精确率-召回率曲线中达到最高 F1 得分的质量得分阈值的精确率、召回率和F1 得分。（RC：召回率，PR：精确率，F1：F1 得分）。

将基于全基因组铂混合物数据训练的模型应用于两个靶标富集集，得到相似的性能，这证实了 NeuSomatic 的鲁棒性（图 18a 和图 18b）。与其他数据集类似，对于肿瘤纯度为 50% 的样本，将正常样本纯度从 100% 降低至 95%，将会最小限度地减少基于全外显子组数据集的 NeuSomatic F1 得分（平均约 0.3%），而 SomaticSeq 和 Strelka2 的 F1得分降低 >5%。

图 18a 示出了使用基于全基因组（铂数据，基因组混合物）和全外显子组（HG003-HG004 外显子组混合物）训练的模型对外显子组混合物数据集进行测试的性能分析。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分指示在图例中，并用实心圆标记在曲线上。

图 18b 示出了使用基于全基因组（铂数据，基因组混合物）和全外显子组（HG003-HG004 外显子组混合物）训练的模型对靶标组混合物数据集进行测试的性能分析。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分指示在图例中，并用实心圆标记在曲线上。

我们进一步评估了 NeuSomatic 在错误率较高的读取上的性能，特别是在来自长读取测序平台的读取上的性能。我们使用基于原始 PacBio 读取模拟的具有 20%、30% 和50% AF 体细胞突变的肿瘤-正常配对样本（图 14a、图 14b 和图 14c）。NeuSomatic 鉴别出的体细胞 SNV 和 INDEL 的 F1 得分分别高达 98.1% 和 86.2%，其性能优于VarDict

值得注意的是，对于不同数据集中的各种 INDEL 尺寸，NeuSomatic 始终优于其他方法（图 11b、图 12b、图 13b、图 14b、图 16c 和图 19）。对于大（大于 50 个碱基）的INDEL，由于大多数具有体细胞 INDEL 的短读取都是软跳过的，因此 INDEL 信息会在堆积计数矩阵中丢失。对于此类情况，NeuSomatic 受益于其他方法的预测，因为诸如 VarDict和 MuTect2 之类的一些方法使用本地程序集进行预测。

图 14a 和图 14b 是不同体细胞突变检测方法在 PacBio 数据集上的性能分析的曲线图。在此数据集中，使用了三种肿瘤和正常样本纯度的情况（70% T:95% N、50% T:95% N 和 25% T:95% N）。图 14a 是用于精确率-召回率分析的曲线图：置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在曲线上，并用实心圆标记。图 14b 是用于针对不同 INDEL 尺寸的 INDEL 准确性（F1 得分）的性能分析的曲线图。图 14c 显示，对于每种方法，精确率-召回率曲线中达到最高 F1 得分的质量得分阈值的精确率、召回率和 F1 得分。（RC：召回率，PR：精确率，F1：F1 得分）。图 19 示出了在Dream 阶段 3、Dream 阶段 4、铂双样本混合物、铂肿瘤掺入、PacBio 和外显子组数据集中的基准真值 (ground truth) INDEL 的尺寸分布。负尺寸对应于缺失。

对于所有讨论的数据集，我们还评定了通过不同体细胞突变检测方法进行 INDEL调用的性能，使用更为宽松的标准，即简单地正确预测体细胞 INDEL 的位置（并忽略确切的 INDEL 序列）。再次，我们观察到 NeuSomatic 相对于其他方案的相似优势，指示主要的改善是由所提出的 CNN 结构而不是后处理 INDEL 解析步骤贡献的（图 20a 和图 20b）。

图 20a 示出了在三种肿瘤纯度（50%、30% 和 20%）以及 95% 正常纯度的情况下，基于针对 PacBio 数据集的所预测的体细胞突变的位置和类型（而忽略了确切的所预测INDEL 序列的准确性）的 INDEL 的性能分析。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在曲线上，并用实心圆标记。

图 20b 示出了针对 Dream 阶段 3、Dream 阶段 4、铂双样本混合物、全外显子组和铂肿瘤掺入数据集的基于所预测的体细胞突变的位置和类型的 INDEL 的性能分析（而忽略了确切的所预测 INDEL 序列的准确性）。对于前四个数据集，使用了三种肿瘤纯度方案（70%、50% 和 25%），而正常样本的纯度为 95%。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在曲线上，并用实心圆标记。

为了评估序列覆盖率对不同技术的影响，我们对全外显子组数据集进行了下采样，以获得序列覆盖率在 20× 和 100× 范围内的样本（图 15a 和图 15b）。对于不同的覆盖率，NeuSomatic 始终优于其他技术。对于低覆盖率样本，随着问题更具挑战性，改善也有所增加。此外，将覆盖范围从 100× 减少到 50× 对 NeuSomatic 的影响很小（对于SNV 约为 1.5%，对于 INDEL 约为 5%），而 SomaticSeq 的 SNV 和 INDEL 的 F1 得分均下降约 20%，并且 Strelka2 的 SNV 的 F1 得分下降了约 13% 且 INDEL 的 F1 得分下降了约 15%。该分析显示了 NeuSomatic 对覆盖率扰动的鲁棒性以及在挑战性场景中的优势，这可以在较低覆盖率下看到。

图 15a 是用于序列覆盖率对全外显子组样本混合物数据集的影响的性能分析的曲线图。在该示例中，肿瘤的纯度为 50%，正常样本的纯度为 95%。Y 轴示出了每种算法在20× 到 100× 范围内进行的样本比对所能获得的最高 F1 得分。

图 15b 是用于序列覆盖率对全外显子组样本混合物数据集的影响的性能分析的曲线图。在该示例中，肿瘤的纯度为 50%，正常样本的纯度为 95%。肿瘤和正常比对的覆盖范围是 20x 至 100x。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分指示在图例中，并用实心圆标记在曲线上。

我们通过基于 DREAM 挑战阶段 3 数据集针对不同纯度进行训练和测试，评定NeuSomatic 针对特定纯度进行训练的鲁棒性。图 21 显示，即使当我们在非常不同的肿瘤纯度上进行训练和测试时，其性能也仅会轻微降低。我们还观察到，使用从多个肿瘤纯度汇总的数据进行的训练与对在靶标肿瘤纯度上进行的训练一样好（图 21）。这表明并入了多个肿瘤纯度的训练集足以获得对肿瘤纯度变化具有鲁棒性的模型。

图 21 示出了针对 Dream 挑战阶段 3 数据集的交叉样本训练的性能分析。我们使用基于不同纯度训练的 NeuSomatic 模型以及基于来自所有不同纯度的集合输入训练的模型，对肿瘤纯度分别为 70%、50% 和 25% 的每个样本进行了测试。置信度或质量得分用于推导出精确率-召回率曲线。通过每种算法获得的最高 F1 得分打印在图例中，并用实心圆标记在曲线上。

在不存在针对体细胞突变

图 22a 显示了不同体细胞突变检测方法在真实数据集 COLO-829 上的性能。图22b 显示了不同体细胞突变检测方法在真实数据集 CLL1 上的性能。图 22c 显示了不同体细胞突变检测方法在真实数据集 TCGA-AZ-6601 上的性能。

为了证明 NeuSomatic 在云上的可扩展性和经济效益，我们还使用集成模式和独立模式在 Microsoft Azure 云平台上处理了来自 TCGA 的 261 个全外显子组测序癌症样本（图 23）。这些样本是从多种癌症类型中采集的，包括结肠直肠腺癌、卵巢浆液性腺癌、宫颈鳞状细胞癌和宫颈内腺癌。尽管云平台使我们能够按需自动启动计算实例 (computeinstances)，但集成模式和独立模式分别平均需要 2.42 小时和 0.72 小时的时间来处理每个样本。使用 Azure 的可抢占型计算实例（使用标准的 H16 实例类型，每个实例具有16 个核心），每个样本的处理成本较低，对于集成模式和独立模式分别为 0.77 USD 和0.23 USD。我们还通过与这些样本中的 44,270 个经过验证的 SNP 进行比较，评定了NeuSomatic 在这些样本上的准确性，这为我们提供了分别为 98.9% 和 97.2% 的整体模式召回率和独立模式召回率。因此，NeuSomatic 不仅可以用于不同的测序技术或测序策略，而且可以在各种计算平台（包括本地 HPC 集群和弹性云计算架构）上运行。

图 23 是用于 Microsoft Azure 实验的 261 个 TCGA 癌症样本的列表。样本取自三种癌症类型：结肠直肠腺癌 (COAD)、卵巢浆液性腺癌 (OV) 以及宫颈鳞状细胞癌和宫颈内腺癌 (CESC)。

NeuSomatic 是第一个基于深度学习的体细胞突变检测框架，它是高性能且通用的。当使用相同的 CNN 架构时，它可以针对改变的肿瘤纯度跨从合成到真实的多个数据集、从全基因组到靶向的多种测序策略以及从短读取测序到高错误率的长读取测序的多种测序技术范围而实现最佳的准确性。具体而言，对于低肿瘤纯度和低等位基因频率，NeuSomatic 明显优于其他最新的体细胞突变检测方法，从而证明了其解决疑难问题的能力。NeuSomatic 利用卷积神经网络的有效实现来快速、准确地解决体细胞突变检测问题。它使用肿瘤/正常比对信息的新颖总结作为一组输入矩阵，可以有效地捕获基因组情境中的主要信号。基于这些矩阵训练所提出的 CNN 架构，可以直接从原始数据中学习特征表示。从观察到的训练数据中学到的深层特征可以准确地鉴别重要的突变标志，这些标志可以将真实的调用结果与由测序错误、交叉污染或覆盖率偏差所引入的伪影区分开来。我们相信，NeuSomatic 通过为体细胞突变检测提供非常广泛适用的方法，可以极大地改进最新技术。

阶段 3 数据由正常样本和肿瘤样本组成，这些样本是通过将 7903 个 SNV 和7604 个 INDEL 突变计算地掺入到具有三个不同 AF（50%、33% 和 20%）的相同正常样本的健康基因组中而构建的，以创建合成但真实的肿瘤正常样本对。阶段 4 的数据具有相似的构造，但是在 30% 和 15% AF 的两个亚克隆中具有 16,268 个 SNV 和 14,194 个INDEL。然后，我们通过将 95% 正常读取与 5% 肿瘤读取混合而构建不纯的正常读取。我们还通过分别以 100:0、70:30、50:50 和 25:75 的比率混合肿瘤和正常读取来构建四种肿瘤混合物。因此，跨这四个肿瘤混合比率的体细胞突变，对于阶段 3 数据，AF 的范围为 5%至 50%，而对于阶段 4 数据集，AF 的范围为 3.75% 至 30%。

我们下载了 200× 铂基因组样本 NA12878 和 NA12877 及其真值种系变体(v2017-1.0)

我们还通过在 NA12878 中的杂合子变体定位和纯合子变体定位（它们是NA12877 中的参考调用结果）处，将来自经下采样（至 50× 覆盖率）NA12878 的读取随机掺入到经下采样（至 50× 覆盖率）NA12877 数据中，从而生成了另一个 50× 虚拟肿瘤样本。对于每个变体，我们用均值 [0.05、0.1、0.2、0.3] 随机分配来自二项式分布的掺入读取的频率。因此，取决于变体的接合性，均值体细胞突变 AF 的范围为 2.5% 至 30%。为了避免在真值集中产生歧义，我们仅使用相关配对末端读取不与任何其他变体重叠的变体（316,050 个 SNV 和 46,978 个 INDEL）。这产生了具有来自 NA12878 的读取的被污染肿瘤样本。我们还使用 NA12877 的另一个独立的经下采样（至 50×）数据用作纯正常样本。

对于这两个实验，均下载了 FastQ 文件和真值种系变体，并先后与 BWA-MEM(v0.7.15)

我们使用 CLL1 慢性淋巴细胞白血病数据集

COLO-829 数据集由 80× 全基因组测序肿瘤样本及其匹配的 60× 正常血液COLO-829BL 样本组成，并具有 454 个经过验证的体细胞 SNV。CLL1 具有覆盖率分别为53× 和 42× 的全基因组测序肿瘤样本和匹配的正常样本，并具有 961 个已公开的体细胞 SNV。

TCGA-AZ-6601

对于真实数据，我们将外推精确率计算为所预测的体细胞突变的百分比，该体细胞突变已经被至少两种独立方法调用，或者已经在 COSMIC 数据库中报告为至少两个相同癌症类型样本中已证实的体细胞突变。然后，我们基于召回率的调和平均数和该外推精确率来计算外推的 F1 得分。

为了用不同靶标富集实验评定 NeuSomatic，我们使用了来自 Ashkenazi Jewishtrio

为了验证在靶标组上的性能，我们将上述比对和真值数据限制为 Illumina 的TruSight 遗传疾病组 BED 文件（216 个 SNV，5 个 INDEL）。由于靶标组区域中真实INDEL 的数量有限，我们仅评估在 SNV 上的性能。

对于长读取分析，我们下载了 HG002 样本的高置信度种系变体（瓶中基因组v3.3.2 版）(ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/)

第一步，我们扫描肿瘤读取比对以找到具有突变证据的候选定位。由于基因组区域的复杂性或测序伪影，这些位置中的许多位置要么具有种系变体，要么被错误调用。我们将一组自由过滤器应用于一组候选定位，以确保此类定位的数量合理。通常，对于 SNV，我们要求 AF≥0.03 或超过两个读取支持 SNV 并且 Phred 缩放的碱基质量得分大于 19（对于真实 WES 数据集，大于 14）作为最低要求。对于 1 个碱基的 INDEL，我们要求 AF≥0.02 或超过一个读取支持。对于大于 1 个碱基的 INDEL，我们要求 AF≥0.03。对于集成方法，我们还包括通过其他体细胞突变检测方法检测到的任何体细胞突变作为输入候选。对于 PacBio 数据集，对于大于 1 个碱基的 SNV 和 INDEL，我们使用 AF≥0.1，对于 1个碱基的 INDEL，AF≥0.15。

对于 DREAM 挑战数据集，我们从输入候选中排除了 dbSNP

对于每个候选位置，我们制备了一个 3 维矩阵 M，其具有尺寸为 5 × 32 的 k个通道（图 10a、图 10c 和图 10d）。每个通道中的五个横行对应于四种 DNA 碱基 A、C、G和 T 以及空位字符（“-”）。矩阵的 32 列中的每一列代表比对的一列。

对于每个候选定位，我们提取了肿瘤和正常读取比对。如图 10a 所示，然后，我们将肿瘤和正常样本与参考的读取比对作为 MSA。为此，当读取中存在插入时，我们通过在参考序列上添加空位而扩充参考序列。必须注意，此过程不需要对输入 BAM 文件的原始读取比对进行任何进一步的重新比对，而只需通过在发生插入的位置处分配其他列，即可将比对重构为 MSA 格式。如果在特定位置之后的多次读取中存在多个不同插入，我们将其视为经左比对序列并将它们放在同一组列中（例如，参见图 10a 中虚拟示例的第九列中的 A碱基插入和 C 碱基插入）。使用此读取表示法，我们在每一列中找到 A/C/G/T/- 字符的频率，并为肿瘤和正常样本（矩阵 M 中的通道 C2 和 C3）记录单独的矩阵。在通道 C1 中，我们在每一列中记录参考碱基（或空位）。通道 Ci (4 ≤ i ≤ k) 记录肿瘤和正常样本中的其他比对信号，诸如覆盖率、碱基质量、映射质量 、链、剪切信息、编辑距离、比对得分和配对末端信息。例如，对于碱基质量通道，每个样本都有一个尺寸为 5 × 32 的矩阵，该矩阵记录在每列中具有给定碱基（对于给定的横行）的读取的平均碱基质量。另举一例，对于编辑距离通道，每个样本都有一个尺寸为 5 × 32 的矩阵，该矩阵记录在每列中具有给定碱基（对于给定的横行）的读取的平均编辑距离。矩阵 M 的一个通道专用于指定候选者所在的列。在当前实现中，我们在独立模式的 NeuSomatic-S 方法中总共使用了 26 个通道。

对于 NeuSomatic 的集成扩展，我们还包括额外的通道来捕获由所使用的六个单独方法中的每个方法报告的特征。在该实现中，我们使用了 93 个额外通道来表示从其他方法中提取的特征以及由 SomaticSeq 报告的比对。因此，NeuSomatic 的集成模式具有针对每个候选矩阵的 119 个输入通道。

对于每个候选定位，我们在候选者周围七个碱基的窗口中报告比对信息。我们保留 32 列，以考虑到具有插入的经扩充比对。在少数情况下，如果插入量较大，则 32 列可能不足以表示比对。对于此类情况，我们将截短插入，以便可以记录候选者附近的至少三个碱基。

所提出的 CNN（图 10b）由如下构造的九个卷积层组成。输入矩阵被馈入具有 64个输出通道、1 × 3 核尺寸和 Relu 激活的第一卷积层中，之后被馈入到批归一化和最大池化层中。然后，这一层的输出被馈入一组四个块中，这些块具有类似于 ResNet 结构的捷径标识连接。这些块由一个具有 3 × 3 个核的卷积层和之后的批归一化和一个具有 5× 5 个核的卷积层组成。在这些捷径块之间，我们使用批归一化层和最大池化层。最终块的输出被馈入尺寸为 240 的全连接层中。然后将所得的特征向量馈入两个 softmax 分类器和一个回归器中。第一个分类器是 4 向分类器，可从非体细胞、SNV、插入和缺失这四类中预测突变类型。第二个分类器根据 0、1、2 和 ≥3 这四个类别来预测所预测突变的长度。非体细胞调用结果标注为零尺寸突变，SNV 和 1 个碱基的 INDEL 标注为尺寸 1，而 2个碱基的且尺寸 ≥3 的 INDEL 分别标注为 2 和 ≥3 尺寸的突变。回归器预测矩阵中突变的列，以确保预测靶向正确的位置，并使用平滑的 L1 损失函数对其进行优化。

CNN 的参数少于 900 K，这使我们能够通过使用大的批次规模来实现高效的实现。在配备 8 台 Tesla K80 Nvidia GPU 的机器上，全基因组训练过程耗时约 8 小时。

对于 DREAM 挑战数据集、铂数据集和靶标富集数据集，我们将基因组区域随机分为 50% 训练集和 50% 测试集。对于 PacBio 数据集，我们基于 1 号染色体上的 120 兆碱基区域训练 NeuSomatic，并在 22 号染色体上对其进行测试。

对于每个数据集，我们合并了来自所有不同肿瘤/正常纯度场景的所生成的训练输入矩阵，并使用合并的集合来训练网络。然后，我们将该统一受训模型应用于每个独立的肿瘤/正常纯度设置下的测试。

DREAM 挑战数据集具有阶段 3 的 15,507 个体细胞突变和阶段 4 的 30,462个体细胞突变。为了更好地进行网络训练，我们使用 BAMSurgeon

我们使用 1000 个的批次规模和 SGD 优化器以 0.01 的学习速率和 0.9 的动量对网络进行训练，并且每 400 代训练就将学习速率乘以 0.1。

一般而言，由于输入候选定位的非体细胞（参考或种系）调用结果比真实体细胞突变多得多，因此在每一代训练中，我们都使用训练集中的所有真实体细胞突变并随机选择的非体细胞候选者，该非体细胞候选者的数目是真实体细胞突变数的两倍。我们使用加权的 softmax 分类损失函数来平衡每个类别中候选者的数量。对于 DREAM 挑战数据，由于我们在训练集中添加了更多的合成突变，因此我们增加了非体细胞类别的权重，以在测试集中实现更高的精确率。

为了评定合成靶标富集数据集，我们使用了全外显子组和全基因组数据作为训练集。

为了测试真实 WGS 样本 CLL1 和 COLO-829，我们使用了基于 DREAM 挑战阶段3 训练的 SomaticSeq 和 NeuSomatic 模型。对于真实 WES 样本 TCGA-AZ-6601，我们使用来自另一个 TCGA WES 数据集 TCGA-AZ-4315

对于超参数调优，我们在 DREAM 挑战阶段 3 实验中使用了 10% 的基因组，并在所有其他实验中使用了推导出的参数。

我们进一步探索了不同的网络架构，诸如具有 4 至 16 个 ResNet 块的预激活ResNet 架构（包括 ResNet-18 和 ResNet-34 架构）（图 24a 至图 24e），以及所提出的残差 NeuSomatic 架构的一些变体（图 24f 至图 24m）。为了评估这些网络，我们将 DREAM阶段 3 数据集中的训练数据分为两个相等的部分，使用其中的一个来训练不同架构，并使用另一个以独立模式评估这些不同架构。图 25 在准确性、网络参数数量、存储器使用和速度方面比较了这些架构。通常，与传统的体细胞突变检测方法相比，所有这些网络都可以获得相对较高的准确性。这一发现揭示了所提出的数据总结方法的重要性，该方法在候选者的基因组情境中捕获主要信号并有助于卷积网络在体细胞突变检测问题上的有效实现。与具有所提出的残差块的那些（图 24f 至图 24m）相比，具有两个 3 × 3 卷积过滤器（图24a 至图 24e）的缺省 ReSNet 架构具有较低的平均准确性。另外，具有跨步卷积的网络（图 24a 至图 24g）具有更多的网络参数和运行时要求。总而言之，尽管每种网络架构在一些比较方面均显示出优势，但我们还是选择了所提出的 NeuSomatic 网络架构（图 10b；图24k）作为缺省网络架构以兼顾所有这些因素，而其他网络则可以轻松实现根据用户的使用情况和时间/计算约束对其进行调整。

图 24a 至图 24m 示出了所测试的不同网络架构。图 24a 至图 24e 示出了具有不同数量的预激活残差块的具有缺省 3x3 卷积层的 ResNet 架构。在这里，跨步卷积与通道扩展一起使用。图 24f 和图 24g 示出了具有 3x3 和 5x5 的卷积层以及一些扩张卷积的多个定制的残差块。在这里，跨步卷积与通道扩展一起使用。图 24h 示出了具有 3x3 和5x5 的卷积层以及一些扩张卷积的四个定制的残差块。这里，不使用跨步卷积。图 24i 至图 24m 示出了具有不同残差块和全连接尺寸的 NeuSomatic 残差架构。

图 25 示出了图 24a 至图 24m 所示的不同网络架构的性能分析。在这里，所有网络都在经过了 600 代训练之后以 1000 个的批次规模进行评定。

在我们的分析中，我们使用了缺省设置的 Strelka2 (v2.8.4)、Mutect2(v4.0.0.0)、SomaticSniper (v1.0.5.0)、MuSE (v1.0rc)、VarDict (v1.5.1)、VarScan2(v2.3.7) 和 SomaticSeq (v2.7.0) 体细胞突变检测算法。

我们使用 VarDict 作为基于 PacBio 数据的 NeuSomatic 的替代方法。为了能够在高错误率的长读取上检测体细胞突变，我们使用了具有“ −m 10000 −Q 1 −q 5 −X 1−c 1 −S 2 −E 3 −g 4 −k 0”参数设置的 VarDict。此外，与 NeuSomatic 一样，对于 SNV，使用 AF≥0.1，对于 INDEL，使用 AF≥0.15。

为了训练 SomaticSeq，我们还遵循了与 NeuSomatic 所用相同的 50% 训练/测试区域划分。此外，与 NeuSomatic 一样，对于每个数据集，我们将来自所有不同肿瘤/正常纯度场景的训练数据组合，以训练 SomaticSeq SNV 分类器和 INDEL 分类器。然后使用这些统一的分类器来预测每个单独的肿瘤/正常纯度设置。

对于精确率-召回率分析，基于每种工具分配的置信度或质量得分对体细胞突变进行分类。对于 MuSE，我们使用层分配作为排序标准。对于 VarDict、VarScan2、MuTect2、Strelka2 和 SomaticSniper，我们分别使用了 VCF 文件中报告的 SSF 值、SSC 值、TLOD值、SomaticEVS 值和 SSC 值进行排序。对于 SomaticSeq 和 NeuSomatic，我们使用了QUAL 领域中的体细胞突变质量得分。NeuSomatic 基于 CNN 预测的概率报告所预测的体细胞突变的质量得分。

为了分析真实样本的性能，我们使用了来自不同方法的 PASS 体细胞调用结果（对于 VarDict，我们仅限于具有 StrongSomatic 状态的调用结果）。对于 NeuSomatic，我们将 WGS 的质量得分阈值设为 0.97，将 WES 的质量得分阈值设为 0.6。

对于全基因组数据，在双 14 核 Intel Xeon CPU E5-2680 v4 2.40 Ghz 机器上，扫描 30× 肿瘤和正常比对以找到候选者、提取特征并制备输入矩阵可能耗费约 3.9小时。在配备 8 台 Tesla K80 Nvidia GPU 的机器上，全基因组训练过程可能耗费约 8小时（规模为 580,000 的每一代训练约 90 秒）。取决于候选体细胞突变的截止值 AF，使用 8 台 Tesla K80 Nvidia GPU 计算对于 30× 全基因组数据上候选突变的网络预测可能耗费约 35 分钟（AF 截止值为 0.05，3.9 M 个候选者）至约 100 分钟（AF 截止值为0.03，11.5 M 个候选者）。对于 125× 全外显子数据，网络预测的整个扫描、制备和计算可能耗费约 30 分钟。将使用 NeuSomatic 集成方法和独立方法（仅 CPU 模式）在 125× 全外显子组数据集和 30× 全基因组数据集上预测体细胞突变的端到端运行时与图 26a 和图 26b 中的其他体细胞突变检测技术进行比较。

图 26a 示出了不同体细胞突变检测算法的运行时比较。显示了 CPU 核心-小时，以预测 125x 全外显子组测序数据集上的体细胞突变。图 26b 示出了不同体细胞突变检测算法的运行时比较。显示了 CPU 核心-小时，以预测 30x 全基因组测序数据集上的体细胞突变。

NeuSomatic 用 Python 和 C++ 编写。它的深度学习框架是使用 PyTorch 0.3.1实现的，以使 GPU 应用程序能够进行训练/测试。根据知识共享署名非商业共享(Creative Commons Attribution-NonCommercial-ShareAlike) 4.0 国际许可，源代码位于 https://github.com/bioinform/neusomatic。本文的结果是基于 NeuSomatic v0.1.3的。

有关实验设计的更多信息，请参见与本文链接的《自然研究报告汇编》(NatureResearch Reporting Summary)。

这项研究的序列数据是从各种来源收集的，即，欧洲核苷酸档案库（登录号：PRJEB3246；https://www.ebi.ac.uk/ena）、序列读取档案库（登录号：SRX1026041；https://www.ncbi.nlm.nih.gov/sra）、国际癌症基因组联合体（项目：ICGC-TCGA DREAM突变调用挑战，受控访问：https://icgc.org/）、癌症基因组图集（登录号：TCGA-AZ-6601、TCGA-AZ-4315；受控访问：https://gdc.cancer.gov/）、欧洲基因组-现象档案库（登录号：EGAS00000000052、EGAD00001000023；受控访问：https://www.ebi.ac.uk/ega/）和瓶中基因组（登录号：HG002、HG003、HG004；ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/）。使用 https://github.com/bioinform/neusomatic/blob/paper/etc/data_scripts.zip 的脚本从上述数据集中生成了合成数据。所有其他相关数据可应要求提供。

1.Cibulskis, K. 等人，Sensitive detection of somatic point mutationsin impure and heterogeneous cancer samples.Nat. Biotechnol. 31, 213 (2013)。

2.Fan, Y. 等人，MuSE: accounting for tumor heterogeneity using asample-specific error model improves sensitivity and specificity in mutationcalling from sequencing data.Genome Biol. 17, 178 (2016)。

3.Lai, Z. 等人，VarDict: a novel and versatile variant caller fornext-generation sequencing in cancer research.Nucleic Acids Res. 44, e108–e108 (2016)。

4.Koboldt, D. C. 等人，VarScan 2: somatic mutation and copy numberalteration discovery in cancer by exome sequencing.Genome Res. 22, 568–576(2012)。

5.Kim, S. 等人，Strelka2: fast and accurate calling of germline andsomatic variants.Nat. Methods 15, 591–594 (2018)。

6.Larson, D. E. 等人，SomaticSniper: identification of somatic pointmutations in whole genome sequencing data.Bioinformatics 28, 311–317 (2011)。

7.Wang, Q. 等人，Detecting somatic point mutations in cancer genomesequencing data: a comparison of mutation callers.Genome Med. 5, 91 (2013)。

8.Alioto, T. S. 等人，A comprehensive assessment of somatic mutationdetection in cancer using whole-genome sequencing.Nat. Commun. 6, 10001(2015)。

9.Roberts, N. D. 等人，A comparative analysis of algorithms forsomatic SNV detection in cancer.Bioinformatics 29, 2223-2230 (2013)。

10.Fang, L. T. 等人，An ensemble approach to accurately detect somaticmutations using SomaticSeq.Genome Biol. 16, 197 (2015)。

11.Poplin, R. 等人，A universal SNP and small-indel variant callerusing deep neural networks.Nat. Biotechnol. 36, 983 (2018)。

12.Luo, R.、Sedlazeck, F. J.、Lam, T.-W. 和 Schatz, M.，Clairvoyante: amulti-task convolutional deep neural network for variant calling in singlemolecule sequencing. https://www.biorxiv.org/content/early/2018/04/28/310458(2018)。

13.Simple convolutional neural network for genomic variant callingwith tensorflow. https://towardsdatascience.com/simple-convolution-neural-network-for-genomic-variant-calling-with-tensorflow-c085dbc2026f (2017)。

14.Esteva, A. 等人，Dermatologist-level classification of skin cancerwith deep neural networks.Nature 542, 115–118 (2017)。

15.Torracinta, R. 等人，Adaptive somatic mutations calls with deeplearning and semi-simulated data. https://www.biorxiv.org/content/early/2016/10/04/079087 (2016)。

16.Running DeepVariant. https://cloud.google.com/genomics/docs/tutorials/deepvariant (2018)。

17.He, K.、Zhang, X.、Ren, S. 和 Sun, J.，Deep residual learning forimage recognition.Proc. IEEE conf. comput. vis. pattern recognit.770–778(2016) http://openaccess.thecvf.com/content_cvpr_2016/html/He_Deep_Residual_Learning_CVPR_2016_paper.html。

18.Eberle, M. A. 等人，A reference data set of 5.4 million phasedhuman variants validated by genetic inheritance from sequencing a three-generation 17-member pedigree.Genome Res. 27, 157-164 (2017)。

19.Ewing, A. D. 等人，Combining tumor genome simulation withcrowdsourcing to benchmark somatic single-nucleotide-variant detection.Nat.Methods 12, 623 (2015)。

20.Zook, J. M. 等人，Extensive sequencing of seven human genomes tocharacterize benchmark reference materials.Sci. Data 3, 160025 (2016)。

21.Xu, C.，A review of somatic single nucleotide variant callingalgorithms for next-generation sequencing data.Comput.Struct. Biotechnol. J.16, 15–24 (2018)。

22.Zook, J. M. 等人，Integrating human sequence data sets provides aresource of benchmark SNP and indel genotype calls.Nat. Biotechnol. 32, 246(2014)。

23.Puente, X. S. 等人，Whole-genome sequencing identifies recurrentmutations in chronic lymphocytic leukaemia.Nature 475, 101 (2011)。

24.Morse, H. G. 和 Moore, G. E.，Cytogenetic homogeneity in eightindependent sites in a case of malignant melanoma.Cancer Genet. Cytogenet.69,108–112 (1993)。

25.Pleasance, E. D. 等人，A comprehensive catalogue of somaticmutations from a human cancer genome.Nature 463, 191 (2010)。

26.Network, C. G. A. 等人，Comprehensive molecular characterization ofhuman colon and rectal cancer.Nature 487, 330 (2012)。

27.Grasso, C. S. 等人，Genetic mechanisms of immune evasion incolorectal cancer.Cancer Discov.8, 730-749 (2018)。

28.Li, H.，Aligning sequence reads, clone sequences and assemblycontigs with BWA-MEM. https://arxiv.org/abs/1303.3997 (2013)。

29.Van der Auwera, G. A. 等人，From FastQ data to high-confidencevariant calls: the genome analysis toolkit best practices pipeline.Curr.Protoc.Bioinforma.43, 10-11 (2013)。

30.Grossman, R. L. 等人，Toward a shared vision for cancer genomicdata.N. Engl. J. Med. 375, 1109–1112 (2016)。

31.Forbes, S. A. 等人，COSMIC: somatic cancer genetics at high-resolution.Nucleic Acids Res. 45, D777–D783 (2016)。

32.Chaisson, M. J. 等人，Resolving the complexity of the human genomeusing single-molecule sequencing.Nature 517, 608 (2015)。

33.Mu, J. C. 等人，VarSim: a high-fidelity simulation and validationframework for high-throughput genome sequencing with cancerapplications.Bioinformatics 31, 1469-1471 (2014)。

34.Lau, B. 等人，LongISLND: in silico sequencing of lengthy and noisydatatypes.Bioinformatics 32, 3829-3832 (2016)。

35.Sedlazeck, F. J. 等人，Accurate detection of complex structuralvariations using single-molecule sequencing.Nat. Methods 15, 461-468 (2018)。

36.Sherry, S. T. 等人，dbSNP: the NCBI database of geneticvariation.Nucleic Acids Res. 29, 308–311 (2001)。

体细胞突变的准确检测具有挑战性，但对于了解癌症的形成、进展和治疗至关重要。我们最近提出了 NeuSomatic，这是第一个基于深度卷积神经网络的体细胞突变检测方法，并展示了在计算机模拟数据上的性能优势。在这项研究中，我们使用了 SEQC-II 联合体的首批全面且充分表征的体细胞参考样本，以探究在癌症突变检测中利用深度学习框架的最佳实践。利用该联合体为这些参考样本建立的高置信度体细胞突变，我们鉴别了用于在源自表示真实场景的样本的多个数据集上建立鲁棒模型的策略。所提出的策略在用于新鲜和 FFPE DNA 输入、不同的肿瘤/正常纯度和不同的覆盖率（范围从 10× 至 2000×）的多种测序技术诸如 WGS、WES、AmpliSeq 靶标测序中实现了高鲁棒性。NeuSomatic 在一般情况下以及在诸如低覆盖率、低突变频率、DNA 损伤和困难基因组区域等挑战性情况下，都大大优于传统检测方法。

体细胞突变是关键的癌症驱动因素。准确的体细胞突变检测可对癌症患者进行精确的诊断、预后和治疗

最近，已经提出了一种基于深度学习的体细胞突变检测方法，称为 NeuSomatic，该方法使用卷积神经网络 (CNN) 直接从原始数据中学习特征表示法

最近，由 FDA 领导的测序质量控制 II 期 (SEQC-II) 联合体的体细胞突变工作组开发了参考匹配的肿瘤-正常样本：人三阴性乳腺癌细胞系 (HCC1395) 和匹配的 B 淋巴细胞来源的正常细胞系 (HCC1395BL)

作为第一个全面且特征明确的经配对肿瘤-正常参考癌症样本，该数据集以及基于多个位点和技术制备的伴随测序数据为实现两个重要目的提供了独特的资源。首先，使用联合体针对这些参考样本开发的高置信度体细胞突变调用结果集，与传统方案相比，我们在真实癌症样本中对基于深度学习的体细胞突变检测实施了深入分析。其次，我们探索了使用 SEQC-II 数据训练 NeuSomatic 中的 CNN 的各种模型构建策略，并鉴别了基于源自表示真实场景的样本的多个数据集的有效训练方法。我们评估了所提出的基于全基因组测序 (WGS)、全外显子组测序 (WES) 和 AmpliSeq 靶向测序数据集的策略，覆盖率范围为10× 至 2000×。WGS 和 WES 数据源自福尔马林固定石蜡包埋 (FFPE) 和新鲜 DNA，使用三种具有不同输入量的文库制备方案，并且来自多个平台/测序位点。使用肿瘤和正常滴定法，我们针对 5% 至 100% 肿瘤纯度，5% 的污染匹配正常样本以及 10× 至 300× WGS覆盖率，评估了不同方法。所提出的用以训练和实现 NeuSomatic 中的深度学习框架的策略在所有上述真实场景中均实现了高鲁棒性，并且显著优于传统的经配对肿瘤-正常体细胞突变检测方法。

我们对 SEQC-II 参考癌症样本的分析表明，在 NeuSomatic 中实现的深度学习方案有助于克服体细胞突变检测中的主要挑战，这些挑战使用传统技术难以解决。因此，源自我们的研究的深度学习模型和策略可以为研究团体提供可行的最佳实践建议，以进行可靠的癌症突变检测。

为了对体细胞突变检测问题进行全谱分析，我们使用了由 SEQC-II 联合体的体细胞突变工作组开发的首个全面的全基因组特征性参考肿瘤-正常配对乳腺癌细胞系（HCC1395 和 HCC1395BL）

图 27 示出了在由 SEQC-II 联合体分类为四个置信度水平（高、中、低和未分类）的针对 HCC1395 的调用结果超集合中，基准真值 SNV 和 INDEL 体细胞突变的 VAF 分布。高置信度和中置信度的调用结果被归为一组体细胞突变的“真值”。

我们使用各种受训网络模型评估了 NeuSomatic 的基于深度学习的体细胞突变检测方法，并将其与七个广泛使用的传统体细胞突变调用算法进行了比较，该传统算法包括 MuTect2

图 28a 至图 28e 显示了 NeuSomatic 在 SEQC-II 数据集中的 123 个复本上的总体性能。图 28a 显示了当应用于 SEQC-II 数据时，基于 DREAM 挑战阶段 3 数据集训练的 NeuSomatic 模型优于其他技术。图 28b 显示了使用 SEQC-II 参考样本训练不同的 NeuSomatic 模型导致平均 F1 得分的额外 4% 至 5% 的绝对改善。图 28c 显示，使用基于 SEQC-II 数据训练的模型，NeuSomatic 跨 WGS、WES、FFPE、AmpliSeq 和不同的文库制备数据集中不同纯度/覆盖率的各种复本集实现了与其他技术相比的一致优势。在该子图中，对于每个复本，最佳 F1 得分是通过不同方法计算得出的。热图示出了体细胞突变检测方法中的任一种方法的 F1 得分与最佳 F1 得分之间的绝对差异。在每张图中，显示了每种方法跨 123 个复本的平均 F1 得分。图 28d 比较了不同技术跨 SNV 和 INDEL 的六个数据集中 123 个复本的性能。对于每个复本，最佳 F1 得分是通过不同方法计算得出的。热图示出了体细胞突变检测方法中的任一种方法的 F1 得分与最佳 F1 得分之间的绝对差异。显示了每种方法跨 123 个复本的平均 F1 得分。

在我们的分析中，我们使用了几种不同的训练模型。首先，我们使用了最近发布的已经可用的模型

为了鉴别使用真实癌症样本建立网络模型的最有效策略，我们使用 SEQC-II 参考样本评估了 NeuSomatic 的另外十种训练方法（图 28e）。

第一模型（SEQC-WGS-掺入）基于一组计算机模拟肿瘤-正常复本的 WGS 对进行训练，其中通过将突变掺入到不同的正常复本品来设计计算机模拟肿瘤。第二模型 (SEQC-WGS-GT-50) 使用 HCC1395 中的基准真值体细胞突变基于 50% 基因组上的一组真实 WGS肿瘤-正常复本对进行训练。第三模型（SEQC-WGS-GT50-掺入WGS10）通过将来自第一模型的训练数据的 10% 添加到用于第二模型的训练数据中来制备，以利用大量的掺入突变和真实的体细胞突变。这三个模型都基于所有数据集上进行了测试。对于诸如 FFPE 和 WES 之类的特定数据集，我们还使用一组合成的计算机模拟肿瘤和正常复本的合成对准备了另外六个专门模型。对于所有模型，我们评估了在基因组 50% 保留区域（未用于 SEQC-WGS-GT-50 模型）上的性能。我们还使用类似于 SEQC-WGS-GT-50 模型的所有基准真值突变但基于全基因组训练了一个模型 (SEQC-WGS-GT-ALL)。SEQC-WGS-GT-ALL 并非直接适用于已经看到所有真实突变的 WGS 数据集，它对于在 HCC1395 以外的其他数据集或样本上的性能分析将是有益的。

使用基于 SEQC-II 样本训练的模型，我们通过将平均 F1 得分额外改善约 4%至 5% 而提高了 DREAM3 模型的平均性能（图 28b）。所提出的模型构建策略跨各种样本类型和测序策略中始终是最佳的，就 SNV 和 INDEL 的平均 F1 得分而言，其性能分别比传统方法高出 5.7% 和 7.8% 以上。类似地，在所有样本中，我们都观察到了比其他传统技术改善了 5.6% 以上的中位数 F1 得分（图 28c 和图 28d）。

我们在使用 HiSeqX10、HiSeq4000 和 NovaSeq S6000 平台在六个测序中心测序的 21 个 WGS 复本上评估了上述体细胞调用技术和网络模型的性能（图 29a 至图 29d）。

NeuSomatic SEQC-WGS-GT50-掺入WGS10 模型的性能始终优于其他方案，但在复本之间存在细微差异，表明了鲁棒性和再现性（图 29a）。NeuSomatic 的 SNV 和 INDEL 的平均 F1 得分分别为 94.6% 和 87.9%，分别比其他传统体细胞突变检测方案的 SNV 和INDEL 的平均 F1 得分高出 5.6% 和 10.2%。精确率-召回率分析表明，与其他技术相比，NeuSomatic 的高精确率推动了这种优势（图 29b）。比较不同的模型训练策略还发现，NeuSomatic-S INDEL 调用从使用基准真值体细胞突变的训练中受益更多，而一般而言，与DREAM3 模型相比，我们使用 SEQC-II 参考样本观察到了高达 11% 的显著改善（图 29c）。

由于肿瘤纯度和匹配的正常样本中的肿瘤细胞的污染会极大地影响突变检测的准确性，因此我们通过 WGS 样本的肿瘤-正常滴定探究了在不同测序深度和样本纯度下的模型鲁棒性。我们首先针对 10× 至 300× 覆盖率研究了与纯正常样本配对的纯度为 5%至 100% 的肿瘤样本（图 30a 至图 30h）。通常，尽管肿瘤纯度和覆盖率有所变化，NeuSomatic 仍保持优于其他方案的优势，这反映了其鲁棒性（图 30a、图 30e 和图 30f）。在更具挑战性的情况比如较低的覆盖率下，它比传统方案具有最大优势（例如，对于具有10× 覆盖率和 100% 纯度的样本，F1 得分为约 20%）。

我们进一步在 80x 覆盖率下，针对 10% 至 100% 范围内的肿瘤样本纯度，分析了正常样本对于 5% 肿瘤污染的鲁棒性（图 30b）。NeuSomatic 对肿瘤-正常的交叉污染具有高鲁棒性，F1 得分的中位数绝对变化少于 5%。在纯正常样本中具有高 F1 得分的其他技术中，Strelka2 也显示出对肿瘤污染的高鲁棒性（F1 得分的中位数变化为 8.4%）。MuTect2、MuSE、Lancet 和 Tnscope 尽管在纯正常样本场景中具有高 F1 得分，但在使用受污染的正常样本时，F1 得分却下降了约 50%。

基于基准真值训练的模型比 DREAM3 模型具有最高的优势，通常主要是由于更高的精确率所致（图 30c 和图 30g）。与基于基准真值训练的模型相比，基于计算机模拟肿瘤训练的 SEQC-WGS-掺入模型也显示出较低的精确率。通常，INDEL、低纯度和低覆盖率样本从基于 SEQC-II 数据训练中受益最大。

分析跨各种覆盖率和纯度设置的不同 INDEL 尺寸的 F1 得分，表明 NeuSomatic对于插入和缺失两者的尺寸变化的鲁棒性（图 30d）。

为了衡量文库制备对预测鲁棒性的影响，我们使用我们的模型来测试使用TruSeq-Nano 和 Nextera Flex 方案制备的六个复本以及三个 DNA 输入量：1ng、10ng 和100ng（图 31a 至图 31f 和图 31h）。在不同的文库制备方法中，NeuSomatic 始终优于其他技术。对于 1 ng TruSeq-Nano 文库，由于去除多余读取后的有效覆盖率受限（约 7x），所有方法的性能均较差。平均而言，与传统的体细胞突变检测技术相比，NeuSomatic 将SNV 的 F1 得分改善了 8.4%。对于 INDEL，Lancet 的基于程序集的 INDEL 调用者比NeuSomatic SEQC-WGS-GT50-掺入WGS10 模型的性能高出约 4%（图 31a）。相反，NeuSomatic 的 SEQC-WGS-GT50 模型在 INDEL 方面的性能与 Lancet 相似（图 31b 和图31g）。与 DREAM3 模型相比，基于 SEQC-II 掺入或基准真值数据的训练导致 SNV 的 F1得分总体改善了约 8.4%。在更具挑战性的情况下（包括输入为 1ng 的 TruSeq-Nano 文库），优势更加明显。我们还观察到，NeuSomatic-S 从这些模型中受益更多。

我们使用在六个测序位点进行测序的 16 个 WES 复本以及 AmpliSeq 数据集中的三个复本，对我们的模型进行了测试（图 31c 至图 31f、图 31i 和图 31j）。尽管基于SEQC-II WGS 样本训练的模型具有不同的覆盖率和平台偏差，但 NeuSomatic 在 2000x覆盖率下在 WES 和 AmpliSeq 数据集上均表现良好。对于 WES，NeuSomatic 的 SNV F1平均得分达到了 95.4%，比替代方案的 F1 平均得分改善了 2.6%。基于 WES 和 WGS 数据训练的模型在 WES 数据集上的表现类似，F1 得分为约 95%。在 AmpliSeq 数据集上，NeuSomatic SEQC-WG-GT50 和 SEQC-WES-掺入模型的平均 F1 得分大于 90%，与其他模型/方案以及 Strelka2 相比是最高的。

为了衡量 NeuSomatic 对经 FFPE 处理样本的预测的鲁棒性，我们使用了 8 个WGS 和 7 个 WES FFPE 复本，分别用 1 小时、2 小时、6 小时和 24 小时的四个不同的甲醛固定时间进行了制备。我们用 FFPE 正常匹配样本和新鲜的正常匹配样本评估了每个FFPE 复本。尽管存在 FFPE 伪像，NeuSomatic 仍然始终具有比其他技术优越的性能，并且其在不同固定时间和所使用的匹配正常样本上的表现在很大程度上保持不变（图 32a、图32e 和图 32f）。基于 WGS FFPE 数据，NeuSomatic 得到的 SNV 和 INDEL 的平均 F1 得分分别为 86.1% 和 76.9%，比替代技术的平均 F1 得分改善了 4% 和 6%。类似地，对于FFPE WES 数据，NeuSomatic 的平均 F1 得分为 78.9%，比传统方案的平均 F1 得分高出4% 以上（图 32c）。

通常，当我们利用 SEQC-II 样本进行训练时，我们观察到了相对于 DREAM3 模型的显著提升（图 32b 和图 32d）。仅对于使用 NeuSomatic 进行 INDEL 预测而言，基于FFPE 样本训练的模型似乎比基于新鲜样本训练的模型有所改善，但对于 SNV，基于新鲜样本训练的模型更为优越。

尽管基于 SEQC-II 训练的通用模型已表现出比其他传统体细胞突变检测方案始终更好的性能，但在这里我们使用跨不同 SEQC-II 数据集的九个复本，探讨了经样本特异性训练的模型是否可以给我们带来额外的准确性提升。对于每个样本，我们使用为该样本制备的计算机模拟肿瘤-正常复本训练了 NeuSomatic 和 NeuSomatic-S 模型。我们还通过将用于 SEQC-WGS-掺入模型的 10% 的训练候选者与针对每个样本推导出的训练数据组合而训练了一个截然独特的模型。我们将这两个样本特异性模型与通用 SEQC-WGS-掺入模型进行了比较（图 33a 至图 33c）。平均而言，与用于 NeuSomatic 的 SEQC-WGS-掺入模型相比，样本特异性模型的 SNV 和 INDEL F1 得分分别改善了约 0.5% 和约 5%。文库制备和 FFPE 样本从样本特异性训练中受益最大。例如，用 Nextera Flex 方案制备的具有1ng DNA 量的文库制备样本的 SNV 和 INDEL 绝对 F1 得分分别获得了 1.6% 和 19.4%的改善。同样，对于 SNV 和 INDEL，具有匹配的新鲜正常样本的 24 小时持续时间 FFPEWGS 样本的 F1 得分分别改善了 1.8 和 14.8 个百分点。对于 NeuSomatic-S，只有INDEL 似乎从样本特异性训练中受益。

我们进一步评估了跨多个数据集检测不同尺寸的 INDEL 的准确性（图 34a 和图34b）。尽管 NeuSomatic 没有明确地将本地程序集用于 INDEL 检测，但是在大范围的INDEL 尺寸、覆盖率和肿瘤纯度下，它仍然始终优于其他方法，包括基于程序集的技术，例如 Lancet（图 34a 和图 30d）。总体而言，对于 2 个以上碱基对 (bps) 的插入和缺失而言，NeuSomatic 得到的 F1 得分比 Lancet 分别高 24.4% 和 6.5%，而 Lancet 是最好的替代传统技术。比较跨不同 NeuSomatic 训练方法的 INDEL 准确性，DREAM3 模型主要受到插入检测准确性低的困扰，而 SEQC-II 受训模型则可靠地鉴别了不同尺寸的插入和缺失两者（图 34b）。

我们分析了跨体细胞突变 VAF 的不同方法的准确性（图 34c 和图 35 至图37c）。在不同的数据集上，我们观察到 NeuSomatic 在 VAF 小于 20% 的情况下显示出优于其他方案的优势（图 34c 和图 35）。我们还观察到，NeuSomatic 对于 VAF 的变化具有高鲁棒性，对于低至 5% 的 VAF 具有一致的预测。对于具有 5% 至 20% VAF 的突变，SEQC-II 受训模型显然比 DREAM3 模型表现更好（图 36a 和图 36b）。对于所有覆盖率和肿瘤纯度设置，观察到相似的鲁棒性（图 37a 至图 37c）。

为了说明对于基因组情境复杂性的鲁棒性，我们评估了不同体细胞调用方案在困难基因组区域上的性能，该困难基因组区域包括不同长度的串联重复序列 (TR) 和片段复制（图 34d 和图 38 至图 41c）。尽管许多其他方案在困难基因组区域上的性能显著下降（图 38），但 NeuSomatic 仍然具有高鲁棒性，与全基因组分析相比，其与其他方案相比的改善甚至更为显著（图 34d 和图 38）。平均而言，对于这些困难基因组区域，NeuSomatic的 SNV 和 INDEL 两者的 F1 得分比最好的替代传统体细胞调用者高出 15% 以上。用SEQC-II 样本训练 NeuSomatic 可以大大改善在这些困难基因组区域上的准确性，其中最大的优势是使用基于基准真值突变训练的模型的串联重复区域（图 39a 至图 40b）。在不同的覆盖率和肿瘤纯度下，观察到了在困难基因组区域上的类似性能（图 41a 至图 41c）。

我们探索了 NeuSomatic 的深度学习框架在跨各种实验设置的集合而检测体细胞突变中的鲁棒性，这些实验设置的集合是在具有不同覆盖率、样本纯度和文库制备的真实 WFP、WES 和 AmpliSeq 样本中见到的，这些检测在使用多个平台在六个中心测序的FFPE 和新鲜 DNA 上进行。实验证实，NeuSomatic 具有从原始数据中捕获真实的体细胞突变标志以将真实调用结果与测序伪影区分开来的潜力。基于 SEQC-II 参考样本训练的NeuSomatic 模型，都使用掺入突变和 SEQC-II 研究推导出的一组基准真值体细胞突变两者，可以提高性能以实现最高的准确性。该分析重点介绍了可用于各种场景中的构建模型的最佳实践策略。与基线 DREAM3 网络模型相比，从 SEQC-II 参考样本推导出的所提出模型显示可通过学习正确的突变信号而减少 SNV 和 INDEL 两者的假阳性，改善插入调用准确性，并增强对困难基因组区域的体细胞检测。

尽管 NeuSomatic 保持对匹配正常样本中肿瘤污染的鲁棒性，但许多体细胞调用方法如 MuTect2、Lancet、MuSE 和 Tnscope 受到了严重影响，因为它们拒绝了正常样本中具有相对应读取的真实体细胞突变。因此，他们的召回率显著下降。基于 WGS 训练的NeuSomatic 模型在靶向测序数据上同样表现良好。另一方面，TNscope 采用机器学习（随机森林）模型，而 Lancet 使用局部彩色的 de Bruijn 图方法来检测突变。两种工具均设计为用于全基因组的体细胞突变检测，因此不适用于靶向测序（通常覆盖少于 1 MB 的基因组）。

分析跨包括突变类型、INDEL 长度、VAF 和基因组区域在内的不同规范的预测准确性，结果表明，NeuSomatic 在挑战性（诸如覆盖率低、VAF 小于 20%、困难基因组区域如串联重复序列、长度超过 2bp 的 INDEL、具有 DNA 损伤的 FFPE 或污染的匹配正常样本）情况下比其他方案具有最大的改善。作为将测序或样本特异性伪影与真实体细胞突变相混淆的结果，传统方案具有许多假阳性（因此精确率较低），这些假阳性通常见于低复杂性基因组区域中（图 34d 和图 38）。但是，通过基于来自多个平台和多个中心的 WGS 样本进行训练，NeuSomatic 学习了易于出错的基因组情境，从而在不同条件下（包括难以调用的低复杂性基因组情境）不断提高了准确性。类似地，通过不同方法对遗漏的调用结果的分析表明，由例如 Strelka2 之类的传统方案产生的私有假阴性调用结果（这些调用结果被NeuSomatic 正确地预测到）具有低 VAF（图 42）。这表明了其他方案在检测低 VAF 突变方面的劣势。另一方面，通过从原始数据中学习正确的标志，NeuSomatic 可以更准确地区别挑战性的低 VAF 调用结果与伪影。

图 42 显示了基于 WGS 数据集的私有 FN 调用结果的 VAF 分布的小提琴图比较。在每个子图中，我们将一种传统的体细胞突变检测方案与 NeuSomatic 进行了比较。对于 X-vs-Y 子图，我们鉴别了 X

在我们的分析中，我们使用了 SEQC-II 参考匹配的样本、人三阴性乳腺癌细胞系(HCC1395) 和匹配的 B 淋巴细胞来源的正常细胞系 (HCC1395BL)。详细的样本信息可在SEQC-II 参考样本手稿

该研究中使用的所有测序数据均可通过 NCBI 的 SRA 数据库 (SRP162370) 公开获得。对于所有样本，首先使用 Trimmomatic15 修剪 FastQ 文件，然后与 BWA-MEM(v0.7.15)16 比对，之后与 Picard MarkDuplicates (https://broadinstitute.github.io/picard) 比对。

使用不同的训练数据集来推导出有效的 NeuSomatic CNN 模型，既使用了用合成的掺入突变制备的计算机模拟肿瘤复本，也使用了具有已知的高置信度体细胞突变真值集的真实肿瘤复本（图 28e）：

作为基线 WGS 模型，我们采用了 DREAM3 模型，该模型是最近通过基于 ICGC-TCGA DREAM 挑战阶段 3 数据训练

作为基线 WES 模型，我们使用了最近通过基于 12 个 TCGA 样本集

为了制备用于构建该模型的训练数据，我们采用 BAMSurgeon 通过将计算机模拟突变掺入到一个 HCC1395BL 复本中并将其与独特的 HCC1395BL 复本配对作为正常匹配来构建合成肿瘤。使用这种方法，我们制备了十个计算机模拟肿瘤-正常样本对。十个复本对中的八个是从 WGS 复本中选择的，这些复本在四个不同位点进行了测序，平均覆盖率为40× 至 95×。通过合并来自 Illumina 的多个 NovaSeq 测序复本而获得其他两个复本对，从而获得约 220× 肿瘤样本覆盖率和约 170× 正常匹配样本覆盖率的计算机模拟。在每一个计算机模拟肿瘤中，我们掺入了约 92 K 个 SNV 和约 22 K 个 INDEL。从 β 分布（形状参数 α= 2 且 β= 5）中随机选择掺入突变的 AF。对于十个重复对中的每一个，我们还通过混合 95% 正常读取和 5% 肿瘤读取来构建不纯正常样本。因此，我们总共使用了20 个计算机模拟肿瘤-正常配对来训练 SEQC-WGS-掺入模型。然后，我们基于在这 20 个复本对中鉴别出的约 272M 个候选突变训练 NeuSomatic 和 NeuSomatic-S，这些候选突变包括约 2M 个具有体细胞 SNV/INDEL 标记的候选者和标记为非体细胞的约 270M 个候选者。

使用真实 WGS 肿瘤正常复本以及 SEQC-II 高置信度真值集作为基准真值体细胞突变来构建该模型。我们使用八个 WGS 肿瘤-正常复本作为基础样本来训练该模型。前七个 WGS 复本对来自 HiSeq 和 NovaSeq 平台上的六个不同测序中心，平均覆盖率为 40× 至 95×，最后一个是通过将 Illumina 的 9 个 NovaSeq 测序复本组合而构建的，得到具有约 390× 覆盖率的复本对。对于这八个复本对中的每一个，我们构建了另外两个纯度变异，一个通过将 95% 正常读取和 5% 肿瘤读取混合而获得 95% 正常纯度，另一个通过将 10% 肿瘤读取和 90% 正常读取混合而获得 10% 肿瘤纯度。因此，对于每个复本对，我们有一个具有 100% 肿瘤和正常的版本、一个 100% 纯肿瘤与 95% 纯正常匹配的版本和一个 10% 纯肿瘤与 100% 纯正常匹配的版本。因此，我们总共使用了 24 个肿瘤-正常复本来训练 SEQC-WGS-GT-50 模型。为了进行公正的评估，我们仅使用了 50% 的基因组来训练该模型，而其余 50% 则用于评估。为了制备训练和评估区域，我们将基因组分为尺寸约 92K bps 的小区域，并随机选择一半的断裂区域进行训练，另一半进行评估。我们将每个灰色区突变周围的 5 个碱基的填充区域排除在训练区域之外，所排除区域包括在SEQC-II 调用结果的超集合中具有低置信度标签的调用结果，以及在 SEQC-II 调用结果的超集合中那些具有超过 30% 的 VAF 的具有未分类标签的调用结果。然后，我们基于在这 24 个复本对中鉴别出的约 137M 个候选突变训练 NeuSomatic 和 NeuSomatic-S，这些候选突变包括约 416K 个具有体细胞 SNV/INDEL 标记的候选者和标记为非体细胞的约136.5M 个候选者。

该模型的制备类似于 SEQC-WGS-GT-50 模型，但是使用了全基因组上的所有高置信度基准真值体细胞突变。因此，它不能用来评估用于训练目的的 WGS 样本，并且鉴于所有样本都已用于训练，因此它可能会偏向于真值集突变。该模型对于其他数据集或非SEQC-II 样本的性能分析仍然是有益的。我们基于在用于 SEQC-WGS-GT-50 模型的这 24个复本对中鉴别出的约 272M 个候选突变来训练用于 NeuSomatic 和 NeuSomatic-S 的模型，这些候选突变包括约 839K 个具有体细胞 SNV/INDEL 标记的候选者和标记为非体细胞的约 271M 个候选者。

该模型的训练数据是通过将 SEQC-WGS-掺入模型中使用的 10% 的训练候选者与SEQC-WGS-GT-50 中的所有候选者合并而成的。这种组合利用了大量的掺入突变和真实肿瘤-正常数据中所见的真实体细胞突变特征。我们基于 164M 个候选突变的组合集训练了NeuSomatic 和 NeuSomatic-S，这些候选突变包括具有体细胞 SNV/INDEL 标签的约 574K个候选者和标记为非体细胞的约 163.5M 个候选者。

与 SEQC-WGS-掺入模型类似，我们使用 BAMSurgeon 构建了 7 个计算机模拟肿瘤-正常 WES 复本，以训练该模型。复本对是从 WES 数据集中选择的，这些复本在四个不同位点进行了测序，平均覆盖率为 60× 至 550×。在每一个计算机模拟肿瘤中，我们掺入了约 97 K 个 SNV 和约 19 K 个 INDEL。从 β 分布（形状参数 α= 2 且 β= 5）中随机选择掺入突变的 AF。然后，我们基于在这 7 个复本对中鉴别出的约 3.7M 个候选突变训练NeuSomatic 和 NeuSomatic-S，这些候选突变包括约 755K 个具有体细胞 SNV/INDEL 标记的候选者和标记为非体细胞的约 3M 个候选者。

与 SEQC-WGS-掺入模型类似，我们使用 BAMSurgeon 构建了 8 个计算机模拟肿瘤-正常 WGS FFPE 复本，以训练该模型。复本对是从 FFPE 数据集中选择的，这些复本在两个不同位点进行了测序，具有四种不同制备时间。在每一个计算机模拟肿瘤中，我们掺入了约 92 K 个 SNV 和约 22 K 个 INDEL。从 β 分布（形状参数 α= 2 且 β= 5）中随机选择掺入突变的 AF。在我们的训练中，我们还将每个计算机模拟肿瘤与一个新鲜的 WGS 复本进行匹配，以包括 FFPE 肿瘤与新鲜的正常场景。因此，我们总共使用了 16 个计算机模拟肿瘤-正常配对来训练 SEQC-FFPE-掺入模型。我们基于在这 7 个复本对中鉴别出的约191M 个候选突变训练 NeuSomatic 和 NeuSomatic-S，这些候选突变包括约 1.7M 个具有体细胞 SNV/INDEL 标记的候选者和标记为非体细胞的约 190M 个候选者。

与其他掺入模型类似，我们使用 BAMSurgeon 构建了 7 个计算机模拟肿瘤-正常WES FFPE 复本，以训练该模型。复本对是从在两个不同位点测序的 WES FFPE 数据集中选择的，并在四个不同的时间间隔内制备。由于我们没有在该数据集中包括两个具有相同FFPE 制备时间和测序位点的正常复本，因此我们将每个复本对的肿瘤比对和正常比对混合在一起，并将混合物分成两个相等的比对。然后，我们将一个比对作为正常复本处理，而将突变掺入到另一个比对中以构建肿瘤复本。在每一个计算机模拟肿瘤中，我们掺入了约90 K 个 SNV 和约 19 K 个 INDEL。从 β 分布（形状参数 α= 2 且 β= 5）中随机选择掺入突变的 AF。在我们的训练中，我们还将每个计算机模拟肿瘤与一个新鲜的 WES 复本进行匹配，以包括 FFPE 肿瘤与新鲜的正常场景。因此，我们总共使用了 14 个计算机模拟肿瘤-正常配对来训练 SEQC-FFPE-WES-掺入模型。我们基于在这 7 个复本对中鉴别出的约9.6M 个候选突变训练 NeuSomatic 和 NeuSomatic-S，这些候选突变包括约 1.4M 个具有体细胞 SNV/INDEL 标记的候选者和标记为非体细胞的约 8.2M 个候选者。

这些模型的训练数据是通过将 10% 的用于 SEQC-WGS-掺入模型的训练候选者与分别来自 SEQC-WES-掺入、SEQC-FFPE-掺入和 SEQC-FFPE-WES-掺入的所有候选者组合而制备的。这种组合既利用了大量的掺入 WGS 突变，又利用了 WES 和 FFPE 样本的样本偏差。

除了上述通用模型外，对于跨多种数据类型的一组九个样本，我们还推导出了样本特异性模型。选择的样本包括通过 NovaSeq 仪器测序的 WGS 样本、WES 样本、使用Nextera Flex 文库制备试剂盒制备的 DNA 量为 1ng 的样本、肿瘤纯度为 50% 的 30×WGS 样本、肿瘤纯度为 10% 的 100× WGS 样本、以及各自用福尔马林处理 24 小时并与新鲜/FFPE 正常样本匹配的两个 WGS 和 两个 WES FFPE 样本。对于每个样本，我们使用随机掺入制备了计算机模拟肿瘤。对于 10% 肿瘤样本，从 β 分布（形状参数 α= 2 且 β=18）中随机选择掺入突变的 AF。对于其他样本，使用 β 分布（形状参数 α= 2 且 β= 5）选择 AF。然后，对于每个样本，我们使用计算机模拟肿瘤-正常复本训练了 NeuSomatic 和NeuSomatic-S 模型。此外，对于每个样本，我们还通过将用于 SEQC-WGS-掺入模型的 10%的训练候选者与针对该样本推导出的训练数据组合而训练了一个独特的模型。

除了 NeuSomatic 以外，我们使用了七种体细胞突变调用者：MuTect2(4.beta.6)

我们以集成和独立两种模式使用了 NeuSomatic (0.1.4)。对于集成模式，除了通过 NeuSomatic 的“比对扫描”步骤鉴别的候选变体之外，我们还合并了来自五个单独的调用者（MuTect2、SomaticSniper、Strelka2、MuSE 和 VarDict）的调用结果，并将它们的出现表示为用于每个候选变体的额外输入通道。对于 NeuSomatic 和 NeuSomatic-S，我们在预处理步骤中使用了“--scan_maf 0.01--min_mapq 10 --snp_min_af 0.03 --snp_min_bq15 --snp_min_ao 3 --ins_min_af 0.02 --del_min_af 0.02 --num_threads 28 --scan_window_size 500 --max_dp 100000”。对于训练，我们使用了“--coverage_thr 300--batch_size 8000”。

我们使用了“瓶中基因组”联合体推导出的一组困难基因组区域

为了更好地理解不同技术的性能差异，我们对 NeuSomatic 的假阴性 SNV 的VAF 分布与其他方案基于 WGS 数据集进行了成对比较（图 42）。在每个 X-vs-Y 分析中，我们鉴别了 XFN，这是在 21 个 WGS 复本的至少 11 个中被算法 X（X 中的假阴性）遗漏的基准真值 SNV 集。同样，我们鉴别了 YFN，这是在 21 个 WGS 复本的至少 11 个中被算法 Y（Y 中的假阴性）遗漏的基准真值 SNV 集。然后，我们计算 X 的私有假阴性，定义为XFN/YFN，以及 Y 的私有假阴性，定义为 YFN/XFN。然后，对于那些集合中每个集合中的体细胞突变，我们计算突变的 VAF 的分布。

为了公平地比较，我们在 50% 保留基因组区域（未用于训练 SEQC-WGS-GT-50 模型的区域）上评估了所有模型和体细胞突变算法。这个约 1.4GB 的区域包含来自用于HCC1395 的 SEQC-II 真值集的约 21K SNV 和约 1.3K INDEL，除用于 SEQC-WGS-GT-ALL模型外，这些 NUV 未用于训练任何 NeuSomatic 模型。

SEQC-II 联合体将针对 HCC1395 标记为高置信度和低置信度的调用结果归为一组，作为本文使用的体细胞突变的“真值集”。我们采用此真值集来计算所有渠道在HCC1395 的所有复本中的真阳性和假阴性。根据 SEQC-II 联合体的建议，我们还将低置信度调用结果列入了黑名单，因为它们的验证率较低。由于该联合体提供的真值集的 VAF 检测极限为 5% 并且深度检测极限为 50×，因此对于覆盖率数据较高的调用结果或 VAF 较低的调用结果，我们无法确定其真实的体细胞状态。因此，对于不属于真值集的任何渠道报告的私有调用结果，我们排除了在评估中被认为模棱两可的调用结果。总而言之，对于手头的平均覆盖率为 C 且肿瘤纯度为 P 的肿瘤复本，如果渠道报告该复本有一个为 d 个读取所支持的私有体细胞突变（该体细胞突变不在真值集中），则仅当其在 100% 纯度和 50x覆盖率下的支持读取推断数量 ≥3 时，我们才将此调用结果标记为假阳性。换句话说，如果 d ≥ 3CP/50，则此调用结果为假阳性；否则，我们将从评估中排除此调用结果。

对于 WES 和 AmpliSeq 数据，由于评估区域中真正 indel 的数量非常有限，因此我们仅报告了 SNV 评估。

1.Xu, C.，A review of somatic single nucleotide variant callingalgorithms for next-generation sequencing data.Comput.Struct. Biotechnol. J.16, 15–24 (2018)。

2.Cibulskis, K. 等人，Sensitive detection of somatic point mutationsin impure and heterogeneous cancer samples.Nat. Biotechnol. 31, 213 (2013)。

3.Fan, Y. 等人，MuSE: accounting for tumor heterogeneity using asample-specific error model improves sensitivity and specificity in mutationcalling from sequencing data.Genome Biol. 17, 178 (2016)。

4.Larson, D. E. 等人，SomaticSniper: identification of somatic pointmutations in whole genome sequencing data.Bioinformatics 28, 311–317 (2011)。

5.Kim, S. 等人，Strelka2: fast and accurate calling of germline andsomatic variants.Nat. Methods 15, 591–594 (2018)。

6.Lai, Z. 等人，VarDict: a novel and versatile variant caller fornext-generation sequencing in cancer research.Nucleic Acids Res. 44, e108–e108 (2016)。

7.Freed, D.、Pan, R. 和 Aldana, R.TNscope: Accurate Detection ofSomatic Mutations with Haplotype-based Variant Candidate Detection andMachine Learning Filtering. bioRxiv 250647

8.Narzisi, G. 等人，Lancet: genome-wide somatic variant calling usinglocalized colored DeBruijn graphs.Communications biology 1, 20 (2018)。

9.Koboldt, D. C. 等人，VarScan 2: somatic mutation and copy numberalteration discovery in cancer by exome sequencing.Genome Res. 22, 568–576(2012)。

10.Fang, L. T. 等人，An ensemble approach to accurately detect somaticmutations using SomaticSeq.Genome Biol. 16, 197 (2015)。

11.Sahraeian, S. M. E.、Liu, R.、Lau, B.、Podesta, K.、Mohiyuddin, M. 和Lam, H. Y.，Deep convolutional neural networks for accurate somatic mutationdetection.Nature communications.10, 1041 (2019)。

12.The Somatic Mutation Working Group of SEQC-II Consortium.Achieving reproducibility and accuracy in cancer mutation detection withwhole-genome and whole-exome sequencing. https://doi.org/10.1101/626440(2019)。

13.Fang, L. T. 等人，Establishing reference samples for detection ofsomatic mutations and germline variants with NGS technologies, https://doi.org/10.1101/625624 (2019)。

14.Ewing, A. D. 等人，Combining tumor genome simulation withcrowdsourcing to benchmark somatic single-nucleotide-variant detection.Nat.Methods 12, 623 (2015)。

15.Bolger, A. M.、Lohse, M. 和 Usadel, B.，Trimmomatic: a flexibletrimmer for Illumina sequence data. Bioinformatics 30, 2114-2120 (2014)。

16.Li, H.，Aligning sequence reads, clone sequences and assemblycontigs with BWA-MEM. https://arxiv.org/abs/1303.3997 (2013)。

17.Ewing, A. D. 等人，Combining tumor genome simulation withcrowdsourcing to benchmark somatic single-nucleotide-variant detection.Nat.Methods 12, 623 (2015)。

18.Grossman, R. L. 等人，Toward a shared vision for cancer genomicdata.N. Engl. J. Med. 375, 1109–1112 (2016)。

19.Krusche, Peter、Len Trigg、Paul C. Boutros、Christopher E. Mason、M.Francisco、Benjamin L. Moore、Mar Gonzalez-Porta 等人，“Best practices forbenchmarking germline small-variant calls in human genomes.”Naturebiotechnology 1 (2019)。

尽管已经参照一些说明性实施例描述了本公开，但应当理解，本领域技术人员可以在本公开原则的精神和范围内设计出许多其他的修改和实施例。更具体地，在不违背本公开的精神的情况下，在上述公开、附图和所附权利要求的范围内，所述主题组合排列的组成部分和/或布置可以进行合理的变化和修改。除了在所述组成部分和/或布置中的变化和修改外，对于本领域技术人员来说，替代用途也将是显而易见的。

本说明书中提到的和/或应用数据表中列出的所有美国专利、美国专利申请公开、美国专利申请、国外专利、国外专利以及非专利公开均通过引用整体并入本文。如果需要采用各种专利、申请和出版物的概念来提供其它实施例，则可以修改实施例的各方面。

当特征或元件在本文中被称为在另一特征或元件“上”时，它可以直接位于另一特征或元件上，或者也可以存在中间特征和/或元件。相反，当特征或元件被称为“直接位于”另一特征或元件“上”时，则不存在中间的特征或元件。还应理解，当特征或元件被称为“连接”、“附接”或“耦接”至另一特征或元件时，它可以直接连接、附接或耦接至另一特征或元件，或可能存在中间特征或元件。相反，当特征或元件被称为“直接连接”、“直接附接”或“直接耦接”到另一特征或元件时，则不存在中间的特征或元件。尽管关于一个实施例进行了描述或显示，但是如此描述或显示的特征和元件可以应用于其他实施例。本领域的技术人员还将意识到，提及与另一特征“相邻”设置的结构或特征可具有与相邻特征重叠或位于其下的部分。

本文使用的术语只是为了描述特定实施例的目的，并非旨在限制本发明。例如，如本文所用，单数形式“一”、“一个”、“该”和“所述”也旨在包括复数形式，除非上下文另外明确指出。将进一步理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指定了所述特征、步骤、操作、元件和/或组件的存在，但并不排除存在或添加一个或多个其他特征、步骤、操作、元件、组件和/或它们的组。如本文所用，术语“和/或”包括一个或多个相关联列出的项目的任意和所有组合，并且可以缩写为“/”。

本文可使用空间相对术语诸如“下方”、“在……下”、“下面”、“上方”、“上面”等来描述如图中所示出的一个元件或特征与另一元件或特征的关系。应理解的是，除了图中所描绘的取向之外，空间相对术语还意图涵盖装置在使用或操作中的不同方位。例如，如果附图中的装置被翻转，则被描述为在其他元件或特征“下方”或“在……下方”的元件将定向为在其他元件或特征“上方”。因此，示例性术语“下方”可涵盖上方和下方两种取向。可以其他方式定向装置（转动 90° 或以其他取向），并且相应地解释本文所使用的空间相对描述语。类似地，除非另外具体指出，否则术语“向上”、“向下”、“竖直”、“水平”等在本文中仅用于解释的目的。

尽管本文可以使用术语“第一”和“第二”来描述各种特征/元件（包括步骤），但是除非上下文另外指出，否则这些特征/元件不应受这些术语的限制。这些术语可以用于将一个特征/元件与另一特征/元件区分开来。因此，在不脱离本发明教导的情况下，下面讨论的第一特征/元件可以被称为第二特征/元件，并且类似地，下面讨论的第二特征/元件可以被称为第一特征/元件。

在整个本说明书和随后的权利要求中，除非上下文另外要求，否则词语“包含”以及诸如“包括”和“含有”的变体意味着可以在方法和物品中共同使用各种组分（例如，组合物以及包括设备和方法的设施）。例如，术语“包含”将被理解为暗示包括任何陈述的元件或步骤，但是不排除任何其他元件或步骤。

如本文在说明书和权利要求书中所用，包括在实施例中所用，并且除非另有明确说明，否则所有数字可以被读作与以词语“约”或“近似地”开头一样，即使该术语没有明确出现。当描述幅度和/或位置时，可以使用短语“约”或“大约”来指示所描述的值和/或位置在值和/或位置的合理预期范围内。例如，数字值的值可以是规定值（或值的范围）的 +/-0.1%、规定值（或值的范围）的 +/- 1%、规定值（或值的范围）的 +/- 2%、规定值（或值的范围）的 +/- 5%、规定值（或值的范围）的 +/- 10% 等。此处给出的任何数值应除非上下文另外指出，否则也应理解为包括约或大约该值。例如，如果公开值“10”，则还公开“约 10”。本文叙述的任何数值范围旨在包括其中包含的所有子范围。还应理解，如本领域技术人员适当理解的那样，当公开的值“小于或等于”该值时，还公开了“大于或等于该值”以及值之间的可能范围。例如，如果公开了值“X”，则还公开了“小于或等于 X”以及“大于或等于 X”（例如，其中 X 是数值）。还应理解，在整个申请中，以多种不同格式提供数据，并且该数据表示端点和起点以及数据点的任意组合的范围。例如，如果公开了特定数据点“10”和特定数据点“15”，则应理解，大于、大于或等于、小于、小于或等于以及等于 10 和 15 视为已公开，且介于 10 和 15 之间也视为已公开。还应理解，还公开了两个特定单位之间的每个单位。例如，如果公开了 10 和 15，则还公开了 11、12、13 和 14。

尽管上面描述了各种说明性实施例，但是在不脱离如权利要求所描述的本发明的范围的情况下，可以对各种实施例进行多种改变中的任何一种。例如，在替代实施例中，可以经常改变所描述的各种方法步骤的顺序，而在其他替代实施例中，可以完全跳过一个或多个方法步骤。在一些实施例中可以包括各种装置和系统实施例的任选特征，而在其他实施例中可以不包括。因此，前面的描述主要是为了示例性目的而提供，并且不应解释为限制本发明的范围，因为本发明的范围在权利要求中所阐述。

本文包括的实例和示例性说明通过示例性说明而非限制的方式示出了可以实践本主题的特定实施例。如所提及的，可以利用其他实施例并从中得出其他实施例，从而可以在不脱离本公开的范围的情况下进行结构和逻辑上的替换和改变。仅出于方便起见，如果事实上公开了多个发明或发明构思，本文中可以单独地或共同地用术语“发明”来指代本发明主题的这些实施例，而无意将本申请的范围限制为任何单个发明或发明构思。因此，尽管本文中已经示例性说明和描述了特定实施例，但是为实现相同目的而计算的任何布置都可以代替所示的特定实施例。本公开意图覆盖各种实施例的任何和所有修改或变化。通过阅读以上描述，以上实施例的组合以及本文中未具体描述的其他实施例对于本领域技术人员将是显而易见的。