Systems and methods for germline and somatic variant calling using neural networks

Abstract

The present disclosure provides systems and methods for utilizing neural networks, such as convolutional neural networks, to analyze genomic sequence data generated by a sequencer and generate accurate predictive data that identifies and describes germline and / or somatic variants within the sequence data.

Description

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims priority to U.S. Provisional Application No. 62 / 718,338, filed on August 13, 2018, and U.S. Provisional Application No. 62 / 855,541, filed on May 31, 2019, each of which is hereby incorporated by reference in its entirety.

[0003] Incorporated by Reference

[0004] All publications and patent applications mentioned in this specification are herein incorporated by reference to the same extent as if each individual publication or patent application was specifically and individually indicated to be incorporated by reference. Technical Field

[0005] Embodiments of the present invention relate generally to systems and methods for variant calling, and more particularly to variant calling using machine learning algorithms. Background Art

[0006] Identifying genetic mutations using sequencing data is a crucial task in bioinformatics that has implications for the diagnosis, prognosis, and treatment of a wide range of diseases, including cancer. This is a nontrivial task, especially when the sequencing technology used to generate the sequencing data has a high error rate (as is common with single-molecule sequencing technologies) or when mutations occur at a low frequency (especially true for cancer mutations) or within complex genomic regions. The low frequency of cancer mutations in sequenced samples may be due to sample contamination (because tumor samples may contain some DNA from normal cells) or tumor heterogeneity.

[0007] Traditional approaches narrow the scope of the problem by limiting specific techniques and types of mutations to be tested. As a result, there are different sets of tools and frameworks for identifying germline and cancer mutations, which are further stratified according to whether the mutation is large or small. In addition, methods used for low-error sequencing technologies (e.g., short-read sequencing) are different from methods used for high-error sequencing technologies (e.g., single-molecule sequencing). The differences are mainly reflected in the data signatures and the statistical models used to distinguish errors from true mutations. The use of explicit data signatures and statistical models limits the generalizability of traditional methods to multiple technologies. In addition, even traditional methods require adjustment of parameters to capture specific characteristics of the sample, which may depend on multiple factors, including the technology used, sequencing protocol, sequencing sample, sample purity, etc.

[0008] Deep learning is a machine learning technique that can implicitly capture data signatures and statistical models using multiple data layers, given the right training data. This makes deep learning a potentially attractive solution to the problem of mutation detection in terms of its generality. Several deep learning solutions for classifying germline mutations are known.

[0009] Poplin, R. et al. (“Creatinga Universal SNP and Small Indel VariantCaller With Deep Neural Networks”, bioRxiv, 092890 (2016)) proposed an approach where candidate variant locations are first identified by genome scanning, and then stacked images of reads aligned around each identified candidate are provided as input to a convolutional neural network (CNN), which then predicts the genotype of the candidate variant. This approach is not very practical because it requires a complex neural network architecture and a lot of memory.

[0010] Luo, R. et al. (“Clairvoyante:a Multi-task Convolutional Deep Neural Network For Variant Calling In Single Molecule Sequencing. bioRxiv, 310458(2018)) used a data summarization approach that involves scanning the genome to identify candidate variant sites and preparing multiple input matrices for each candidate position, one matrix for each variant type. Similar techniques are discussed in Chin, J. (“Simple Convolutional Neural Network for Genomic Variant Calling with TensorFlow”, https: / / towardsdatascience.com / simple-convolution-neural-network-for-genomic-variant-calling-with-tensorflow-c085dbc2026f, published on July 16, 2017). These methods are also suboptimal because they do not clearly represent all insertion sequences in the entire window. For example, Clairvoyante Insertion events are assigned a separate input matrix that captures only insertions at the variant locations, while ignoring true / false insertions at other locations within the window. Full-map information about the precise location of the inserted sequence throughout the window is critical to distinguish true insertion variants from background noise.

[0011] Torracinta, R. et al. (“Adaptive somatic mutations calls with deeplearning and semi-simulated data. bioRxiv, 79087 (2016)) describe applying a six-layer fully connected neural network to a set of manually extracted features. But a system with only fully connected layers cannot take advantage of the capabilities provided by convolutional neural networks, namely, learning feature representations directly from raw sequence data using patterns found in the local genomic context of candidate variants. Furthermore, because the fully connected network is more complex, the approach is less general and less scalable than what can be achieved using CNNs. Summary of the invention

[0012] The present invention relates generally to systems and methods for variant calling, and more particularly to variant calling using machine learning algorithms.

[0013] In some embodiments, a method for germline variant calling is provided. The method may include obtaining a reference sequence, a plurality of sequence reads, and positions of candidate variants within the sequence reads; obtaining an extended sequence read by inserting one or more spaces into the one or more sequence reads; obtaining an extended reference sequence by inserting one or more spaces into the reference sequence; converting segments of the extended sequence reads around the candidate variants into a sample matrix; converting segments of the extended reference sequence around the candidate variants into a reference matrix; providing the sample matrix and the reference matrix to a trained neural network; and, at the output of the trained neural network, obtaining prediction data about variants within the plurality of sequence reads.

[0014] In some embodiments, the method further includes detecting one or more inserted bases within the plurality of sequence reads, wherein expanding the sequence reads and the reference sequence includes: for each inserted base detected in any of the sequence reads, inserting a space into the reference sample at the position of the inserted base.

[0015] In some embodiments, the method further comprises: for each inserted base detected in any of the sequence reads, inserting a space at the position of the inserted base into each sequence read in which no insertion was detected at the position of the inserted base.

[0016] In some embodiments, the sample matrix includes at least four rows representing four types of nucleotide bases, each row representing the base number of the corresponding nucleotide base type at different positions within the fragments read by the extended sequence; and at least one row representing the number of inserted spaces at different positions within the fragments read by the extended sequence.

[0017] In some embodiments, the reference matrix has the same dimensionality as the sample matrix, and wherein the reference matrix provides a complete representation of the specific locations of different nucleotide bases and spaces within the augmented reference sequence.

[0018] In some embodiments, the trained neural network comprises a trained convolutional neural network.

[0019] In some embodiments, the method further comprises providing at least one of the following to the trained neural network: a variant position matrix representing the positions of candidate variants within the fragments of the expanded sequence reads; a coverage matrix representing the coverage or depth of the fragments of the expanded sequence reads; an alignment feature matrix representing alignment features of the expanded sequence reads; and a knowledgeable base matrix representing information about publicly known information about one or more variants.

[0020] In some embodiments, the predicted data about the variant includes at least one of: a predicted type of the variant; a predicted position of the variant; a predicted length of the variant; and a predicted genotype of the variant.

[0021] In some embodiments, the prediction data about the variant includes a predicted type of the variant, and wherein the neural network is configured to generate one of a plurality of values ​​for the predicted type of the variant, the plurality of values ​​comprising: a first value indicating a probability that the variant is a false positive; a second value indicating a probability that the variant is a single nucleotide polymorphism variant; a third value indicating a probability that the variant is a deletion variant; and a fourth value indicating a probability that the variant is an insertion variant.

[0022] In some embodiments, a method for somatic variant calling is provided. The method may include obtaining a plurality of normal sequence reads and a plurality of tumor sequence reads; converting the fragments of the normal sequence reads and the fragments of the tumor sequence reads into a normal sample matrix and a tumor sample matrix, respectively; feeding the normal sample matrix and the tumor sample matrix into a trained convolutional neural network; and obtaining, at the output of the trained convolutional neural network, a predicted type of somatic variants within the plurality of tumor sequence reads.

[0023] In some embodiments, the plurality of tumor sequence reads represent genetic information of a tumor sample of the patient, and the plurality of normal sequence reads represent genetic information of a normal sample of the patient.

[0024] In some embodiments, converting normal sequence reads to a normal sample matrix includes expanding segments of normal sequence reads by inserting one or more spaces into one or more normal sequence reads; and converting segments of tumor sequence reads to a tumor sample matrix includes expanding segments of tumor sequence reads by inserting one or more spaces into one or more tumor sequence reads.

[0025] In some embodiments, the tumor sample matrix includes at least one row for each nucleotide base type, each row representing the number of occurrences of the corresponding nucleotide base type at each position within the fragments read by the tumor sequence; and at least one row representing the number of inserted spaces at each position within the fragments read by the tumor sequence.

[0026] In some embodiments, the method further comprises providing the trained convolutional neural network with one or more matrices representing features obtained from one or more other variant callers that have analyzed the plurality of tumor sequence reads and / or the plurality of normal sequence reads.

[0027] In some embodiments, the method further comprises obtaining a reference sequence; converting the reference sequence into a reference matrix; and feeding the reference matrix into the trained convolution matrix together with the normal sample matrix and the tumor sample matrix.

[0028] In some embodiments, a non-transitory computer-readable medium is provided, comprising instructions that, when executed by one or more processors of a computer system, cause the computer system to perform operations including: obtaining a plurality of normal sequence reads and a plurality of tumor sequence reads; converting fragments of the normal sequence reads and fragments of the tumor sequence reads into a normal sample matrix and a tumor sample matrix, respectively; feeding the normal sample matrix and the tumor sample matrix into a trained convolutional neural network; and obtaining, at the output of the trained convolutional neural network, predicted types of somatic variants within the plurality of normal sequence reads.

[0029] In some embodiments, a computer system is provided, comprising one or more processors and coupled to one or more non-transitory computer-readable memories storing instructions, which, when executed by the computer system, cause the computer system to perform operations including: obtaining a plurality of tumor sequence reads; obtaining augmented tumor sequence reads by inserting one or more spaces into the one or more tumor sequence reads; converting segments of the tumor sequence reads into a tumor sample matrix; feeding the normal sample matrix and the tumor sample matrix into a trained neural network; and obtaining, at the output of the trained neural network, predicted types of somatic variants within the plurality of tumor sequence reads.

[0030] In some embodiments, a method for variant calling is provided. The method may include: obtaining a reference sequence and a plurality of sequence reads; optionally performing a first alignment of the plurality of sequence reads with the reference sequence, unless the plurality of sequence reads and the reference sequence obtained are obtained in an already aligned configuration; identifying candidate variant positions from the aligned sequence reads and the reference sequence; expanding sequence reads and / or reference sequences around the candidate variant positions to achieve a second alignment of the plurality of sequence reads with the reference sequence; generating a reference matrix for the candidate variant positions from the expanded reference sequence, and generating a sample matrix for the candidate variant positions from the plurality of expanded sequence reads; inputting the reference matrix and the sample matrix into a neural network; and determining with the neural network whether a variant type exists at the candidate variant position.

[0031] In some embodiments, the step of augmenting the sequence reads and / or reference sequences comprises introducing one or more spaces into the sequence reads and / or reference sequences to account for insertions and / or deletions in the sequence reads.

[0032] In some embodiments, the method further includes generating a plurality of training matrices from a training data set, wherein the training matrix has a structure corresponding to the sample matrix and the reference matrix, wherein the training data set comprises sequence data, the sequence data comprises a plurality of mutations, the plurality of mutations comprises single nucleotide variants, insertions and deletions; and training the neural network using the plurality of training matrices.

[0033] In some embodiments, the training data set comprises a plurality of subsets, wherein each subset comprises tumor purity levels ranging from 0% to 100%, wherein at least two of the subsets each have a different tumor purity level.

[0034] In some embodiments, at least three of the subsets each have a different level of tumor purity.

[0035] In some embodiments, the plurality of subsets include a first subset having a tumor purity level of less than about 30%, a second subset having a tumor purity level between about 30% and 70%, and a third subset having a third tumor purity level of at least about 70%.

[0036] In some embodiments, the plurality of subsets include a first subset having a tumor purity level of less than about 40%, a second subset having a tumor purity level between about 40% and 60%, and a third subset having a tumor purity level of at least about 60%.

[0037] In some embodiments, the plurality of subsets includes a subset having a tumor purity level of less than about 10%.

[0038] In some embodiments, the plurality of subsets includes a subset having a tumor purity level of less than about 5%.

[0039] In some embodiments, the training data set includes synthetic data.

[0040] In some embodiments, the synthetic data comprises artificially generated mutations, wherein the artificially generated mutations comprise single nucleotide variants, insertions, and deletions.

[0041] In some embodiments, the training data set comprises real data, wherein the real data comprises real mutations, wherein the real mutations comprise single nucleotide variants, insertions, and deletions.

[0042] In some embodiments, the training data set comprises a plurality of subsets, wherein each subset comprises variant allele frequencies ranging from 0% to 100%, wherein at least two of the subsets each have a different variant allele frequency.

[0043] In some embodiments, at least three of the subsets each have a different variant allele frequency level.

[0044] In some embodiments, at least one of the subsets has a variant allele frequency of at least 2.5%.

[0045] In some embodiments, at least one of the subsets has a variant allele frequency of at least 5%.

[0046] In some embodiments, at least one of the subsets has a variant allele frequency of at least 10%.

[0047] In some embodiments, the method further comprises inputting at least one prediction from at least one mutation calling algorithm into a neural network.

[0048] In some embodiments, the at least one prediction comprises at least three predictions from at least three independent mutation calling algorithms.

[0049] In some embodiments, the at least one prediction comprises at least five predictions from at least five independent mutation calling algorithms.

[0050] In some embodiments, the training data set includes a collection of synthetic data and real data.

[0051] In some embodiments, the training data set includes at least 5% synthetic data.

[0052] In some embodiments, the training data set includes at least 10% synthetic data.

[0053] In some embodiments, the training dataset comprises whole genome sequencing data.

[0054] In some embodiments, the training dataset comprises whole exome sequencing data.

[0055] In some embodiments, the training dataset comprises targeted sequencing data.

[0056] In some embodiments, the training data set includes data obtained from formalin-fixed paraffin-embedded samples.

[0057] In some embodiments, the training dataset includes at least two of whole genome sequencing data, whole exome sequencing data, targeted sequencing data, and data obtained from formalin-fixed paraffin-embedded samples.

[0058] In some embodiments, the training dataset includes at least three of whole genome sequencing data, whole exome sequencing data, targeted sequencing data, and data obtained from formalin-fixed paraffin-embedded samples.

[0059] In some embodiments, the training dataset includes whole genome sequencing data, whole exome sequencing data, targeted sequencing data, and data obtained from formalin-fixed paraffin-embedded samples.

[0060] In some embodiments, a method for variant calling is provided. The method may include obtaining a reference sequence, a plurality of tumor sequence reads, and a plurality of normal sequence reads; optionally performing a first alignment of the plurality of tumor sequence reads and the plurality of normal sequence reads with the reference sequence, unless the plurality of tumor sequence reads and the plurality of normal sequence reads and the reference sequence obtained are obtained in an already aligned configuration; identifying a candidate variant position from the aligned tumor sequence reads, the normal sequence reads, and the reference sequence; augmenting the tumor sequence reads and / or the normal sequence reads and / or the reference sequence around the candidate variant position to achieve a second alignment of the plurality of tumor sequence reads and the plurality of normal sequence reads with the reference sequence; generating a reference matrix for the candidate variant position from the augmented reference sequence, generating a tumor matrix for the candidate variant position from the plurality of augmented tumor sequence reads, and generating a normal matrix for the candidate variant position from the plurality of augmented normal sequence reads; inputting the reference matrix, the tumor matrix, and the normal matrix into a neural network; and determining with the neural network whether a variant type exists at the candidate variant position.

[0061] In some embodiments, the method further includes generating a plurality of training matrices from a training data set, wherein the training matrices have structures corresponding to a tumor matrix, a normal matrix, and a reference matrix, wherein the training data set includes tumor sequence data and normal sequence data; and training the neural network using the plurality of training matrices.

[0062] In some embodiments, both the tumor sequence data and the normal sequence data include multiple mutations including single nucleotide variants, insertions, and deletions.

[0063] In some embodiments, the normal sequence data includes up to 5% of the tumor sequence data.

[0064] In some embodiments, the normal sequence data includes up to 10% of the tumor sequence data.

[0065] In some embodiments, the tumor sequence data comprises a tumor purity level between about 10% and 100%.

[0066] In some embodiments, the training data set comprises a plurality of tumor sequence data subsets, wherein each tumor sequence data subset comprises a tumor purity level ranging from 10% to 100%, wherein at least two of the tumor sequence data subsets each have a different tumor purity level.

[0067] In some embodiments, at least three of the tumor sequence data subsets each have a different level of tumor purity.

[0068] In some embodiments, the plurality of tumor sequence data subsets include a first tumor sequence data subset having a tumor purity level of less than about 30%, a second tumor sequence data subset having a tumor purity level between about 30% and 70%, and a third tumor sequence data subset having a tumor purity level of at least about 70%.

[0069] In some embodiments, the plurality of tumor sequence data subsets include a first tumor sequence data subset having a tumor purity level of less than about 40%, a second tumor sequence data subset having a tumor purity level between about 40% and 60%, and a third tumor sequence data subset having a tumor purity level of at least about 60%.

[0070] In some embodiments, the training data set includes synthetic data.

[0071] In some embodiments, the synthetic data comprises artificially generated mutations, wherein the artificially generated mutations comprise single nucleotide variants, insertions, and deletions.

[0072] In some embodiments, the training data set comprises real data, wherein the real data comprises real mutations, wherein the real mutations comprise single nucleotide variants, insertions, and deletions.

[0073] In some embodiments, the training dataset comprises whole genome sequencing data.

[0074] In some embodiments, the training dataset comprises whole exome sequencing data.

[0075] In some embodiments, the training dataset comprises targeted sequencing data.

[0076] In some embodiments, the training data set includes data obtained from formalin-fixed paraffin-embedded samples.

[0077] In some embodiments, a system is provided. The system may include a processor configured to perform the steps described in any one of claims 18 to 64.

[0078] BRIEF DESCRIPTION OF THE DRAWINGS

[0079] The novel features of the present invention are particularly set forth in the appended claims. The features and advantages of the present invention may be better understood by reference to the following detailed description in conjunction with the accompanying drawings, which set forth illustrative embodiments in which the principles of the present invention are utilized, wherein:

[0080] FIG1 illustrates an exemplary system including a sequencing apparatus and a computer system according to some embodiments;

[0081] Figures 2a and 2b illustrate exemplary methods for germline variant calling according to some embodiments;

[0082] FIG3 illustrates an expanded step of the exemplary method of FIG2 according to some embodiments;

[0083] FIG4 illustrates exemplary data that may be provided to and obtained from a neural network according to some embodiments;

[0084] FIG5 illustrates an exemplary architecture of a neural network according to some embodiments;

[0085] FIG6 illustrates an exemplary method for somatic variant calling according to some embodiments;

[0086] FIG. 7 illustrates an expansion step that may be implemented as part of the exemplary method of FIG. 6 according to some embodiments;

[0087] FIG8 illustrates another example of data that may be provided to and obtained from a neural network in accordance with some embodiments;

[0088] FIG9 illustrates an exemplary architecture of a neural network according to some embodiments;

[0089] Figure 10a shows the generation of the input matrix for a given candidate somatic SNV; Figure 10b shows the input matrix and network architecture; Figure 10c shows the generation of the input matrix for a given candidate somatic deletion; Fig.10dThe generation of an input matrix for a given candidate somatic insertion is shown;

[0090] Figures 11a to 11c show the performance of different detection methods on the platinum two-sample mixture dataset;

[0091] Figures 12a to 12c show the performance of different detection methods on the Dream phase 3 dataset;

[0092] Figures 13a to 13c show the performance of different detection methods on the Dream stage 4 dataset;

[0093] Figures 14a to 14c show the performance of two different detection methods on the PacBio dataset;

[0094] Figures 15a-15b show the effect of sequence coverage on the performance of different detection methods on a whole exome sample mixture dataset;

[0095] Figures 16a to 16d show the performance of different detection methods on the platinum tumor incorporation dataset;

[0096] Figures 17a and 17c show the performance of different somatic mutation detection methods on whole exome datasets; Figures 17b and 17d show the performance of somatic mutation detection methods on targeted group datasets;

[0097] FIG18a shows the performance analysis of testing an exome mixture dataset using a model trained on the whole genome and the whole exome; FIG18b shows the performance analysis of testing a targeted group mixture dataset using a model trained on the whole genome and the whole exome;

[0098] FIG19 shows the size distribution of ground truth INDELs in Dream stage 3, Dream stage 4, platinum two-sample mixture, platinum tumor spike-in, PacBio, and exome datasets;

[0099] Figures 20a and 20b. Figure 20a shows the performance analysis of various detection methods testing INDELs based on predicted somatic mutation locations and types for various datasets;

[0100] Figure 21 shows the performance analysis of cross-sample training on the Dream Challenge Phase 3 dataset;

[0101] Figure 22a shows the performance of different somatic mutation detection methods on the real dataset COLO-829; Figure 22b FIG22 c shows the performance of different somatic mutation detection methods on the real dataset CLL1; FIG22 c shows the performance of different somatic mutation detection methods on the real dataset TCGA-AZ-6601;

[0102] Figure 23 is a list of 261 TCGA cancer samples used for Microsoft Azure experiments;

[0103] Figures 24a to 24m show the different network architectures tested;

[0104] FIG25 shows the performance analysis of the different network architectures shown in FIG24a to FIG24m;

[0105] Figures 26a and 26b show runtime comparisons of different somatic mutation detection algorithms;

[0106] FIG27 shows the VAF distribution of ground truth SNVs and INDEL somatic mutations in a superset of call results for HCC1395 classified into four confidence levels (high, medium, low, and unclassified) by the SEQC-II consortium;

[0107] Figures 28a to 28e show the overall performance of NeuSomatic on the 123 replicates in the SEQC-II dataset;

[0108] Figures 29a to 29d show the performance comparison of various models on the 21 WGS replica datasets. Fig.29a The F1 score (%) is shown, and Fig. 29b shows the precision-recall comparison across different somatic mutation callers. Here, the SEQC-WGS-GT50-spiked WGS10 trained model was used for NeuSomatic and NeuSomatic-S. Fig.29c The comparison of F1 scores (%) of various NeuSomatic trained models is shown. Figure 29d shows the F1 score (%) performance of different somatic mutation detection methods and NeuSomatic network trained models based on WGS datasets;

[0109] Figures 30a to 30d show the performance comparison on the tumor purity dataset. Figure 30a shows the comparison of SNV and INDEL F1 scores (%) for different somatic mutation callers across different coverages (10× to 300×) and tumor purity (5% to 100%). Figure 30b shows the robustness of matched normal samples to tumor contamination: at 80× coverage and 5% to 100% tumor purity, the F1 score (%) changes when the matched normal samples are contaminated with 5% tumor compared to the F1 score (%) of pure normal samples. For Figures 30a and 30b, the SEQC-WGS-GT50-spiked WGS10 trained model was used for NeuSomatic and NeuSomatic-S. Figure 30c shows the comparison of F1 scores (%) for various NeuSomatic trained models. Figure 30d shows the performance analysis for different INDEL sizes. Negative INDEL size reflects deletion;

[0110] Figures 30e and 30f show the precision-recall analysis based on the tumor purity dataset. Figure 30e shows the SNVs, and Figure 30f shows the INDEL accuracy comparison of different somatic mutation callers across different coverages (10× to 300×) and tumor purities (5% to 100%);

[0111] Figure 30g shows the precision-recall comparison of various NeuSomatic-S trained models based on the tumor-normal titration dataset;

[0112] Figure 30h shows the F1 score (%) performance of different somatic mutation detection methods and the NeuSomatic network trained model based on the tumor-normal titration dataset;

[0113] Figures 31a to 31f show the performance comparisons for library preparation, WES, and AmpliSeq datasets. Fig.31aFigure 31a shows the comparison of SNV and INDEL F1 scores (%) for different somatic mutation callers across different library-kits and DNA amounts. Figure 31b shows the comparison of precision-recall for various NeuSomatic trained models for different library-kits and DNA amounts. Figure 31c and Figure 31e show the comparison of SNV F1 scores (%) for different somatic mutation callers. Figure 31d and Figure 31f show the comparison of F1 scores (%) for various NeuSomatic trained models. For Figure 31a, Figure 31c, and Figure 31e, the SEQC-WGS-GT50-spiked-in WGS10 trained model was used for NeuSomatic and NeuSomatic-S;

[0114] Figure 31g shows the precision-recall comparison of various NeuSomatic-S trained models based on the library preparation dataset;

[0115] Figure 31h shows the F1 score (%) performance of different somatic mutation detection methods and the NeuSomatic network trained model based on the library preparation dataset;

[0116] Figure 31i shows the F1 score (%) performance of different somatic mutation detection methods and the NeuSomatic network trained model based on the WES dataset;

[0117] Figure 31j shows the F1 score (%) performance of different somatic mutation detection methods and the NeuSomatic network trained model based on the AmpliSeq dataset;

[0118] Figures 32a to 32d show the performance comparison for the FFPE datasets. Figures 32a and 32b show the performance based on 16 FFPE WGS replicates. Figures 32c and 32d show the performance based on 14 FFPE WGS replicates. Figures 32a and 32c show the F1 score (%) comparison of different somatic mutation callers across FFPE and fresh normal samples. Figures 32b and 32d show the F1 score (%) comparison of various NeuSomatic trained models. Figure 32b shows a precision-recall plot for 16 FFPE WGS replicates using different somatic mutation callers for FFPE and fresh normal samples. For Figures 32a and 32c, the SEQC-WGS-GT50-spiked-in WGS10 trained model was used for NeuSomatic and NeuSomatic-S;

[0119] Figure 32e shows the F1 score (%) performance of different somatic mutation detection methods and the NeuSomatic network trained model based on the FFPE WGS dataset;

[0120] Figure 32f shows the F1 score (%) performance of different somatic mutation detection methods and the NeuSomatic network trained model based on the FFPE WES dataset;

[0121] Figures 33a and 33b are graphs analyzing the impact of training on target sample based training using NeuSomatic (Figure 33a) and NeuSomatic-S (Figure 33b) on nine replicates from different datasets. For each sample, two sample-specific models were used: one trained on the target sample only, and the other trained using the additional 10% of the training data from the SEQC-WGS spike-in model;

[0122] Figure 33c shows the F1 score (%) performance of the SEQC-WGS-spiking model and two sample-specific NeuSomatic network trained models;

[0123] Figures 34a to 34d show performance analysis of different INDEL sizes, different VAF distributions, and on difficult regions. Figures 34a and 34b show performance analysis of different INDEL sizes across the SEQC-II dataset using different callers and different training methods. Negative INDEL sizes reflect deletions. Figure 34c shows performance analysis of mutations with different VAF ranges across the SEQC dataset using different callers. Figure 34d shows performance analysis of different somatic callers on difficult regions, which include tandem repeats (TRs) and segmental duplications of different sizes;

[0124] Figure 35 shows the performance analysis of mutations with different VAF ranges across the SEQC dataset using different callers.

[0125] Figures 36a and 36b show the performance analysis of mutations with different VAF ranges across the SEQC dataset using different training methods for NeuSomatic (Figure 36a) and NeuSomatic-S (Figure 36b);

[0126] Figures 37a to 37c show the results of using different callers (Figure 37a) and targeting NeuSomatic ( Figure 37b ) and NeuSomatic-S ( Figure 37c ) training methods ( Figure 37b and Figure 37c ) on mutations with different VAF ranges for the tumor / normal titration dataset;

[0127] Figure 38 shows the performance comparison of using different callers on the whole genome and difficult genomic regions including tandem repeats (TRs) and segmental duplications of different sizes;

[0128] Figures 39a and 39b show different responses to NeuSomatic (Figure 39a) and NeuSomatic-S ( Fig.39b ) on a difficult region consisting of tandem repeats (TRs) of varying sizes and segmental duplications;

[0129] Figures 40a and 40b show the performance comparisons using different training methods for NeuSomatic and NeuSomatic-S on the whole genome and difficult genomic regions, which include tandem repeats (TRs) and segmental duplications of different sizes;

[0130] Figures 41a to 41c show the results of using different callers (Figure 41a) and targeting NeuSomatic ( Fig.41b ) and NeuSomatic-S (Figure 41c) training methods (Figure 41b and Figure 41c) on difficult regions on the tumor / normal titration dataset, which include tandem repeats (TRs) of different sizes and segmental duplications; and

[0131] Figure 42 shows a violin plot comparison of the VAF distribution of private FN calling results based on the WGS dataset. DETAILED DESCRIPTION

[0132] In addition, the present disclosure (among other things) describes a method for germline variant calling, which may include obtaining a reference sequence, a plurality of sequence reads, and positions of candidate variants within the sequence reads; obtaining an augmented sequence read by inserting one or more spaces into the one or more sequence reads; obtaining an augmented reference sequence by inserting one or more spaces into the reference sequence; converting fragments of the augmented sequence reads around the candidate variant into a sample matrix; converting fragments of the augmented reference sequence around the candidate variant into a reference matrix; providing the sample matrix and the reference matrix to a trained neural network; and, at the output of the trained neural network, obtaining prediction data about variants within the plurality of sequence reads. The disclosed systems and methods allow for capturing important variant signals directly from raw data, and consistently achieve high accuracy rates for different sequencing technologies, sample purities, and sequencing strategies (such as whole genome versus target enrichment).

[0133] FIG. 1 illustrates an illustrative system 100 that includes a sequencing device 110 communicatively coupled to a computer system 102. The sequencing device 110 may be coupled to the computer system 102 directly (e.g., via one or more communication cables) or via a network 130, which may be the Internet or any other combination of a wide area network, a local area network, a wired network, and / or a wireless network. In some embodiments, the computer system 102 may be included in or integrated into the sequencing device 110. In some embodiments, the sequencing device 110 may sequence a sample comprising genetic material and generate resulting sequencing data. The sequencing data may be sent to the computer system 102 (e.g., via the network 130), or stored on a storage device and sent to the computer system 102 at a later stage (e.g., via the network 130). In some embodiments, the computer system 102 may or may not include a display 108 and one or more input devices (not shown) for receiving commands from a user or operator (e.g., a technician or geneticist). In some embodiments, computer system 102 and / or sequencing device 110 may be remotely accessible by a user or other device via network 130. Thus, in some embodiments, various methods discussed herein may be remotely executed on computer system 102.

[0134] The computer system 102 may include one or a combination of multiple computer devices of any type, such as a personal computer, a laptop computer, a network server (e.g., a local server or a server included in a public / private / hybrid cloud), a mobile device, etc., some or all of which may be interconnected. The computer system 102 may include one or more processors (not shown), each of which may have one or more cores. In some embodiments, the computer system 102 may include one or more general-purpose processors (e.g., CPUs), special-purpose processors such as graphics processors (GPUs), digital signal processors, or any combination of these and other types of processors. In some embodiments, some or all of the processors in the computer system may be implemented using custom or customizable circuits such as application-specific integrated circuits (ASICs) or field-programmable gate arrays (FPGAs). In some embodiments, the computer system 102 may also retrieve and execute non-transitory computer-readable instructions stored in one or more memories or storage devices (not shown) that are integrated into or otherwise communicatively coupled to the computer system 102. The memory / storage device may include any combination of non-transitory computer-readable storage media, including various types of semiconductor memory chips (DRAM, SRAM, SDRAM, flash memory, programmable read-only memory), etc. Disks and / or optical disks may also be used. The memory / storage device may also include a removable storage medium that is readable and / or writable; examples of such media include compact disks (CDs), read-only digital versatile disks (e.g., DVD-ROMs, dual-layer DVD-ROMs), read-only and readable Blu-ray® disks, ultra-density optical disks, flash memory cards (e.g., SD cards, mini SD cards, micro SD cards, etc.), etc. In some embodiments, data and other information (e.g., sequencing data) may be stored in one or more remote locations (e.g., cloud storage) and synchronized with other components of the system 100.

[0135] Germline variant calling

[0136] 2a illustrates an exemplary method 200 for germline variant calling according to some embodiments. The method 200 can be implemented, for example, in the form of software (i.e., a set of instructions stored in one or more non-transitory media that can be accessed and executed by, for example, one or more processors of the computer system 102) or in the form of firmware, hardware, or any combination thereof.

[0137] Method 200 may begin at step 210, where a reference sequence may be obtained. The reference sequence may be obtained, for example, from one or more private or public repositories, such as RefSeq, an open access, annotated, and curated collection of nucleotide base sequences established by the National Center for Biotechnology Information (NCBI) or the NCBIGenomes FTP website, which stores a set of complete genomes of different organisms. It should be understood that in some embodiments, a particular copy of the reference sequence may be stored locally (e.g., stored in the memory of computer system 102), while in other embodiments, the reference sequence may be obtained from a remote server, for example, via network 130. Furthermore, in some embodiments, the entire reference sequence may be obtained, while in other embodiments, one or more segments of the reference sequence may be obtained, for example, only one or more segments associated with a specific assay. Therefore, as used herein, "reference sequence" generally refers to one or more segments of a reference sequence, which may or may not include the entire reference sequence.

[0138] In step 220, a plurality of sequence reads corresponding to a genetic sample (e.g., a sample containing a patient's DNA or RNA material) sequenced by the sequencing device 110 may be obtained. As described above, the sequence reads may be obtained directly from the sequencing device 110 or from one or more local or remote volatile or non-volatile memories, storage devices, or databases communicatively coupled to the computer system 102. The sequence reads obtained may have been pre-processed (e.g., pre-aligned), or they may be "raw", in which case the method 200 may also include a pre-processing (e.g., pre-alignment) step (not shown). Moreover, while in some embodiments, a complete sequence read (generated by the sequencing device 110) may be obtained, in other embodiments, only a segment of a sequence read may be obtained. Therefore, as used herein, "obtaining a sequence read" generally refers to obtaining one or more segments of one or more (e.g., adjacent) sequence reads.

[0139] At step 230, multiple candidate variant positions within the sequence read can be obtained. In some embodiments, this includes scanning the multiple sequence reads obtained, comparing them with the reference sequence obtained, and determining a set of one or more positions within the sequence read that appear to include some types of variants (e.g., insertion or deletion variants (INDELs), single nucleotide variants (SNVs), or structural variants (SVs)). Finding candidate variant positions can include, for example, checking all positions in the reference and determining the position where at least one sequence read is different from the reference. When determining the candidate positions, some filters can also be used, for example, for a specific position to be considered as a candidate variant position, a filter for the number or percentage of reads that must be different from the reference. In other embodiments, the sequence reads may have been scanned and analyzed (e.g., by independent software and / or hardware modules), in which case multiple candidate positions can be obtained from a memory or database coupled to the computer system 102. At step 240, the method can advance from multiple identified candidate variant positions to the next candidate variant position.

[0140] At step 250, the sequence reads and reference sequence are expanded around the current candidate variant position to achieve an accurate multiple sequence alignment (MSA). This step can be better illustrated in conjunction with FIG3, which shows a plurality of exemplary sequence reads 311, an exemplary reference sequence 310, and an exemplary candidate variant position 350. As shown in FIG3, expanding the reference sequence 310 and the sequence reads 311 may include inserting one or more gaps or spaces (e.g., spaces 325) into the reference sequence 310 and / or into one or more sequence reads 311, thereby generating an expanded reference sequence 320 and a plurality of sequence reads 321 that are accurately aligned with the expanded reference sequence 320.

[0141] Specifically, as shown in FIG3 , in some embodiments, the expansion may include: determining or detecting an inserted base in a specific sequence read; inserting a space into a reference sample at the position of the inserted base; inserting a space at the same position into each sequence read in which no insertion is detected at the same position (in some cases, this corresponds to all sequence reads except the specific sequence read in which the insertion is detected); and repeating the process for each inserted base detected in each sequence read. It should be understood that two or more spaces may be inserted respectively in positions where potential insertions of two bases or longer are detected. It should be understood that other suitable techniques may be used to insert spaces into sequence reads and reference sequences, so that an accurate multiple sequence alignment (MSA) may be ultimately achieved, meaning that all bases within the expanded sequence reads (or at least those positions where no deletions or insertions occur) may be accurately aligned with each other and with the expanded reference sequence.

[0142] In some embodiments, augmentation may also include detecting missing bases in a particular sequence read and inserting spaces into the sequence read at potential deletion variant positions. It should be understood that in the event of a missing base being detected, the reference sequence does not need to be augmented.

[0143] It should also be understood that "inserting spaces" may or may not include actual modifications to the memory in which the sequence reads and reference sequences are stored, and may be replaced by inherent or implicit portions of other operations (e.g., other steps of method 200) to achieve the same end result. For example, instead of modifying the original locations of the sequence reads and reference sequences, the sequence reads and reference sequences may be selectively copied from the original locations to independent memory locations in a manner that will cause new copies to store spaces at all correct locations. As used herein, "storing spaces" may mean storing any random or predetermined value (e.g., 0) other than the four values ​​corresponding to four types of nucleotide bases (A, C, T, and G). As another example, as described below, when the sequence reads and reference sequences are converted into a count matrix, spaces may not be physically inserted in any step, but may be accounted for algorithmically.

[0144] In some embodiments, expanding the sequence reads and reference sequence around the candidate variant position may mean expanding only a portion of the variable length or predetermined length that includes the candidate variant position. However, in other embodiments, the entire length of the available sequence reads and reference sequence may be expanded at once, in which case step 250 may be performed only once (e.g., before step 240 or even before step 230), rather than repeating step 250 for each candidate variant position.

[0145] At step 260, the segments (i.e., windows) of the expanded sequence reads may be converted (or "summed up") into a sample matrix. The segments may have a fixed length and may or may not place the candidate variant position in the middle. For example, in the example shown in FIG3 , the expanded segment is 16 columns (positions) long and includes 7 columns to the left and 8 columns to the right of the candidate variant position 350, thereby placing the candidate variant position 350 substantially in the middle of the expanded segment.

[0146] In some embodiments, the sample matrix can be a count matrix having the same number of columns as the fragments being converted. For each type of nucleotide base (A, C, G, and T), the matrix can include a row, each row representing the number of times the corresponding nucleotide base appears at each position (column) of the extended sequence read. Referring to the example of FIG. 3 , row “A” in the sample matrix 331 includes the value “6000000306520000”, indicating that base “A” appears 6 times at the first position within the extended read fragment; 3 times at the eighth position; 6 times at the tenth position; 5 times at the eleventh position; twice at the twelfth position; and 0 times (no occurrence) at any other position within the fragment.

[0147] In some embodiments, the sample matrix may also include a row that indicates the number of times a space (denoted as "-" in FIG. 3) occurs at each position within the expanded sequence read. For example, in FIG. 3, the first row (labeled "-") of the sample matrix 331 includes the value "0000014000140000", indicating that the expanded sequence read includes 1, 4, 1, and 4 spaces at positions 6, 7, 11, and 12 within the segment, respectively. It should be understood that, although the row of spaces "-" is shown as the first row of the matrix, in other embodiments, the row of spaces may be placed elsewhere within the matrix (e.g., at the end). It should also be understood that the row / column notation is arbitrary, and in other embodiments, the notation may be reversed, meaning that the sample matrix may have 5 columns (A, C, T, G, -), and the number of rows corresponds to the length of the segment. Accordingly, the matrix can generally be described as having at least five "rows" (horizontal rows or columns), one row representing counts of each nucleotide type, and at least one row representing counts of spaces within the expanded sequence reads.

[0148] In step 270, the segments of the augmented reference sequence corresponding to the segments of the sequence reads can be similarly converted into a reference matrix having the same dimensions as the sample matrix. For example, in FIG3, augmented reference 320 is converted into reference matrix 330. Similar to sample matrix 331, reference matrix 330 can have four rows indicating the number of times each type of nucleotide occurs at each position within the augmented reference sequence, and a row of spaces indicating the number of spaces at each position within the augmented reference sequence. However, it should be understood that in other embodiments, reference matrix 330 can have a different summary / representation of different bases of the reference sequence than that shown in FIG3.

[0149] In the embodiment shown in FIG. 3 , the reference matrix 320 is also normalized to have the same range of values ​​as the sample matrix 331. As shown in FIG. 3 , normalization may include multiplying each count by a normalization factor, such as the total number of sequence reads (in this example, multiplied by 6). It will be appreciated that in other embodiments, instead of normalizing the reference matrix, normalization may be performed on the sample matrix or on both matrices. In other embodiments, the matrices may not need to be normalized at all. It will also be appreciated that step 270 need not be performed after step 260, and in some embodiments may be performed before or in parallel with step 260.

[0150] At step 280, the reference matrix and the sample matrix can be provided as inputs to a trained deep neural network, and at step 290, an output of the trained deep neural network can be obtained, wherein the output can include various predictions regarding variants included within the sequence reads represented by the sample matrix, as discussed in more detail below.

[0151] Those skilled in the art should also understand that, unless explicitly indicated to the contrary, the number and order of steps in method 200 are not necessarily limited to the number and order of steps shown in FIG. 2a. Therefore, in different embodiments, some steps of method 200 may be reordered or performed in parallel; certain steps may be set as optional or omitted; and some additional steps may be added without departing from the spirit and scope of the present disclosure.

[0152] For example, FIG. 2 b shows another embodiment, wherein at step 210 ', a reference sequence and a plurality of corresponding sequence reads are obtained. The reference sequence and the plurality of corresponding sequence reads may be obtained in any order (i.e., sequentially or simultaneously). Next, at step 220 ', the reference sequence and the plurality of corresponding sequence reads are compared. In some embodiments, the reference sequence and the plurality of corresponding sequence reads may be obtained in an compared format, which may eliminate the need to perform a first comparison step. At step 230 ', candidate variant positions are identified from the compared sequence reads and the reference sequence. At step 240 ', multiple sequence reads and / or reference sequences are expanded around the candidate variant positions in the same manner as described above (i.e., by inserting spaces to account for insertions and / or deletions) to achieve a second comparison of the multiple sequence reads with the reference sequence.

[0153] In some embodiments, a window of about 2 to 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, or 100 bases at any candidate variant position may be expanded and aligned. In some embodiments, a window of at least about 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, or 100 bases on either side of a candidate variant position may be expanded and aligned. In some embodiments, a window of no more than about 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, or 100 bases on either side of a candidate variant position may be expanded and aligned.

[0154] At step 250', a reference matrix for candidate variant positions is generated from the augmented reference sequence, and a sample matrix for candidate variant positions is generated from the plurality of augmented sequence reads. As described above, the matrix includes the frequency of each base or gap (i.e., A, C, G, T, or gap) at each sequence position. At step 260', the reference matrix and the sample matrix are input into the trained neural network.

[0155] At step 270', the neural network determines whether a variant type exists at the candidate variant position. For example, the neural network may include multiple convolutional layers that process the reference matrix and the sample matrix and generate outputs processed by one or more classifiers and regressors (i.e., a classifier of mutation type and a length classifier, and a position regressor) to determine the variant type, size, and position. In some embodiments, the neural network may include up to about 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, or 100 layers. In some embodiments, the convolutional layers may be structured into blocks. In some embodiments, the number of blocks is about half the number of layers.

[0156] If the reference sequence and corresponding sequence reads include more than one variant candidate position, steps 230' to 260' may be repeated for each variant candidate position. In addition, steps 210' to 260' may be repeated by obtaining additional reference sequences and additional corresponding sequence reads.

[0157] Although the examples shown in Figures 2a and 2b are described in the section for germline variant calling, these examples can be applicable to somatic variant calling. On the contrary, the embodiments described in the section for somatic variant calling can also be applicable to germline variant calling. For example, the multiple sequence reads described in Figure 2b can be subdivided into two groups: normal sequence reads and tumor sequence reads. Normal sequence reads can be obtained from a normal tissue sample of a subject, and tumor sequence reads can be obtained from a tumor tissue sample of a subject. Instead of generating a single sample matrix, a normal matrix and a tumor matrix will be generated together with a reference matrix. The rest of the method will be performed in a similar manner.

[0158] FIG4 illustrates an exemplary embodiment in which a reference matrix 330 and a sample matrix 331 (along with other optional inputs discussed below) are fed into a trained deep neural network 350 , which then outputs one or more predictions about the sequence reads represented by the sample matrix 331 .

[0159] In some embodiments, network 350 may be a convolutional neural network (CNN), an illustrative example of which is shown in FIG5. As shown in FIG5, network 350 may include a plurality of convolutional layers (510A-510I) connected in series, which ultimately feed into a first fully connected layer 520A, the output of which feeds into four independent fully connected layers 520A-520C that output various predictions. Figure 5 In the example of FIG. 5 , different convolutional layers 510 may have different filter sizes (e.g., 1x3, 3x3, or 5x5) and may be interconnected through additional processing layers (not explicitly shown), such as rectified linear units (ReLU), pooling layers, batch normalization (BN) layers, etc. In addition, as shown in FIG. 5 , in some embodiments, the input into a particular layer may sometimes be connected to the output of a subsequent layer via an identification shortcut connection 515 to help keep the signal in deeper layers of the network.

[0160] In some embodiments, the network 350 may output one or more predictions associated with candidate variants within the fragments of the augmented sequence reads. As shown in the example of FIG. 5 , one of such predictions may be obtained at the output of the fully connected layer 520B and may indicate whether the candidate variant position corresponds to a variant, and if so, the type of variant. More specifically, in some embodiments, the prediction may include at least four probability values ​​associated with four possible types of variants: NONE (no variant, i.e., a false positive call result); SNP / SNV (single nucleotide polymorphism / variant); INS (insertion variant); and DEL (deletion variant). Based on these probability values, it may be determined which type of variation is most likely, as well as the confidence of the specific prediction.

[0161] In some embodiments, the neural network 350 may also output (e.g., at the output of the fully connected layer 520C) the predicted position of the variant within the segment of the augmented sequence read represented by the sample matrix 331. In other embodiments, this output may be omitted, and the position of the variant may be assumed to be known (e.g., at the center of the segment) based on the manner in which the segment of the sequence read is selected for a given candidate variant, as described above.

[0162] In some embodiments, the neural network 350 may also output (e.g., at the output of the fully connected layer 520D) the predicted length of the variant. For example, if the predicted variant is one base long (e.g., a SNP variant or a one-base DEL or INS variant), a length of "1" may be output; if the predicted variant is a two-base long DEL or INS variant, a length of "2" may be output; and so on. In some embodiments, if the output indicates that the variant is 2 bases long or longer, a post-processing step may resolve the exact sequence that is being inserted or deleted.

[0163] In some embodiments, neural network 350 may also output (e.g., at the output of fully connected layer 520E) a predicted genotype associated with the variant, e.g., indicating the probability that the variant is (1) homozygous reference (non-variant); (2) heterozygous variant (only the maternal or paternal copy has the variant); (3) homozygous variant (where both copies have the same variant); or (4) other (each copy has a different variant). It should be understood that in other embodiments, some of these outputs may be omitted and / or additional outputs may be added.

[0164] Referring now to the inputs to the network, in addition to the reference matrix 330 and sample matrix 331 discussed above, in some embodiments, other inputs may be provided to further improve the accuracy of the predictions and / or expand the amount of information that can be obtained at the output of the network. As shown in the example of FIG. 4 , the additional inputs may include, for example, a variant position matrix 340, a coverage matrix 343, and one or more alignment feature matrices 344.

[0165] In some embodiments, all input matrices may be provided to the network 350 as one large three-dimensional matrix. The dimensions of the matrix may be, for example, 5 xs xk, where 5 corresponds to the number of rows / columns (-, A, C, T, G); s corresponds to the length of the segment (e.g., 16); and k corresponds to the number of different two-dimensional (e.g., 5 xs) matrices described above. In some embodiments, k may be as high as 30, or even higher.

[0166] Variant position matrix 340 can include a two-dimensional representation of candidate variant positions within a reference matrix and a sample matrix. For example, the matrix can have five rows (-, A, C, T, G), and each column can represent a position within a fragment. For example, variant position matrix 340 can include one value (e.g., 1) in all rows of a column corresponding to a candidate variant position, and another value (e.g., 0) in all rows of all other columns.

[0167] The coverage matrix 343 can represent the coverage or depth of a specific fragment of a sequence read. For example, in some embodiments, the coverage matrix 343 can include the same value in all its elements, which represents the coverage / depth (e.g., average coverage / depth) of different reads within the fragment. In other embodiments, the coverage matrix 343 can include different values ​​at different columns, each value representing the corresponding coverage / depth of each column. The alignment feature matrix 344 can represent various metrics associated with the quality of sequence reads and their alignments. Such metrics can include, for example, base quality, mapping quality, chain bias, shearing information, etc. These alignment features can help improve prediction accuracy by providing more context for sequence reads, for example, by informing the network of the quality of sequence bases and / or their mappings. Additional inputs to the network 350 (not shown for brevity) can include various data about known variants. Such data can be obtained from public and / or private knowledge bases such as dbSNP, COSMIC, ExAC, etc.

[0168] It should be understood that FIG. 5 shows only one exemplary configuration of network 350, and neural networks with other suitable configurations / architectures may be used to analyze the augmented sequence reads (represented by the sample matrix) and predict the type of variant and its other characteristics without departing from the scope and spirit of the present disclosure. Such alternative configurations may include fewer layers, additional layers, layers with different parameters, additional inputs or outputs, fewer inputs or outputs, and the like. In addition, in some embodiments, network 350 may not be a convolutional neural network (CNN), but may be another type of deep neural network (DNN), i.e., another type of artificial neural network (ANN) having multiple layers between the input layer and the output layer, without departing from the scope and spirit of the present disclosure.

[0169] Somatic variant calling

[0170] In some embodiments, the neural network 350 can also be trained to perform somatic variant calling. FIG6 illustrates an exemplary method 600 for performing somatic variant calling using the trained neural network 350.

[0171] At step 610, method 600 includes obtaining a reference sequence, for example, in a manner similar to that described above in conjunction with step 210 of method 200. At step 620, method 600 includes obtaining a plurality of tumor sequence reads and a plurality of corresponding normal (non-tumor) sequence reads. For example, in some embodiments, the tumor sequence reads may include sequencing results of a patient's tumor (cancerous) tissue, while the normal sequence reads may include sequencing results of normal (non-cancerous) tissue taken from the same patient. In other embodiments, the normal sequence reads may include sequencing results of normal tissue taken from a different patient. Both types of sequence reads may be obtained, for example, in a manner similar to that described above in conjunction with step 220 of method 200.

[0172] At step 630, multiple candidate somatic variant positions can be obtained by obtaining predetermined variant positions, i.e., a set of whitelist candidate positions, from a memory or a public / private database, or by performing a comparison between tumor sequence reads and a reference sequence (e.g., in the manner described above in conjunction with step 230 of method 200).

[0173] At step 640, the method proceeds to one of a plurality of candidate somatic variant positions. At step 660, segments of normal sequence reads are converted to a normal sample matrix around the candidate somatic variant position. Similarly, at step 665, segments of tumor sequence reads (substantially corresponding in length and location to segments of normal sequence reads) are converted to a tumor sample matrix around the candidate somatic variant position. At step 670, segments of reference sequence (substantially corresponding in length and location to segments of normal sequence reads) are converted to a reference matrix around the candidate somatic variant position.

[0174] In some embodiments, before converting the above three fragments into their corresponding reference matrices, each of them can be augmented using the augmentation technique described above in conjunction with step 250 of method 200. In other embodiments, different types of augmentation techniques can be applied to the fragments before converting them into matrices; for example, a k-mer collapse method can be used, in which homopolymers of length L>k are collapsed to length k, where k is any number greater than or equal to 1. In other embodiments, the augmentation step can be omitted entirely, and the original, unprocessed sequence reads and reference sequences can be directly converted into their corresponding matrices.

[0175] In the example shown in FIG7 , an augmentation similar to that performed in step 250 of method 200 is performed. In this example, the fragments of the reference sequence 310, the fragments of the tumor sequence reads 311-A, and the fragments of the normal sequence reads 311-B are first augmented to generate augmented reference sequences 320, augmented tumor sequence reads 321-A, and augmented normal sequence reads 321-B, respectively, and then converted into reference matrices 330, tumor sample matrices 331-A, and normal sample matrices 331-B, respectively.

[0176] Now returning to method 600, at step 680, the three matrices are fed into (provided as input) a trained neural network, and at step 690, the output of the trained neural network is obtained, wherein the output includes at least a prediction (probability estimate) of the type of somatic variant contained in the tumor sequence read (e.g., at or near the candidate somatic variant position).

[0177] An exemplary illustration of a trained neural network 350 for somatic variant calling is shown in Figure 8. In this example, the network 350 receives as its inputs a reference matrix 330, a tumor sample matrix 331-A, and a normal sample matrix 331-B.

[0178] It should be understood that although Figure 8 shows all three types of matrices as input to the network 350, in some embodiments, the reference matrix 330 and / or the normal sample matrix 331-B can be omitted and not provided to the network 350, and the network 350 can still perform somatic variant calling. In embodiments where the normal sample matrix is ​​omitted, obtaining normal sequence reads can be omitted from step 620 of method 600, and step 660 of converting the normal sequence reads into a matrix can be omitted entirely.

[0179] In some embodiments, to further improve the accuracy of prediction and / or expand the amount of information that can be obtained at the output of the network, the network may also receive the following additional matrices as input: such as a somatic variant position matrix 340 representing the positions of candidate somatic variants; a tumor coverage matrix 343-A and a normal coverage matrix 343-B representing the corresponding coverage of fragments of tumor sequence reads and fragments of normal sequence reads; and one or more alignment feature matrices 344, which represent the features of sequence reads in a manner similar to that described above in conjunction with FIG4. In addition, although not shown in FIG8 for simplicity, in some embodiments, data on known variants may be obtained from public and / or private knowledge bases such as dbSNP, COSMIC, ExAC, etc., and provided as input to the network 350.

[0180] In some embodiments, to further improve the accuracy of somatic variant calling, one or more other caller feature matrices 345 may also be provided to the neural network 350 at its input. These matrices may represent (e.g., summarize) one or more features obtained from one or more other (e.g., third-party) variant calling applications or algorithms that have processed the same fragment of the sequence read. Such features may include the type, length, and position of the variant predicted by the other application / algorithm, as well as quality scores, mapping scores, variant significance scores, and other alignment features. As described above in conjunction with Figure 4 As discussed, in some embodiments, all input matrices feeding into network 350 may be "combined" and provided as one large three-dimensional matrix.

[0181] FIG9 shows an exemplary trained neural network 350 that receives the various inputs described above and outputs predictions that indicate at least the variant type, position, and length of somatic variants contained in the fragments of somatic sequence reads represented by the somatic sample matrix 331-A at its input. In some embodiments, the network 350 may also include a fully connected layer 520E (as shown in FIG5 ) that outputs the genotype of the variant, but in the case of somatic variant calling, the genotype is not important and can be omitted from the network.

[0182] Based on the examples of FIG. 9 and FIG. 5 , it should be understood that in some embodiments, a network 350 having the same or substantially the same architecture can be used to predict germline variants and somatic variants, depending on the type of data it has been trained on (as discussed below), and depending on the type of input provided to it. However, in some embodiments, the architecture of the network 350 can be modified and optimized for either germline variant calling or somatic variant calling.

[0183] Neural network training

[0184] Those skilled in the art will appreciate that before the neural network 350 can begin to implement accurate germline variant calls or somatic variant calls, it first needs to be trained based on germline or somatic training data (i.e., training sequences), respectively. Training can include, for example, implementing all the steps of method 200 or 600 on a large number of training sequences, but then also providing the network with "ground truth" data (e.g., actual known variant types and their position lengths, genotypes, etc.) so that after each run, its trainable filters and other trainable parameters can be adjusted through a process known as "backpropagation" so that the network gradually minimizes its output error.

[0185] In some embodiments, the network 350 can be trained based on a genome such as the NA12878 genome with a well-characterized benchmark true value variant. In some embodiments, instead of or in addition to the real genome, various simulation-based strategies can be used to train the network. For example, in order to train a network for germline variant calls, a synthetic sample with a given set of variants can be simulated using the method described in Mu, JC et al. ("VarSim: a high-fidelity simulation and validation framework for high-throughput genome sequencing with cancer applications." Bioinformatics 31, 1469–1471 (2015)). As another example, in order to train a network for somatic variant calls, a normal sample in which a random variant has been incorporated can be fed into the network, for example, using the method described in Eving, A. et al. ("Combining tumor genome simulation with crowdsourcing to benchmark somatic single-nucleotide-variant detection," Nature methods (2015)). Alternatively or additionally, two different normal germline samples can be mixed in different proportions to generate synthetic tumor / normal training samples. As another example, a virtual tumor / normal training sample with a desired allele frequency distribution can be generated by switching reads between variant mapping of two germline samples.

[0186] Example 1

[0187] Preface

[0188] Somatic mutations are key hallmarks in cancer development, progression, and treatment. Accurate detection of somatic mutations is challenging due to cross-contamination of tumor-normal samples, tumor heterogeneity, sequencing artifacts, and coverage. In general, it is critical for accurate somatic mutation detection algorithms to effectively filter out false-positive calls introduced by the above issues and accurately retain elusive true-positive calls that may occur at low allele frequencies (AFs) or in regions of low complexity.

[0189] To date, a number of tools have been developed to address the problem of somatic mutation detection, including MuTect2 1 、MuSE2 、VarDict 3 、VarScan2 4 、Strelka2 5 and SomaticSniper 6 These tools employ different statistical and algorithmic approaches and perform well in certain cancers or sample types for which they were designed; however, they are limited in their generalizability across a wider range of sample types and sequencing technologies and therefore may show suboptimal accuracy in these scenarios. 7,8,9 In our early work, SomaticSeq 10 In , we used an orthogonal approach to maximize sensitivity by integrating algorithms. It also used machine learning to integrate nearly a hundred features to maintain high precision, resulting in improved precision over all individual methods. However, the machine learning backbone used in SomaticSeq relies on a set of extracted features for mutation localization. As a result, it cannot fully capture the original information in the genomic context of somatic mutations to further distinguish true somatic mutations from background errors, limiting its performance in challenging situations such as low-complexity regions and low tumor purity.

[0190] Here, we address the limitations of generalizability and complexity of statistical models for tumor sequencing data by leveraging deep convolutional neural networks (CNNs). CNNs have recently demonstrated remarkable performance in classification problems across different domains, including germline variant calling. 11、12、13 and skin cancer classification 14 Even so, the application of CNNs to the challenging problem of somatic mutation detection remains unexplored. The only previous attempt based on deep learning 15 A six-layer fully connected neural network is applied to a set of manually extracted features. This approach lacks the capabilities provided by the CNN architecture, which uses patterns found in local regions to learn feature representations directly from raw data. In addition, due to the complexity of the fully connected network, it is not as general and scalable as CNN.

[0191] We introduce NeuSomatic, the first CNN-based somatic mutation detection method that can effectively leverage signals derived from sequence alignment and other methods to accurately identify somatic mutations. Unlike other deep learning-based methods that focus on germline variants, NeuSomatic is addressing a larger unmet need in accuracy due to the complexity of tumor samples. It can effectively capture important mutation signals directly from raw data and consistently achieve high accuracy rates for different sequencing technologies, sample purities, and sequencing strategies such as whole genome vs. target enrichment.

[0192] result

[0193] NeuSomatic Overview

[0194] The input to the NeuSomatic network is candidate somatic mutations identified by scanning sequence alignments of tumor samples and matched normal samples (Figures 10a to 10d). Somatic mutations reported by other methods can also be included in this candidate list. For each candidate locus, we construct a 3D feature matrix M (of size k × 5 × 32) consisting of k channels, each of size 5 × 32, to capture the signal from the region centered on the locus. Each channel of the matrix M has five rows representing the four nucleotide bases and the gap base (“-”), and 32 columns representing the alignment columns around the candidate position.

[0195] Figure 10a shows a virtual example of input matrix preparation for a given candidate somatic SNV. Sequence alignment information in a window of seven bases around the candidate somatic mutation is extracted. The reference sequence is then augmented by adding gaps to account for insertions in the reads. The augmented alignments are then summarized into a reference matrix, a tumor count matrix, and a normal count matrix. The count matrix records the number of A / C / G / T and gap (“-”) characters in each column of the alignment, while the reference matrix records the reference base in each column. The count matrix is ​​then normalized by coverage to reflect the base frequency in each column. Independent channels are retained to record tumor and normal coverage.

[0196] Figure 10b shows the input 3D matrix and the proposed NeuSomatic network architecture. The input matrix consists of a reference channel, tumor and normal frequency channels, coverage and location channels, followed by several channels summarizing the alignment features. When used in ensemble mode, NeuSomatic also includes additional channels for other individual method features. The NeuSomatic network architecture consists of nine convolutional layers consisting of four blocks with shortcut identity connections. We use two softmax classifiers and a regressor in the last layer to predict mutation type, size, and location.

[0197] Figure 10c shows a virtual example of input matrix preparation for a given candidate somatic deletion. Sequence alignment information in a window of 7 bases around the candidate somatic mutation is extracted. The reference sequence is then augmented by adding gaps to account for insertions in the reads. The augmented alignments are then summarized into a reference matrix, a tumor count matrix, and a normal count matrix. The count matrix records the number of A / C / G / T and gap (“-”) characters in each column of the alignment, while the reference matrix records the reference base in each column. The count matrix is ​​then normalized by coverage to reflect the base frequency in each column. Independent channels are retained to record tumor and normal coverage.

[0198] Figure 10d shows a virtual example of input matrix preparation for a given candidate somatic insertion. Sequence alignment information in a window of 7 bases around the candidate somatic mutation is extracted. The reference sequence is then augmented by adding gaps to account for the insertion in the reads. The augmented alignments are then summarized into a reference matrix, a tumor count matrix, and a normal count matrix. The count matrix records the number of A / C / G / T and gap (“-”) characters in each column of the alignment, while the reference matrix records the reference base in each column. The count matrix is ​​then normalized by coverage to reflect the base frequency in each column. Independent channels are retained to record tumor and normal coverage.

[0199] The first three channels are the reference channel, tumor frequency channel, and normal frequency channel, which summarize the reference bases around the candidate loci and the frequencies of different bases in this region. To capture insertions, we augmented the reference sequence with gaps around the candidate loci corresponding to the insertions in the read alignment (Figure 10a, Figure 10c, and Figure 10d). Therefore, each column of the tumor frequency matrix and the normal frequency matrix represents the frequency of A / C / G / T / gap bases in the corresponding multiple sequence alignment (MSA) column of the tumor sample and the normal sample, respectively. The remaining channels summarize other features of the reads supporting different bases, such as coverage, base quality, mapping quality, strand bias, and shearing information. If NeuSomatic is used in integrated mode, we also use additional channels for features reported by separate somatic mutation detection methods. With this concise yet comprehensive structured representation, NeuSomatic can use the necessary information in tumor, normal, and reference to distinguish hard-to-capture somatic mutations with low AF from germline variants and sequencing errors. This design also enables the use of convolutional filters in CNNs to capture contextual patterns in matrix sub-blocks.

[0200] To compare with other CNN methods used in genomic problems, DeepVariant 11Use read pileups as input for germline variant calling. Instead, we use the base frequency summary of each column as input to our network. This simplifies the CNN structure, allowing for a substantially more efficient implementation. For example, DeepVariant takes about 1000 CPU core hours to call germline variants for a 30× whole genome sample. 16 , while the standalone version of NeuSomatic can detect somatic mutations from 30× tumor-normal paired samples in about 156 CPU core hours, despite processing two (tumor-normal) samples instead of one and looking for candidates at lower somatic AFs than 50% or 100% AF of germline. Another germline variant calling method, Clairvoyante 12 , using three channels to summarize the base counts of allele counts, deletions, and insertions in the center of the window. Instead, we summarized all of these events in a single base frequency matrix using the reference augmentation method described earlier, which clearly represents all insertion and deletion (INDEL) events in the entire window.

[0201] NeuSomatic adopts a novel CNN structure that can predict the type and length of candidate somatic mutations given a feature matrix M (Figure 10b). The proposed CNN consists of nine convolutional layers, which consist of four blocks with shortcut signature connections, inspired by ResNet. 17 , but with different formations to suit the proposed input structure. We use two softmax classifiers and a regressor in the last layer. The first classifier identifies whether the candidate is a non-somatic call result, SNV, insertion or deletion. The second classifier predicts the length of the somatic mutation through four categories (0 indicates non-somatic, or the length is 1, 2 or greater than 2), and the regressor predicts the location of the somatic mutation. Using the output of these classifiers, we identify a set of somatic mutations. If the length of the predicted INDEL is greater than 2, we implement a simple post-processing step on the reads overlapping with this position to resolve the INDEL sequence from the aligned CIGAR string of the read. This has been shown to perform well for data generated by the Illumina sequencer. For sequencing data with higher error rates, more complex local rearrangement post-processing is performed to resolve the INDEL sequence.

[0202] Since NeuSomatic can be used in both standalone and integrated modes, we use NeuSomatic-S to represent the standalone mode and keep NeuSomatic to represent the integrated mode. We compare NeuSomatic and NeuSomatic-S with state-of-the-art somatic mutation detection methods, including MuTect2 1 、MuSE 2 、SomaticSniper 6 、Strelka2 5 、VarDict 3 and VarScan2 4 ) and the ensemble method SomaticSeq 10 We compare and contrast the performance using multiple synthetic datasets and real datasets. We report below that the difficulty of somatic mutation detection on synthetic datasets increases in turn, taking into account the AF of somatic mutations in the dataset.

[0203] Comparison based on a platinum sample mixture dataset

[0204] For the first synthetic dataset, as in previous studies 5,10 Similarly, we combined two normal platinum genomes 18 Samples, NA12877 and NA12878, were mixed at tumor purity ratios of 70:30, 50:50, and 25:75 to create three tumor-contamination curves, and at a ratio of 5:95 to create contaminated normal samples. We also included testing on 100% pure normal samples and 50% pure tumor samples. We used germline variants in NA12878 (reference call results in NA12877) as the truth set for evaluation. Both NeuSomatic-S and NeuSomatic significantly outperformed all other methods (Figure 11a, Figure 11b, and Figure 11c). The improved performance of NeuSomatic over other methods is due to the lower and more challenging tumor purity (25:75 mixture). Overall, NeuSomatic achieved F1 scores of up to 99.6% and 97.2% for SNVs and INDELs, respectively, which is an improvement of up to 7.2% relative to the best method for this dataset at the lowest sample purity. For samples with 50% tumor purity, reducing the purity of normal samples from 100% to 95% had little effect on NeuSomatic performance (<0.3%), but it caused a decrease in SomaticSeq accuracy of approximately 3%.

[0205] Figures 11a and 11b are graphs for performance analysis of the platinum two-sample mixture dataset. In this dataset, four cases of tumor and normal sample purity were used (50% T:100% N, 70% T:95% N, 50% T:95% N, and 25% T:95% N). Figure 11a is a graph for precision-recall analysis: Confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. Fig.11b is a graph for performance analysis of INDEL accuracy (F1 score) for different INDEL sizes. Fig.11c is a table showing the performance of different somatic mutation detection methods on the platinum two-sample mixed dataset. For each method, we report the precision, recall, and F1 score for the quality score threshold that achieves the highest F1 score in the precision-recall curve. (RC: recall, PR: precision, F1: F1 score).

[0206] Comparison based on ICGC-TCGA DREAM challenge dataset

[0207] For the second synthetic dataset, we use the ICGC-TCGA DREAM challenge phase 3 and phase 4 datasets 19 , which were constructed by computationally incorporating mutations into the healthy genomes of paired normal samples with different AFs (see Methods). We mixed tumor and normal samples to create five tumor / normal purity cases. NeuSomatic-S outperformed all other independent methods for both Phase 3 and Phase 4 datasets, with SNVs averaging over 8% and INDELs averaging over 22% (Fig. 12a, Fig. 12b, Fig. 12c, Fig. 13a, Fig.13b, Fig. 13c). This improvement in performance increases as tumor purity decreases. We further observe that NeuSomatic (in ensemble mode) significantly outperforms both SomaticSeq and NeuSomatic-S, even though NeuSomatic-S still outperforms SomaticSeq in more challenging cases such as 25:75 mixtures of SNVs and 25:75 and 50:50 mixtures of INDELs. Overall, NeuSomatic achieves F1 scores up to 96.2% and 93.5% for SNVs and INDELs, respectively, which is up to 34.6% improvement over the best method at the lowest sample purity. For samples with 50% tumor purity, reducing the normal sample purity from 100% to 95% has little effect on NeuSomatic’s performance (~1.2% on average), while SomaticSeq and Strelka2 experience a >3% decrease in F1 score.

[0208] Figures 12a and 12b are graphs for performance analysis on the DREAM stage 3 dataset. In this dataset, five cases of purity of tumor and normal samples (100% T:100% N, 50% T:100% N, 70% T:95%N, 50% T:95%N, and 25% T:95%N) were used. Figure 12a is a graph for precision-recall analysis: Confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. Figure 12b is a graph for performance analysis of INDEL accuracy (F1 score) for different INDEL sizes. Figure 12c is a table showing the performance of different somatic mutation detection methods on the Dream challenge stage 3 dataset. For each method, we report the precision, recall, and F1 score for the quality score threshold that achieves the highest F1 score in the precision-recall curve. (RC: Recall, PR: Precision, F1: F1 score).

[0209] Figures 13a and 13b are graphs for performance analysis on the DREAM stage 4 dataset. In this dataset, five cases of purity of tumor and normal samples (100% T:100% N, 50% T:100% N, 70% T:95%N, 50% T:95%N, and 25% T:95%N) were used. Figure 13a is a graph for precision-recall analysis: Confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. Figure 13b is a graph for performance analysis of INDEL accuracy (F1 score) for different INDEL sizes. Figure 13c is a table showing the performance of different somatic mutation detection methods on the Dream challenge stage 4 dataset. For each method, we report the precision, recall, and F1 score for the quality score threshold that achieves the highest F1 score in the precision-recall curve. (RC: Recall, PR: Precision, F1: F1 score).

[0210] Comparison based on platinum tumor incorporation datasets

[0211] For the third synthetic dataset, as in previous studies 1,10 Similarly, we constructed tumor samples by spiking reads from NA12878 into NA12877 at variant positions in NA12878 with spiking frequencies taken from a binomial distribution with mean [0.05, 0.1, 0.2, 0.3], and used independent NA12877 reads as pure normal samples. Note that unlike earlier strategies that mixed samples at a fixed ratio to obtain somatic mutations with fixed AFs, this mixing approach generates somatic mutations with different AFs ranging from 0.025 to 0.3. NeuSomatic achieved F1 scores of up to 80.9% and 66.7% for SNVs and INDELs, respectively, which is up to 4% improvement over the best method (Figures 16a to Fig.16d ). For low AF somatic mutations, the performance improvement is even higher (AF = 0.025, 11% improvement; AF = 0.05, 8% improvement) (Figure 16b). Figure 16a shows the precision-recall analysis, where the confidence or quality score is used to derive the precision-recall curve. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. Fig.16bFigure 16c shows the performance analysis of INDEL accuracy (F1 score) for different INDEL sizes. Figure 16d shows the precision, recall, and F1 score of the quality score threshold that achieves the highest F1 score in the precision-recall curve for each method. (RC: recall, PR: precision, F1: F1 score).

[0212] Comparison between whole exome and targeted panels

[0213] To evaluate the performance of NeuSomatic for enrichment of different targets, we used the AshkenaziJewish trio 20 We trained NeuSomatic and SomaticSeq on the whole exome dataset and applied the trained models to the whole exome and this set. For the whole exome, NeuSomatic achieved F1 scores of 99.3% and 88.6% for SNVs and INDELs, respectively. For the targeted set, NeuSomatic and NeuSomatic-S consistently outperformed other methods with SNVF1 scores of 99.2%.

[0214] Figure 17a shows the performance analysis of various methods on a mixture of exome samples. In this dataset, four cases of tumor and normal sample purity were used (50% T:100% N, 70% T:95% N, 50% T:95% N, and 25% T:95% N). The confidence or quality score is used to derive the precision-recall curve. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. The training here is NeuSomatic, NeuSomatic-S, and SomaticSeq based on exome data.

[0215] Figure 17b shows the performance analysis of various methods on a target panel sample mixture. In this dataset, four cases of tumor and normal sample purity were used (50% T:100% N, 70% T:95% N, 50% T:95% N, and 25% T:95% N). The confidence or quality score is used to derive the precision-recall curve. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. The training here is NeuSomatic, NeuSomatic-S, and SomaticSeq based on exome data.

[0216] Figure 17c shows the performance of different somatic mutation detection methods on a mixed dataset of whole exome samples. For each method, we report the precision, recall, and F1 score of the quality score threshold that achieves the highest F1 score in the precision-recall curve. (RC: recall, PR: precision, F1: F1 score).

[0217] Figure 17d shows the performance of different somatic mutation detection methods on the targeted panel dataset. For each method, we report the precision, recall, and F1 score for the quality score threshold that achieves the highest F1 score in the precision-recall curve. (RC: recall, PR: precision, F1: F1 score).

[0218] The robustness of NeuSomatic was confirmed by applying the model trained on the whole genome platinum mixture data to both target enrichment sets, yielding similar performance (Figures 18a and 18b). Similar to the other datasets, for samples with 50% tumor purity, reducing the purity of normal samples from 100% to 95% resulted in minimal reduction of the NeuSomatic F1 score based on the whole exome dataset (on average, approximately 0.3%), while the F1 scores of SomaticSeq and Strelka2 decreased by >5%.

[0219] Figure 18a shows the performance analysis of the models trained on the whole genome (Platinum data, genome mixture) and the whole exome (HG003-HG004 exome mixture) tested on the exome mixture dataset. The confidence or quality scores are used to derive the precision-recall curves. The highest F1 score obtained by each algorithm is indicated in the legend and marked on the curves with a solid circle.

[0220] Figure 18b shows the performance analysis of the target group mixture datasets tested using models trained on the whole genome (Platinum data, genome mixture) and the whole exome (HG003-HG004 exome mixture). The confidence or quality scores are used to derive the precision-recall curves. The highest F1 score obtained by each algorithm is indicated in the legend and marked on the curves with a solid circle.

[0221] Comparison based on PacBio dataset

[0222] We further evaluated the performance of NeuSomatic on reads with higher error rates, especially reads from long-read sequencing platforms. We used tumor-normal paired samples with 20%, 30%, and 50% AF somatic mutations simulated based on raw PacBio reads (Figure 14a, Figure 14b, and Figure 14c). NeuSomatic identified somatic SNVs and INDELs with F1 scores as high as 98.1% and 86.2%, respectively, outperforming VarDict. 3 , with win rates up to 34.4% for SNVs and 53.2% for INDELs. This analysis demonstrates the power of NeuSocomic in detecting somatic mutations even when the error rate of sequence reads is as high as with PacBio long raw reads.

[0223] Comparison of different INDEL sizes

[0224] It is noteworthy that NeuSomatic consistently outperforms other methods for various INDEL sizes in different datasets (Fig. 11b, Fig. 12b, Fig. 13b, Fig. 14b, Fig. 16c, and Fig. 19). For large (greater than 50 bases) INDELs, the INDEL information is lost in the piled count matrix because most of the short reads with somatic INDELs are soft-skipped. For such cases, NeuSomatic benefits from the predictions of other methods because some methods such as VarDict and MuTect2 use local assemblies for predictions.

[0225] Figures 14a and 14b are graphs for performance analysis of different somatic mutation detection methods on the PacBio dataset. In this dataset, three cases of purity of tumor and normal samples were used (70% T:95% N, 50% T:95% N, and 25% T:95% N). Figure 14a is a graph for precision-recall analysis: Confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle. Figure 14b is a graph for performance analysis of INDEL accuracy (F1 score) for different INDEL sizes. Figure 14c shows the precision, recall, and F1 score for the quality score threshold that achieves the highest F1 score in the precision-recall curve for each method. (RC: recall, PR: precision, F1: F1 score). Figure 19 shows the size distribution of ground truth INDELs in Dream stage 3, Dream stage 4, Platinum two-sample mixture, Platinum tumor spike-in, PacBio, and Exome datasets. Negative sizes correspond to deletions.

[0226] INDEL Types and Positional Accuracy

[0227] For all discussed datasets, we also evaluated the performance of INDEL calling by different somatic mutation detection methods, using a more relaxed criterion of simply correctly predicting the location of somatic INDELs (and ignoring the exact INDEL sequence). Again, we observed a similar advantage of NeuSomatic over the other schemes, indicating that the main improvement was contributed by the proposed CNN structure rather than the post-processing INDEL resolution step (Figure 20a and Figure 20b).

[0228] Figure 20a shows the performance analysis of INDELs based on the location and type of somatic mutations predicted for the PacBio dataset (ignoring the accuracy of the exact predicted INDEL sequence) at three tumor purities (50%, 30%, and 20%) and 95% normal purity. The confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle.

[0229] Figure 20b shows the performance analysis of INDELs based on the location and type of predicted somatic mutations (ignoring the accuracy of the exact predicted INDEL sequences) for Dream stage 3, Dream stage 4, Platinum two-sample mixture, whole exome, and Platinum tumor spike-in datasets. For the first four datasets, three tumor purity schemes (70%, 50%, and 25%) were used, while the purity of normal samples was 95%. Confidence or quality scores were used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed on the curve and marked with a solid circle.

[0230] Read coverage analysis

[0231] To evaluate the impact of sequence coverage on different techniques, we downsampled the whole exome dataset to obtain samples with sequence coverage ranging from 20× to 100× (Figure 15a and Figure 15b). For different coverages, NeuSomatic consistently outperformed the other techniques. For low coverage samples, the improvement increased as the problem became more challenging. In addition, reducing the coverage from 100× to 50× had a small impact on NeuSomatic (~1.5% for SNVs and ~5% for INDELs), while SomaticSeq’s F1 score for both SNVs and INDELs dropped by ~20%, and Strelka2’s F1 score dropped by ~13% for SNVs and ~15% for INDELs. This analysis shows the robustness of NeuSomatic to coverage perturbations and its advantage in challenging scenarios, which can be seen at lower coverage.

[0232] Figure 15a is a graph of the performance analysis of the effect of sequence coverage on a whole exome sample mixture dataset. In this example, the purity of the tumor is 50% and the purity of the normal sample is 95%. The Y-axis shows the highest F1 score that each algorithm can achieve for sample alignments ranging from 20× to 100×.

[0233] FIG15b is a graph of a performance analysis of the effect of sequence coverage on a whole exome sample mixture data set. In this example, the purity of the tumor is 50% and the purity of the normal sample is 95%. The coverage of the tumor and normal alignments ranges from 20x to 100x. The confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is indicated in the legend and marked on the curves with a solid circle.

[0234] Training robustness

[0235] We evaluated the robustness of NeuSomatic trained on a specific purity by training and testing on different purities on the DREAM Challenge Phase 3 dataset. Figure 21 shows that performance only slightly degrades when we train and test on very different tumor purities. We also observed that training with data aggregated from multiple tumor purities performs just as well as training on the target tumor purity (Figure 21). This suggests that training sets incorporating multiple tumor purities are sufficient to obtain a model that is robust to changes in tumor purity.

[0236] Figure 21 shows the performance analysis of cross-sample training on the Dream Challenge Phase 3 dataset. We tested each sample with tumor purity of 70%, 50%, and 25%, using NeuSomatic models trained on different purities and models trained on pooled inputs from all different purities. Confidence or quality scores are used to derive precision-recall curves. The highest F1 score obtained by each algorithm is printed in the legend and marked on the curve with a solid circle.

[0237] Comparison based on real data

[0238] In the absence of somatic mutations 21 High-quality, comprehensive ground truth datasets such as the Genome-in-a-Bottle gold set for germline variants 22 In this case, it would be impossible to calculate F1 accuracy outside of synthetic data. Fortunately, we can employ existing datasets with validated somatic mutations to evaluate the accuracy performance of NeuSomatic on real data (see “Methods” for more details on how to extrapolate F1 scores on real data). We used two datasets: CLL1 23 The COLO-829 genome-wide dataset contains 961 validated somatic SNVs from patients with chronic lymphocytic leukemia. 24,25is a whole genome dataset derived from a persistent metastatic melanoma cell line with 454 verified somatic SNVs. To evaluate NeuSomatic with these two real WGS samples, we used a model trained on the DREAM challenge phase 3. As shown in Figure 22a and Figure 22b, NeuSomatic achieved the highest extrapolated F1 scores of 99.7% and 93.2% for the COLO-829 melanoma sample and CLL1 chronic lymphocytic leukemia sample, respectively. We also used the TCGA colorectal adenocarcinoma 26,27 NeuSomatic was evaluated on a whole exome sequencing (WES) sample (TCGA-AZ-6601), achieving the highest extrapolated F1 score of over 99.6% ( Figure 22 c).

[0239] Figure 22a shows the performance of different somatic mutation detection methods on the real dataset COLO-829. Figure 22b Figure 22c shows the performance of different somatic mutation detection methods on the real dataset CLL1. Figure 22c shows the performance of different somatic mutation detection methods on the real dataset TCGA-AZ-6601.

[0240] To demonstrate the scalability and cost-effectiveness of NeuSomatic on the cloud, we also processed 261 whole exome sequencing cancer samples from TCGA on the Microsoft Azure cloud platform using both integrated and standalone modes (Figure 23). These samples were collected from a variety of cancer types, including colorectal adenocarcinoma, ovarian serous adenocarcinoma, cervical squamous cell carcinoma, and endocervical adenocarcinoma. Although the cloud platform enables us to automatically launch compute instances on demand, it takes an average of 2.42 hours and 0.72 hours for the integrated and standalone modes to process each sample, respectively. Using Azure's preemptible compute instances (using the standard H16 instance type, each with 16 cores), the processing cost per sample is lower, at $0.77 for the integrated mode and $0.23 for the standalone mode. We also evaluated the accuracy of NeuSomatic on these samples by comparing with 44,270 validated SNPs in these samples, which provided us with an overall pattern recall of 98.9% and an individual pattern recall of 97.2%, respectively. Therefore, NeuSomatic can not only be used for different sequencing technologies or sequencing strategies, but also run on various computing platforms, including local HPC clusters and elastic cloud computing architectures.

[0241] Figure 23 is a list of 261 TCGA cancer samples used for the Microsoft Azure experiment. The samples were taken from three cancer types: colorectal adenocarcinoma (COAD), ovarian serous adenocarcinoma (OV), and cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC).

[0242] discuss

[0243] NeuSomatic is the first deep learning-based somatic mutation detection framework that is both high-performance and general. It achieves best-in-class accuracy for altered tumor purity across multiple datasets from synthetic to real, multiple sequencing strategies from whole genome to targeted, and multiple sequencing technologies ranging from short-read sequencing to high-error-rate long-read sequencing, when using the same CNN architecture. Specifically, for low tumor purity and low allele frequency, NeuSomatic significantly outperforms other state-of-the-art somatic mutation detection methods, demonstrating its ability to solve difficult problems. NeuSomatic leverages an efficient implementation of convolutional neural networks to solve the somatic mutation detection problem quickly and accurately. It uses a novel summary of tumor / normal alignment information as a set of input matrices that can effectively capture the main signals in the genomic context. Training the proposed CNN architecture based on these matrices allows feature representations to be learned directly from the raw data. The deep features learned from the observed training data can accurately identify important mutational signatures that distinguish true calls from artifacts introduced by sequencing errors, cross-contamination, or coverage bias. We believe that NeuSomatic can greatly improve upon the state-of-the-art by providing a very broadly applicable approach for somatic mutation detection.

[0244] method

[0245] ICGC-TCGA DREAM Challenge Data

[0246] Phase 3 data consisted of normal and tumor samples constructed by computationally spiking 7903 SNVs and 7604 INDEL mutations into the healthy genome of the same normal sample with three different AFs (50%, 33%, and 20%) to create synthetic but realistic tumor-normal sample pairs. Phase 4 data had a similar construction but with 16,268 SNVs and 14,194 INDELs in two subclones at 30% and 15% AF. We then constructed impure normal reads by mixing 95% normal reads with 5% tumor reads. We also constructed four tumor mixtures by mixing tumor and normal reads at ratios of 100:0, 70:30, 50:50, and 25:75, respectively. Therefore, somatic mutations across these four tumor mixing ratios ranged from 5% to 50% AF for the phase 3 data, while 3.75% to 30% AF for the phase 4 dataset.

[0247] Platinum Synthesis Tumor Data

[0248] We downloaded 200× Platinum genome samples NA12878 and NA12877 and their true germline variants (v2017-1.0) 18 , to construct virtual tumor and normal pairs (ENA accession number PRJEB3246). For normal samples, we downsampled NA12877 to 50×. For tumor samples, we constructed three 50× in silico mixture samples with tumor purity of 70%, 50%, and 25%, respectively, by independently downsampling NA12877 at 15×, 25×, and 37.5×, and mixing each with 35×, 25×, and 12.5× downsampled NA12878. We used heterozygous variants in NA12878 and homozygous variants (which were reference calls in NA12877, at least five bases away from NA12877 variants and 300 bases away from each other) as the truth set (1,103,285 SNVs and 174,754 INDELs) for the training and evaluation steps. Thus, the AF of somatic mutations across these three tumor mixture ratios ranged from 12.5% ​​to 70%, depending on the zygosity of the germline variant in NA12878.

[0249] We also generated another 50× virtual tumor sample by randomly spiking reads from downsampled (to 50× coverage) NA12878 into downsampled (to 50× coverage) NA12877 data at heterozygous variant locations in NA12878 and homozygous variant locations that were reference calls in NA12877. For each variant, we randomly assigned the frequency of spiked reads from a binomial distribution with mean [0.05, 0.1, 0.2, 0.3]. Thus, the mean somatic mutation AF ranged from 2.5% to 30%, depending on the zygosity of the variant. To avoid ambiguity in the truth set, we only used variants whose associated paired-end reads did not overlap with any other variant (316,050 SNVs and 46,978 INDELs). This generated a contaminated tumor sample with reads from NA12878. We also use another independent downsampled (to 50×) data of NA12877 as pure normal samples.

[0250] For both experiments, FastQ files and true germline variants were downloaded and compared with BWA-MEM (v0.7.15) 28 and Picard MarkDuplicates (v2.10.10) (https: / / broadinstitute.github.io / picard) and GATK IndelRealigner and Base QualityScore Recalibration (v3.7) 29 Compare.

[0251] Real tumor-normal pair data

[0252] We use the CLL1 chronic lymphocytic leukemia dataset 23 (Access: https: / / www.ebi.ac.uk / ega / datasets / EGAD00001000023) and COLO-829 Persistent Metastatic Melanoma Cell Line Dataset 24,25 (Access: https: / / www.ebi.ac.uk / ega / studies / EGAS00000000052) and evaluate our approach based on real tumor-normal paired data and published validated somatic mutation lists.

[0253] The COLO-829 dataset consists of 80× whole-genome sequenced tumor samples and their matched 60× normal blood COLO-829BL samples, with 454 validated somatic SNVs. CLL1 has whole-genome sequenced tumor samples and matched normal samples with coverage of 53× and 42×, respectively, with 961 published somatic SNVs.

[0254] TCGA-AZ-6601 26,27 The dataset is whole exome sequencing of colon adenocarcinoma tumor samples and their matched normal tissues from TCGA. Tumor and normal samples were sequenced at depths of 145× and 165×, respectively. We used TCGA 30 and COSMIC 31 The 952 validated SNVs in the database were used as the benchmark truth somatic mutations for this sample.

[0255] For real data, we calculated the extrapolated precision as the percentage of predicted somatic mutations that have been called by at least two independent methods or have been reported in the COSMIC database as confirmed somatic mutations in at least two samples of the same cancer type. We then calculated the extrapolated F1 score based on the harmonic mean of the recall and the extrapolated precision.

[0256] Whole-exome and targeted-group data

[0257] To assess NeuSomatic using different target enrichment assays, we used a 20We downloaded deep sequencing (200× coverage) whole exome alignment files for HG003 and HG004 (ftp: / / ftp-trace.ncbi.nlm.nih.gov / giab / ftp / ), as well as high-confidence germline variants (Genomes in a Bottle v3.3.2). We then constructed 70%, 50%, and 25% pure tumor samples using a mixture of random 70×, 50×, and 25× downsampling of HG004 and 30×, 50×, and 75× downsampling of HG003, respectively. We also constructed a 95% pure normal sample by mixing the 95× HG003 downsampled alignment and the 5× HG004 downsampled alignment. We used Agilent SureSelect Human All Exon V5 BED files for our analysis. Similar to the platinum composite tumor data (11,720 SNVs, 878 INDELs), ground truth somatic mutations were identified. The AF of somatic mutations across these three tumor mixture ratios ranged from 12.5% ​​to 70%, depending on the zygosity of the germline variants in HG004.

[0258] To validate the performance on the target panel, we restricted the alignment and ground truth data to the Illumina TruSight genetic disease panel BED file (216 SNVs, 5 INDELs). Due to the limited number of true INDELs in the target panel region, we only evaluated the performance on SNVs.

[0259] PacBio Data

[0260] For long read analysis, we downloaded high confidence germline variants for the HG002 sample (Genome in a Bottle v3.3.2) (ftp: / / ftp-trace.ncbi.nlm.nih.gov / giab / ftp / ) 20 We use the CHM1 dataset 32 (SRA accession number SRX533609) to establish the long read error distribution. We then used VarSim to simulate the 33 In silico long-read sequencer simulator with LongISLND 34 We simulated 100× pure normal samples using a combination of random somatic mutations. We also simulated 100× pure tumor samples with the same error distribution using a set of random somatic mutations. We used NGMLR (v0.2.6) 35to align sequences. We then mixed a 47.5× downsample of the pure normal alignment with a 2.5× downsample of the pure tumor alignment to form a 50× normal pair of 95% purity, and mixed 40×, 35×, and 25× independent downsamples of the normal alignment with 10×, 15×, and 25× downsamples of the pure tumor alignment to construct a 50× tumor mixture of 20%, 30%, and 50% purity, respectively. We restricted the training set to a 120-megabase region on chromosome 1 (with 39,818 true somatic SNVs and 38,804 true somatic INDELs) and the test set to the entire chromosome 22 (with 12,201 true somatic SNVs and 12,185 true somatic INDELs). The AF of somatic mutations across the three tumor mixture ratios ranged from 20% to 50%.

[0261] Candidate mutation preparation

[0262] In the first step, we scanned the tumor read alignments to find candidate positions with evidence of mutations. Many of these positions either have germline variants or are miscalled due to the complexity of the genomic region or sequencing artifacts. We applied a set of liberal filters to the set of candidate positions to ensure that the number of such positions was reasonable. In general, for SNVs, we required AF ≥ 0.03 or more than two reads supporting the SNV and a Phred-scaled base quality score greater than 19 (greater than 14 for real WES datasets) as a minimum requirement. For 1-base INDELs, we required AF ≥ 0.02 or more than one read supporting it. For INDELs greater than 1 base, we required AF ≥ 0.03. For ensemble methods, we also included any somatic mutations detected by other somatic mutation detection methods as input candidates. For the PacBio dataset, we used AF ≥ 0.1 for SNVs and INDELs greater than 1 base and AF ≥ 0.15 for 1-base INDELs.

[0263] For the DREAM challenge dataset, we excluded dbSNP from the input candidates. 36 For fair comparison, we also filtered the dbSNP call results of all other somatic mutation detection tools.

[0264] Input mutation matrix

[0265] For each candidate position, we prepared a 3D matrix M with k channels of size 5 × 32 (Figure 10a, Figure 10c, and Figure 10d). The five rows in each channel correspond to the four DNA bases A, C, G, and T and the space character (“-”). Each of the 32 columns of the matrix represents a column of the alignment.

[0266] For each candidate position, we extracted the tumor and normal read alignments. As shown in Figure 10a, we then took the read alignments of tumor and normal samples to the reference as MSA. To do this, we augmented the reference sequence by adding gaps to the reference sequence when there were insertions in the reads. It is important to note that this process does not require any further realignment of the original read alignments of the input BAM file, but simply reconstructs the alignment into MSA format by allocating additional columns at the locations where the insertions occurred. If there are multiple different insertions in multiple reads after a specific position, we treat them as left-aligned sequences and put them in the same set of columns (e.g., see the A base insertion and C base insertion in the ninth column of the virtual example in Figure 10a). Using this read representation, we find the frequency of A / C / G / T / - characters in each column and record separate matrices for tumor and normal samples (channels C2 and C3 in matrix M). In channel C1, we record the reference base (or gap) in each column. Channels Ci (4 ≤ i ≤ k) record other alignment signals in tumor and normal samples, such as coverage, base quality, mapping quality, strand, shear information, edit distance, alignment score, and paired end information. For example, for the base quality channel, each sample has a matrix of size 5 × 32 that records the average base quality of reads with a given base in each column (for a given row). As another example, for the edit distance channel, each sample has a matrix of size 5 × 32 that records the average edit distance of reads with a given base in each column (for a given row). One channel of the matrix M is dedicated to the column where the candidate is located. In the current implementation, we use a total of 26 channels in the NeuSomatic-S method in standalone mode.

[0267] For the ensemble extension of NeuSomatic, we also included additional channels to capture features reported by each of the six individual methods used. In this implementation, we used 93 additional channels to represent features extracted from other methods as well as alignments reported by SomaticSeq. Thus, the ensemble mode of NeuSomatic has 119 input channels for each candidate matrix.

[0268] For each candidate alignment, we report the alignment information in a window of seven bases around the candidate. We keep 32 columns to account for extended alignments with insertions. In a few cases, if the insertion is large, 32 columns may not be enough to represent the alignment. For such cases, we truncate the insertion so that at least three bases around the candidate can be recorded.

[0269] CNN Architecture

[0270] The proposed CNN (Figure 10b) consists of nine convolutional layers constructed as follows. The input matrix is ​​fed into the first convolutional layer with 64 output channels, 1 × 3 kernel size, and Relu activation, followed by batch normalization and max pooling layers. The output of this layer is then fed into a set of four blocks with shortcut identity connections similar to the ResNet structure. These blocks consist of a convolutional layer with 3 × 3 kernels followed by batch normalization and a convolutional layer with 5 × 5 kernels. Between these shortcut blocks, we use batch normalization layers and max pooling layers. The output of the final block is fed into a fully connected layer of size 240. The resulting feature vector is then fed into two softmax classifiers and a regressor. The first classifier is a 4-way classifier that predicts the mutation type from four categories: non-somatic, SNV, insertion, and deletion. The second classifier predicts the length of the predicted mutation based on four categories: 0, 1, 2, and ≥3. Non-somatic calls were annotated as zero-size mutations, SNVs and 1-base INDELs were annotated as size 1, and 2-base INDELs of size ≥ 3 were annotated as size 2 and ≥ 3 mutations, respectively. The regressor predicted the columns of the mutation in the matrix to ensure that the predictions targeted the correct positions, and was optimized using a smooth L1 loss function.

[0271] The CNN has less than 900K parameters, which enables us to achieve an efficient implementation by using a large batch size. The whole genome training process takes about 8 hours on a machine equipped with 8 Tesla K80 Nvidia GPUs.

[0272] CNN Training

[0273] For the DREAM Challenge dataset, the Platinum dataset, and the Target Enrichment dataset, we randomly split the genomic regions into 50% training and 50% test sets. For the PacBio dataset, we trained NeuSomatic on a 120-megabase region on chromosome 1 and tested it on chromosome 22.

[0274] For each dataset, we merged the generated training input matrices from all different tumor / normal purity scenarios and used the merged set to train the network. We then applied this unified trained model to testing under each independent tumor / normal purity setting.

[0275] The DREAM challenge dataset has 15,507 somatic mutations in stage 3 and 30,462 somatic mutations in stage 4. For better network training, we use BAMSurgeon 19 About 95 K more SNVs and about 95 K more INDELs with AF distribution similar to the original DREAM data were incorporated into stage 3 and stage 4 tumor samples.

[0276] We trained the network using a batch size of 1000 and the SGD optimizer with a learning rate of 0.01 and a momentum of 0.9, and multiplied the learning rate by 0.1 every 400 epochs.

[0277] In general, since there are many more non-somatic (reference or germline) calls than true somatic mutations for input candidate localizations, in each generation of training, we use all true somatic mutations in the training set and randomly select non-somatic candidates, which is twice the number of true somatic mutations. We use a weighted softmax classification loss function to balance the number of candidates in each class. For the DREAM challenge data, since we added more synthetic mutations in the training set, we increased the weight of the non-somatic class to achieve higher precision in the test set.

[0278] To assess synthetic target enrichment datasets, we used whole exome and whole genome data as training sets.

[0279] To test on real WGS samples CLL1 and COLO-829, we used SomaticSeq and NeuSomatic models trained on DREAM challenge phase 3. For real WES sample TCGA-AZ-6601, we used data from another TCGA WES dataset TCGA-AZ-4315. 30We prepared a training set based on the data from the synthetic tumor-normal WES dataset. We mixed tumor and normal alignments from this dataset and split the mixture into two equal alignments. We then used one alignment as a pure normal sample and spiked ~91 K random SNVs and ~9 K random INDELs into the other alignment using BAMSurgeon to generate a synthetic tumor sample for training. We used the model trained on this synthetic tumor-normal WES dataset to test NeuSomatic and SomaticSeq on the real WES dataset TCGA-AZ-6601. For experiments with 261 real TCGA samples, we used a similar approach to prepare a training set using 12 TCGA samples. The model trained on this synthetic dataset was used for testing on 261 TCGA samples.

[0280] Hyperparameter Tuning

[0281] For hyperparameter tuning, we used 10% of the genome in the DREAM challenge phase 3 experiments and used the derived parameters in all other experiments.

[0282] We further explored different network architectures, such as pre-activation ResNet architectures with 4 to 16 ResNet blocks (including ResNet-18 and ResNet-34 architectures) (Figures 24a to 24e), and some variants of the proposed residual NeuSomatic architecture (Figures 24f to 24m). To evaluate these networks, we divided the training data in the DREAM Phase 3 dataset into two equal parts, used one of them to train different architectures, and used the other to evaluate these different architectures in an independent mode. Figure 25 compares these architectures in terms of accuracy, number of network parameters, memory usage, and speed. In general, all of these networks can achieve relatively high accuracy compared to traditional somatic mutation detection methods. This finding reveals the importance of the proposed data summarization method, which captures the main signals in the genomic context of the candidates and contributes to the efficient implementation of convolutional networks on the somatic mutation detection problem. Compared with those with the proposed residual blocks (Figures 24f to 24m), the networks with two 3 × 3 convolutional filters ( Fig.24aThe default ReSNet architecture with strided convolutions (Fig. 24a to Fig. 24g) has lower average accuracy. In addition, the networks with strided convolutions (Fig. 24a to Fig. 24g) have more network parameters and runtime requirements. In summary, although each network architecture shows advantages in some comparison aspects, we choose the proposed NeuSomatic network architecture (Fig. 10b; Figure 24k ) as the default network architecture to take all these factors into account, while other networks can be easily adjusted based on user usage and time / computation constraints.

[0283] Figures 24a to 24m show the different network architectures tested. Figures 24a to 24e show ResNet architectures with default 3x3 convolutional layers with different numbers of pre-activation residual blocks. Here, strided convolutions are used with channel expansion. Figures 24f and 24g show multiple customized residual blocks with 3x3 and 5x5 convolutional layers and some dilated convolutions. Here, strided convolutions are used with channel expansion. Figure 24h shows four customized residual blocks with 3x3 and 5x5 convolutional layers and some dilated convolutions. Here, strided convolutions are not used. Figures 24i to 24m show NeuSomatic residual architectures with different residual blocks and fully connected sizes.

[0284] Figure 25 shows the performance analysis of the different network architectures shown in Figures 24a to 24m. Here, all networks are evaluated with a batch size of 1000 after 600 epochs of training.

[0285] Other somatic mutation detection algorithms

[0286] In our analysis, we used 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), and SomaticSeq (v2.7.0) somatic mutation detection algorithms with default settings.

[0287] We used VarDict as an alternative to NeuSomatic based on PacBio data. To be able to detect somatic mutations on long reads with high error rates, we used VarDict with the parameter setting of “ −m 10000 −Q 1 −q 5 −X 1−c 1 −S 2 −E 3 −g 4 −k 0 ”. In addition, as with NeuSomatic, AF ≥ 0.1 was used for SNVs and AF ≥ 0.15 for INDELs.

[0288] To train SomaticSeq, we also followed the same 50% training / testing region split used by NeuSomatic. Additionally, as with NeuSomatic, for each dataset, we combined the training data from all different tumor / normal purity scenarios to train the SomaticSeq SNV classifier and the INDEL classifier. These unified classifiers were then used to predict each individual tumor / normal purity setting.

[0289] For precision-recall analysis, somatic mutations were classified based on the confidence or quality score assigned by each tool. For MuSE, we used stratum assignment as the ranking criterion. For VarDict, VarScan2, MuTect2, Strelka2, and SomaticSniper, we used the SSF value, SSC value, TLOD value, SomaticEVS value, and SSC value reported in the VCF file for ranking, respectively. For SomaticSeq and NeuSomatic, we used the somatic mutation quality score from the QUAL field. NeuSomatic reports the quality score of the predicted somatic mutations based on the probability predicted by the CNN.

[0290] To analyze the performance on real samples, we used PASS somatic calls from different methods (for VarDict, we restricted ourselves to calls with StrongSomatic status). For NeuSomatic, we set the quality score threshold to 0.97 for WGS and 0.6 for WES.

[0291] Computational complexity

[0292] For whole genome data, scanning 30× tumor and normal alignments to find candidates, extract features, and prepare input matrices can take about 3.9 hours on a dual 14-core Intel Xeon CPU E5-2680 v4 2.40 Ghz machine. The whole genome training process can take about 8 hours (about 90 seconds per generation for a size of 580,000) on a machine equipped with 8 Tesla K80 Nvidia GPUs. Depending on the cutoff AF for candidate somatic mutations, network predictions for candidate mutations on 30× whole genome data can take about 35 minutes (AF cutoff 0.05, 3.9 M candidates) to about 100 minutes (AF cutoff 0.03, 11.5 M candidates) to compute using 8 Tesla K80 Nvidia GPUs. For 125× whole exome data, the entire scan, preparation, and computation of network predictions can take about 30 minutes. The end-to-end runtimes for predicting somatic mutations using the NeuSomatic ensemble approach and the standalone approach (CPU-only mode) on the 125× whole-exome dataset and the 30× whole-genome dataset were compared with other somatic mutation detection techniques in Figures 26a and 26b.

[0293] Figure 26a shows a runtime comparison of different somatic mutation detection algorithms. CPU core-hours are shown to predict somatic mutations on a 125x whole exome sequencing dataset. Figure 26b shows a runtime comparison of different somatic mutation detection algorithms. CPU core-hours are shown to predict somatic mutations on a 30x whole genome sequencing dataset.

[0294] Code availability

[0295] NeuSomatic is written in Python and C++. Its deep learning framework is implemented using PyTorch 0.3.1 to enable GPU applications for training / testing. Source code is available at https: / / github.com / bioinform / neusomatic under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. The results in this article are based on NeuSomatic v0.1.3.

[0296] Report Summary

[0297] For more information on the experimental design, see the Nature Research Reporting Summary linked to this article.

[0298] Data availability

[0299] Sequence data for this study were collected from various sources, namely, the European Nucleotide Archive (accession number: PRJEB3246; https: / / www.ebi.ac.uk / ena), the Sequence Read Archive (accession number: SRX1026041; https: / / www.ncbi.nlm.nih.gov / sra), the International Cancer Genome Consortium (project: ICGC-TCGA The DREAM mutation calling challenge (controlled access: https: / / icgc.org / ), The Cancer Genome Atlas (accession numbers: TCGA-AZ-6601, TCGA-AZ-4315; controlled access: https: / / gdc.cancer.gov / ), the European Genome-Phenome Archive (accession numbers: EGAS00000000052, EGAD00001000023; controlled access: https: / / www.ebi.ac.uk / ega / ), and Genomes in a Bottle (accession numbers: HG002, HG003, HG004; ftp: / / ftp-trace.ncbi.nlm.nih.gov / giab / ftp / ). Synthetic data were generated from the above datasets using scripts at https: / / github.com / bioinform / neusomatic / blob / paper / etc / data_scripts.zip. All other relevant data are available upon request.

[0300] References

[0301] 1.Cibulskis, K. et al., Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31,213 (2013).

[0302] 2. Fan, Y. et al., MuSE: accounting for tumor heterogeneity using a sample-specific error model improves sensitivity and specificity in mutation calling from sequencing data. Genome Biol. 17, 178 (2016).

[0303] 3. Lai, Z. et al., VarDict: a novel and versatile variant caller for next-generation sequencing in cancer research. Nucleic Acids Res. 44, e108–e108 (2016).

[0304] 4. Koboldt, D. C. et al., VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).

[0305] 5. Kim, S. et al., Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods 15, 591–594 (2018).

[0306] 6. Larson, D. E. et al., SomaticSniper: identification of somatic point mutations in whole genome sequencing data. Bioinformatics 28, 311–317 (2011).

[0307] 7. Wang, Q. et al., Detecting somatic point mutations in cancer genome sequencing data: a comparison of mutation callers. Genome Med. 5, 91 (2013).

[0308] 8. Alioto, T. S. et al., A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing. Nat. Commun. 6, 10001 (2015).

[0309] 9. Roberts, N. D. et al., A comparative analysis of algorithms for somatic SNV detection in cancer. Bioinformatics 29, 2223-2230 (2013).

[0310] 10. Fang, L. T. et al., An ensemble approach to accurately detect somatic mutations using SomaticSeq. Genome Biol. 16, 197 (2015).

[0311] 11. Poplin, R. et al., A universal SNP and small-indel variant caller using deep neural networks. Nat. Biotechnol. 36, 983 (2018).

[0312] 12. Luo, R., Sedlazeck, F. J., Lam, T.-W., & Schatz, M. Clairvoyante: A multi-task convolutional deep neural network for variant calling in single molecule sequencing. https: / / www.biorxiv.org / content / early / 2018 / 04 / 28 / 310458 (2018).

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

[0314] 14. Esteva, A., et al. Dermatologist-level classification of skin cancer with deep neural networks. Nature 542, 115–118 (2017).

[0315] 15. Torracinta, R., et al. Adaptive somatic mutations calls with deep learning and semi-simulated data. https: / / www.biorxiv.org / content / early / 2016 / 10 / 04 / 079087 (2016).

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

[0317] 17. He, K., Zhang, X., Ren, S., & Sun, J. Deep residual learning for image 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.

[0318] 18. Eberle, M. A., et al. A reference data set of 5.4 million phased human variants validated by genetic inheritance from sequencing a three-generation 17-member pedigree. Genome Res. 27, 157 - 164 (2017).

[0319] 19. Ewing, A. D., et al. Combining tumor genome simulation with crowdsourcing to benchmark somatic single-nucleotide-variant detection. Nat. Methods 12, 623 (2015).

[0320] 20. Zook, J. M., et al. Extensive sequencing of seven human genomes to characterize benchmark reference materials. Sci. Data 3, 160025 (2016).

[0321] 21. Xu, C. A review of somatic single nucleotide variant calling algorithms for next-generation sequencing data. Comput. Struct. Biotechnol. J. 16, 15–24 (2018).

[0322] 22. Zook, J. M. et al., Integrating human sequence data sets provides a resource of benchmark SNP and indel genotype calls. Nat. Biotechnol. 32, 246 (2014).

[0323] 23. Puente, X. S. et al., Whole-genome sequencing identifies recurrent mutations in chronic lymphocytic leukaemia. Nature 475, 101 (2011).

[0324] 24. Morse, H. G. and Moore, G. E., Cytogenetic homogeneity in eight independent sites in a case of malignant melanoma. Cancer Genet. Cytogenet. 69, 108–112 (1993).

[0325] 25. Pleasance, E. D. et al., A comprehensive catalogue of somatic mutations from a human cancer genome. Nature 463, 191 (2010).

[0326] 26. Network, C. G. A. et al., Comprehensive molecular characterization of human colon and rectal cancer. Nature 487, 330 (2012).

[0327] 27. Grasso, C. S. et al., Genetic mechanisms of immune evasion in colorectal cancer. Cancer Discov. 8, 730-749 (2018).

[0328] 28. Li, H., Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. https: / / arxiv.org / abs / 1303.3997 (2013).

[0329] 29. Van der Auwera, G. A. et al., From FastQ data to high-confidence variant calls: the genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinforma. 43, 10-11 (2013).

[0330] 30. Grossman, R. L. et al., Toward a shared vision for cancer genomic data. N. Engl. J. Med. 375, 1109–1112 (2016).

[0331] 31. Forbes, S. A. et al., COSMIC: somatic cancer genetics at high-resolution. Nucleic Acids Res. 45, D777–D783 (2016).

[0332] 32. Chaisson, M. J. et al., Resolving the complexity of the human genome using single-molecule sequencing. Nature 517, 608 (2015).

[0333] 33. Mu, J. C. et al., VarSim: a high-fidelity simulation and validation framework for high-throughput genome sequencing with cancer applications. Bioinformatics 31, 1469-1471 (2014).

[0334] 34. Lau, B. et al., LongISLND: in silico sequencing of lengthy and noisy data types. Bioinformatics 32, 3829-3832 (2016).

[0335] 35. Sedlazeck, FJ et al. Accurate detection of complex structural variations using single-molecule sequencing. Nat. Methods 15, 461-468 (2018).

[0336] 36. Sherry, ST et al., dbSNP: the NCBI database of genetic variation. Nucleic Acids Res. 29, 308–311 (2001).

[0337] Example 2

[0338] Accurate detection of somatic mutations is challenging but critical for understanding cancer development, progression, and treatment. We recently proposed NeuSomatic, the first deep convolutional neural network-based somatic mutation detection method, and demonstrated superior performance on in silico data. In this study, we used the first comprehensive and well-characterized somatic reference samples from the SEQC-II consortium to explore best practices for leveraging deep learning frameworks in cancer mutation detection. Using the high-confidence somatic mutations established by the consortium for these reference samples, we identified strategies for building robust models on multiple datasets derived from samples representing real-world scenarios. The proposed strategy achieved high robustness across multiple sequencing technologies such as WGS, WES, AmpliSeq targeted sequencing for fresh and FFPE DNA inputs, different tumor / normal purities, and different coverages ranging from 10× to 2000×. NeuSomatic significantly outperformed conventional detection methods both in general and in challenging cases such as low coverage, low mutation frequency, DNA damage, and difficult genomic regions.

[0339] Preface

[0340] Somatic mutations are key cancer drivers. Accurate somatic mutation detection enables precise diagnosis, prognosis, and treatment of cancer patients. 1Several tools have been developed to identify somatic mutations from next-generation sequencing technologies. 2-11 In general, traditional techniques use various statistical / algorithmic modeling methods to obtain tumor-normal paired sequencing data. 2-10 A set of handcrafted sequencing features extracted from a sample accurately pinpoints somatic mutations from background noise, germline variants, and / or cross-contamination. These manually designed methods are not easily generalizable beyond the specific cancer type, sample type, or sequencing strategy developed for it.

[0341] Recently, a deep learning-based somatic mutation detection method called NeuSomatic has been proposed, which uses a convolutional neural network (CNN) to learn feature representations directly from raw data. 11 NeuSomatic uses a novel summary of tumor / normal alignment information as a set of input matrices that can be employed to efficiently train models to learn how to effectively distinguish true somatic mutations from artifacts. Network models trained by NeuSomatic can capture important mutational signals directly from read alignments and genomic contexts without manual intervention. Thus, NeuSomatic provides a framework that can be easily applied to different problem statements, including various sequencing technologies, cancer types, tumor and normal purity, and mutant allele frequencies, by training on appropriate data. Furthermore, it can be implemented as a standalone somatic mutation detection method or combined with existing methods to achieve the highest accuracy. NeuSomatic has been demonstrated to significantly outperform conventional techniques based on in silico datasets. Performance evaluation on real samples has been limited to partial sensitivity analysis of a small number of validated somatic variants, given the lack of thoroughly characterized benchmark samples with known “ground truth” somatic mutations. Therefore, despite the advantages seen in implementing NeuSomatic’s CNN-based framework on in silico data, the accuracy and reproducibility of NeuSomatic on real cancer samples has not been fully evaluated to date.

[0342] Recently, the Somatic Mutation Working Group of the Sequencing Quality Control Phase II (SEQC-II) consortium, led by the FDA, developed reference matched tumor-normal samples: a human triple-negative breast cancer cell line (HCC1395) and a matched B lymphocyte-derived normal cell line (HCC1395BL). 12,13 Using orthogonal sequencing technology, multiple sequencing replicates, and multiple bioinformatics analysis pipelines, the SEQC-II consortium has developed a “gold set” of well-defined somatic single nucleotide variants (SNVs) and insertions / deletions (INDELs) for HCC1395.

[0343] As the first comprehensive and well-characterized paired tumor-normal reference cancer samples, this dataset and accompanying sequencing data prepared from multiple loci and technologies provide a unique resource to achieve two important goals. First, using the set of high-confidence somatic mutation calls developed by the consortium for these reference samples, we provide in-depth analysis of deep learning-based somatic mutation detection in real cancer samples compared to traditional approaches. Second, we explored various model building strategies for training CNNs in NeuSomatic using SEQC-II data and identified effective training methods based on multiple datasets derived from samples representing real-world scenarios. We evaluated the proposed strategy on whole genome sequencing (WGS), whole exome sequencing (WES), and AmpliSeq targeted sequencing datasets, ranging from 10× to 2000× coverage. The WGS and WES data were derived from formalin-fixed paraffin-embedded (FFPE) and fresh DNA, using three library preparation protocols with different input amounts, and from multiple platforms / sequencing sites. Using tumor and normal titration, we evaluated different methods for 5% to 100% tumor purity, 5% contamination matched normal samples, and 10× to 300× WGS coverage. The proposed strategy for training and implementing the deep learning framework in NeuSomatic achieved high robustness in all of the above real-world scenarios and significantly outperformed traditional paired tumor-normal somatic mutation detection methods.

[0344] Our analysis of SEQC-II reference cancer samples demonstrates that the deep learning scheme implemented in NeuSomatic helps overcome major challenges in somatic mutation detection that are difficult to address using traditional techniques. Therefore, the deep learning models and strategies derived from our study can provide the research community with actionable best practice recommendations for reliable cancer mutation detection.

[0345] result

[0346] Reference samples and datasets

[0347] To fully characterize the somatic mutation detection problem, we used the first comprehensive genome-wide characterized reference tumor-normal paired breast cancer cell lines (HCC1395 and HCC1395BL) developed by the Somatic Mutation Working Group of the SEQC-II consortium. 12,13. We utilized the high-confidence somatic mutations derived by the consortium (39,536 SNVs and 2,020 INDELs) as our baseline truth set, as shown in Figure 27, which shows the VAF distribution of the baseline truth SNV and INDEL somatic mutations in the HCC1395 supercall result set classified into four confidence levels (high, medium, low, and unclassified) by the SEQC-II consortium. High-confidence and medium-confidence calls were grouped as a set of “true values” for somatic mutations. To extensively assess the consistency and reproducibility of the predictions, we used 123 replicates from a variety of datasets representing real-world cancer detection applications, including real WGS, WES, and AmpliSeq target sequencing with different coverage, tumor purity, and library preparation, which were applied to FFPE and fresh DNA sequenced on multiple platforms at six centers (see “Methods”).

[0348] Figure 27 shows the VAF distribution of ground truth SNVs and INDEL somatic mutations in a superset of calls for HCC1395 classified into four confidence levels (high, medium, low, and unclassified) by the SEQC-II consortium. High and medium confidence calls are grouped as a set of “ground truth” somatic mutations.

[0349] Analysis Overview

[0350] We evaluated NeuSomatic's deep learning-based somatic mutation detection approach using various trained network models and compared it with seven widely used traditional somatic mutation calling algorithms, including MuTect2 2 、MuSE 3 、SomaticSniper 4 、Strelka2 5 、VarDict 6 TNscope 7 and Lancet 8 We evaluated NeuSomatic in its standalone mode (shown as NeuSomatic-S) and in its integrated mode, which includes, in addition to the raw data, predictions reported by MuTect2, SomaticSniper, VarDict, MuSE, and Strelka2 as input channels.

[0351] Figures 28a to 28e show the overall performance of NeuSomatic on 123 replicates in the SEQC-II dataset. Figure 28a shows that the NeuSomatic model trained on the DREAM Challenge Phase 3 dataset outperforms other techniques when applied to the SEQC-II data. Figure 28b shows that training different NeuSomatic models using the SEQC-II reference sample resulted in an additional 4% to 5% absolute improvement in the average F1 score. Figure 28c shows that NeuSomatic achieves consistent advantages over other techniques across various sets of replicates of different purity / coverage in WGS, WES, FFPE, AmpliSeq, and different library preparation datasets using models trained on SEQC-II data. In this sub-figure, for each replicate, the best F1 score is calculated by different methods. The heatmap shows the absolute difference between the F1 score of any one of the somatic mutation detection methods and the best F1 score. In each figure, the average F1 score across 123 replicates for each method is shown. Figure 28d compares the performance of different techniques across 123 replicates in six datasets for SNVs and INDELs. For each replicate, the best F1 score was calculated by different methods. The heat map shows the absolute difference between the F1 score of any method and the best F1 score among the somatic mutation detection methods. The average F1 score across 123 replicates is shown for each method.

[0352] In our analysis, we used several different trained models. First, we used a recently released model that is already available. 11 , the model is based on the DREAM Challenge Stage 3 dataset 13 Despite the huge differences in sample types, sequencing platforms, coverage, spike-in mutation frequencies, and sample heterogeneity used to train the DREAM3 model, the model outperformed other traditional techniques on real cancer datasets of various features, with the average F1 score of the model across different samples being about 4% higher than other traditional techniques for both SNVs and INDELs (Figure 28a). While this superiority supports the robustness of NeuSomatic to the assigned variants, it also suggests that the deep learning framework used in NeuSomatic can perform even better by learning sequencing features and mutational signatures from real cancer samples, especially for predicting INDELs and for PCR-enriched datasets with higher coverage such as AmpliSeq.

[0353] To identify the most effective strategy for building network models using real cancer samples, we evaluated ten additional training methods of NeuSomatic using the SEQC-II reference samples (Figure 28e).

[0354] The first model (SEQC-WGS-Spike-In) was trained on a set of WGS pairs of in silico tumor-normal replicates, where the in silico tumors were designed by spiking mutations into different normal replicates. The second model (SEQC-WGS-GT-50) was trained on a set of real WGS tumor-normal replicate pairs on 50% of the genome using the ground truth somatic mutations in HCC1395. The third model (SEQC-WGS-GT50-Spike-In WGS10) was prepared by adding 10% of the training data from the first model to the training data used for the second model to take advantage of the large number of spike-in mutations and real somatic mutations. All three models were tested on all datasets. For specific datasets such as FFPE and WES, we also prepared six additional specialized models using a set of synthetic pairs of synthetic in silico tumor and normal replicates. For all models, we evaluated the performance on the 50% conserved regions of the genome (not used for the SEQC-WGS-GT-50 model). We also trained a model (SEQC-WGS-GT-ALL) using all ground truth mutations similar to the SEQC-WGS-GT-50 model but based on the whole genome. SEQC-WGS-GT-ALL is not directly applicable to WGS datasets where all true mutations have been seen, and it would be beneficial to analyze its performance on other datasets or samples besides HCC1395.

[0355] Using models trained on SEQC-II samples, we improved the average performance of the DREAM3 model by an additional improvement of about 4% to 5% in the average F1 score (Figure 28b). The proposed model building strategy is consistently the best across a variety of sample types and sequencing strategies, outperforming traditional methods by more than 5.7% and 7.8% in average F1 scores for SNVs and INDELs, respectively. Similarly, across all samples, we observed a median F1 score improvement of more than 5.6% over other traditional techniques (Figures 28c and 28d).

[0356] WGS datasets

[0357] We evaluated the performance of the above somatic cell calling techniques and network models on 21 WGS replicates sequenced at six sequencing centers using HiSeqX10, HiSeq4000, and NovaSeq S6000 platforms (Figure 29a to Figure 29d).

[0358] The NeuSomatic SEQC-WGS-GT50-incorporated WGS10 model consistently outperformed the other options, with slight variations between replicates, demonstrating robustness and reproducibility (Figure 29a). NeuSomatic achieved an average F1 score of 94.6% for SNVs and 87.9% for INDELs, 5.6% and 10.2% higher than the average F1 scores of other traditional somatic mutation detection options for SNVs and INDELs, respectively. Precision-recall analysis showed that this advantage was driven by NeuSomatic’s high precision compared to the other techniques (Figure 29b). Comparing different model training strategies also revealed that NeuSomatic-S INDEL calls benefited more from training with ground truth somatic mutations, and in general, we observed significant improvements of up to 11% using the SEQC-II reference sample compared to the DREAM3 model (Figure 29c).

[0359] Tumor purity and contaminated normal samples

[0360] Since tumor purity and contamination of tumor cells in matched normal samples can greatly affect the accuracy of mutation detection, we explored the robustness of the model at different sequencing depths and sample purities by tumor-normal titration of WGS samples. We first studied tumor samples with purities ranging from 5% to 100% paired with pure normal samples for coverages ranging from 10× to 300× (Figure 30a to Figure 30h). In general, NeuSomatic maintained its advantage over other schemes despite changes in tumor purity and coverage, reflecting its robustness (Figure 30a, Figure 30e, and Figure 30f). It had the greatest advantage over traditional schemes in more challenging situations such as lower coverage (e.g., F1 score of about 20% for samples with 10× coverage and 100% purity).

[0361] We further analyzed the robustness of normal samples to 5% tumor contamination for tumor sample purity ranging from 10% to 100% at 80x coverage (Figure 30b). NeuSomatic was highly robust to tumor-normal cross-contamination, with a median absolute change in F1 score of less than 5%. Among other techniques with high F1 scores in pure normal samples, Strelka2 also showed high robustness to tumor contamination (median change in F1 score of 8.4%). MuTect2, MuSE, Lancet, and Tnscope, despite having high F1 scores in pure normal sample scenarios, showed a drop in F1 scores of about 50% when using contaminated normal samples.

[0362] Models trained on ground truth had the highest advantage over the DREAM3 model, generally due primarily to higher precision (Figure 30c and Figure 30g). SEQC-WGS-spiking models trained on in silico tumors also showed lower precision compared to models trained on ground truth. In general, INDEL, low purity, and low coverage samples benefited the most from training on SEQC-II data.

[0363] Analysis of the F1 scores for different INDEL sizes across various coverage and purity settings demonstrated the robustness of NeuSomatic to size variations for both insertions and deletions ( FIG. 30 d ).

[0364] Library preparation and DNA input

[0365] To measure the impact of library preparation on prediction robustness, we used our model to test six replicates prepared using TruSeq-Nano and Nextera Flex protocols and three DNA input amounts: 1 ng, 10 ng, and 100 ng (Figure 31a-Figure 31f and Figure 31h). Across different library preparation methods, NeuSomatic consistently outperformed other technologies. For the 1 ng TruSeq-Nano library, all methods performed poorly due to limited effective coverage (~7x) after removing redundant reads. On average, NeuSomatic improved the F1 score for SNVs by 8.4% compared to traditional somatic mutation detection technologies. For INDELs, Lancet's assembly-based INDEL caller outperformed the NeuSomatic SEQC-WGS-GT50-incorporated WGS10 model by ~4% (Figure 31a). In contrast, NeuSomatic’s SEQC-WGS-GT50 model performs similarly to Lancet on INDELs (Figures 31b and 31g). Training on SEQC-II spike-in or ground truth data leads to an overall improvement of ~8.4% in F1 scores for SNVs compared to the DREAM3 model. The advantage is more pronounced in more challenging cases, including TruSeq-Nano libraries with 1 ng input. We also observed that NeuSomatic-S benefits more from these models.

[0366] Capture (WES) and targeted (AmpliSeq) panels

[0367] We tested our models using 16 WES replicates sequenced at six sequencing sites and three replicates from the AmpliSeq dataset (Figures 31c-31f, 31i, and 31j). Despite the different coverage and platform bias of the models trained on SEQC-II WGS samples, NeuSomatic performed well on both WES and AmpliSeq datasets at 2000x coverage. For WES, NeuSomatic achieved an average SNVF1 score of 95.4%, a 2.6% improvement over the average F1 scores of the alternatives. Models trained on WES and WGS data performed similarly on the WES dataset, with F1 scores of approximately 95%. On the AmpliSeq dataset, the average F1 scores of NeuSomatic SEQC-WG-GT50 and SEQC-WES-Spike models were greater than 90% on average, the highest compared to other models / schemes and Strelka2.

[0368] Effects of FFPE Processing

[0369] To measure the robustness of NeuSomatic predictions on FFPE-processed samples, we used 8 WGS and 7 WES FFPE replicates prepared with four different formaldehyde fixation times: 1, 2, 6, and 24 hours. We evaluated each FFPE replicate with both a FFPE normal matched sample and a fresh normal matched sample. Despite the presence of FFPE artifacts, NeuSomatic consistently outperformed the other techniques, and its performance remained largely constant across different fixation times and matched normal samples used (Fig. 32a, Fig.32e Based on WGS FFPE data, NeuSomatic achieved an average F1 score of 86.1% and 76.9% for SNV and INDEL, respectively, which was 4% and 6% better than the average F1 score of alternative technologies. Similarly, for FFPE WES data, NeuSomatic achieved an average F1 score of 78.9%, which was more than 4% higher than the average F1 score of the traditional method (Figure 32c).

[0370] In general, we observed significant improvements over the DREAM3 model when training with SEQC-II samples (Figure 32b and Figure 32d). Only for INDEL prediction using NeuSomatic did the model trained on FFPE samples appear to improve over the model trained on fresh samples, but for SNVs, the model trained on fresh samples was superior.

[0371] Sample-specific models

[0372] Although generic models trained on SEQC-II have shown consistently better performance than other traditional somatic mutation detection schemes, here we explored whether sample-specific trained models could give us additional accuracy gains using nine replicates across different SEQC-II datasets. For each sample, we trained NeuSomatic and NeuSomatic-S models using in silico tumor-normal replicates prepared for that sample. We also trained a distinct model by combining 10% of the training candidates used for the SEQC-WGS-spiking model with the training data derived for each sample. We compared these two sample-specific models with the generic SEQC-WGS-spiking model (Figures 33a to 33c). On average, the sample-specific models improved their SNV and INDEL F1 scores by approximately 0.5% and approximately 5%, respectively, compared to the SEQC-WGS-spiking model used for NeuSomatic. Library preparation and FFPE samples benefited the most from sample-specific training. For example, the library preparation samples with 1 ng DNA amount prepared with the Nextera Flex protocol achieved 1.6% and 19.4% improvement in absolute F1 scores for SNVs and INDELs, respectively. Similarly, the F1 scores of 24-hour duration FFPE WGS samples with matched fresh normal samples improved by 1.8 and 14.8 percentage points for SNVs and INDELs, respectively. For NeuSomatic-S, only INDELs seemed to benefit from sample-specific training.

[0373] INDEL Performance

[0374] We further evaluated the accuracy of detecting INDELs of different sizes across multiple datasets (Figures 34a and Fig.34b). Although NeuSomatic does not explicitly use native assemblies for INDEL detection, it still consistently outperforms other methods, including assembly-based techniques such as Lancet, across a wide range of INDEL sizes, coverage, and tumor purity (Figure 34a and Figure 30d). Overall, for insertions and deletions larger than 2 base pairs (bps), NeuSomatic achieved 24.4% and 6.5% higher F1 scores than Lancet, which was the best alternative traditional technique. Comparing INDEL accuracy across different NeuSomatic training methods, the DREAM3 model was primarily plagued by low insertion detection accuracy, while the SEQC-II trained model reliably identified both insertions and deletions of different sizes (Figure 34b).

[0375] Performance analysis of different variant allele frequencies (VAFs)

[0376] We analyzed the accuracy of different methods for VAF across somatic mutations (Fig. 34c and Fig. 35 to Fig.37c ). Across different datasets, we observed that NeuSomatic showed an advantage over other schemes in cases where VAF was less than 20% (Figure 34c and Figure 35). We also observed that NeuSomatic was highly robust to changes in VAF, with consistent predictions for VAFs as low as 5%. For mutations with VAFs between 5% and 20%, the SEQC-II trained model clearly performed better than the DREAM3 model (Figure 36a and Figure 36b). Similar robustness was observed for all coverage and tumor purity settings (Figure 37a to Figure 37c).

[0377] Performance on difficult genomic regions

[0378] To illustrate the robustness to the complexity of genomic context, we evaluated the performance of different somatic calling schemes on difficult genomic regions, including tandem repeats (TRs) of varying lengths and segmental duplications (Figure 34d and Figure 38 to Figure 41c). While many other schemes showed significant performance degradation on difficult genomic regions (Figure 38), NeuSomatic remained highly robust, with even more significant improvements over other schemes compared to whole-genome analysis (Figure 34d and Figure 38). On average, NeuSomatic's F1 scores for both SNVs and INDELs for these difficult genomic regions were more than 15% higher than the best alternative traditional somatic caller. Training NeuSomatic with SEQC-II samples resulted in substantial improvements in accuracy on these difficult genomic regions, with the largest gains being for tandem repeat regions using models trained on ground truth mutations (Figure 39a to Figure 40b). Similar performance on difficult genomic regions was observed at different coverages and tumor purities (Figure 41a to Figure 41c).

[0379] discuss

[0380] We explored the robustness of NeuSomatic's deep learning framework in detecting somatic mutations across a diverse set of experimental settings seen in real WFP, WES, and AmpliSeq samples with varying coverage, sample purity, and library preparation, performed on FFPE and fresh DNA sequenced at six centers using multiple platforms. The experiments demonstrated NeuSomatic's potential to capture true somatic mutation signatures from raw data to distinguish true calls from sequencing artifacts. NeuSomatic models trained on SEQC-II reference samples, using both spiked-in mutations and a set of benchmark truth somatic mutations derived from the SEQC-II study, improved performance to achieve the highest accuracy. The analysis highlights best practice strategies for building models that can be used in a variety of scenarios. Compared to the baseline DREAM3 network model, the proposed model derived from the SEQC-II reference samples was shown to reduce false positives for both SNVs and INDELs by learning the correct mutation signatures, improve insertion calling accuracy, and enhance somatic detection of difficult genomic regions.

[0381] While NeuSomatic remains robust to tumor contamination in matched normal samples, many somatic calling methods such as MuTect2, Lancet, MuSE, and Tnscope suffer from this, as they reject true somatic mutations with corresponding reads in normal samples. As a result, their recall rates drop significantly. NeuSomatic models trained on WGS perform equally well on targeted sequencing data. On the other hand, TNscope employs a machine learning (random forest) model, while Lancet uses a locally colored de Bruijn graph approach to detect mutations. Both tools are designed for whole-genome somatic mutation detection and are therefore not suitable for targeted sequencing (which typically covers genomes less than 1 MB).

[0382] Analyzing the prediction accuracy across different specifications including mutation type, INDEL length, VAF, and genomic region, the results showed that NeuSomatic had the greatest improvement over other protocols in challenging cases such as low coverage, VAF less than 20%, difficult genomic regions such as tandem repeats, INDELs longer than 2 bp, FFPE with DNA damage, or contaminated matched normal samples. As a result of confusing sequencing or sample-specific artifacts with true somatic mutations, traditional protocols have many false positives (and therefore low precision), which are often seen in low-complexity genomic regions (Figure 34d and Figure 38). However, by training on WGS samples from multiple platforms and multiple centers, NeuSomatic learned error-prone genomic contexts, thereby consistently improving accuracy across different conditions, including difficult-to-call low-complexity genomic contexts. Similarly, analysis of missed calls by different methods showed that private false-negative calls generated by traditional protocols such as Strelka2, which were correctly predicted by NeuSomatic, had low VAFs (Figure 42). This demonstrates the disadvantage of other schemes in detecting low VAF mutations. On the other hand, by learning the correct landmarks from the raw data, NeuSomatic can more accurately distinguish challenging low VAF calls from artifacts.

[0383] Figure 42 shows a violin plot comparison of the VAF distribution of private FN calling results based on WGS datasets. In each subplot, we compared a traditional somatic mutation detection scheme with NeuSomatic. For the X-vs-Y subplot, we identified FN, which is the set of ground truth SNVs missed by algorithm X (FN in X) in at least 11 of the 21 WGS replicates. Similarly, we identified Y FN , which is the set of ground truth SNVs missed by algorithm Y (FN in Y) in at least 11 of the 21 WGS replicates. The figure then shows the VAF distribution of the private FN call results for X and Y. In other words, the violin plot shows in blue the VAF distribution of X. FN / Y FN VAF distribution of call results in the set, with Y shown in red FN / X FN VAF distribution of call results in the set. For most traditional schemes such as Strelka2, private FNs correctly predicted by NeuSomatic have lower VAFs, which indicates the disadvantage of traditional schemes in detecting low-VAF mutations. In each X-vs-Y subplot, the number of private FNs of X and Y is reported at the top.

[0384] method

[0385] SEQC-II Tumor-Normal Sequencing Data and Ground Truth

[0386] In our analysis, we used SEQC-II reference matched samples, a human triple-negative breast cancer cell line (HCC1395), and a matched B lymphocyte-derived normal cell line (HCC1395BL). Detailed sample information is available in the SEQC-II reference sample manuscript. 12、13 The SEQC-II Somatic Mutation Working Group has established a gold standard truth set of somatic mutations for these samples. 13 (Figure 27). The truth set was defined using multiple tumor-normal sequencing replicates from different sequencing centers and orthogonal mutation detection bioinformatics pipelines. In order to assign confidence scores to mutation call results and minimize platform, center, and pipeline specific biases, SomaticSeq was implemented 10The machine learning framework is based on a set of classifiers trained on computer-simulated tumor-normal copy pairs, which are generated by incorporating synthetic mutations into HCC1395BL alignment files, and each copy pair is matched to another HCC1395BL. Using these classifiers, mutation call results are divided into four confidence levels (high confidence, medium confidence, low confidence, and unclassified) based on the reproducibility of cross-aligners and cross-sequencing centers. High-confidence and medium-confidence call results are classified into a set of "truth sets" (v1.0) of somatic mutations, which contain a total of 39,536 SNVs and 2,020 INDELs. The truth set of somatic mutations in HCC1395 is available on the NCBI ftp website (ftp: / / ftp-trace.ncbi.nlm.nih.gov / seqc / ftp / Somatic_Mutation_WG / ) and is publicly available.

[0387] All sequencing data used in this study are publicly available through NCBI's SRA database (SRP162370). For all samples, FastQ files were first trimmed using Trimmomatic15 and then aligned with BWA-MEM (v0.7.15)16 and then with Picard MarkDuplicates (https: / / broadinstitute.github.io / picard).

[0388] Training dataset and model

[0389] Different training datasets were used to derive effective NeuSomatic CNN models, using both in silico tumor replicas prepared with synthetic spike-in mutations and real tumor replicas with a known high-confidence ground truth set of somatic mutations (Figure 28e):

[0390] DREAM3 Model

[0391] As a baseline WGS model, we adopt the DREAM3 model, which was recently trained on the ICGC-TCGA DREAM challenge phase 3 data. 14 Developed by 11The Phase 3 dataset consists of normal and tumor samples constructed by computationally incorporating 7903 SNVs and 7604 INDEL mutations into the healthy genome of the same normal sample with three different AFs (50%, 33%, and 20%) to create synthetic but realistic tumor-normal sample pairs. The BAMSurgeon with a similar AF distribution as the original DREAM data was used. 17 , an additional ~95K SNVs and ~95K INDELs were incorporated into the tumor samples 11 , for better network training. The model was trained by combining training data (in 50% of the genome) from five different tumor-normal purity settings (100T:100N, 50T:100N, 70T:95N, 50T:95N, and 25T:95N11). NeuSomatic and NeuSomatic-S were trained on ~29.2M candidate mutations identified in these five replicate pairs, including ~450K candidates with somatic SNV / INDEL markers and ~28.7M candidates marked as non-somatic.

[0392] TCGA Model

[0393] As a baseline WES model, we used a recently developed 18 TCGA model developed for training 11 . The tumor and normal alignments for each of these samples were mixed and split into two equal alignments. One alignment was treated as the normal replicate, and mutations were incorporated into the other alignment to construct a tumor replicate. For each sample, approximately 88K SNVs and approximately 44K INDELs were incorporated to generate synthetic tumor samples for training. NeuSomatic and NeuSomatic-S were trained based on approximately 5.9M candidate mutations identified in these 12 replicate pairs, including approximately 1.5M candidates with somatic SNV / INDEL markers and approximately 4.4M candidates marked as non-somatic.

[0394] SEQC-WGS-Incorporation Model

[0395] To prepare training data for building the model, we used BAMSurgeon to construct synthetic tumors by spiking in silico mutations into one HCC1395BL replicate and pairing it with a unique HCC1395BL replicate as a normal match. Using this approach, we generated ten in silico tumor-normal sample pairs. Eight of the ten replicate pairs were selected from WGS replicates that were sequenced at four different sites with an average coverage of 40× to 95×. The other two replicate pairs were obtained by merging multiple NovaSeq sequencing replicates from Illumina, resulting in in silico simulations with ~220× tumor sample coverage and ~170× normal matched sample coverage. In each in silico tumor, we spiked ~92 K SNVs and ~22 K INDELs. AFs into which mutations were spiked were randomly selected from a beta distribution (shape parameters α = 2 and β = 5). For each of the ten replicate pairs, we also constructed an impure normal sample by mixing 95% normal reads and 5% tumor reads. Therefore, we used a total of 20 in silico tumor-normal pairs to train the SEQC-WGS-spiking model. We then trained NeuSomatic and NeuSomatic-S based on the ~272M candidate mutations identified in these 20 replicate pairs, including ~2M candidates with somatic SNV / INDEL markers and ~270M candidates marked as non-somatic.

[0396] SEQC-WGS-GT-50 Model

[0397] The model was constructed using real WGS tumor-normal replicates and the SEQC-II high-confidence truth set as the baseline truth somatic mutations. We used eight WGS tumor-normal replicates as the base sample to train the model. The first seven WGS replicate pairs were from six different sequencing centers on HiSeq and NovaSeq platforms with an average coverage of 40× to 95×, and the last one was constructed by combining 9 NovaSeq sequencing replicates from Illumina, resulting in a replicate pair with approximately 390× coverage. For each of these eight replicate pairs, we constructed two additional purity variants, one with 95% normal purity obtained by mixing 95% normal reads and 5% tumor reads, and the other with 10% tumor purity obtained by mixing 10% tumor reads and 90% normal reads. Therefore, for each replicate pair, we have a version with 100% tumor and normal, a version with 100% pure tumor and 95% pure normal matching, and a version with 10% pure tumor and 100% pure normal matching. Therefore, we used a total of 24 tumor-normal replicates to train the SEQC-WGS-GT-50 model. For fair evaluation, we used only 50% of the genome to train the model, while the remaining 50% was used for evaluation. To prepare the training and evaluation regions, we divided the genome into small regions of size about 92K bps and randomly selected half of the break regions for training and the other half for evaluation. We excluded the 5-base padding region around each gray zone mutation from the training region, the excluded regions included the calls with low confidence labels in the superset of SEQC-II calls, and those with unclassified labels with a VAF of more than 30% in the superset of SEQC-II calls. We then trained NeuSomatic and NeuSomatic-S based on ~137M candidate mutations identified in these 24 replicate pairs, including ~416K candidates with somatic SNV / INDEL markers and ~136.5M candidates marked as non-somatic.

[0398] SEQC-WGS-GT-ALL model

[0399] This model was prepared similarly to the SEQC-WGS-GT-50 model, but uses all high-confidence ground-truth somatic mutations on the whole genome. Therefore, it cannot be used to evaluate WGS samples used for training purposes, and given that all samples were used for training, it may be biased towards the truth set mutations. The model is still beneficial for performance analysis on other datasets or non-SEQC-II samples. We trained models for NeuSomatic and NeuSomatic-S based on ~272M candidate mutations identified in these 24 replicate pairs used for the SEQC-WGS-GT-50 model, including ~839K candidates with somatic SNV / INDEL markers and ~271M candidates marked as non-somatic.

[0400] SEQC-WGS-GT50-incorporated WGS10 model

[0401] The training data for the model was created by merging 10% of the training candidates used in the SEQC-WGS-spiked-in model with all candidates in SEQC-WGS-GT-50. This combination exploits the large number of spike-in mutations and the signature of true somatic mutations seen in real tumor-normal data. We trained NeuSomatic and NeuSomatic-S on a combined set of 164M candidate mutations, including ~574K candidates with somatic SNV / INDEL labels and ~163.5M candidates labeled as non-somatic.

[0402] SEQC-WES-Incorporation Model

[0403] Similar to the SEQC-WGS-incorporation model, we constructed 7 in silico tumor-normal WES replicates using BAMSurgeon to train the model. Replicate pairs were selected from WES datasets that were sequenced at four different loci with average coverage ranging from 60× to 550×. In each in silico tumor, we spiked ~97 K SNVs and ~19 K INDELs. AFs for incorporation mutations were randomly selected from a beta distribution with shape parameters α = 2 and β = 5. We then trained NeuSomatic and NeuSomatic-S on ~3.7M candidate mutations identified in these 7 replicate pairs, including ~755K candidates with somatic SNV / INDEL markers and ~3M candidates labeled as non-somatic.

[0404] SEQC-FFPE-spiking model

[0405] Similar to the SEQC-WGS-spiking model, we constructed 8 in silico tumor-normal WGSFFPE replicates using BAMSurgeon to train the model. Replicate pairs were selected from FFPE datasets sequenced at two different sites with four different preparation times. In each in silico tumor, we spiked ~92 K SNVs and ~22 K INDELs. AFs with spiked mutations were randomly selected from a beta distribution with shape parameters α = 2 and β = 5. In our training, we also matched each in silico tumor with a fresh WGS replicate to include FFPE tumor with fresh normal scenarios. Therefore, we used a total of 16 in silico tumor-normal pairs to train the SEQC-FFPE-spiking model. We trained NeuSomatic and NeuSomatic-S on ~191M candidate mutations identified in these 7 replicate pairs, including ~1.7M candidates with somatic SNV / INDEL markers and ~190M candidates marked as non-somatic.

[0406] SEQC-FFPE-WES-spiking model

[0407] Similar to other spike-in models, we constructed 7 in silico tumor-normal WES FFPE replicates using BAMSurgeon to train this model. Replicate pairs were selected from WES FFPE datasets sequenced at two different sites and prepared at four different time intervals. Since we did not include two normal replicates with the same FFPE preparation time and sequencing sites in this dataset, we mixed the tumor and normal alignments for each replicate pair and split the mixture into two equal alignments. We then treated one alignment as a normal replicate and spiked mutations into the other alignment to construct the tumor replicate. In each in silico tumor, we spiked ~90 K SNVs and ~19 K INDELs. AFs for spiked mutations were randomly selected from a beta distribution (shape parameters α = 2 and β = 5). In our training, we also matched each in silico tumor with a fresh WES replicate to include FFPE tumor vs. fresh normal scenarios. Therefore, we used a total of 14 in silico tumor-normal pairs to train the SEQC-FFPE-WES-spiking model. We trained NeuSomatic and NeuSomatic-S on approximately 9.6M candidate mutations identified in these 7 replicate pairs, including approximately 1.4M candidates with somatic SNV / INDEL markers and approximately 8.2M candidates marked as non-somatic.

[0408] SEQC-WES-spiked into WGS10, SEQC-FFPE-spiked into WGS10, SEQC-FFPE-WES-spiked into WGS10 models

[0409] Training data for these models was prepared by combining 10% of the training candidates for the SEQC-WGS-spike-in model with all candidates from SEQC-WES-spike-in, SEQC-FFPE-spike-in, and SEQC-FFPE-WES-spike-in, respectively. This combination exploits both the large number of spike-in WGS mutations and the sample bias of WES and FFPE samples.

[0410] Sample-specific models

[0411] In addition to the general models described above, we derived sample-specific models for a set of nine samples across multiple data types. The samples selected included WGS samples sequenced by the NovaSeq instrument, WES samples, samples with 1ng of DNA prepared using the Nextera Flex Library Prep Kit, 30× WGS samples with 50% tumor purity, 100× WGS samples with 10% tumor purity, and two WGS and two WES FFPE samples each treated with formalin for 24 hours and matched to fresh / FFPE normal samples. For each sample, we prepared an in silico tumor using random spike-in. For the 10% tumor sample, AFs with spike-in mutations were randomly selected from a beta distribution (shape parameters α= 2 and β=18). For the other samples, AFs were selected using a beta distribution (shape parameters α= 2 and β=5). Then, for each sample, we trained the NeuSomatic and NeuSomatic-S models using in silico tumor-normal replicates. In addition, for each sample, we also trained a unique model by combining 10% of the training candidates for the SEQC-WGS-spiking model with the training data derived for that sample.

[0412] Somatic mutation detection algorithm

[0413] In addition to NeuSomatic, we used seven somatic mutation callers: MuTect2 (4.beta.6) 2 、SomaticSniper (1.0.5.0) 4 , Lancet (1.0.7) 8 、Strelka2 (2.8.4) 5 、TNscope(201711.03) 7 、MuSE (v1.0rc) 3 and VarDict (v1.5.1) 6, and ran each of them with default parameters or parameters recommended by the user manual. For SomaticSniper, we used “-q 1 -Q 15 -s 1e-05”. For TNscope (201711.03), we used the version implemented in Seven Bridges's CGC with the following command: “sentieondriver -i $tumor_bam -i $normal_bam -r $ref --algo TNscope --tumor_sample $tumor_sample_name --normal_sample $normal_sample_name -d $dbsnp$output_vcf”. For Lancet, we used “--cov-thr 10 --cov-ratio 0.005 --max-indel-len 50 -e 0.005”. The high confidence outputs marked as “PASS” in the final VCF file will be used for our comparative analysis. For VarDict, we also restricted the call results with “Somatic” status. The results from each caller used for comparison are all candidate mutants that would otherwise be considered by the user as “real” mutations detected by this caller.

[0414] We used NeuSomatic (0.1.4) in both ensemble and standalone modes. For ensemble mode, in addition to candidate variants identified by the “alignment scan” step of NeuSomatic, we merged the calls from five separate callers (MuTect2, SomaticSniper, Strelka2, MuSE, and VarDict) and represented their occurrences as additional input channels for each candidate variant. For NeuSomatic and NeuSomatic-S, we used “--scan_maf0.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” in the preprocessing step. For training, we used "--coverage_thr 300 --batch_size 8000".

[0415] Difficult genomic regions

[0416] We used a set of difficult genomic regions derived from the Genomes in a Bottle consortium 19 These regions include tandem repeats (divided into two different categories of repeats less than 50 bp and greater than 50 bp) as well as segmental duplication regions. We evaluated different techniques on these regions to show robustness to the complexity of the genomic context.

[0417] Analysis of false negative call results

[0418] To better understand the performance differences between the different techniques, we performed a pairwise comparison of the VAF distribution of false negative SNVs for NeuSomatic with other schemes on WGS datasets (Figure 42). In each X-vs-Y analysis, we identified XFN, which is the set of ground truth SNVs that were missed by algorithm X in at least 11 of the 21 WGS replicates (false negatives in X). Similarly, we identified YFN, which is the set of ground truth SNVs that were missed by algorithm Y in at least 11 of the 21 WGS replicates (false negatives in Y). We then calculated the private false negatives of X, defined as XFN / YFN, and the private false negatives of Y, defined as YFN / XFN. Then, for the somatic mutations in each of those sets, we calculated the distribution of the VAF of the mutations.

[0419] Evaluation Process

[0420] For fair comparison, we evaluated all models and somatic mutation algorithms on the 50% held-out genomic regions (regions not used to train the SEQC-WGS-GT-50 model). This ~1.4GB region contains ~21K SNVs and ~1.3K INDELs from the SEQC-II truth set for HCC1395, which were not used to train any NeuSomatic model except for the SEQC-WGS-GT-ALL model.

[0421] The SEQC-II consortium grouped the calls labeled as high-confidence and low-confidence for HCC1395 as the “truth set” for somatic mutations used in this paper. We adopted this truth set to calculate the true positives and false negatives for all channels in all copies of HCC1395. According to the recommendation of the SEQC-II consortium, we also blacklisted the low-confidence calls because of their low validation rate. Since the truth set provided by the consortium has a VAF detection limit of 5% and a deep detection limit of 50×, we cannot determine the true somatic status of calls with higher coverage data or calls with lower VAF. Therefore, for private calls reported by any channel that is not part of the truth set, we excluded calls that were considered ambiguous in the evaluation. In summary, for a tumor replica with average coverage C and tumor purity P at hand, if the pipeline reports that the replica has a private somatic mutation supported by d reads (which is not in the truth set), we mark this call as a false positive only if its inferred number of supporting reads at 100% purity and 50x coverage is ≥3. In other words, if d ≥3CP / 50, the call is a false positive; otherwise, we exclude the call from the evaluation.

[0422] For WES and AmpliSeq data, we only reported SNV assessments due to the very limited number of true indels in the assessed regions.

[0423] References

[0424] 1. Xu, C., A review of somatic single nucleotide variant calling algorithms for next-generation sequencing data. Comput. Struct. Biotechnol. J. 16, 15–24 (2018).

[0425] 2.Cibulskis, K. et al., Sensitive detection of somatic point mutations in impure and heterogeneous cancersamples. Nat. Biotechnol. 31, 213 (2013).

[0426] 3. Fan, Y. et al., MuSE: accounting for tumor heterogeneity using a sample-specific error model improves sensitivity and specificity in mutation calling from sequencing data. Genome Biol. 17, 178 (2016).

[0427] 4. Larson, D. E. et al., SomaticSniper: identification of somatic point mutations in whole genome sequencing data. Bioinformatics 28, 311–317 (2011).

[0428] 5. Kim, S. et al., Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods 15, 591–594 (2018).

[0429] 6. Lai, Z. et al., VarDict: a novel and versatile variant caller for next-generation sequencing in cancer research. Nucleic Acids Res. 44, e108–e108(2016).

[0430] 7. Freed, D., Pan, R. and Aldana, R. TNscope: Accurate Detection of Somatic Mutations with Haplotype-based Variant Candidate Detection and Machine Learning Filtering. bioRxiv 250647

[0431] 8. Narzisi, G. et al., Lancet: genome-wide somatic variant calling using localized colored DeBruijn graphs. Communications biology 1, 20 (2018).

[0432] 9. Koboldt, D. C. et al., VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).

[0433] 10. Fang, L. T. et al., An ensemble approach to accurately detect somatic mutations using SomaticSeq. Genome Biol. 16, 197 (2015).

[0434] 11. Sahraeian, S. M. E., Liu, R., Lau, B., Podesta, K., Mohiyuddin, M. and Lam, H. Y., Deep convolutional neural networks for accurate somatic mutation detection. Nature communications. 10, 1041 (2019).

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

[0436] 13. Fang, L. T. et al., Establishing reference samples for detection of somatic mutations and germline variants with NGS technologies, https: / / doi.org / 10.1101 / 625624 (2019).

[0437] 14. Ewing, A. D. et al., Combining tumor genome simulation with crowdsourcing to benchmark somatic single - nucleotide - variant detection. Nat. Methods 12, 623(2015).

[0438] 15. Bolger, A. M., Lohse, M. and Usadel, B., Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114 - 2120 (2014).

[0439] 16. Li, H., Aligning sequence reads, clone sequences and assembly contigs with BWA - MEM.https: / / arxiv.org / abs / 1303.3997 (2013).

[0440] 17. Ewing, A. D. et al., Combining tumor genome simulation with crowdsourcing to benchmark somatic single - nucleotide - variant detection. Nat. Methods 12, 623 (2015).

[0441] 18. Grossman, R. L. et al., Toward a shared vision for cancer genomic data. N. Engl. J. Med. 375, 1109–1112 (2016).

[0442] 19. Krusche, Peter, Len Trigg, Paul C. Boutros, Christopher E. Mason, M. Francisco, Benjamin L. Moore, Mar Gonzalez-Porta, et al., “Best practices for benchmarking germline small-variantcalls in human genomes.” Naturebiotechnology 1 (2019).

[0443] Although the present disclosure has been described with reference to some illustrative embodiments, it should be understood that those skilled in the art can design many other modifications and embodiments within the spirit and scope of the principles of the present disclosure. More specifically, without violating the spirit of the present disclosure, within the scope of the above disclosure, the drawings and the appended claims, the components and / or arrangements of the subject combination arrangement may be reasonably varied and modified. In addition to the variations and modifications in the components and / or arrangements, alternative uses will also be apparent to those skilled in the art.

[0444] All U.S. patents, U.S. patent application publications, U.S. patent applications, foreign patents, foreign patents, and non-patent publications mentioned in this specification and / or listed in the application data sheet are incorporated herein by reference in their entirety. If it is necessary to adopt the concepts of the various patents, applications, and publications to provide other embodiments, various aspects of the embodiments may be modified.

[0445] When a feature or element is referred to as being "on" another feature or element in this article, it can be directly located on another feature or element, or there may also be intermediate features and / or elements. On the contrary, when a feature or element is referred to as being "directly located" on another feature or element, there are no intermediate features or elements. It should also be understood that when a feature or element is referred to as being "connected", "attached" or "coupled" to another feature or element, it can be directly connected, attached or coupled to another feature or element, or there may be intermediate features or elements. On the contrary, when a feature or element is referred to as being "directly connected", "directly attached" or "directly coupled" to another feature or element, there are no intermediate features or elements. Although described or shown with respect to one embodiment, the features and elements described or shown in this way can be applied to other embodiments. Those skilled in the art will also appreciate that the structure or feature mentioned to be "adjacent" to another feature may have a part overlapping or located below the adjacent feature.

[0446] The terms used herein are only for the purpose of describing specific embodiments and are not intended to limit the present invention. For example, as used herein, the singular forms "one", "the", "the" and "said" are also intended to include plural forms, unless the context clearly indicates otherwise. It will be further understood that when the terms "comprise" and / or "include" are used in this specification, they specify the presence of the features, steps, operations, elements and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, elements, components and / or their groups. As used herein, the term "and / or" includes any and all combinations of one or more associated listed items and can be abbreviated as " / ".

[0447] Spatial relative terms such as "below", "under ... ", "below", "above", "above", etc. may be used herein to describe the relationship between an element or feature shown in the figure and another element or feature. It should be understood that, in addition to the orientation depicted in the figure, spatial relative terms are also intended to cover different orientations of the device in use or operation. For example, if the device in the accompanying drawings is turned over, the element described as "below" or "under ... " of other elements or features will be oriented to be "above" other elements or features. Therefore, the exemplary term "below" can cover both upper and lower orientations. The device can be oriented in other ways (rotated 90 ° or in other orientations), and the spatial relative descriptors used herein are interpreted accordingly. Similarly, unless otherwise specifically noted, the terms "upward", "downward", "vertical", "horizontal", etc. are only used for the purpose of explanation in this article.

[0448] Although the terms "first" and "second" may be used herein to describe various features / elements (including steps), these features / elements should not be limited by these terms unless the context indicates otherwise. These terms may be used to distinguish one feature / element from another feature / element. Thus, the first feature / element discussed below may be referred to as the second feature / element, and similarly, the second feature / element discussed below may be referred to as the first feature / element without departing from the teachings of the present invention.

[0449] Throughout this specification and the claims that follow, unless the context requires otherwise, the word "comprise" and variations such as "include" and "comprising" mean that various components (e.g., compositions and apparatus including devices and methods) can be used together in methods and articles. For example, the term "comprising" will be understood to imply the inclusion of any stated elements or steps, but not the exclusion of any other elements or steps.

[0450] As used herein in the specification and claims, including in the examples, and unless otherwise expressly stated, all numbers may be read as beginning with the word "about" or "approximately", even if the term does not appear explicitly. When describing amplitude and / or position, the phrase "about" or "approximately" may be used to indicate that the value and / or position described is within the reasonable expected range of the value and / or position. For example, the value of a numerical value may be + / -0.1% of a specified value (or range of values), + / -1% of a specified value (or range of values), + / -2% of a specified value (or range of values), + / -5% of a specified value (or range of values), + / -10% of a specified value (or range of values), etc. Any numerical value given herein should also be understood to include about or approximately that value unless the context otherwise indicates. For example, if the value "10" is disclosed, "about 10" is also disclosed. Any numerical range described herein is intended to include all sub-ranges contained therein. It is also understood that when a value is disclosed as "less than or equal to" the value, "greater than or equal to the value" and possible ranges between values ​​are also disclosed, as appropriately understood by one skilled in the art. For example, if a value "X" is disclosed, "less than or equal to X" and "greater than or equal to X" are also disclosed (e.g., where X is a numerical value). It is also understood that throughout the application, data is provided in a variety of different formats and that the data represents endpoints and starting points and ranges for any combination of data points. For example, if a specific data point "10" and a specific data point "15" are disclosed, it is understood that greater than, greater than or equal to, less than, less than or equal to, and equal to 10 and 15 are considered disclosed, and between 10 and 15 are also considered disclosed. It is also understood that every unit between two specific units is also disclosed. For example, if 10 and 15 are disclosed, 11, 12, 13, and 14 are also disclosed.

[0451] Although various illustrative embodiments are described above, any of a number of changes may be made to the various embodiments without departing from the scope of the invention as described in the claims. For example, in alternative embodiments, the order of the various method steps described may often be changed, while in other alternative embodiments, one or more method steps may be skipped entirely. Optional features of various device and system embodiments may be included in some embodiments and may not be included in other embodiments. Therefore, the foregoing description is provided primarily for exemplary purposes and should not be construed as limiting the scope of the invention, as the scope of the invention is set forth in the claims.

[0452] The examples and illustrative descriptions included herein show specific embodiments in which the present subject matter can be practiced by way of illustrative description and not limitation. As mentioned, other embodiments can be utilized and derived therefrom, so that structural and logical replacements and changes can be made without departing from the scope of the present disclosure. For convenience only, if multiple inventions or inventive concepts are in fact disclosed, the term "invention" may be used herein to refer to these embodiments of the present subject matter individually or collectively, without intending to limit the scope of the present application to any single invention or inventive concept. Therefore, although specific embodiments have been exemplified and described herein, any arrangement calculated to achieve the same purpose may replace the specific embodiments shown. The present disclosure is intended to cover any and all modifications or changes of various embodiments. By reading the above description, the combination of the above embodiments and other embodiments not specifically described herein will be apparent to those skilled in the art.

Claims

1. A method for germline variant calling, the method comprising: obtaining a reference sequence, a plurality of sequence reads, and positions of candidate variants within the sequence reads; Obtaining an expanded sequence read by inserting one or more spaces into one or more sequence reads; Obtaining an augmented reference sequence by inserting one or more spaces into the reference sequence; converting fragments of the expanded sequence reads surrounding the candidate variant into a sample matrix; converting segments of the augmented reference sequence surrounding the candidate variant into a reference matrix; providing the sample matrix and the reference matrix to a trained neural network; as well as At the output of the trained neural network, prediction data regarding variants within the plurality of sequence reads is obtained.

2. The method of claim 1 , further comprising detecting one or more inserted bases within the plurality of sequence reads, wherein augmenting the sequence reads and the reference sequence comprises: For each inserted base detected in any of the sequence reads, a space is inserted into the reference sample at the position of the inserted base.

3. The method according to claim 2, further comprising: For each inserted base detected in any of the sequence reads, a space is inserted at the position of the inserted base into each of the sequence reads in which no insertion was detected at the position of the inserted base.

4. The method of claim 1, wherein the sample matrix comprises: at least four rows representing four types of nucleotide bases, each row representing the base number of the corresponding nucleotide base type at a different position within the fragment of the expanded sequence read; and At least one row representing the number of inserted spaces at different positions within the fragments of the expanded sequence read.

5. The method of claim 4, wherein the reference matrix has the same dimensionality as the sample matrix, and wherein the reference matrix provides a complete representation of the locations of different nucleotide bases and spaces within the augmented reference sequence.

6. The method of claim 1, wherein the trained neural network comprises a trained convolutional neural network.

7. The method of claim 1 , further comprising providing the trained neural network with at least one of: a variant position matrix representing the positions of the candidate variants within the segments of the expanded sequence reads; a coverage matrix representing the coverage or depth of the fragments read by the expanded sequence; an alignment feature matrix representing alignment features of the expanded sequence reads; A knowledgeable base matrix that represents information about publicly known information about one or more variants.

8. The method of claim 1, wherein the predictive data about the variant comprises at least one of the following: the predicted type of the variant; the predicted position of the variant; the predicted length of the variant; and The predicted genotype of the variant.

9. The method of claim 1, wherein the prediction data about a variant comprises a predicted type for the variant, and wherein the neural network is configured to generate one of a plurality of values ​​for the predicted type for the variant, the plurality of values ​​comprising: a first value indicating the probability that the variant is a false positive; a second value indicating the probability that the variant is a single nucleotide polymorphism variant; a third value indicating the probability that the variant is a deletion variant; and A fourth value indicates the probability that the variant is an insertion variant.

10. A method for somatic variant calling, the method comprising: obtaining multiple normal sequence reads and multiple tumor sequence reads; Converting the normal sequence read segments and the tumor sequence read segments into a normal sample matrix and a tumor sample matrix, respectively; wherein converting the normal sequence read segments into the normal sample matrix comprises extending the normal sequence read segments by inserting one or more spaces into one or more normal sequence reads, and converting the tumor sequence read segments into the tumor sample matrix comprises extending the tumor sequence read segments by inserting one or more spaces into one or more tumor sequence reads; Feeding the normal sample matrix and the tumor sample matrix into a trained convolutional neural network; as well as At the output of the trained convolutional neural network, predicted types of somatic variants within the plurality of tumor sequence reads are obtained. 11 . The method of claim 10 , wherein the plurality of tumor sequence reads represent genetic information of a tumor sample of a patient, and the plurality of normal sequence reads represent genetic information of a normal sample of the patient.

12. The method of claim 10, wherein the tumor sample matrix comprises: at least one row for each nucleotide base type, each row representing the number of times the corresponding nucleotide base type occurs at each position within the fragment of the tumor sequence read; and At least one row representing the number of inserted spaces at each position within the fragments of the tumor sequence read.

13. The method according to claim 10, further comprising: The trained convolutional neural network is provided with one or more matrices representing one or more features obtained from one or more other variant callers that have analyzed the plurality of tumor sequence reads and / or the plurality of normal sequence reads.

14. The method according to claim 10, further comprising: Obtain reference sequences; converting the reference sequence into a reference matrix; as well as The reference matrix is ​​fed into a trained convolution matrix together with the normal sample matrix and the tumor sample matrix.

15. A non-transitory computer-readable medium comprising instructions that, when executed by one or more processors of a computing system, cause the computing system to perform operations comprising: obtaining multiple normal sequence reads and multiple tumor sequence reads; The fragments read from the normal sequence and the fragments read from the tumor sequence are converted into a normal sample matrix and a tumor sample matrix respectively, wherein: Converting the segments of the normal sequence reads to the normal sample matrix comprises augmenting the segments of the normal sequence reads by inserting one or more spaces into one or more normal sequence reads, and converting the segments of the tumor sequence reads to the tumor sample matrix comprises augmenting the segments of the tumor sequence reads by inserting one or more spaces into one or more tumor sequence reads; Feeding the normal sample matrix and the tumor sample matrix into a trained convolutional neural network; as well as At the output of the trained convolutional neural network, predicted types of somatic variants within the plurality of normal sequence reads are obtained.

16. A computing system comprising one or more processors and coupled to one or more non-transitory computer-readable memories storing instructions that, when executed by the computing system, cause the computing system to perform operations comprising: obtaining multiple normal sequence reads and multiple tumor sequence reads; The fragments read from the normal sequence are converted into a normal sample matrix, where: Converting the segments of the normal sequence reads into the normal sample matrix includes expanding the segments of the normal sequence reads by inserting one or more spaces into one or more normal sequence reads; Obtaining an augmented tumor sequence read by inserting one or more spaces into one or more tumor sequence reads; Converting the fragments read from the tumor sequence into a tumor sample matrix; Feeding the normal sample matrix and the tumor sample matrix into a trained neural network; as well as At the output of the trained neural network, predicted types of somatic variants within the plurality of tumor sequence reads are obtained.

17. A method for variant calling, the method comprising: Obtain reference sequences and multiple sequence reads; Optionally performing a first alignment of the plurality of sequence reads to the reference sequence, unless the plurality of sequence reads and the reference sequence are obtained in an already aligned configuration; identifying candidate variant positions from the aligned sequence reads and a reference sequence; augmenting the sequence reads and / or the reference sequence around the candidate variant position to achieve a second alignment of the plurality of sequence reads to the reference sequence, wherein the step of augmenting the sequence reads and / or the reference sequence comprises introducing one or more spaces into the sequence reads and / or the reference sequence to account for insertions and / or deletions in the sequence reads; generating a reference matrix for the candidate variant positions from the augmented reference sequence and generating a sample matrix for the candidate variant positions from the plurality of augmented sequence reads; Inputting the reference matrix and the sample matrix into a neural network; as well as The neural network is used to determine whether a variant type exists at the candidate variant position.

18. The method according to claim 17, further comprising: generating a plurality of training matrices from a training data set, wherein the training matrices have structures corresponding to the sample matrix and the reference matrix, wherein the training data set comprises sequence data, the sequence data comprises a plurality of mutations, the mutations comprise single nucleotide variants, insertions, and deletions; and The neural network is trained using the plurality of training matrices.

19. The method of claim 18, wherein the training data set comprises a plurality of subsets, wherein each subset comprises a tumor purity level ranging from 0% to 100%, wherein at least two of the subsets each have a different tumor purity level.

20. The method of claim 19, wherein at least three of the subsets each have a different level of tumor purity.

21. The method of claim 19, wherein the plurality of subsets comprises a first subset having a tumor purity level of less than about 30%, a second subset having a tumor purity level between about 30% and 70%, and a third subset having a third tumor purity level of at least about 70%.

22. The method of claim 19, wherein the plurality of subsets comprises a first subset having a tumor purity level of less than about 40%, a second subset having a tumor purity level between about 40% and 60%, and a third subset having a tumor purity level of at least about 60%.

23. The method of any one of claims 19 to 22, wherein the plurality of subsets comprises subsets having a tumor purity level of less than about 10%.

24. The method of any one of claims 19 to 22, wherein the plurality of subsets comprises subsets having a tumor purity level of less than about 5%.

25. The method of claim 18, wherein the training data set comprises synthetic data.

26. The method of claim 25, wherein the synthetic data comprises artificially generated mutations, wherein the artificially generated mutations comprise single nucleotide variants, insertions, and deletions.

27. The method of claim 18, wherein the training data set comprises real data, wherein the real data comprises real mutations, wherein the real mutations comprise single nucleotide variants, insertions, and deletions.

28. The method of claim 18, wherein the training data set comprises a plurality of subsets, wherein each subset comprises variant allele frequencies ranging from 0% to 100%, wherein at least two of the subsets each have a different variant allele frequency level.

29. The method of claim 28, wherein at least three of the subsets each have a different variant allele frequency level.

30. The method of claim 28, wherein at least one of the subsets has a variant allele frequency of at least 2.5%.

31. The method of claim 28, wherein at least one of the subsets has a variant allele frequency of at least 5%.

32. The method of claim 28, wherein at least one of the subsets has a variant allele frequency of at least 10%.

33. The method of claim 17, further comprising inputting at least one prediction from at least one mutation calling algorithm into the neural network.

34. The method of claim 33, wherein the at least one prediction comprises at least three predictions from at least three independent mutation calling algorithms.

35. The method of claim 33, wherein the at least one prediction comprises at least five predictions from at least five independent mutation calling algorithms.

36. The method of claim 18, wherein the training data set comprises a mixture of synthetic data and real data.

37. The method of claim 36, wherein the training data set comprises at least 5% synthetic data.

38. The method of claim 36, wherein the training data set comprises at least 10% synthetic data.

39. The method of claim 18, wherein the training dataset comprises whole genome sequencing data.

40. The method of claim 18, wherein the training dataset comprises whole exome sequencing data.

41. The method of claim 18, wherein the training dataset comprises targeted sequencing data.

42. The method of claim 18, wherein the training data set comprises data obtained from formalin-fixed paraffin-embedded samples.

43. The method of claim 18, wherein the training dataset comprises at least two of whole genome sequencing data, whole exome sequencing data, targeted sequencing data, and data obtained from formalin-fixed paraffin-embedded samples.

44. The method of claim 18, wherein the training dataset comprises at least three of whole genome sequencing data, whole exome sequencing data, targeted sequencing data, and data obtained from formalin-fixed paraffin-embedded samples.

45. The method of claim 18, wherein the training dataset comprises whole genome sequencing data, whole exome sequencing data, targeted sequencing data, and data obtained from formalin-fixed paraffin-embedded samples.

46. ​​A method for variant calling, the method comprising: obtaining a reference sequence, multiple tumor sequence reads, and multiple normal sequence reads; Optionally performing a first alignment of the plurality of tumor sequence reads and the plurality of normal sequence reads with the reference sequence, unless the plurality of tumor sequence reads and the plurality of normal sequence reads and the reference sequence are obtained in an already aligned configuration; identifying candidate variant positions from aligned tumor sequence reads, normal sequence reads, and a reference sequence; augmenting the tumor sequence reads and / or the normal sequence reads and / or the reference sequence around the candidate variant position to achieve a second alignment of the plurality of tumor sequence reads and the plurality of normal sequence reads with the reference sequence, wherein the step of augmenting the tumor sequence reads and / or the normal sequence reads and / or the reference sequence comprises introducing one or more spaces into the tumor sequence reads and / or the normal sequence reads and / or the reference sequence to account for insertions and / or deletions in the sequence reads; generating a reference matrix for the candidate variant positions from the augmented reference sequence, and generating a tumor matrix for the candidate variant positions from the plurality of augmented tumor sequence reads, and generating a normal matrix for the candidate variant positions from the plurality of augmented normal sequence reads; Inputting the reference matrix, the tumor matrix and the normal matrix into a neural network; as well as The neural network is used to determine whether a variant type exists at the candidate variant position.

47. The method of claim 46, further comprising: generating a plurality of training matrices from a training data set, wherein the training matrices have structures corresponding to the tumor matrix, the normal matrix, and the reference matrix, wherein the training data set comprises tumor sequence data and normal sequence data; as well as The neural network is trained using the plurality of training matrices.

48. The method of claim 47, wherein both the tumor sequence data and the normal sequence data comprise a plurality of mutations comprising single nucleotide variants, insertions, and deletions.

49. The method of claim 47, wherein the normal sequence data comprises at most 5% of the tumor sequence data.

50. The method of claim 47, wherein the normal sequence data comprises at most 10% of the tumor sequence data.

51. The method of claim 47, wherein the tumor sequence data comprises a tumor purity level between about 10% and 100%.

52. The method of claim 47, wherein the training data set comprises a plurality of tumor sequence data subsets, wherein each tumor sequence data subset comprises a tumor purity level ranging from 10% to 100%, wherein at least two of the tumor sequence data subsets each have a different tumor purity level.

53. The method of claim 52, wherein at least three of the tumor sequence data subsets each have a different level of tumor purity.

54. The method of claim 52, wherein the plurality of tumor sequence data subsets comprise a first tumor sequence data subset having a tumor purity level of less than about 30%, a second tumor sequence data subset having a tumor purity level between about 30% and 70%, and a third tumor sequence data subset having a tumor purity level of at least about 70%.

55. The method of claim 52, wherein the plurality of tumor sequence data subsets comprise a first tumor sequence data subset having a tumor purity level of less than about 40%, a second tumor sequence data subset having a tumor purity level between about 40% and 60%, and a third tumor sequence data subset having a tumor purity level of at least about 60%.

56. The method of claim 47, wherein the training data set comprises synthetic data.

57. The method of claim 56, wherein the synthetic data comprises artificially generated mutations, wherein the artificially generated mutations comprise single nucleotide variants, insertions, and deletions.

58. The method of claim 47, wherein the training data set comprises real data, wherein the real data comprises real mutations, wherein the real mutations comprise single nucleotide variants, insertions, and deletions.

59. The method of claim 47, wherein the training dataset comprises whole genome sequencing data.

60. The method of claim 47, wherein the training dataset comprises whole exome sequencing data.

61. The method of claim 47, wherein the training dataset comprises targeted sequencing data.

62. The method of claim 47, wherein the training data set comprises data obtained from formalin-fixed paraffin-embedded samples.

63. A system, comprising: A processor configured to implement the steps recited in any one of claims 17 to 62.

Images