Semi-supervised learning methods and systems for training ensembles of deep convolutional neural networks

Abstract

The disclosed technology relates to constructing a convolutional neural network based classifier for variant classification. More specifically, it relates to training a convolutional neural network based classifier on training data using a back-propagation gradient update method, which gradually matches the output of the convolutional neural network based classifier with the corresponding true value label. The convolutional neural network based classifier includes a residual block group, each residual block group is parameterized by the number of convolution filters in the residual block, the convolution window size of the residual block, and the hole convolution rate of the residual block, the convolution window size varies between the residual block groups, and the hole convolution rate varies between the residual block groups. The training data includes benign training examples and pathogenic training examples of translation sequence pairs generated from benign variants and pathogenic variants.

Description

[0001] appendix

[0002] The Appendix includes a bibliography of potentially relevant references listed in papers written by the inventors. Subject matter herein is covered in the U.S. Provisional Patent Applications to which this application claims priority / benefit. These references are available upon request from legal counsel and are also available through the Global Archives.

[0003] Priority application

[0004] This application claims priority or the benefit of U.S. Provisional Patent Application No. 62 / 573,144, filed on October 16, 2017, entitled “Training a Deep Pathogenicity Classifier Using Large-Scale Benign Training Data,” by Hong Gao, Kai-How Farh, Laksshman Sundaram, and Jeremy Francis McRae (Attorney Docket No. ILLM1000-1 / IP-1611-PRV); U.S. Provisional Patent Application No. 62 / 573,149, filed on October 16, 2017, entitled “Pathogenicity Classifier Based on Deep Convolutional Neural Networks (CNNS),” by Kai-How Farh, Laksshman Sundaram, Samskruthi ReddyPadigepati, and Jeremy Francis McRae (Attorney Docket No. ILLM1000-2 / IP-1612-PRV); U.S. Provisional Patent Application No. 62 / 573,149, filed on October 16, 2017, entitled “Pathogenicity Classifier Based on Deep Convolutional Neural Networks (CNNS),” by Kai-How Farh, Laksshman Sundaram, Samskruthi ReddyPadigepati, and Jeremy Francis McRae (Attorney Docket No. ILLM1000-3 / IP-1614-PRV); Sundaram and Jeremy Francis McRae, entitled “Deep Semi-Supervised Learning for Generating Large-Scale Pathogenic Training Data,” filed on October 16, 2017 (Attorney Docket No. ILLM 1000-3 / IP-1613-PRV); and U.S. Provisional Patent Application No. 62 / 582,898, entitled “Pathogenicity Classification of Genomic Data Using Deep Convolutional Neural Networks (CNNs),” filed on November 7, 2017 by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM1000-4 / IP-1618-PRV). These provisional applications are hereby incorporated herein by reference for all purposes.

[0005] Incorporation

[0006] The following are incorporated herein by reference for all purposes as if fully set forth herein:

[0007] PCT Patent Application No. PCT / US18 / 55840, entitled “Deep Convolutional Neural Network Training Method Based on Deep Learning,” filed by Hong Gao, Kai-How Farh, Laksshman Sundaram and Jeremy Francis McRae on October 15, 2018 (Attorney Docket No. ILLM 1000-8 / IP-1611-PCT) and subsequently published as PCT Publication No. WO______.

[0008] PCT Patent Application No. PCT / US18 / 55878, entitled “Deep Convolutional Neural Networks for Variant Classification,” filed on October 15, 2018 by Laksshman Sundaram, Kai-How Farh, Hong Gao, Samskruthi ReddyPadigepati, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-9 / IP-1612-PCT) and subsequently published as PCT Publication No. WO______.

[0009] U.S. nonprovisional patent application titled “Deep Convolutional Neural Network Training Technique Based on Deep Learning” filed concurrently by Hong Gao, Kai-How Farh, Laksshman Sundaram, and Jeremy Francis McRae (Attorney Docket No. ILLM1000-5 / IP-1611-US).

[0010] U.S. nonprovisional patent application titled “DEEP CONVOLUTIONAL NEURAL NETWORKS FOR VARIANT CLASSIFICATION” filed concurrently by Laksshman Sundaram, Kai-How Farh, Hong Gao, Samskruthi ReddyPadigepati, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-6 / IP-1612-US).

[0011] U.S. non-patent application titled “Semi-supervised learning for training a set of deep convolutional neural networks,” filed concurrently by Laksshman Sundaram, Kai-How Farh, Hong Gao, and Jeremy Francis McRae (Attorney Docket No. ILLM1000-7 / IP-1613-US).

[0012] Reference 1 – S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, “WAVENET: A GENERATIVE MODEL FOR RAW AUDIO,” arXiv:1609.03499, 2016;

[0013] Reference 2 – M. Chrzanowski, A. Coates, G. Diamos, A. Gibiansky, Y. Kang, X. Li, J. Miller, A. Ng, J. Raiman, S. Sengupta and M. Shoeybi, “DEEP VOICE: REAL-TIME NEURAL TEXT-TO-SPEECH,” arXiv:1702.07825, 2017;

[0014] Reference 3 – F. Yu and V. Koltun, “MULTI-SCALE CONTEXT AGGREGATION BY DILATED CONVOLUTIONS,” arXiv:1511.07122, 2016;

[0015] Reference 4 – K. He, X. Zhang, S. Ren, and J. Sun, “DEEP RESIDUAL LEARNING FOR IMAGE RECOGNITION,” arXiv:1512.03385, 2015;

[0016] Reference 5 – R. K. Srivastava, K. Greff, and J. Schmidhuber, “HIGHWAY NETWORKS,” arXiv:1505.00387, 2015;

[0017] Reference 6 – G. Huang, Z. Liu, L. van der Maaten and K. Q. Weinberger, “DENSELY CONNECTED CONVOLUTIONAL NETWORKS,” arXiv:1608.06993, 2017;

[0018] Reference 7–C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “GOING DEEPER WITH CONVOLUTIONS,” arXiv:1409.4842, 2014;

[0019] Reference 8–S. Ioffe and C. Szegedy, “BATCH NORMALIZATION: ACCELERATING DEEPNETWORK TRAINING BY REDUCING INTERNAL COVARIATE SHIFT,” arXiv:1502.03167, 2015;

[0020] Reference 9–J. M. Wolterink, T. Leiner, M. A. Viergever, and “DILATED CONVOLUTIONAL NEURAL NETWORKS FOR CARDIOVASCULAR MR SEGMENTATION IN CONGENITAL HEART DISEASE,” arXiv:1704.03669, 2017;

[0021] Reference 10–L. C. Piqueras, “AUTOREGRESSIVE MODEL BASED ON A DEEP CONVOLUTIONAL NEURAL NETWORK FOR AUDIO GENERATION,” Tampere University of Technology, 2016;

[0022] Reference 11–J. Wu, “Introduction to Convolutional Neural Networks,” Nanjing University, 2017;

[0023] Reference 12–I. J. Goodfellow, D. Warde-Farley, M. Mirza, A. Courville, and Y. Bengio, “CONVOLUTIONAL NETWORKS”, Deep Learning, MIT Press, 2016; and

[0024] Document 13 – J. Gu, Z. Wang, J. Kuen, L. Ma, A. Shahroudy, B. Shuai, T. Liu, X. Wang, and G. Wang, “RECENT ADVANCES IN CONVOLUTIONAL NEURAL NETWORKS,” arXiv:1512.07108, 2017.

[0025] Document 1 describes a deep convolutional neural network architecture that uses a group of residual blocks of convolution filters with the same convolution window size, a batch normalization layer, a rectified linear unit (abbreviated ReLU) layer, a dimension change layer, a dilated convolution layer with an exponentially increasing dilated convolution rate, a skip connection, and a softmax classification layer to accept an input sequence and produce an output sequence that scores the entries in the input sequence. The disclosed technology uses the neural network components and parameters described in Document 1. In one implementation, the disclosed technology modifies the parameters of the neural network components described in Document 1. For example, unlike Document 1, the dilated convolution rate in the disclosed technology evolves non-exponentially from a lower residual block group to a higher residual block group. In another example, unlike Document 1, the convolution window size in the disclosed technology varies between groups of residual blocks.

[0026] Reference 2 describes the details of the deep convolutional neural network architecture described in Reference 1.

[0027] Document 3 describes the dilated convolution used by the disclosed technology. As used herein, dilated convolution is also referred to as "expanded convolution". Dilated / expanded convolution allows for a large receptive field with few trainable parameters. Dilated / expanded convolution refers to a convolution in which a kernel is applied to an area greater than its length by skipping input values ​​with a certain step size (also referred to as a dilated convolution rate or dilation factor). Dilated / expanded convolution increases the spacing between elements of the convolution filter / kernel so that adjacent input entries (e.g., nucleotides, amino acids) are considered at a larger spacing when performing a convolution operation. This enables the incorporation of long-range contextual dependencies in the input. Dilated convolution retains part of the convolution calculation for reuse when processing adjacent nucleotides.

[0028] Document 4 describes the residual blocks and residual connections used by the disclosed technique.

[0029] Reference 5 describes the skip connections used by the disclosed technology. As used herein, the skip connections are also referred to as "highway networks".

[0030] Reference 6 describes a densely connected convolutional network architecture used by the disclosed technology.

[0031] Reference 7 describes a dimension-changing convolution layer and a block-based processing pipeline used by the disclosed technique. An example of a dimension-changing convolution is a 1×1 convolution.

[0032] Reference 8 describes a batch normalization layer used by the disclosed technique.

[0033] Reference 9 also describes the dilated / atrous convolutions used by the disclosed technique.

[0034] Document 10 describes various architectures of deep neural networks that can be used by the disclosed technology, including convolutional neural networks, deep convolutional neural networks, and deep convolutional neural networks with void / dilated convolutions.

[0035] Document 11 describes details of a convolutional neural network that can be used by the disclosed technology, including an algorithm for training a convolutional neural network with subsampling layers (e.g., pooling) and fully connected layers.

[0036] Reference 12 describes the details of various convolution operations that can be used by the disclosed technology.

[0037] Reference 13 describes various architectures of convolutional neural networks that can be used by the disclosed technology.

[0038] Incorporation by Reference Form for Electronic Submission of Applications

[0039] The following table file in ASCII text format is submitted with this application and is incorporated by reference. The name, creation date and size of the file are:

[0040]

[0041]

[0042] Supplementary Table 1: Details of variants derived from each species used in the analysis. This table includes intermediate results from each pipeline of these data sources. Note: This table is listed in Supplementary Table 1.txt.

[0043] Supplementary Table 2: Attrition of missense variants present in other species at common human allele frequency. The attrition was calculated based on the missense:synonymous ratio in common variants (>0.1%) compared to rare variants (<0.1%) using variants with the same status between humans and other species. Note: This table is listed in Supplementary Table 2.txt.

[0044] Supplementary Table 3: Depletion of missense variants present in other species at common human allele frequency, restricted to genes with >50% average nucleotide conservation between humans and other mammals. The depletion was calculated based on the missense:synonymous ratio in common variants (>0.1%) compared to rare variants (<0.1%) using variants with the same status between humans and other species. Note: This table is listed in Supplementary Table 3.txt.

[0045] Supplementary Table 4: Attrition of missense variants present as fixed alternatives in related species pairs at common human allele frequency. The attrition was calculated based on the missense:synonymous ratio in common variants (>0.1%) compared to rare variants (<0.1%) using variants with the same status between human and related species pairs. Note: This table is listed in Supplementary Table 3.txt.

[0046] Supplementary Table 6: Domain-specific annotations for the SCN2A gene. Wilcoxon rank sum p-values ​​indicate the divergence of PrimateAI scores for specific domains compared to the entire protein. The domains highlighted in bold cover approximately 7% of the protein but have the majority of ClinVar pathogenic annotations. This correlates well with the average PrimateAI score for each domain and is the top 3 pathogenic domains based on the PrimateAI model. Note: This table is listed in Supplementary Table 6.txt.

[0047] Supplementary Table 7: Raw counts used in calculating the effect of allele frequency on expected missense:synonymous ratios. Expected counts for synonymous and missense variants are based on variants in intronic regions, calculated using a trinucleotide background to control for mutation rate and gene conversion. Note: This table is listed in Supplementary Tables.xlsx.

[0048] Supplementary Table 13: List of protein names used to train deep learning models for 3-state secondary structure and 3-state solvent accessibility prediction. The label column indicates whether each protein was used in the training / validation / testing phase of model training. Note: This table is listed in Supplementary Table 13.txt.

[0049] Supplementary Table 18: List of 605 genes with nominally significant disease association (p<0.05) in the DDD study when calculated only by protein truncating variants Note: This table is listed in Supplementary Table 18.txt.

[0050] Supplementary Table 20: Results of DNM enrichment test for each gene among all genes with at least one observed de novo mutation (DNM). P values ​​are provided when all DNMs are included and after removing missense DNMs with PrimateAI scores < 0.803. FDR-corrected P values ​​are provided in a similar manner. Counts of protein truncations (PTVs) and missense DNMs observed by the DDD cohort only and by the full meta-analysis cohort are included. Similar counts of observed and expected missense DNMs are also included in the following cases: first, when all missense DNMs are included, and second, when all missense DNMs with PrimateAI scores < 0.803 are removed. Note: This table is listed in Supplementary Table 20.txt and Supplementary Table 20 Summary.txt.

[0051] Supplementary Table 21: Enrichment test results for de novo mutations in genes with FDR < 0.1. Includes counts of observed protein truncating (PTV) de novo mutations and counts of other protein altering de novo mutations, once with all missense de novo mutations and once with only damaging missense mutations. P values ​​when all missense sites are included vs. P values ​​after excluding low-scoring missense sites are provided. Note: This table is listed in Supplementary Table 21.txt.

[0052] Data File S1: List of all variants present in other species. The ClinVar Significance column lists the existing non-conflicting ClinVar annotations. Note: This table is included in Data File S1.txt.

[0053] Data File S2: List of all fixed substitutions originating from related species pairs. Note: This table is listed in Data File S2.txt.

[0054] Data File S3: List of hidden benign test variants associated with IBS in primates. Benign test variants are uncommon human variants associated with IBS in more than 1 primate species. Note: This table is listed in Data File S3.txt.

[0055] Data File S4: List of unmarked variants IBS and primates matched to hidden benign test variants. Unmarked variants were matched to benign test variants for mutation rate, coverage bias, and alignment to primate species. Note: This table is listed in Data File S4.txt.

[0056] Pathogenicity_prediction_model: Code written in the Python programming language that implements the disclosed technology according to one implementation. Note: This table is listed in pathogenicity_prediction_model.txt. Technical Field

[0057] The disclosed technology relates to artificial intelligence computers and digital data processing systems and corresponding data processing methods and products for intelligent simulation (i.e., knowledge-based systems, reasoning systems, and knowledge acquisition systems); and includes systems for uncertainty reasoning (e.g., fuzzy logic systems), adaptive systems, machine learning systems, and artificial neural networks. In particular, the disclosed technology relates to using deep learning-based techniques to train deep convolutional neural networks. Background Art

[0058] The subject matter discussed in this section should not be considered to be prior art simply because it is mentioned in this section. Similarly, issues mentioned in this section or associated with the subject matter provided as background should not be assumed to have been previously recognized in the prior art. The subject matter in this section merely represents different approaches, which themselves may also correspond to implementations of the claimed technology.

[0059] Machine Learning

[0060] In machine learning, input variables are used to predict output variables. Input variables are often called features and are represented by X = (X 1 ,X 2 ,...,X k ) means that each X i ,i∈1,...,k are features. The output variable is usually called the response variable or dependent variable and is represented by the variable Yi. The relationship between Y and the corresponding X can be written in general form:

[0061] Y=f(X)+∈

[0062] In the above equation, f is the feature (X 1 ,X 2 ,...,X k ) and ∈ is a random error term. The error term is independent of X and has a mean of zero.

[0063] In practice, feature X is available without Y or without knowing the exact relationship between X and Y. Since the error term has a mean of zero, the goal is to estimate the value f.

[0064]

[0065] In the above equation, is an estimate of ∈ and is usually considered a black box, meaning that only The relationship between the input and output was known, but the question of why it worked remained unanswered.

[0066] Use learning to find this function Supervised learning and unsupervised learning are two approaches used in machine learning for this task. In supervised learning, labeled data is used for training. By showing the input and the corresponding output (= label), the function is optimized so that it is close to the output. In unsupervised learning, the goal is to find hidden structures in unlabeled data. The algorithm does not have any measure of accuracy on the input data, which makes it different from supervised learning.

[0067] Neural Networks

[0068] Figure 1A An implementation of a fully connected neural network with multiple layers is shown. A neural network is a system of interconnected artificial neurons (e.g., a1, a2, a3) that exchange messages between each other. The neural network shown has three inputs, two neurons in a hidden layer and two neurons in an output layer. The hidden layer has an activation function f(·) and the output layer has an activation function g(·). The connection has digital weights (e.g., w11, w21, w12, w31, w22, w32, v11, v22) that are tuned during the training process so that the correctly trained network responds correctly to recognition when fed an image. The input layer processes the original input, and the hidden layer processes the output from the input layer based on the weights of the connection between the input layer and the hidden layer. The output layer obtains the output from the hidden layer and processes the output based on the weights of the connection between the hidden layer and the output layer. The network includes multiple layers of feature detection neurons. Each layer has many neurons that respond to different combinations of inputs from the previous layer. The layers are constructed so that the first layer detects a set of original patterns in the input image data, the second layer detects patterns of patterns, and the third layer detects patterns of these patterns.

[0069] A survey of deep learning applications in genomics can be found in the following publications:

[0070] T.Ching et al., Opportunities And Obstacles For Deep Learning InBiology And Medicine, www.biorxiv.org:142760, 2017;

[0071] Angermueller C, T,Parts L,Stegle O.Deep Learning For Computational Biology.Mol Syst Biol.2016;12:878;

[0072] Park Y,Kellis M.2015Deep Learning For RegulatoryGenomics.Nat.Biotechnol.33,825–826.(doi:10.1038 / nbt.3313);

[0073] Min,S.,Lee,B.&Yoon,S.Deep Learning In Bioinformatics.Brief.Bioinform.bbw068(2016);

[0074] Leung MK, Delong A, Alipanahi B et al.Machine Learning In GenomicMedicine: A Review of Computational Problems and Data Sets2016; and

[0075] Libbrecht MW, Noble WS. Machine Learning Applications In Genetics and Genomics. Nature Reviews Genetics 2015;16(6):321-32. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] In the drawings, like reference numerals generally refer to like parts in different views. Moreover, the drawings are not necessarily drawn to scale, but emphasis is generally placed on illustrating the principles of the disclosed technology. In the following description, various implementations of the disclosed technology are described with reference to the following drawings, in which:

[0077] Figure 1A One implementation of a feed-forward neural network with multiple layers is shown.

[0078] Figure 1B Describes one implementation of how convolutional neural networks work.

[0079] Figure 1C A block diagram for training a convolutional neural network according to one implementation of the disclosed technology is described.

[0080] Figure 1D is an implementation of a subsampling layer (averaging / max pooling) according to an implementation of the disclosed technology.

[0081] Figure 1E An implementation of a ReLU nonlinear layer in accordance with one implementation of the disclosed technology is shown.

[0082] Figure 1F Describes an implementation of double-layer convolution of the convolutional layer.

[0083] Figure 1G Describes adding residual connections that reinject prior information downstream via feature maps.

[0084] Figure 1H An implementation of residual blocks and skip connections is described.

[0085] Fig. 1I The batch normalization forward pass is shown.

[0086] Figure 1J Illustrates the batch normalization transformation at test time.

[0087] Figure 1K The batch normalization backward pass is shown.

[0088] Figure 1L Describes the use of batch normalization layers before and after convolutional or fully connected layers.

[0089] Figure 1M An implementation of 1D convolution is shown.

[0090] Figure 1N Explains how Global Average Pooling (GAP) works.

[0091] Fig.1O Illustrated dilated convolution.

[0092] Figure 1P One implementation of stacked atrous convolutions is shown.

[0093] Figure 1Q An example computing environment capable of operating the disclosed technology is shown.

[0094] Figure 2 An example structure of a deep residual network for pathogenicity prediction (referred to herein as “PrimateAI”) is shown.

[0095] Figure 3 A schematic diagram depicting the architecture of PrimateAI’s deep learning network for pathogenicity classification.

[0096] Figure 4A , 4B and 4C is Supplementary Table 16, showing example model structure details of the pathogenicity prediction deep learning model PrimateAI.

[0097] Figure 5 and 6 Deep learning networks for predicting protein secondary structure and solvent accessibility are illustrated.

[0098] Fig. 7A and 7B is Supplementary Table 11, showing example model structure details of the 3-state secondary structure prediction deep learning (DL) model.

[0099] Fig. 8A and 8B is Supplementary Table 12, showing example model structure details of the 3-state solvent accessibility prediction deep learning model.

[0100] Fig. 9 One implementation for generating reference sequences and alternative protein sequences from benign and pathogenic variants is described.

[0101] Fig.10 One implementation of aligning a reference sequence and an alternative protein sequence is shown.

[0102] Fig.11 It is an implementation to generate a position frequency matrix (abbreviated PFM), where PFM is also called position weighted matrix (abbreviated PWM) or position specific scoring matrix (abbreviated PSSM).

[0103] Fig.12 , 13 , 14 and 15 show the treatment of the secondary structure and solvent accessibility subnetworks.

[0104] Fig.16 Operation of the variant pathogenicity classifier. As used herein, the term variant is also referred to as single nucleotide polymorphism (abbreviation SNP), also commonly referred to as single nucleotide variant (abbreviation SNV).

[0105] Fig.17 Residual block plot is illustrated.

[0106] Fig.18 The neural network architecture of the secondary structure and solvent accessibility subnetworks is described.

[0107] Fig.19 The neural network architecture of the variant pathogenicity classifier is shown.

[0108] Fig. 20 The predicted pathogenicity scores at each amino acid position in the SCN2A gene are shown, with annotations for major functional domains.

[0109] Fig.21D Shown is a comparison of classifiers predicting a benign outcome for a test set containing 10,000 common primate variants held back from training.

[0110] Fig.21EIllustrated are the distributions of PrimateAI prediction scores for de novo missense variants arising in patients with developmental disorders (DDD) compared with their unaffected siblings, and the corresponding Wilcoxon rank-sum test P values.

[0111] Fig.21F Depicted is a comparison of classifiers when isolating de novo missense variants in DDD cases compared to controls. Wilcoxon rank sum test P values ​​for each classifier are shown.

[0112] Fig.22A Among the 605 genes associated with significant (P<0.05) de novo protein truncating variants, the enrichment of unexpected de novo missense mutations in affected individuals of the DDD cohort is shown.

[0113] Fig. 22B Depicted are the distributions of PrimateAI prediction scores for de novo missense variants occurring in patients with DDD compared with unaffected siblings in 605 relevant genes, and the corresponding Wilcoxon rank-sum test P values.

[0114] Fig. 22C Shown is a comparison of each classifier in separating de novo missense variants in cases compared to controls in 605 genes.

[0115] Fig.22D Comparison of classifiers by area under the curve (AUC) shown for each classifier (as shown on the receiver operating characteristic curve) is depicted.

[0116] Fig.22E The separation accuracy and area under the curve (AUC) of each classifier are depicted.

[0117] Fig.23A , 23B , 23C and 23D show the effect of the data used for training on the classification accuracy.

[0118] Fig.24 Correction for the effect of sequence coverage on confirmation of common primate variants is shown.

[0119] Fig.25A , 25B , 25C and 26 describe the recognition of protein motifs by the disclosed neural network. Fig.26 Included are line plots showing the impact of each position in and near the perturbed variant on the predicted deep learning score for that variant.

[0120] Fig. 27 The weight correlation patterns of simulated BLOSUM62 and Grantham score matrices are illustrated.

[0121] Fig.28A ,28B and 28C show the performance evaluation of the deep learning network PrimateAI and other classifiers.

[0122] Fig.29A and 29B The distribution of prediction scores of the four classifiers is illustrated.

[0123] Fig. 30A , 30B and 30C compared the accuracy of the PrimateAI network and other classifiers in separating pathogenic and benign variants in 605 disease-associated genes.

[0124] Fig.31A and 31B The correlation between the classifier's performance on ClinVar variants selected by human experts and the performance in an empirical implementation is shown.

[0125] Fig.32 is Supplementary Table 14, showing the performance of the 3-state secondary structure and 3-state solvent accessibility prediction models on annotated samples of the Protein Data Bank.

[0126] Fig.33 is Supplementary Table 15, showing the performance comparison of deep learning networks using annotated secondary structure tags of human proteins from the DSSP database.

[0127] Fig.34 Supplementary Table 17 shows the accuracy values ​​and de novo variant P values ​​for each of the 20 classifiers we evaluated on 10,000 hidden primate variants in DDD cases vs controls.

[0128] Fig.35 is Supplementary Table 19, showing the performance comparison of different classifiers for de novo variants in the DDD Cases Vs Controls database (limited to 605 disease-associated genes).

[0129] Fig.36 The computing environment of the disclosed semi-supervised learner is shown.

[0130] Fig.37 , 38 , 39, 40 and 41 show the various cycles of the disclosed semi-supervised learning.

[0131] Fig.42 is an illustration of the iterative balanced sampling process.

[0132] Fig.43 An implementation of a computing environment for generating a benign implementation is described.

[0133] Fig.44One implementation for generating benign human missense SNPs is described.

[0134] Fig.45 One implementation of a human homologous missense SNP is shown. A missense SNP in a non-human species with reference and alternative codons that match humans.

[0135] Fig.46 One implementation is described for classifying SNPs of non-human primate species (eg, chimpanzees) that have reference codons that match humans as benign.

[0136] Fig.47 An implementation of calculating enrichment scores and comparing them is described.

[0137] Fig.48 An implementation of one implementation of a benign SNP is described.

[0138] Fig.49A , 49B , 49C, 49D, and 49E describe missense / synonymous ratios in the human allele frequency spectrum.

[0139] Fig.50A , 50B , 50C and 50D show purifying selection of missense variants that are identical to the other species states.

[0140] Fig.51 Shown are the expected missense:synonymous ratios in the human allele frequency spectrum in the absence of purifying selection.

[0141] Fig.52A , 52B , 52C, and 52D depict the missense:synonymous ratios of CpG and non-CpG variants.

[0142] Fig.53 , 54 and 55 illustrate missense:synonymous ratios for human variants that are identical to the six primate states.

[0143] Fig.56 is a simulation showing the saturation of new common missense variants discovered by increasing the size of the human cohort investigated.

[0144] Fig.57 The accuracy of PrimateAI in different conservation profiles of the genome is shown.

[0145] Fig.58 is Supplementary Table 5, showing the contribution of common human variants and variants present in non-human primates to labeling one implementation of benign training.

[0146] Fig.59is Supplementary Table 8, showing the effect of allele frequency on the expected missense:synonymous ratio.

[0147] Fig.60 is Supplementary Table 9, showing the ClinVar analysis.

[0148] Fig.61 is Supplementary Table 10, showing the number of missense variants found in ClinVar for other species according to one implementation.

[0149] Fig.62 Table 1 shows one implementation for discovering 14 additional candidate genes for intellectual disability.

[0150] Fig.63 is Table 2, showing one implementation of the mean difference in Grantham scores between pathogenic and benign variants in ClinVar.

[0151] Fig.64 One implementation of the enrichment analysis for each gene is shown.

[0152] Fig.65 One implementation of gene set enrichment analysis is shown.

[0153] Fig.66 is a simplified block diagram of a computer system that can be used to implement the disclosed techniques. DETAILED DESCRIPTION

[0154] The following discussion is given to enable any person skilled in the art to make and use the disclosed technology, and is provided in the context of a specific application and its requirements. Various modifications to the disclosed implementations will be apparent to those skilled in the art, and the general principles defined herein may be applied to other implementations and applications without departing from the spirit and scope of the disclosed technology. Thus, the disclosed technology is not intended to be limited to the implementations shown, but rather should be given the widest scope consistent with the principles and features disclosed herein.

[0155] introduce

[0156] Convolutional Neural Networks

[0157] Convolutional neural networks are a special type of neural network. The fundamental difference between densely connected layers and convolutional layers is that dense layers learn global patterns in their input feature space, while convolutional layers learn local patterns: in the case of images, patterns found in a small 2D window of the input. This key feature gives convolutional neural networks two interesting properties: (1) the patterns they learn are translation invariant, and (2) they can learn spatial hierarchies of patterns.

[0158] Regarding the first one, after learning a certain pattern in the bottom right corner of an image, the convolutional layer can recognize it anywhere: for example, in the top left corner. If the densely connected network is presented at a new location, it will have to relearn this pattern. This makes convolutional neural networks data efficient, as they require fewer training samples to learn their representations that have the ability to generalize.

[0159] Regarding the second, the first convolutional layer can learn small local patterns such as edges, the second convolutional layer will learn larger patterns made up of features from the first layer, etc. This allows convolutional neural networks to effectively learn increasingly complex and abstract visual concepts.

[0160] A convolutional neural network learns highly nonlinear mappings by interconnecting layers of artificial neurons arranged in many different layers with activation functions that relate these layers. It includes one or more convolutional layers interspersed with one or more subsampling layers and nonlinear layers, which are typically followed by one or more fully connected layers. Each element of a convolutional neural network receives input from a set of features in the previous layer. Convolutional neural networks learn simultaneously because neurons in the same feature map have the same weights. These local shared weights reduce the complexity of the network, so that when multi-dimensional input data enters the network, the convolutional neural network avoids the complexity of data reconstruction in the feature extraction and regression or classification process.

[0161] Convolutions operate on a 3D tensor (called a feature map) with two spatial axes (height and width) and a depth axis (also called a channels axis). For RGB images, the dimension of the depth axis is 3, because images have three color channels; red, green, and blue. For black and white images, the depth is 1 (grayscale). The convolution operation extracts sub-blocks from its input feature map and applies the same transformation to all of them, producing an output feature map. This output feature map is still a 3D tensor: it has a width and a height. Its depth can be arbitrary, because the output depth is a parameter of the layer, and the different channels on this depth axis no longer represent specific colors as in the RGB input; instead, they represent filters. Filters encode specific aspects of the input data: for example, at the height level, a single filter can encode the concept of "there is a face in the input."

[0162] For example, the first convolutional layer takes a feature map of size (28, 28, 1) and outputs a feature map of size (26, 26, 32): it computes 32 filters on its input. Each of these 32 output channels contains a grid of 26x26 values, which is a map of the response of the filter on the input, indicating the response of this filter pattern at different locations in the input. This is what the term "feature map" means: each dimension in the depth axis is a feature (or filter), and the 2D tensor output [:,:,n] is a 2D spatial map of the response of this filter on the input.

[0163] A convolution is defined by two key parameters: (1) the size of the sub-blocks extracted from the input - these are typically 1x1, 3x3, or 5x5; and (2) the depth of the output feature map, which is the number of filters computed by the convolution. Typically, these start at a depth of 32, continue to a depth of 64, and end at a depth of 128 or 256.

[0164] Convolutions work by sliding these windows of size 3x3 or 5x5 over the 3D input feature map, stopping at each location, and extracting a 3D sub-block of shape (window_height, window_width, input_depth) of the surrounding features. Each such 3D sub-block is then transformed (via tensor product of matrices with the same learned weights, called a convolution kernel) into a 1D vector of shape (output_depth). All of these vectors are then spatially recombined into a 3D output map of shape (height, width, output_depth). Each spatial position in the output feature map corresponds to the same position in the input feature map (e.g., the bottom right corner of the output contains information about the bottom right corner of the input). For example, for a 3×3 window, the vector output [i,j,:] comes from the 3D sub-block input [i-1:i+1,j-1:j+1,:]. Figure 1B The whole process is described in detail in .

[0165] A convolutional neural network consists of convolutional layers that perform convolution operations between input values ​​and convolutional filters (matrices of weights) that are iteratively learned during training through multiple gradient updates. Let (m,n) be the filter size, and w be the matrix of weights, then the convolutional layer performs the convolution of W with the input X by computing the dot product W x+b, where x is an instance of X, and b is the bias. The step size that the convolutional filter slides over the input is called the stride, and the filter area (m×n) is called the receptive field. The same convolutional filter is applied to different positions of the input, which reduces the number of weights learned. It also allows for position-invariant learning, i.e., if there is an important pattern in the input, the convolutional filter can learn it no matter where it is in the sequence.

[0166] Training a Convolutional Neural Network

[0167] Figure 1C A block diagram of training a convolutional neural network according to one implementation of the disclosed technology is depicted. The convolutional neural network is adjusted or trained so that input data results in a specific output estimate. The convolutional neural network is adjusted using backpropagation based on a comparison of the output estimate with a ground truth until the output estimate asymptotically matches or approaches the ground truth.

[0168] Convolutional neural networks are trained by adjusting the weights between neurons based on the difference between the true value and the actual output. This is described mathematically as:

[0169]

[0170] Where δ = (true value) - (actual output)

[0171] In one implementation, the training rule is defined as:

[0172]

[0173] In the above equation: the arrow indicates the update of the value; t m is the target value of neuron m; is the current output of the neuron m being calculated; a n is the input n; and α is the learning rate.

[0174] The intermediate steps in training involve using convolutional layers to generate feature vectors from the input data. Gradients are calculated with respect to the weights in each layer (starting at the output). This is called a backward pass or backpropagation. The weights in the network are updated using a combination of negative gradients and the previous weights.

[0175] In one implementation, the convolutional neural network uses a stochastic gradient update algorithm (such as ADAM) that performs backpropagation of errors by means of gradient descent. An example of a backpropagation algorithm based on a sigmoid function is described below:

[0176]

[0177] In the sigmoid function above, is the weighted sum computed by the neuron. This sigmoid function has the following derivatives:

[0178]

[0179] The algorithm consists of computing the activations of all neurons in the network, producing the output for the forward pass. The activation of a neuron m in a hidden layer is described as:

[0180]

[0181]

[0182] This operation is performed for all hidden layers to obtain the activations, as described below:

[0183]

[0184]

[0185] Then, the error and correct weights of each layer are calculated. The error of the output is calculated as:

[0186]

[0187] The error in the hidden layer is calculated as:

[0188]

[0189] The weights of the output layer are updated as:

[0190]

[0191] The weights of the hidden layer are updated using the following learning rate α:

[0192] v nm ←w nm +αδhman

[0193] In one implementation, the convolutional neural network uses gradient descent optimization to calculate the error across all layers. In such optimization, for an input feature vector x and a predicted output For y prediction when the target is The cost of the loss function is defined as l, that is, l Use function f to transform the predicted output from the input feature vector x The function f is parameterized by the weights of the convolutional neural network, i.e. The loss function is described as Or Q(z,w)=l(f w (x), y), where z is the input and output data pair (x, y). Gradient descent optimization is performed by updating the weights according to the following condition:

[0194]

[0195] w t+1 =w t +v t+1

[0196] In the equation above, α is the learning rate. In addition, the loss is calculated as the average over a set of n data pairs. At linear convergence, when the learning rate α is small enough, the calculation terminates. In other implementations, the gradients are calculated using only selected data pairs fed to the Nesterov accelerated gradient and adaptive gradient to inject computational efficiency.

[0197] In one implementation, the convolutional neural network uses stochastic gradient descent (SGD) to compute the cost function. SGD works by learning the cost function from only one, random, data pair, z t , calculate the weights in the loss function to approximate the gradient, as described below:

[0198]

[0199] w t+1 =w t +v t+1

[0200] In the above equation: α is the learning rate; μ is the momentum; and t is the current weight state before the update. When the learning rate α is both fast enough and slow enough, the convergence rate of SGD is close to O(1 / t). In other implementations, the convolutional neural network uses different loss functions, such as Euclidean loss and softmax loss. In another implementation, the Adam random optimizer is used by the convolutional neural network.

[0201] Convolutional Layer

[0202] The convolutional layers of a convolutional neural network act as feature extractors. The convolutional layers act as adaptive feature extractors that are able to learn and decompose input data into hierarchical features. In one implementation, the convolutional layer takes two images as input and produces a third image as output. In such implementations, convolution operates on two images in two dimensions (2D), where one image is the input image and the other image (called the "kernel") is applied as a filter on the input image, thereby producing an output image. Thus, for an input vector f of length n and a kernel g of length m, the convolution f*g of f and g is

[0203] is defined as:

[0204]

[0205] The convolution operation consists of sliding a kernel over the input image. For each location of the kernel, the overlapping values ​​of the kernel and the input image are multiplied and the results are added. The sum of the products is the value of the output image at the point in the input image where the kernel is centered. The different outputs produced by many kernels are called feature maps.

[0206] Once the convolutional layers are trained, they can be applied to perform recognition tasks on new inference data. Since convolutional layers learn from the training data, they avoid explicit feature extraction and learn implicitly from the training data. Convolutional layers use convolutional filter kernel weights that are determined and updated as part of the training process. Convolutional layers extract different features of the input, which are combined at higher layers. Convolutional neural networks use different numbers of convolutional layers, each with different convolutional parameters such as kernel size, stride, padding, feature maps, and number of weights.

[0207] Subsampling Layer

[0208] Figure 1D An implementation of a subsampling layer according to an implementation of the disclosed technology. The subsampling layer reduces the resolution of features extracted by the convolutional layer, making the extracted features or feature maps more robust to noise and deformation. In one implementation, the subsampling layer uses two types of pooling operations: average pooling and maximum pooling. The pooling operation divides the input into non-overlapping two-dimensional spaces. For average pooling, the average of the four values ​​in the area is calculated. For maximum pooling, the maximum of the four values ​​is selected.

[0209] In one implementation, the subsampling layer includes a pooling operation of a group of neurons in the previous layer: their output is mapped to only one input in max pooling, and their output is mapped to the average of the inputs in average pooling. In max pooling, the output of the pooling neuron is the maximum value residing in the input, as described below:

[0210]

[0211] In the above equation, N is the total number of elements in the neuron group.

[0212] In average pooling, the output of the pooling neuron is the average of the input values ​​residing with the group of input neurons as described below:

[0213]

[0214] In the above equation, N is the total number of elements in the input neuron group.

[0215] exist Figure 1D In , the input size is 4×4. For 2×2 subsampling, the 4×4 image is divided into four non-overlapping 2×2 matrices. For average pooling, the average of the four values ​​is the integer output. For max pooling, the maximum of the four values ​​in the 2×2 matrix is ​​the integer output.

[0216] Non-linear layers

[0217] Figure 1E An implementation of a nonlinear layer according to an implementation of the disclosed technology is shown. The nonlinear layer uses different nonlinear trigger functions to send different recognition signals to possible features on each hidden layer. The nonlinear layer uses various specific functions to implement nonlinear triggering, including rectified linear unit (ReLU), hyperbolic tangent, absolute value of hyperbolic tangent, sigmoid, and continuous trigger (nonlinear) functions. In one implementation, the ReLU activation implements the function y=max(x,0) and keeps the input and output sizes of the layer the same. The advantage of using ReLU is that the training speed of convolutional neural networks is many times faster. ReLU is a non-continuous, non-saturated activation function that is linear with respect to the input if the input value is greater than zero, and zero otherwise. Mathematically, the ReLU activation function is described as:

[0218]

[0219]

[0220] In other implementations, convolutional neural networks use a power unit activation function, which is a continuous, non-saturating function described as follows:

[0221]

[0222] In the above equations, a, b and c are parameters that control shift, scale and power, respectively. If c is an odd number, the power activation function can produce x and y-antisymmetric activations, and if c is an even number, it can produce y-symmetric activations. In some implementations, the unit produces non-rectified linear activations.

[0223] In other implementations, convolutional neural networks use the sigmoid unit activation function, which is a continuous saturating function described by the following logistic function:

[0224]

[0225] In the equation above, β = 1. The sigmoid unit activation function does not produce negative activations and is only antisymmetric with respect to the y-axis.

[0226] Convolution Example

[0227] Figure 1F Depicts an implementation of two convolution layers in a convolutional layer. Figure 1FIn , an input of size 2048 dimensions is convolved. At Convolution 1, the input is convolved by a convolutional layer comprising two channels of sixteen kernels of size 3×3. The resulting 16 feature maps are then rectified with the help of a ReLU activation function at ReLU1 and then pooled in Pooling 1 with the help of average pooling using a 16 channel pooling layer with a kernel of size 3×3. At Convolution 2, the output of Pooling 1 is then convolved by another convolutional layer comprising sixteen channels of thirty kernels of size 3×3. This is followed by another ReLU2 and average pooling in Pooling 2 with a kernel size of 2×2. The convolutional layers use different strides and paddings, such as zero, one, two, and three. According to one implementation, the resulting feature vector is five hundred and twelve (512) dimensions.

[0228] In other implementations, the convolutional neural network uses different numbers of convolutional layers, subsampling layers, nonlinear layers, and fully connected layers. In one implementation, the convolutional neural network is a shallow network with fewer layers and more neurons per layer, such as one, two, or three fully connected layers with one hundred (100) to two hundred (200) neurons per layer. In another implementation, the convolutional neural network is a deep network with more layers and fewer neurons per layer, such as five (5), six (6), or eight (8) fully connected layers per layer, with thirty (30) to fifty (50) neurons per layer.

[0229] Forward pass

[0230] The output of the neuron at row x, column y in the l-th convolutional layer and the k-th feature map in the feature map is determined by the following equation:

[0231]

[0232] The output of the neuron at the xth row, yth column, and kth feature map in the lth subsampling layer is determined by the following equation:

[0233]

[0234] The output of the i-th neuron in the l-th output layer is determined by the following equation:

[0235]

[0236] Backward Propagation

[0237] The output bias of the kth neuron in the output layer is determined by the following equation:

[0238]

[0239] The input bias of the kth neuron in the output layer is determined by the following equation:

[0240]

[0241] The weight and bias changes for the kth neuron in the output layer are determined by the following equation:

[0242]

[0243]

[0244] The output bias of the kth neuron in the hidden layer is determined by the following equation:

[0245]

[0246] The input bias of the kth neuron in the hidden layer is determined by the following equation:

[0247]

[0248] The weight and bias changes in row x, column y of the mth feature map of the initial layer receiving input from k neurons in the hidden layer are determined by the following equation:

[0249]

[0250]

[0251] The output bias of row x, column y in the m-th feature map of the subsampling layer S is determined by the following equation:

[0252]

[0253] The input bias of row x and column y in the m-th feature map of the subsampling layer S is determined by the following equation:

[0254]

[0255] The weight and bias changes in row x, column y of the m-th feature map of the subsampling layer S and convolutional layer C are determined by the following equations:

[0256]

[0257]

[0258] The output bias of row x, column y in the k-th feature map of convolutional layer C is determined by the following equation:

[0259]

[0260] The input bias of row x and column y in the k-th feature map of convolutional layer C is determined by the following equation:

[0261]

[0262] The weight and bias changes in row r and column c of the m-th convolution kernel of the k-th feature map of the l-th convolution layer C are:

[0263]

[0264]

[0265] Residual Connection

[0266] Figure 1G Depicted is the addition of residual connections that reinject initial information downstream via feature maps. Residual connections involve reinjecting previous representations into the downstream data stream by adding past output tensors to later output tensors, which helps prevent information loss along the data processing stream. Residual connections address two common problems that plague any large-scale deep learning model: vanishing gradients and representational bottlenecks. In general, adding residual connections to any model with more than 10 layers is likely to be beneficial. As discussed above, residual connections involve making the output of an earlier layer available as an input to a later layer, effectively creating shortcuts in a sequential network. Instead of being concatenated to a later activation, the earlier output is added to the later activation, assuming that both activations are the same size. If they are of different sizes, a linear transformation that reshapes the earlier activation to the target shape can be used. Additional information about residual connections can be found in K. He, X. Zhang, S. Ren, and J. Sun, “Deep Residual Learning for Image Recognition,” arXiv:1512.03385, 2015, which is incorporated herein by reference for all purposes as if fully set forth herein.

[0267] Residual Learning and Skip Connections

[0268] Figure 1HDepicts one implementation of residual blocks and skip connections. The main idea of ​​residual learning is that residual mappings are easier to learn than the original mappings. Residual networks stack multiple residual units to mitigate the drop in training accuracy. Residual blocks utilize special additive skip connections to combat vanishing gradients in deep neural networks. At the beginning of a residual block, the data flow is split into two streams: the first stream carries the unchanged input of the block, while the second stream applies weights and nonlinearities. At the end of the block, the two streams are merged using an element-wise sum. The main advantage of this type of construction is that it allows gradients to flow through the network more easily. Additional information about residual block diagrams and skip connections can be found in AVD Ord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, "WAVENET: A Generative Model for Raw Rates," arXiv:1609.03499, 2016.

[0269] Thanks to residual networks, deep convolutional neural networks (CNNs) can be easily trained and have achieved improved accuracy for image classification and object detection. A convolutional feed-forward network connects the output of layer l as input to layer l+1, which results in the following layer transition: l =H l (x l-1 ). The residual block adds a skip connection that uses the identity function: l =H l (x l-1 )+x l-1 Bypassing nonlinear transformations. The advantage of the residual block is that the gradient can flow directly from the later layers to the previous layers through the identity function. However, the identity function and H are combined by summing l output, which may hinder the flow of information in the network.

[0270] Dilated Convolution

[0271] Fig.1O illustrates a dilated convolution. A dilated convolution, sometimes called a atrous convolution, literally means with a hole. The French name has its origin in the algorithm trous, which computes the fast binary wavelet transform. In these types of convolutional layers, the inputs corresponding to the receptive field of the filter are not adjacent points. This is in Fig.1O The input spacing depends on the stretch factor.

[0272] WaveNet

[0273] Wavenet is a deep neural network for generating raw audio waveforms.

[0274] WaveNet is different from other convolutional networks because it is able to obtain a relatively large "field of view" at a low cost. In addition, it is able to add conditioning of the signal both locally and globally, which allows Wavenet to be used as a text-to-speech (TTS) engine with multiple voices, where TTS gives local conditioning and voice-specific to global conditioning.

[0275] The main building block of WaveNet is the causal dilated convolution. As a dilation on the causal dilated convolution, WaveNet also allows stacking of these convolutions, such as Figure 1P As shown in the figure. To get the same receptive field with dilated convolutions, another dilated layer is needed. A stack is a repetition of dilated convolutions, connecting the outputs of the dilated convolutional layers to a single output. This allows WaveNet to have a large "field of view" for one output node at a relatively low computational cost. For comparison, to get a field of view of 512 inputs, a fully convolutional network (FCN) would require 511 layers. In the case of a dilated convolutional network, we would need eight layers. Stacked dilated convolutions only require seven layers with two stacks or six layers with four stacks. To get an idea of ​​the difference in computational power required to cover the same field of view, the table below shows the number of weights required for the network assuming one filter per layer and a filter width of two. Furthermore, it is assumed that the network is using 8-bit binary encoding.

[0276]

[0277] WaveNet adds skip connections before making residual connections, which bypass all the following residual blocks. Each of these skip connections is summed before passing through a series of activation functions and convolutions. Intuitively, this is the sum of the information extracted at each layer.

[0278] Batch Normalization

[0279] Batch Normalization is a method used to speed up the training of deep networks by making data normalization an integral part of the network architecture. Batch Normalization can adaptively normalize the data even if the mean and variance change over time during training. It works by internally maintaining an exponential moving average of the batch mean and variance of the data seen during training. The main effect of Batch Normalization is that it helps with gradient propagation (much like residual connections) and therefore allows deep networks. Some very deep networks can only be trained when they contain multiple Batch Normalization layers.

[0280] Batch Normalization can be thought of as another layer that can be inserted into the model architecture, just like a fully connected or convolutional layer. Batch Normalization layers are usually used after convolutional or densely connected layers. It can also be used before convolutional or densely connected layers. Both implementations can be used by the disclosed techniques and are Figure 1L . Batch Normalization layers take an axis argument that specifies the feature axis that should be normalized. This argument defaults to -1, the last axis in the input tensor. This is the correct value when using Dense, Conv1D, RNN, and Conv2D layers with the data format set to channels_last. However, in the niche use case of Conv2D layers with the data format set to channels_first, the feature axis is axis 1; the axis argument in Batch Normalization can be set to 1.

[0281] Batch normalization provides a definition for feeding forward the input and computing the gradients with respect to the parameters and its own input via a backward pass. In practice, a batch normalization layer is inserted after a convolutional or fully connected layer, but before the output is fed to the activation function. For convolutional layers, different elements of the same feature map (i.e., activations at different locations) are normalized in the same way so as to obey the convolutional properties. Therefore, all activations in a mini-batch are normalized at all locations rather than per activation.

[0282] Internal covariate shift is the main reason why deep architectures are notoriously slow to train. This stems from the fact that deep networks not only need to learn new representations at each layer, but also need to account for changes in their distribution.

[0283] Covariate shift is generally a known problem in the field of deep learning and often occurs in real-world problems. A common covariate shift problem is the difference in the distribution of the training and test sets, which can lead to suboptimal generalization performance. This problem is usually handled by a normalization or whitening preprocessing step. However, especially the whitening operation is computationally expensive and is therefore impractical in an online setting, especially if the covariate shift occurs throughout different layers.

[0284] Internal covariate shift refers to the phenomenon that the distribution of network activations changes across layers due to changes in network parameters during training. Ideally, each layer should be transformed into a space where they have the same distribution, but the functional relationship remains unchanged. To avoid expensive computation of the covariance matrix at each layer and each step to decorrelate and whiten the data, we normalize the distribution of each input feature in each layer across each mini-batch to have zero mean and one standard deviation.

[0285] Forward pass

[0286] During the forward pass, the mini-batch mean and variance are calculated. Using these mini-batch statistics, the data is normalized by subtracting the mean and dividing by the standard deviation. Finally, the data is scaled and shifted using the learned scaling and shifting parameters. Fig. 1I The batch normalization forward pass f is depicted in BN .

[0287] exist Fig. 1I In the corresponding β is the batch mean and is the batch variance. The learned scaling and shifting parameters are denoted by γ and β respectively. For clarity, in this paper, the batch normalization process is described for each activation and the corresponding index is omitted.

[0288] Since normalization is a differentiable transformation, the error is propagated into these learned parameters and it is therefore possible to recover the representative power of the network by learning the identity transformation. Conversely, by learning the same scaling and shifting parameters as the corresponding batch statistics, the batch normalization transformation will have no effect on the network if it is the best operation to perform. At test time, the batch mean and variance are replaced by the respective population statistics since the input does not depend on other samples from the mini-batch. Another approach is to keep a running average of the batch statistics during training and use these averages to compute the network output at test time. At test time, one can Figure 1J The following shows the batch normalization transformation. Figure 1J Medium, μ D and Accordingly, the population mean and variance are expressed rather than batch statistics.

[0289] Backward pass

[0290] Since normalization is a differentiable operation, we can Figure 1K The backward pass is computed as depicted in .

[0291] 1D Convolution

[0292] 1D convolution extracts local 1D sub-blocks or sub-sequences from a sequence, such as Figure 1M As shown, a 1D convolution obtains each output time step from a temporal sub-block in the input sequence. A 1D convolution layer recognizes local patterns in the sequence. Because the same input transformation is performed on each sub-block, a pattern learned at one position in the input sequence can be recognized later at a different position, making the translations of the 1D convolution layer invariant to temporal translations. For example, a 1D convolution layer that processes a sequence of bases using a convolution window of size 5 should be able to learn bases or base sequences of length 5 or less, and should be able to recognize base motifs in any context in the input sequence. Therefore, a basic level 1D convolution is able to learn basic morphologies.

[0293] Global Average Pooling

[0294] Figure 1NIt shows how Global Average Pooling (GAP) works. Global Average Pooling can be used to replace the Fully Connected (FC) layers for classification by taking the spatial average of the last layer elements for scoring. This reduces the training load and bypasses the problem of overfitting. Global Average Pooling applies a structure before the model and it is equivalent to a linear transformation with predefined weights. Global Average Pooling reduces the number of parameters and eliminates the fully connected layers. Fully connected layers are usually the layers with the largest parameters and dense connections, and Global Average Pooling provides a lower cost way to achieve similar results. The main idea of ​​Global Average Pooling is to generate an average value from each last layer feature map as a confidence factor for the score, which is directly fed into the softmax layer.

[0295] Global average pooling has three benefits: (1) there are no extra parameters in the global average pooling layer, so overfitting is avoided at the global average pooling layer; (2) since the output of global average pooling is the average of the entire feature map, global average pooling will be more robust to spatial translation; and (3) due to the large number of parameters in the fully connected layers, which usually require more than 50% of all the parameters in the entire network, replacing them with global average pooling layers can significantly reduce the size of the model, and this makes global average pooling very useful in model compression.

[0296] Global average pooling makes sense because stronger features in the last layer are expected to have higher averages. In some implementations, global average pooling can be used as a proxy for the classification score. Feature maps under global average pooling can be interpreted as confidence maps and enforce correspondence between feature maps and categories. Global average pooling can be particularly effective if the features of the last layer are abstract enough for direct classification; however, if multi-level features should be combined into groups such as part models, global average pooling alone is not enough, which is best performed by adding a simple fully connected layer or other classifier after global average pooling.

[0297] Deep Learning in Genomics

[0298] Genetic variation can help explain many diseases. Every person has a specific genetic code, and many genetic variants exist within a group of individuals. Most of the bad genetic variants have been removed from the genome by natural selection. It is important to determine which genetic variants may be pathogenic or harmful. This will help researchers focus on genetic variants that may cause disease and speed up the diagnosis and treatment of many diseases.

[0299] Modeling the properties and functional roles of variants, such as pathogenicity, is an important but challenging task in the field of genomics. Despite the rapid development of functional genomic sequencing technology, the interpretation of the functional consequences of variants remains a huge challenge due to the complexity of cell type-specific transcriptional regulatory systems.

[0300] Advances in biochemical technology over the past few decades have produced next-generation sequencing (NGS) platforms that can quickly generate genomic data at a much lower cost than before. Such extremely large amounts of sequenced DNA remain difficult to annotate. Supervised machine learning algorithms generally work well when large amounts of labeled data are available. In bioinformatics and many other data-rich disciplines, the process of labeling instances is expensive; however, unlabeled instances are cheap and readily available. For situations where the amount of labeled data is relatively small and the amount of unlabeled data is significantly larger, semi-supervised learning is a cost-effective alternative to manual labeling.

[0301] There is an opportunity to use semi-supervised algorithms to build deep learning-based pathogenicity classifiers that can accurately predict the pathogenicity of variants. It may be possible to obtain a database of pathogenic variants that is free of human confirmation bias.

[0302] Regarding pathogenicity classifiers, deep neural networks are a class of artificial neural networks that use multiple layers of complex transformations with nonlinearities to successfully model high-level features. Deep neural networks provide feedback via backpropagation, which carries the difference between observed and predicted outputs to adjust parameters. Deep neural networks have evolved with the availability of large-scale training datasets, the power of parallel and distributed computing, and sophisticated training algorithms. Deep neural networks have led to significant advances in many fields such as computer vision, speech recognition, and natural language processing.

[0303] Convolutional Neural Networks (CNNs) and Recurrent Neural Networks (RNNs) are components of deep neural networks. Convolutional neural networks have been particularly successful in image recognition and their architecture includes convolutional layers, nonlinear layers, and pooling layers. Recurrent neural networks aim to exploit the sequential information of input data through circular connections between building blocks such as perceptrons, long short-term memory units, and gated recurrent units. In addition, many other emerging deep neural networks have been proposed in limited contexts, such as deep spatiotemporal neural networks, multidimensional recurrent neural networks, and convolutional autoencoders.

[0304] The goal of training a deep neural network is to optimize the weight parameters of each layer, gradually combining simpler features into complex features, so that the most suitable hierarchical representation can be learned from the data. A single loop of the optimization process is organized as follows. First, given a training dataset, the forward pass sequentially computes the output of each layer and propagates the function signal forward through the network. In the final output layer, the objective loss function measures the error between the inferred output and the given label. To minimize the training error, the backward pass backpropagates the error signal using the chain rule and calculates the gradient of all weights in the neural network. Finally, the weight parameters are updated using an optimization algorithm based on stochastic gradient descent. While batch gradient descent performs parameter updates for each complete dataset, stochastic gradient descent provides a stochastic approximation by performing updates for each small set of data instances. Several optimization algorithms are derived from stochastic gradient descent. For example, the Adagrad and Adam training algorithms perform stochastic gradient descent while adaptively changing the learning rate according to the update frequency and gradient moment of each parameter, respectively.

[0305] Another core element of deep neural network training is regularization, which refers to strategies designed to avoid overfitting and thus achieve good generalization performance. For example, weight decay adds a penalty term to the objective loss function, causing the weight parameters to converge to smaller absolute values. Dropout randomly removes hidden units from a neural network during training and can be thought of as a collection of possible subnetworks. To improve dropout capabilities, a new activation function, Maxout, and a dropout variant for recurrent neural networks called rnnDrop have been proposed. In addition, batch normalization methods provide a new regularization method by normalizing the scalar features of each activation within a mini-batch and learning with each mean and variance as parameters.

[0306] Given that sequencing data is multidimensional and high-dimensional, deep neural networks have great potential for application in bioinformatics research due to their wide applicability and high predictive power. Convolutional neural networks have been used to solve sequence-based problems in genomics, such as motif discovery, pathogenic variant identification, and gene expression inference. Convolutional neural networks use a weight-sharing strategy that is particularly useful for studying DNA because it can capture sequence motifs, which are short, recurring local patterns in DNA that are presumed to have significant biological functions. The hallmark of convolutional neural networks is the use of convolutional filters. Unlike traditional classification methods that are based on carefully designed and artificially crafted features, convolutional filters perform adaptive feature learning, similar to the process of mapping raw input data into an information representation of knowledge. In this sense, convolutional filters act as a series of motif scanners, as a set of such filters are able to identify relevant patterns in the input and update themselves during the training procedure. Recurrent neural networks can capture long-range dependencies in sequence data of varying lengths, such as protein or DNA sequences.

[0307] Therefore, powerful computational models that predict pathogenicity of variants would be of great benefit to both basic science and translational research.

[0308] Common polymorphisms represent natural experiments whose fitness has been tested by natural selection over many generations. By comparing the allele frequency distributions of missense and synonymous substitutions in humans, we found that the presence of a missense variant at high allele frequency in nonhuman primate species reliably predicts that the variant is also under neutral selection in human populations. In contrast, common variants in more distant species are under negative selection as evolutionary distance increases.

[0309] We used common variation from six nonhuman primate species to train a semisupervised deep learning network that accurately classifies clinically de novo missense mutations using sequence alone. Based on more than 500 known species, the primate lineage contains enough common variation to systematically model the effects of most human variants of unknown significance.

[0310] The human reference genome contains more than 70 million potential protein-altering missense substitutions, the vast majority of which are rare mutations whose effects on human health have not yet been characterized. These variants of unknown significance pose a challenge to genomic interpretation in clinical applications and have long been an obstacle to population screening sequencing and personalized medicine.

[0311] Categorizing common variation across human populations is an effective strategy for identifying clinically benign variants, but the common variation available in modern humans is limited by ancient bottleneck events in our species. Humans and chimpanzees share 99% sequence identity, suggesting that natural selection on chimpanzee variants has the potential to model the effects of identical variants in humans. The average coalescence time of allelic action polymorphisms in human populations is a fraction of the time since species divergence, so naturally occurring chimpanzee variation largely explores mutational space that does not overlap with human variation, except in rare cases where haplotypes are maintained by balancing selection.

[0312] Recently available pooled exome data from 60,706 individuals allowed us to test this hypothesis by comparing the allele frequency spectra of missense and synonymous variants. Singleton variants in ExAc closely matched the expected 2.2:1 missense:synonymous ratio predicted from de novo mutations after adjusting for mutation rates using a trinucleotide background, but at higher allele frequencies the number of missense variants observed was reduced by filtering out deleterious variants through natural selection. The pattern of missense:synonymous ratios across the spectrum of alleles suggests that the majority of missense variants with a population frequency <0.1% are mildly deleterious, that is, neither pathogenic enough to warrant immediate elimination from the population nor neutral enough to allow presence at high allele frequencies, consistent with previous observations from more limited population data. These findings support the widespread empirical practice of filtering out alleles greater than 0.1%–1% in diagnostic laboratories, despite well-documented expectations arising from balancing selection and founder effects, because of the potential benign nature of an overtly inherited disorder.

[0313] Repeating the analysis with a subset of human variants that have the same status as common chimpanzee variants (observed more than once in chimpanzee population sequencing), we found that missense:synonymous ratios are roughly constant across the spectrum of allele frequencies. The high allele frequencies of these variants in chimpanzee populations suggest that they have been subjected to natural selection in chimpanzees, and their neutral effects on human population fitness provide compelling evidence that the selection pressure for missense variation is highly consistent in both species. The lower missense:synonymous ratios observed in chimpanzees are consistent with the larger effective population size in ancestral chimpanzee populations, which could more effectively filter out mildly deleterious variants.

[0314] In contrast, rare chimpanzee variants (observed only once in sequenced chimpanzee populations) showed a moderate decrease in missense:synonymous ratios at higher allele frequencies. Simulating a cohort of the same size from human variation data, we estimate that only 64% of variants observed once in a cohort of this size have an allele frequency greater than 0.1% in the general population, compared to 99.8% of variants observed multiple times in this cohort, suggesting that not all rare chimpanzee variants are subject to selection. In summary, we estimate that 16% of the confirmed chimpanzee missense variants have an allele frequency less than 0.1% in the general population and are negatively selected at higher allele frequencies.

[0315] We next characterized human variants with variant states identical to those seen in other nonhuman primate species (bonobos, gorillas, orangutans, rhesus monkeys, and marmosets). Similar to chimpanzees, we observed missense:synonymous ratios to be roughly equal across the range of allele frequencies, with the exception of a slight attrition of missense variants at high allele frequencies, which is expected due to the inclusion of a small number of rare variants (approximately 5-15%). These results imply that selection for missense variants is largely uniform within the primate lineage, at least in New World macaques, which are estimated to have diverged from the human ancestral lineage approximately 35 million years ago.

[0316] Human missense variants with the same status as variants in other primates were strongly enriched for benign results in ClinVar. After excluding variants with unknown or conflicting annotations, we found that human variants with primate homologs were approximately 95% likely to be annotated as benign or likely benign in ClinVar, compared with 45% of missense variants in general. The small fraction of ClinVar variants from nonhuman primates that were classified as pathogenic was similar to the fraction of pathogenic ClinVar variants observed by confirming rare variants from a similarly sized healthy human cohort. A large fraction of these variants were annotated as pathogenic or likely pathogenic, suggesting that they received their classification before the advent of large allele frequency databases and may be curated differently currently.

[0317] The field of human genetics has long relied on model organisms to infer the clinical impact of human mutations, but the long evolutionary distance from most genetically traceable models has raised concerns about the extent to which these findings generalize to humans. To examine the consistency of natural selection on missense variants in humans and more distant species, we extended our analyses beyond the primate lineage to include common variants derived from four additional mammalian species (mouse, pig, goat, cattle) and two more distant vertebrate species (chicken, zebrafish) in general. In contrast to previous primate analyses, we found that missense variants were significantly depleted at common allele frequencies compared with rare allele frequencies, especially at larger evolutionary distances, suggesting that a large fraction of common missense variants in more distant species would be under negative selection in human populations. Nonetheless, the missense variants observed in more distant vertebrates still raise the possibility of a benign outcome, as the fraction of common missense variants attributable to natural selection is much lower than the ~50% attrition of human missense variants at baseline. Consistent with these results, we found that approximately 85% of human missense variants observed in mouse, dog, pig, and cow were annotated as benign or likely benign in ClinVar, compared with 95% for primate variants and 45% for the ClinVar database as a whole.

[0318] The presence of closely related species pairs at different evolutionary distances also provides an opportunity to assess the functional consequences of fixed missense substitutions in the human population. In closely related species pairs on the mammal family tree (branch length < 0.1), we found that fixed missense variants were depleted at common allele frequencies compared to rare allele frequencies, indicating that in humans and even within the primate lineage, a large fraction of interspecies fixed substitutions will be non-neutral. Comparison of the extent of missense depletion showed that interspecies fixed substitutions were significantly less neutral than intraspecies polymorphisms. Interestingly, interspecies variants between closely related mammals were not substantially more pathogenic in ClinVar compared to common intraspecies polymorphisms (83% were annotated as likely benign or probably benign), suggesting that these changes do not confer loss of protein function but rather reflect modulation of protein function, conferring species-specific adaptive advantages.

[0319] The large number of possible variants of unknown significance and the critical importance of accurate variant classification for clinical applications have prompted many attempts to use machine learning to address the problem, but these efforts have been largely limited by the insufficient number of common human variants and the questionable quality of annotations in curated databases. Variants from six nonhuman primates contribute more than 300,000 unique missense variants that do not overlap with common human variation and are largely benign, greatly expanding the size of the training dataset that can be used for machine learning methods.

[0320] Unlike earlier models that employed a large number of human-engineered features and meta-classifiers, we applied a simple deep learning residual network that takes as input only the amino acid sequences flanking the target variant and alignments of orthologous sequences from other species. To feed the network with information relevant to protein structure, we trained two separate networks to learn secondary structure and solvent accessibility from sequence alone, and used these as subnetworks of a larger deep learning network to predict effects on protein structure. Starting with sequence avoids potential biases in annotation of protein structure and functional domains, which may not be fully confirmed or consistently applied.

[0321] We use semi-supervised learning to overcome the problem of the training set containing only variants with benign labels by initially training an ensemble of networks to distinguish likely benign primate variants from random unknown variants matched in mutation rate and sequencing coverage. This ensemble of networks is used to score the entire set of unknown variants and influence the selection of unknown variants for the next iteration of the seeding classifier by biasing them towards those with greater disease predictive consequence, taking incremental steps at each iteration to prevent the model from converging prematurely to suboptimal results.

[0322] Common primate variants also provide a clean validation dataset for evaluating existing methods that is completely independent of previously used training data (which is difficult to evaluate objectively due to the proliferation of meta-classifiers). We evaluate the performance of our model and four other popular classification algorithms (Sift, Polyphen2, CADD, M-CAP) using 10,000 fixed primate common variants. Since approximately 50% of human missense mutations can be removed by natural selection at common allele frequency, we calculated the 50th percentile score of each classifier on a set of randomly selected missense variants (matched in mutation rate to the 10,000 fixed primate common variants) and used this threshold to evaluate primate common variants. Our deep learning models significantly outperformed other classifiers in accuracy on this independent validation dataset, using either deep learning networks trained only on human common variants or using both human common variants and primate variants.

[0323] Recent trio sequencing studies have catalogued thousands of de novo mutations in patients with neurodevelopmental disorders and their healthy siblings, allowing us to assess the strength of various classification algorithms at separating de novo missense mutations in cases vs controls. For each of the four classification algorithms, we scored de novo missense variants in cases compared with controls and report Wilcoxon rank-sum test p-values ​​for the difference between the two distributions, demonstrating that a deep learning approach trained on primate variants (p~10-33) performed significantly better than other classifiers (p~10-13 to 10-19) in this clinical setting. Based on the ~1.3-fold enrichment of de novo missense variants reported previously for this cohort and previous estimates that ~20% of missense variants produce loss-of-function effects, we would expect a perfect classifier to separate the two classes with a p-value of 10-40, indicating that our classifiers still have room for improvement.

[0324] The accuracy of the deep learning classifier scaled with the size of the training dataset, with variant data from each of the six primate species independently contributing to the improvement of the classifier's accuracy. The large number and diversity of extant nonhuman primate species, and evidence that selective pressures on protein-altering variants are largely consistent across the primate lineage, suggest that systematic primate population sequencing is an effective strategy for cataloging the millions of human variants of unknown significance that currently limit clinical genomic interpretation. Of the 504 known nonhuman primate species, approximately 60% are facing extinction due to hunting and habitat loss, prompting an urgent need for conservation efforts worldwide to benefit these unique, irreplaceable species and ourselves.

[0325] While not as much aggregate whole genome data is available as exome data, limiting detection of the effects of natural selection in deep intronic regions, we were also able to calculate observed vs expected counts of cryptic splicing variants far from exonic regions. Overall, we observed a 60% depletion in cryptic splicing variants at a distance >50 nt from exon-intron boundaries. The attenuated signal is likely a combination of smaller sample sizes with whole genome data compared to exomes and greater difficulty in predicting the effects of deep intronic variants.

[0326] the term

[0327] All documents and similar materials cited in this application, including but not limited to patents, patent applications, articles, books, papers, and web pages, regardless of the format of these documents and similar materials, are expressly incorporated by reference in their entirety. If one or more of the incorporated documents and similar materials differs from or contradicts this application, including but not limited to defined terms, term usage, described techniques, etc., this application shall prevail.

[0328] As used herein, the following terms have the indicated meanings.

[0329] Base refers to a nucleotide base or nucleotide, A (adenine), C (cytosine), T (thymine), or G (guanine).

[0330] The terms "protein" and "translational sequence" are used interchangeably herein.

[0331] The terms "codon" and "base triplet" are used interchangeably herein.

[0332] The terms "amino acid" and "translation unit" are used interchangeably herein.

[0333] This application uses the phrases "variant pathogenicity classifier," "convolutional neural network based classifier for variant classification," and "deep convolutional neural network based classifier for variant classification" interchangeably.

[0334] The term "chromosome" refers to the genetic carrier of living cells, which is derived from chromatin chains containing DNA and protein components (especially histones). This article adopts the conventional internationally recognized chromosome coding system for human individual genomes.

[0335] The term "site" refers to a unique position on a reference genome (e.g., chromosome ID, chromosome position and orientation). In some implementations, a site can be the position of a fragment on a residue, sequence tag, or sequence. The term "locus" can be used to refer to a specific position of a nucleic acid sequence or polymorphism on a reference chromosome.

[0336] The term "sample" herein refers to a sample of a biological fluid, cell, tissue, organ or organism that is usually derived from a mixture of nucleic acids containing nucleic acids or containing at least one nucleic acid sequence to be sequenced and / or performed in stages. Such samples include, but are not limited to, sputum / oral fluid, amniotic fluid, blood, blood fractions, fine needle aspiration biopsy samples (e.g., surgical biopsy, fine needle aspiration biopsy, etc.), urine, peritoneal fluid, pleural fluid, tissue explants, organ cultures and any other tissue or cell preparations, or fractions or derivatives thereof or separated therefrom. Although samples are usually taken from human subjects (e.g., patients), samples can be taken from any organism with chromosomes, including but not limited to dogs, cats, horses, goats, sheep, cattle, pigs, etc. The sample can be used directly as it is obtained from a biological source or after a pretreatment that modifies the sample properties. For example, such a pretreatment may include preparing plasma from blood, diluting viscous fluids, etc. The pretreatment method may also include, but is not limited to, filtration, precipitation, dilution, distillation, mixing, centrifugation, freezing, freeze-drying, concentration, amplification, nucleic acid fragmentation, inactivation of interfering components, reagent addition, cracking, etc.

[0337] The term "sequence" includes or represents a nucleotide chain coupled to each other. Nucleotides can be based on DNA or RNA. It should be understood that a sequence can include multiple subsequences. For example, a single sequence (for example, a PCR amplicon) can have 350 nucleotides. Sample interpretation can include multiple subsequences within these 350 nucleotides. For example, sample interpretation can include a first and a second flanking subsequence with, for example, 20-50 nucleotides. The first and second flanking subsequences can be located on either side of a repeat fragment with a corresponding subsequence (for example, 40-100 nucleotides). Each flanking subsequence can include (or include part of) a primer subsequence (for example, 10-30 nucleotides). For ease of reading, the term "subsequence" will be referred to as a "sequence", but it should be understood that the two sequences are not necessarily separated from each other on a public chain. In order to distinguish the various sequences described herein, different labels (for example, target sequence, primer sequence, flanking sequence, reference sequence, etc.) can be given to the sequence. Other terms such as "allele" can be given different labels to distinguish similar objects.

[0338] The term "paired end sequencing" refers to a sequencing method that sequences both ends of a target fragment. Paired end sequencing can facilitate the detection of genomic rearrangements and duplicate fragments, as well as the detection of gene fusions and novel transcripts. Methods for paired end sequencing are described in PCT Publication WO07010252, PCT Application Serial No. PCTGB2007 / 003798, and U.S. Patent Application Publication US 2009 / 0088327, each of which is incorporated herein by introduction. In one example, a series of operations can be performed as follows; (a) generating a nucleic acid cluster; (b) linearizing the nucleic acid; (c) hybridizing a first sequencing primer, and performing repeated cycles of extension, scanning, and deblocking as described above; (d) "reversing" the target nucleic acid on the surface of the flow cell by synthesizing a complementary copy; (e) linearizing the resynthesized chain; (f) hybridizing a second sequencing primer, and performing repeated cycles of extension, scanning, and deblocking as described above. The conversion operation can be performed by delivering reagents for a single cycle of bridge amplification as described above.

[0339] The term "reference genome" or "reference sequence" refers to any specific known genome sequence of any organism that can be used to reference an identified sequence from a subject, whether partial or complete. For example, reference genomes for human subjects and many other organisms can be found at the National Center for Biotechnology Information (ncbi.nlm.nih.gov). "Genome" refers to the complete genetic information of an organism or virus expressed as a nucleic acid sequence. A genome includes both genes and non-coding sequences of DNA. A reference sequence may be larger than the reading compared to it. For example, it can be at least about 100 times larger, or at least about 1000 times larger, or at least about 10,000 times larger, or at least about 105 times larger, or at least about 106 times larger, or at least about 107 times larger. In one example, a reference genome sequence is a sequence of a full-length human genome. In another example, a reference genome sequence is limited to a specific human chromosome, such as chromosome 13. In some embodiments, a reference chromosome is a chromosome sequence from a human genome version hg19. Such a sequence may be referred to as a chromosome reference sequence, although the term reference genome is intended to cover such a sequence. Other examples of reference sequences include genomes of other species, and chromosomes, subchromosome regions (e.g., chains) of any species, etc. In various embodiments, the reference genome is a consensus sequence or other combination derived from multiple individuals. However, in some applications, the reference sequence can be taken from a specific individual.

[0340] The term "interpretation" refers to a collection of sequence data describing a fragment of a nucleotide sample or reference. The term "interpretation" may refer to sample interpretation and / or reference interpretation. Typically, although not necessarily, interpretation represents a short sequence of continuous base pairs in a sample or reference. Interpretation may be symbolically represented by a base pair sequence (in ATCG) of a sample or reference fragment. It may be stored in a storage device and appropriately processed to determine whether the interpretation matches the reference sequence or meets other criteria. Interpretation may be obtained directly from a sequencing device, or may be obtained indirectly from stored sequence information associated with the sample. In some cases, interpretation is a DNA sequence with sufficient length (e.g., at least about 25bp) that can be used to identify a larger sequence or region, such as a larger sequence or region that can be aligned and specifically assigned to a chromosome or genomic region or gene.

[0341] Next generation sequencing methods include, for example, sequencing by synthesis technology (Illumina), pyrophosphate sequencing (454), ion semiconductor technology (ion torrent sequencing), single molecule real-time sequencing (Pacific Biosciences), and sequencing by ligation (SOLiD sequencing). Depending on the sequencing method, the length of each read may be about 30bp to more than 10,000bp. For example, the Illumina sequencing method using the SOLiD sequencer produces a nucleic acid read of about 50bp. For another example, ion torrent sequencing produces a nucleic acid read of up to 400bp, while 454 pyrophosphate sequencing produces a nucleic acid read of about 700bp. For another example, a single molecule real-time sequencing method can produce a read of 10,000bp to 15,000bp. Therefore, in certain embodiments, the length of the nucleic acid sequence read is 30-100bp, 50-200bp, or 50-400bp.

[0342] The term "sample interpretation", "sample sequence" or "sample fragment" refers to sequence data of a genomic sequence of interest from a sample. For example, a sample interpretation includes sequence data from a PCR amplicon having forward and reverse primer sequences. The sequence data can be obtained from any selected sequence method. The sample interpretation can be, for example, from a sequencing by synthesis (SBS) reaction, sequencing by ligation, or any other suitable sequencing method that requires determination of the length and / or identity of a repeated sequence. The sample interpretation can be a common (e.g., average or weighted) sequence derived from a plurality of sample interpretations. In certain embodiments, providing a reference sequence includes identifying a locus of interest based on the primer sequence of the PCR amplicon.

[0343] The term "original fragment" refers to sequence data of a portion of a genomic sequence of interest that at least partially overlaps a designated position of interest or a secondary position within a sample interpretation or sample fragment. Non-limiting examples of original fragments include duplex stitched fragments, simplex stitched fragments, duplex non-stitched fragments, and simplex non-stitched fragments. The term "original" is used to indicate that the original fragment includes sequence data that has a certain relationship with the sequence data in the sample interpretation, regardless of whether the original fragment shows supporting variants corresponding to and identifying or confirming the potential variation in the sample interpretation. The term "original fragment" does not mean that the fragment necessarily contains supporting variants that can verify the variant identification in the sample interpretation. For example, when the variant identification application determines that the sample interpretation shows a first variant, the variant identification application can determine that one or more original fragments lack a "supporting" variant of the corresponding type, which can otherwise be expected to occur under the conditions of a given variant in the sample interpretation.

[0344] The term "mapping", "aligned", "alignment" or "aligning" refers to the process of comparing a read or label to a reference sequence to determine whether the reference sequence contains the read sequence. If the reference sequence contains an interpretation, the interpretation can be mapped to the reference sequence, or in some implementations, can be mapped to a specific position in the reference sequence. In some cases, the alignment simply tells whether the interpretation is a member of a specific reference sequence (i.e., whether the interpretation exists or does not exist in the reference sequence). For example, an alignment of an interpretation with a reference sequence of human chromosome 13 will indicate whether the interpretation exists in the reference sequence of chromosome 13. A tool that provides this information can be referred to as a set membership tester. In some cases, the alignment also indicates the position in the reference sequence to which the interpretation or label is mapped. For example, if the reference sequence is the entire human genome sequence, the alignment can indicate that an interpretation exists on chromosome 13, and can further indicate that the interpretation is on a specific chain and / or site of chromosome 13.

[0345] The term "indel" refers to the insertion and / or deletion of bases in the DNA of an organism. Micro-indels refer to indels that result in a net change of 1 to 50 nucleotides. In the coding regions of the genome, unless the length of the indel is a multiple of 3, it will produce a frameshift mutation. Indels can be contrasted with point mutations. Indels insert and delete nucleotides from a sequence, while point mutations are a form of substitution that replaces one nucleotide without changing the total number of DNA. Indels can also be contrasted with tandem base mutations (TBMs), which can be defined as substitutions at adjacent nucleotides (mostly substitutions at two adjacent nucleotides, but substitutions at three adjacent nucleotides have been observed).

[0346] The term "variant" refers to a nucleic acid sequence different from a nucleic acid reference. Typical nucleic acid sequence variants include, but are not limited to, single nucleotide polymorphisms (SNPs), short deletions and insertion polymorphisms (Indels), copy number variations (CNVs), microsatellite markers or short tandem repeats and structural variations. Somatic variation recognition is the work of identifying variations that are present at low frequencies in DNA samples. In the context of cancer treatment, somatic variation recognition is of interest. Cancer is caused by the accumulation of DNA mutations. DNA samples from tumors are usually heterogeneous, including some normal cells, some cells in the early stages of cancer progression (less mutations) and some late cells (more mutations). Due to this heterogeneity, when tumors are sequenced (e.g., from FFPE samples), somatic mutations will usually occur at a lower frequency. For example, SNVs may only be visible in 10% of the readings covering a given base. Variants classified as somatic or germline by variant classifiers are also referred to herein as "tested variants."

[0347] The term “noise” refers to erroneous variant calls resulting from one or more errors in the sequencing process and / or the variant calling application.

[0348] The term "variant frequency" refers to the relative frequency of alleles (genetic variation) at a specific locus in a population, expressed as a fraction or percentage. For example, the fraction or percentage can be the fraction of all chromosomes in the population carrying the allele. For example, the sample variant frequency represents the relative frequency of alleles / variants at a specific locus / position along the genomic sequence of interest on a "population", and the "population" corresponds to the number of interpretations and / or samples obtained for the genomic sequence of interest from an individual. As another example, the baseline variant frequency represents the relative frequency of alleles / variants at a specific locus / position along one or more baseline genomic sequences, wherein the "population" corresponds to the number of interpretations and / or samples obtained from one or more baseline genomic sequences of a normal individual population.

[0349] The term "variant allele frequency (VAF)" refers to the observed percentage of sequencing reads matching a variant divided by the total coverage of the target position. VAF is a measure of the proportion of sequencing reads that carry a variant.

[0350] The terms "position", "designated position" and "locus" refer to the position or coordinates of one or more nucleotides within a nucleotide sequence. The terms "position", "designated position" and "locus" also refer to the position or coordinates of one or more base pairs in a nucleotide sequence.

[0351] The term "haplotype" refers to a combination of alleles at adjacent sites on a chromosome that are inherited together. A haplotype can be a locus, multiple loci, or an entire chromosome, depending on the number of recombination events that occur between a given set of loci (if any).

[0352] The term "threshold" refers herein to a numerical or non-numerical value used as a critical value to characterize a sample, nucleic acid, or a portion thereof (e.g., interpretation). The threshold value may vary based on empirical analysis. The threshold value may be compared to a measured value or a calculated value to determine whether the source generating the value should be classified in a particular manner. The threshold value may be identified empirically or analytically. The selection of the threshold value depends on the confidence that the user wishes to classify. The threshold value may be selected for a specific purpose (e.g., to balance sensitivity and selectivity). As used herein, the term "threshold value" represents a point at which an analysis process may be changed and / or a point at which an action may be triggered. The threshold value does not need to be a predetermined number. Instead, the threshold value may be, for example, a function based on multiple factors. The threshold value may be adapted to the environment. In addition, the threshold value may indicate a range between an upper limit, a lower limit, or a limit.

[0353] In some implementations, a metric or score based on sequencing data can be compared to a threshold. As used herein, the term "metric" or "score" can include a value or result determined from sequencing data, or can include a function based on a value or result determined from sequencing data. Like a threshold, a metric or score can adapt to the environment. For example, a metric or score can be a normalized value. As an example of a score or metric, one or more implementations can use a count score when analyzing data. The count score can be based on the number of sample interpretations. The sample interpretation may have undergone one or more filtering stages so that the sample interpretation has at least one common feature or quality. For example, each sample interpretation for determining the count score may have been compared to a reference sequence, or may be designated as a potential allele. The number of sample interpretations with common features can be counted to determine the interpretation count. The count score can be based on the interpretation count. In some embodiments, the count score can be a value equal to the interpretation count. In other embodiments, the count score can be based on the interpretation count and other information. For example, the count score can be based on the interpretation count of a specific allele of a genetic locus and the total number of interpretations of a genetic locus. In some implementations, the count score can be based on the data previously obtained from the interpretation count and the genetic locus. In some implementations, the count score can be a normalized score between predetermined values. The count score can also be a function of the reading count from other loci of the sample, or a function of the reading count of other samples that are run simultaneously with the sample of interest. For example, the count score can be a function of the reading count of other loci in the reading count and sample and / or the reading count from other samples of a specific allele. As an example, the reading count from other loci and / or the reading count from other samples can be used for the count score of normalization specific allele.

[0354] The term "coverage" or "fragment coverage" refers to a count or other measure of multiple sample interpretations of the same fragment of a sequence. The reading count can represent a count of the number of readings covering the corresponding fragment. Alternatively, the coverage can be determined by multiplying the reading count by a specified factor based on historical knowledge, sample knowledge, locus knowledge, etc.

[0355] The term "read depth" (usually a number followed by an "×") refers to the number of sequencing reads with overlapping alignments at a target position. This is usually expressed as an average or percentage exceeding a cutoff value over a set of time intervals (e.g., exons, genes, or panels). For example, a clinical report might say that the panel average coverage is 1,105-fold, with 98% of target bases covered more than 100-fold.

[0356] The term "base call quality score" or "Q score" refers to a PHRED-scaled probability ranging from 0-20, with the range being inversely proportional to the probability of a single sequenced base being correct. For example, a T base call with a Q of 20 might be considered correct with a confidence P value of 0.01. Any base call with a Q < 20 should be considered low quality, and any variant that identifies a low quality that supports a significant portion of the sequencing reads for that variant should be considered a potential false positive.

[0357] The term "variant call" or "number of variant call" refers to the number of sequencing reads that support the presence of a variant.

[0358] Sequencing process

[0359] The embodiment set forth herein may be applicable to analyzing nucleic acid sequences to identify sequence variations. Implementations may be used to analyze potential variants / alleles of genetic position / locus, and determine the genotype of the genetic locus, or in other words, provide genotype identification for the locus. For example, nucleic acid sequences may be analyzed according to the methods and systems described in U.S. Patent Application Publication No. 2016 / 0085910 and U.S. Patent Application Publication No. 2013 / 0296175, all of which are expressly incorporated herein by reference in their entirety.

[0360] In one implementation, the sequencing process includes receiving a sample that contains or is suspected of containing a nucleic acid, such as DNA. The sample may be from a known or unknown source, such as an animal (e.g., human), a plant, a bacterium, or a fungus. The sample may be obtained directly from a source. For example, blood or saliva may be taken directly from an individual. Alternatively, the sample may not be obtained directly from a source. Then, one or more processors direct the system to prepare the sample for sequencing. Preparation may include removing foreign matter and / or separating certain substances (e.g., DNA). Biological samples may be prepared to include features for a specific assay. For example, biological samples may be prepared for synthetic sequencing (SBS). In certain embodiments, preparation may include amplifying certain regions of the genome. For example, preparation may include amplifying known predetermined genetic loci including STRs and / or SNPs. The genetic loci may be amplified using predetermined primer sequences.

[0361] Next, one or more processors guide the system to sequence the sample. Sequencing can be performed by a variety of known sequencing schemes. In a specific implementation, sequencing includes SBS. In SBS, a plurality of fluorescently labeled nucleotides are used to sequence the DNA clusters (possibly millions of clusters) present in a plurality of amplifications on an optical substrate surface (e.g., a surface that at least partially defines a channel in a flow cell). The flow cell may include a nucleic acid sample for sequencing, wherein the flow cell is placed in an appropriate flow cell support.

[0362] Nucleic acids can be prepared so that they contain known primer sequences adjacent to unknown target sequences. In order to start the first SBS sequencing cycle, one or more differently labeled nucleotides and DNA polymerases, etc. can flow into / flow through the flow cell through a fluid flow subsystem. One type of nucleotide can be added at a time, or the nucleotides used in the sequencing process can be specially designed to have reversible termination characteristics, so that each cycle of the sequencing reaction occurs simultaneously in the presence of several types of labeled nucleotides (e.g., A, C, T, G). Nucleotides can include detectable labeling parts, such as fluorophores. In the case of four nucleotides mixed together, polymerase can select the correct base for integration, and each sequence can extend one base. Unbound nucleotides can be washed away by flowing washing solution through the flow cell. One or more lasers can excite nucleic acids and induce fluorescence. The fluorescence emitted from nucleic acids is based on the fluorophores incorporated into the bases, and different fluorophores can emit light of different wavelengths. A deblocking agent can be added to the flow cell to remove the reversible terminator group from the DNA chain extended and detected. Then, the deblocking agent can be washed away by flowing washing solution through the flow cell. Then, the flow cell is ready for further sequencing cycles, starting with the introduction of the above-mentioned labeled nucleotides. The fluid and detection operations can be repeated several times to complete the sequencing run. Exemplary sequencing methods are described in, for example, Bentley et al., Nature 456:53-59 (2008), International Publication No. WO 04 / 018497; and U.S. Patent No. 7,057,026; International Publication No. WO 91 / 06678; International Publication No. WO 07 / 123744; U.S. Patent No. 7,329,492; U.S. Patent No. 7,211,414; U.S. Patent No. 7,315,019; U.S. Patent No. 7,405,281 and U.S. Patent Application Publication No. 2008 / 0108082, each of which is incorporated herein by reference.

[0363] In some embodiments, nucleic acids can be attached to a surface and amplified before or during sequencing. For example, bridge amplification can be used to amplify to form nucleic acid clusters on a surface. Useful bridge amplification methods are described in, for example, U.S. Patent No. 5,641,658; U.S. Patent Application Publication No. 2002 / 0055100; U.S. Patent No. 7,115,400; U.S. Patent Application Publication No. 2004 / 0096853; U.S. Patent Application Publication No. 2004 / 0002090; U.S. Patent Application Publication No. 2007 / 0128624; and U.S. Patent Application Publication No. 2008 / 0009420, each of which is incorporated herein by reference in its entirety. Another useful method for amplifying nucleic acids on a surface is rolling circle amplification (RCA), for example, as described in Lizardi et al., Nat. Genet. 19: 225-232 (1998), and U.S. Patent Application Publication No. 2007 / 0099208A1, each of which is incorporated herein by reference.

[0364] An exemplary SBS protocol utilizes modified nucleotides with removable 3' blocks, for example, as described in International Publication No. WO 04 / 018497, U.S. Patent Application Publication No. 2007 / 0166705A1, and U.S. Patent No. 7,057,026, each of which is incorporated herein by reference. For example, as a result of a bridge amplification scheme, for example, a repeatedly cycled SBS reagent can be delivered to a flow cell with a target nucleic acid attached thereto. A linearization solution can be used to convert the nucleic acid clusters into a single-stranded form. The linearization solution can include, for example, a restriction endonuclease capable of cutting one strand of each cluster. Other methods of cleavage can be used as an alternative to restriction endonucleases or nicking enzymes, including, among others, chemical cleavage (e.g., cleavage of diol bonds with periodate), cleavage of abasic sites by cleavage with endonucleases (e.g., "USER", such as "USER" supplied by NEB, Ipswich, Mass., USA, part number M5505S), cleavage of ribonucleotides incorporated into the amplification product by exposure to heat or alkali, which otherwise consists of deoxyribonucleotides, photochemical cleavage, or peptide linker cleavage. Following the linearization operation, the sequencing primer can be delivered to a flow cell under conditions for hybridizing the sequencing primer to the target nucleic acid to be sequenced.

[0365] Then, the flow cell can be contacted with an SBS extension reagent having a modified nucleotide (with a removable 3' block) and a fluorescent label, provided that the primer hybridized with each target nucleic acid can be extended by adding a single nucleotide. Only one nucleotide is added to each primer, because once the modified nucleotide is incorporated into the growing polynucleotide chain complementary to the template region being sequenced, there is no free 3'-OH group available to guide further sequence extension, so the polymerase cannot add other nucleotides. The SBS extension reagent can be removed and replaced with a scanning reagent, which is included in a component that protects the sample under radiation excitation. Exemplary components of scanning reagents are described in U.S. Patent Application Publication No. 2008 / 0280773A1 and U.S. Patent Application No. 13 / 018,255, each of which is incorporated herein by reference. The extended nucleic acid can then be fluorescently detected in the presence of a scanning reagent. Once fluorescence is detected, the 3' block can be removed using a decapping reagent suitable for the used blocking group. Exemplary deblocking reagents that can be used for each blocking group are described in WO004018497, US 2007 / 0166705A1 and U.S. Patent No. 7,057,026, which are each incorporated herein by reference. The deblocking agent can be washed off, and the remaining target nucleic acid is hybridized with the extension primer with 3'-OH group, and other nucleotides can be added now. Therefore, the cycle of adding extension reagent, scanning reagent and deblocking reagent can be repeated, and optional washing is carried out between one or more operations until the desired sequence is obtained. When each modified nucleotide has different marks (known to correspond to specific bases) connected thereto, a single extension reagent delivery operation can be used in each cycle to carry out the above-mentioned cycle. Different marks promote the distinction between the nucleotides added during each incorporation operation. Alternatively, each cycle can include a separate operation of extending the reagent delivery, followed by a separate operation of scanning reagent delivery and detection, in which case two or more nucleotides can have the same mark and can be distinguished based on a known delivery order.

[0366] Although sequencing operations have been discussed above with respect to a particular SBS protocol, it will be appreciated that other protocols for sequencing any of a variety of other molecular analyses may be performed as desired.

[0367] Then, one or more processors of the system receive the sequencing data for subsequent analysis. The sequencing data can be formatted in various ways, such as in .BAM file format. The sequencing data can include, for example, many sample interpretations. The sequencing data can include multiple sample interpretations of corresponding sample sequences with nucleotides. Although only one sample interpretation has been discussed, it should be understood that the sequencing data can include, for example, hundreds, thousands, hundreds of thousands or millions of sample interpretations. Different sample interpretations may have different numbers of nucleotides. For example, the range of sample interpretations may be between 10 nucleotides and about 500 nucleotides or more. The sample interpretation may span the entire genome of the source (multiple). As an example, the sample interpretation is directly for a predetermined genetic locus, such as those genetic loci with suspicious STRs or suspicious SNPs.

[0368] Each sample interpretation can include nucleotide sequence, which can be called sample sequence, sample fragment or target sequence. Sample sequence can include, for example, primer sequence, flanking sequence and target sequence. The number of nucleotides in the sample sequence can include 30, 40, 50, 60, 70, 80, 90, 100 or more. In some embodiments, one or more samples are interpreted (or sample sequence) and include at least 150 nucleotides, 200 nucleotides, 300 nucleotides, 400 nucleotides, 500 nucleotides or more. In some embodiments, sample interpretation can include more than 1000 nucleotides, 2000 nucleotides or more. Sample interpretation (or sample sequence) can include primer sequence at one end or both ends.

[0369] Next, one or more processors analyze the sequencing data to obtain potential variant identification (multiple) and sample variant identification (multiple) sample variation frequencies. This operation may also be referred to as a variant identification application or a variant identifier. Therefore, a variant identifier identifies or detects a variant, and a variant classifier classifies the detected variant as a somatic cell or germline. Alternative variant identifiers may be utilized according to the implementation of this article, wherein different variant identifiers may be used based on the type of sequencing operation performed, based on sample features of interest, etc. A non-limiting example of a variant identification application, such as the PiscesTM application of Illumina Inc. (San Diego, CA), is hosted at https: / / github.com / Illumina / Pisces and described in the text Dunn, Tamsen & Berry, Gwenn & Emig-Agius, Dorothea & Jiang, Yu & Iyer, Anita & Udar, Nitin & Michael. (2017). Pisces: an accurate and versatile single-sample somatic and germline variant caller 595-595.10.1145 / 3107411.3108203, the entire contents of which are expressly incorporated herein by reference.

[0370] Such a variant identification application may include four modules that are executed sequentially:

[0371] (1) Pisces read stitcher: Reduces noise by stitching paired reads in a BAM (read one and read two of the same molecule) into a common read. The output is a stitched BAM.

[0372] (2) Pisces variant caller: identifies small SNVs, insertions, and deletions. Pisces includes a variant folding algorithm to merge variants that are broken down by interpretation boundaries, a basic filtering algorithm, and a simple Poisson-based variant confidence scoring algorithm. The output is VCF.

[0373] (3) Pisces Variant Quality Recalibrator (VQR): The VQR step will reduce the variant Q score of a suspected variant call if the variant call irresistibly follows a pattern associated with thermal damage or FFPE deamination. The output is an adjusted VCF.

[0374] (4) Pisces variant phaser (Scylla) uses a post-read greedy clustering approach to pack small variants from clonal subpopulations into complex alleles. This allows downstream tools to more accurately determine functional consequences. The output is an adjusted VCF.

[0375] Additionally or alternatively, the operation can utilize the variant calling application Strelka hosted by Illumina Inc. at https: / / github.com / Illumina / strelka and described in the article T Saunders, Christopher & Wong, Wendy & Swamy, Sajani & Becq, Jennifer & J Murray, Lisa & Cheetham, Keira. (2012) TMApplications. Strelka: Accurate somatic small variant identification from sequenced tumor-normal sample pairs. Bioinformatics (Oxford, England). 28.1811-7.10.1093 / bioinformatics / bts271, the entire subject matter of which is expressly incorporated herein by reference in its entirety. In addition or alternatively, the operation can utilize Illumina Inc. hosted at https: / / github.com / Illumina / strelka and in the articles Kim, S., Scheffler, K., Halpern, AL, Bekritsky, MA, Noh, E., M., Chen, X., Beyter, D., Krusche, P., and Saunders, C. T. (2017) TM Strelka2: Fast and Accurate Variant Calling for Clinical Sequencing Applications, the entire subject matter of which is expressly incorporated herein by reference in its entirety. In addition or alternatively, the operation can utilize variant annotation / calling tools such as Nirvana hosted at https: / / github.com / Illumina / Nirvana / wiki by Illumina Inc. and described in the article Stromberg, Michael & Roy, Rajat & Lajugie, Julien & Jiang, Yu & Li, Haochen & Margulies, Elliott. (2017) TM Appl. Nirvana: Clinical Grade Variant Annotator. 596-596.10.1145 / 3107411.3108204, the entire contents of which are expressly incorporated herein by reference.

[0376] Such variant annotation / calling tools may apply different algorithmic techniques, such as those disclosed in Nirvana:

[0377] a. Identify all overlapping transcripts using interval arrays: For functional annotation, we can identify all transcripts that overlap with a variant and can use an interval tree. However, since the set of intervals can be static, we are able to further optimize this into an interval array. An interval tree returns all overlapping transcripts in O(min(n,k lg n)) time, where n is the number of intervals in the tree and k is the number of overlapping intervals. In practice, since k is really small compared to n for most variants, the effective runtime on the interval tree is O(k lg n). We improve this to O(lg n+k) by creating an interval array that stores all intervals in a sorted array, so we only need to find the first overlapping interval and then enumerate the remaining (k-1).

[0378] b.CNVs / SVs (Yu): Annotations for “copy number variations” and “structural variants” can be provided. Similar to the annotation of small variants, transcripts overlapping with SVs and previously reported structural variants can be annotated in online databases. Unlike small variants, not all overlapping transcripts need to be annotated, because too many transcripts will overlap with large SVs. Instead, all overlapping transcripts belonging to partially overlapping genes can be annotated. Specifically, for these transcripts, the affected introns, exons, and the consequences caused by the structural variant can be reported. An option is provided to allow the output of all overlapping transcripts, but basic information of these transcripts can be reported, such as gene symbols, whether the marker overlaps with the transcript canonical or partially overlaps. For each SV / CNV, it is also interesting to know whether these variants have been studied and their frequency in different populations. Therefore, we report overlapping SVs in external databases such as 1000 genomes, DGV, and ClinGen. In order to avoid using arbitrary cutoffs to determine which SV overlaps, all overlapping transcripts can be used instead, and the mutual overlap can be calculated, that is, the overlap length divided by the minimum of the two SV lengths.

[0379] c. Reporting Supplemental Annotations: There are two types of supplemental annotations: small and structural variants (SVs). SVs can be modeled as intervals, and the interval array discussed above is used to identify overlapping SVs. Small variants are modeled as points and matched by position and (optionally) allele. As such, they are searched using an algorithm similar to binary search. Since the supplemental annotation database can be large, a smaller index is created to map chromosomal positions to file locations where the supplemental annotations reside. The index is a sorted array of objects (consisting of chromosomal positions and file positions) that can be binary searched using the positions. To keep the index size small, multiple positions (up to some maximum count) are compressed into a single object that stores the value of the first position and only the increments for subsequent positions. Since we use binary search, the running time is O(lg n), where n is the number of items in the database.

[0380] d.VEP cache file

[0381] e. Transcript Database: The transcript cache (cache) and supplementary database (SAdb) files are serialized dumps of data objects (e.g., transcripts and supplementary annotations). We use the Ensembl VEP cache as the data source for the cache. To create the cache, all transcripts are inserted into an interval array and the final state of the array is stored in the cache file. Therefore, during annotation, we only need to load a pre-computed interval array and perform searches on it. Since the cache is loaded into memory and searches are very fast (as described above), finding overlapping transcripts in Nirvana is very fast (configured to be less than 1% of the total run time?).

[0382] f. Supplementary Database: The data sources for SAdb are listed under Supplementary Materials. A small variant of SAdb is produced by performing a k-way merge of all data sources so that each object in the database (identified by a reference name and location) has all relevant supplementary annotations. The problems encountered when parsing the data source files are documented in detail in the Nirvana homepage. To limit memory usage, only the SA index is loaded into memory. This index allows fast lookups of file locations to obtain supplementary annotations. However, since the data must be fetched from disk, adding supplementary annotations has been identified as the largest bottleneck in Nirvana (configured to be 30% of the total runtime).

[0383] g. Consequences and Sequence Ontology: Nirvana's functional annotations (when provided) follow the Sequence Ontology (SO) (http: / / www.sequenceontology.org / ) guidelines. Occasionally, we have the opportunity to identify issues in the current SO and work with the SO team to improve the annotation state.

[0384] Such variant annotation tools can include preprocessing. For example, Nirvana includes a large number of annotations from external data sources, such as ExAC, EVS, 1000 Genomes Project, dbSNP, ClinVar, Cosmic, DGV, and ClinGen. To fully utilize these databases, we had to clean the information from them. We implemented different strategies to deal with different conflicts from different data sources. For example, in the case of multiple dbSNP entries for the same position and alternating alleles, we merged all IDs into a comma-separated list of IDs; if there were multiple entries with different CAF values ​​for the same allele, the first CAF value was used. For conflicting ExAC and EVS entries, we considered the number of sample counts and used the entry with higher sample counts. In the 1000 Genomes Project, we removed the allele frequencies of conflicting alleles. Another issue is incorrect information. We mainly extracted allele frequency information from the 1000 Genomes Project, however, we noticed that for GRCh38, the allele frequencies reported in the info field did not exclude samples for which genotypes were not available, resulting in a reduction in the frequency of variants that were not available for all samples. To ensure the accuracy of the annotation, we used all individual-level genotypes to calculate the true allele frequency. It is well known that the same variant can have different representations based on different alignments. To ensure that we can accurately report the information of the identified variants, we must pre-process the variants from different resources so that they have a consistent representation. For all external data sources, we trimmed the alleles to remove repeated nucleotides in the reference allele and the alternative allele. For ClinVar, we parsed the xml file directly and performed a five-prime alignment for all variants, which is often used in vcf files. Different databases can contain the same set of information. To avoid unnecessary duplication, we removed some duplicate information. For example, we removed variants with the 1000 Genomes Project data source in DGV because we have already reported these variants in 1000 Genomes and provided more detailed information.

[0385] According to at least some implementations, variant identification applications provide identification of low-frequency variants, germline identification, etc. As a non-limiting example, variant identification applications can be run on tumor samples only and / or tumor normal paired samples. Variant identification applications can search for single nucleotide variations (SNVs), multinucleotide variations (MNVs), insertions and deletions, etc. Variant identification applications identify variations while filtering mismatches caused by sequencing or sample preparation errors. For each variant, the variant identifier identifies the reference sequence, the position of the variant, and potential one or more variant sequences (e.g., A to C SNV or AG to A deletion). Variant identification applications identify sample sequences (or sample fragments), reference sequences / fragments, and variant identification to indicate the presence of variants. Variant identification applications can identify the original fragments, and output the names of the original fragments, verify the count of the number of original fragments identified by potential variants, the positions of the original fragments in which the supporting variants occur, and other relevant information. Non-limiting examples of original fragments include duplex stitched fragments, simplex stitched fragments, duplex non-stitched fragments, and simplex non-stitched fragments.

[0386] The variant calling application can output the calls in various formats, such as .VCF or .GVCF files. By way of example only, the variant calling application can be included in the MiSeqReporter pipeline (e.g., when The analysis may include a single protocol or a combination of protocols for analyzing sample interpretation in a specified manner to obtain the desired information.

[0387] Then, one or more processors perform the verification operation related to potential variant identification. The verification operation can be based on the hierarchy of quality score and / or layered test, as described below. When the verification operation verifies or verifies potential variant identification, the verification operation passes the variant identification information (from the variant identification application) to the sample report generator. Alternatively, when the verification operation invalidates or disqualifies potential variant identification, the verification operation passes corresponding instructions (for example, negative indicator, no identification indicator, invalid identification indicator) to the sample report generator. The verification operation can also pass the confidence score related to the correct or invalid identification of variant identification.

[0388] Next, one or more processors generate and store sample reports.Sample reports may include, for example, information about multiple genetic loci about a sample.For example, for each genetic locus of a predetermined set of genetic loci, a sample report may provide at least one of the following: provide genotype identification; indicate that genotype identification cannot be performed; provide a confidence score for the certainty of genotype identification; or indicate possible problems with the assay method for one or more genetic loci.Sample reports may also indicate the individual's sex and / or indication sample providing a sample including multiple sources.As used herein, "sample report" may include digital data (e.g., data files) and / or a printed report of a set of genetic loci or a predetermined set of genetic loci or genetic loci or a set of genetic loci.Therefore, generating or providing may include creating a data file and / or printing a sample report, or displaying a sample report.

[0389] The sample report may indicate that a variant identification has been determined, but not verified. When a variant identification is determined to be invalid, the sample report may indicate other information about the basis for determining not to verify the variant identification. For example, the additional information in the report may include a description of the original fragments and the extent to which the original fragments support or conflict with the variant identification (e.g., count). Additionally or alternatively, the additional information in the report may include a quality score obtained according to the implementation described herein.

[0390] Variant Identification Application

[0391] Implementation disclosed herein includes analyzing sequencing data to identify potential variation identification. Variant identification can be performed on stored data to perform previously performed sorting operations. Additionally or alternatively, it can be performed in real time while performing sequencing operations. Each sample interpretation is assigned to a corresponding genetic locus. The sequence of nucleotides that can be interpreted based on the sample, or in other words, the order of nucleotides (e.g., A, C, G, T) in the sample interpretation, the sample interpretation is assigned to the corresponding genetic locus. Based on this analysis, the sample interpretation can be designated as a possible variant / allele including a specific genetic locus. The sample interpretation can be collected (or aggregated or classified) together with other sample interpretations that have been designated as possible variants / alleles including genetic loci. The allocation operation can also be referred to as a recognition operation, in which the sample interpretation is identified as being associated with a specific genetic position / locus. The sample interpretation can be analyzed to locate one or more identification sequences (e.g., primer sequences) that distinguish the nucleotides that the sample interpretation is interpreted from other samples. More specifically, one or more recognition sequences can identify the sample interpretation from other sample interpretations as being associated with a specific genetic locus.

[0392] The distribution operation can include analyzing the n nucleotide sequences of the recognition sequence to determine whether the n nucleotide sequences of the recognition sequence are effectively matched with one or more selection sequences. In a specific implementation, the distribution operation can include analyzing the first n nucleotides of the sample sequence to determine whether the first n nucleotides of the sample sequence are effectively matched with one or more selection sequences. The number n can have various values, which can be programmed into the protocol or input by the user. For example, quantity n can be defined as the number of nucleotides of the shortest selection sequence in the database. Quantity n can be a predetermined number. The predetermined number can be, for example, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29 or 30 nucleotides. However, less or more nucleotides can be used in other implementations. Number n can also be selected by an individual (such as a user of the system). Quantity n can be based on one or more conditions. For example, quantity n can be defined as the number of nucleotides or a specified number of the shortest primer sequence in the database, whichever is smaller. In some implementations, a minimum value for n can be used, such as 15, so that any primer sequence that is less than 15 nucleotides can be designated as an exception.

[0393] In some cases, n nucleotide sequences of recognition sequence may not completely match the nucleotide of selection sequence. However, if recognition sequence is almost identical with selection sequence, recognition sequence can effectively match with selection sequence. For example, if the sequence of n nucleotide sequences of recognition sequence (for example, first n nucleotide) has the mispairing (for example 3) and / or the selection sequence matching of the displacement quantity (for example 2) specified that is no more than a specified number and / or the selection sequence matching of specified number of shifts of selection sequence, then can be genetic locus identification sample interpretation. Can set up rule, make each mispairing or displacement can be counted as the difference between sample interpretation and primer sequence. If the number of difference is less than specified number, then can be corresponding genetic locus (that is, be assigned to corresponding genetic locus) identification sample interpretation. In some implementations, can determine matching score based on the recognition sequence of sample interpretation and the difference number between the selection sequence related to locus. If matching score exceeds the matching threshold specified, then can be designated as the potential locus of sample interpretation by the genetic locus corresponding to the selection sequence. In some implementations, can carry out subsequent analysis to determine whether sample interpretation is called genetic locus.

[0394] If the sample is interpreted and one of the selected sequences in the effective matching database is obtained (i.e., fully matched or almost matched as described above), the sample is interpreted and assigned to the genetic locus relevant to the selected sequence. This can be referred to as locus identification or temporary locus identification, wherein the sample is interpreted for the genetic locus identification sample relevant to the selected sequence. However, as described above, more than one locus may be required to interpret the sample. In such an implementation, further analysis can be performed to identify or assign a sample interpretation only for one of the potential genetic loci. In some implementations, the sample interpretation compared with the reference sequence database is the first interpretation from paired end sequencing. When paired end sequencing is performed, the second interpretation (representing the original fragment) relevant to the sample interpretation will be obtained. After the allocation, the subsequent analysis performed using the interpretation assigned can be based on the type of the genetic locus identified for the interpretation assigned.

[0395] Next, the sample interpretation is analyzed to identify potential variant identifications. Wherein, the analysis results determine the potential variant identifications, sample variant frequencies, reference sequences, and locations where variations occur in the genomic sequence of interest. For example, if the known genetic locus contains SNPs, the specified interpretation of the genetic locus identification can be analyzed to identify the SNPs of the specified interpretation. If the known genetic locus contains polymorphic repetitive DNA elements, the specified interpretation can be analyzed to identify or characterize the polymorphic repetitive DNA elements in the sample interpretation. In some implementations, if the assigned interpretation effectively matches the STR locus and the SNP locus, a warning or flag can be assigned to the sample interpretation. The sample interpretation can be designated as the STR locus and the SNP locus. The analysis can include aligning the assigned interpretation according to the alignment protocol to determine the sequence and / or length of the assigned interpretation. The alignment scheme can include the method described in International Patent Application No. PCT / US2013 / 030867 (Publication No. WO 2014 / 142831) filed on March 15, 2013, which is incorporated herein by reference in its entirety.

[0396] Then, one or more processors analyze the original fragments to determine whether there is a supporting variant in the corresponding position in the original fragment. Various types of original fragments can be identified. For example, the variant identifier can identify the type of the original fragment, and the original fragment shows a variant that verifies the original variant identification. For example, the type of the original fragment can represent a duplex stitching fragment, a simplex stitching fragment, a duplex unstitched fragment or a simplex unstitched fragment. Replace the above examples or in addition to the above examples, other original fragments can also be optionally identified. Relevant to identifying each type of original fragment, the variant identifier also identifies the position where the original fragment supports the variant, and the number of original fragments that show the supporting variant. For example, the variant identifier can output an indication, which identifies the interpretation of 10 original fragments to represent the duplex stitching fragment with a supporting variant at a specific position X. The variant identifier can also output an indication, and the interpretation of its five original fragments is identified as representing the simplex unstitched fragment with a supporting variant at a specific position Y. The variant identifier can also output many original fragments corresponding to the reference sequence, so the supporting variant is not included, otherwise it provides evidence to verify the potential variant call at the genomic sequence of interest.

[0397] Next, counts of original fragments including supported variants are maintained, as well as locations where supported variants occur. Additionally or alternatively, counts of original fragments that do not include supported variants at positions of interest (relative to the position of potential variant identification in the sample interpretation or sample fragments) can be maintained. Additionally or alternatively, counts of original fragments that correspond to reference sequences and do not authenticate or confirm potential variant identification can be maintained. The determined information will be output to the variant identification verification application, including counts and types of original fragments that support potential variant identification, locations in the original fragments that support variants, counts of original fragments that do not support potential variant identification, etc.

[0398] When potential variant identification is identified, the process outputs the indication of potential variant identification, variant sequence, variant position and the reference sequence associated therewith. Variant identification is designated as representing "potential" variants, because errors may cause the identification process to identify the wrong variable. According to the implementation of this article, potential variant identification is analyzed to reduce and eliminate erroneous variants or false positives. Additionally or alternatively, the process analyzes one or more original fragments associated with sample interpretation, and outputs the corresponding variant identification associated with the original fragment.

[0399] Benign training set generation

[0400] Millions of human genomes and exomes have been sequenced, but their clinical utility remains limited by the difficulty of distinguishing disease-causing mutations from benign genetic variation. Here, we demonstrate that common missense variants in other primate species are largely clinically benign in humans, allowing the systematic identification of disease-causing variants by elimination. Using tens of thousands of common variants sequenced from populations of six nonhuman primate species, we trained a deep neural network to identify disease-causing variants in patients with rare diseases with 88% accuracy and discovered 14 new candidate genes for intellectual disability in a gene-wide sense. Cataloging common variants from other primate species could improve the interpretation of millions of variants of uncertain significance, further enhancing the clinical impact of human genome sequencing.

[0401] The clinical feasibility of diagnostic sequencing is limited by the difficulty in interpreting rare genetic variants in human populations and inferring their impact on disease risk. Due to the adverse health effects of rare genetic variants, clinically significant genetic variants are extremely rare in human populations, and for the vast majority of rare genetic variants, their impact on human health is uncertain. These large numbers of rare variants of uncertain clinical significance severely hamper the use of sequencing for personalized medicine and population-wide health screening.

[0402] Most penetrant Mendelian diseases occur at low prevalence in a population, so observation of a variant at high frequency in a population is strong evidence for a favorable outcome. Assessing common variation across diverse human populations is an effective strategy for cataloguing benign variants, but the total number of common variants in present-day humans is limited due to a bottleneck in our species' recent history during which much of the original diversity was lost. The effective population size (Ne) of present-day human populations has increased dramatically from less than 10,000 individuals over the past 15,000-65,000 years, and a small set of common polymorphisms has been traced to the limited capacity for variation in populations of this size. Of the more than 700 million potential protein-altering missense substitutions in the reference genome, only approximately one in a thousand is greater than 0.1% of the allele frequency in the entire population.

[0403] Outside of the human population, chimpanzees comprise the next closest living species, with 99.4% amino acid sequence identity. The close identity of protein-coding sequences between humans and chimpanzees suggests that purifying selection on chimpanzee protein-coding variants could also mimic the fitness consequences of state-identical human mutations.

[0404] Since the average age of neutral polymorphisms in the human ancestral lineage (-4N generations) is a fraction of the time since the species diverged (6 million years ago), the mutational space explored by naturally occurring chimpanzee variants is largely non-overlapping, except for rare instances of haploids maintained by non-balancing selection, excluding chance. If state-consistent polymorphisms similarly affect fitness in both species, the presence of high-allele frequency variants in chimpanzee populations may indicate a benign outcome in humans, expanding the catalog of known variants with benign outcomes established by purifying selection.

[0405] Results - Common variants in other primates are largely benign in humans

[0406] Recently available aggregated exome data from the Exome Aggregation Consortium (ExAC) and the Genome Aggregation Database (gnomAD) for 123,136 humans have allowed us to measure the effects of natural selection on missense and synonymous mutations across the spectrum of allele frequencies. Fig.49A , 51 , 52A, 52B, 52C, and 52D), the rare single variants observed once in the cohort closely matched the expected 2.2 / 1 missense / synonymous mutation rate predicted by de novo mutations, but the number of missense variants observed at higher allele frequencies was reduced due to purification of deleterious mutations by natural selection. Although these variants are also observed in healthy individuals, the gradual decrease in missense / synonymous mutation rates with increasing allele frequencies is consistent with a considerable proportion of missense variants with a population frequency of less than 0.1%, which have mildly toxic consequences. In addition to a few well-documented anomalies caused by balancing selection and founder effects, these findings support the common empirical practice of diagnostic laboratories by filtering out variants with allele frequencies greater than 0.1% to 1% may be benign for explicit genetic diseases.

[0407] We identified common chimpanzee variants sampled twice or more from cohorts of 24 unrelated individuals. We estimated that 99.8% of the variants were common in the common chimpanzee population (allele frequency (AF) > 0.1%), suggesting that they have survived purifying selection. For state-consistent human variants ( Fig.49B), we tested human allele frequencies, excluding the extended major histocompatibility complex region as a known region of balancing selection and variants lacking one-to-one mapping in multiple sequence alignments. For human variants consistent with common chimpanzee variant status, missense / synonymous mutation rates were largely constant across the human allele frequency spectrum (P > 0.5 as tested by chi-square (χ2)), consistent with attrition of common chimpanzee variants in the human population under negative selection and consistent selection coefficients for missense variants in both species. The low missense / synonymous mutation rates observed in human variants consistent with common chimpanzee variant status are consistent with the large effective population size in chimpanzees (Ne-73000), resulting in more effective filtering of mildly deleterious variants.

[0408] In contrast, for individual chimpanzee variants (sampled only once from the cohort), we observed a significant decrease in missense / synonymous mutation rates at common allele frequencies (P less than 5.8 × 10-6), as shown in Figure 2. Fig.49C This suggests that 24% of the individual chimpanzee variants have been filtered through purifying selection in the human population at allele frequencies greater than 0.1%. This loss implies that a significant fraction of chimpanzee single variants are rare deleterious mutations whose fitness-damaging effects prevent them from reaching common allele frequencies in both species. We estimate that only 69% of the individual variants are common (AF greater than 0.1%) in the common chimpanzee population.

[0409] We next identified human variants with consistent variant status observed in at least one of the six nonhuman primate species. Variants in each of the six species were identified through the Great Ape Genome Project (chimpanzee, bonobo, gorilla, and orangutan) or submitted to the database of single nucleotide polymorphisms (bdSNP) from the Primate Genome Project (rhesus macaque, marmoset) and largely represent common variants based on the limited number of sequenced singletons and the low missense:synonymous mutation rates observed in each species (see Supplementary Table 1). Similar to chimpanzees, we found that the missense / synonymous mutation rates of variants from the six nonhuman primate species were roughly the same across the human allele frequency spectrum, which is different from the mild attrition of missense variants at the common allele frequency (e.g., Fig.49D , 53 , 54 and 55, Supplementary data file 1), due to the inclusion of a small number of rare variants (16% of variants in chimpanzees are below an allele frequency of 0.1%, and lower in other species due to the small number of sequenced individuals). These results suggest that selection coefficients for state-consistent missense variants are consistent with primate lineages at least related to New World monkeys, which are estimated to have diverged from the human ancestral lineage 35 million years ago.

[0410] We found that human missense variants consistent with observed primate variant status were significantly enriched in the ClinVar database, which had a benign outcome. After excluding variants of uncertain significance and variants with conflicting annotations, ClinVar variants occurring in at least one nonhuman primate species were annotated as benign or likely benign an average of 90% of the time, compared with ClinVar missense variants that were usually benign 35% of the time (P < 10-40, e.g., 90% confidence interval [CI], 1.1-2.0) and 1.2-1.0 (P < 0.05). Fig.49E ). Except for human variants with allele frequencies greater than 1% to reduce selection bias, the pathogenicity of ClinVar annotations for primate variants was slightly greater than that observed in healthy humans sampled from a similar sized cohort (95% benign or probably benign, p=0.07).

[0411] The field of human genetics has long relied on model organisms to infer the clinical impact of human mutations, and the long evolutionary distances to most genetically tractable animal models raise concerns about the extent to which findings in model organisms can be generalized to humans. We extended our analysis beyond primate lineages to include highly common variants from four other mammalian species (mouse, pig, goat, and cattle) and two more distant vertebrates (chicken and zebrafish). We selected species for which gene-wide variation in dbSNP was well defined, and we determined these variants to be highly common, as the missense / synonymous mutation rate was much less than 2.2 / 1. Compared to primate analyses, there were significant differences in common allele frequencies (e.g., 95% confidence intervals) between the two groups. Fig.50A ), human missense mutations consistent with variant states are significantly lost in more distant species, and the loss increases over longer evolutionary distances (e.g. Fig.50B and Supplementary Tables 2 and 3).

[0412] The presence of missense variants that are toxic to humans and tolerate high allele frequencies in more distant species suggests that the selection coefficients for missense variants that are consistent with the state are well differentiated between humans and more distant species. Nevertheless, since the fraction of missense variants attrition by natural selection at common allele frequencies is usually less than the 50% attrition observed for human missense variants (e.g. Fig.49A ), the presence of missense variants in more distant mammals still increases the likelihood of a benign outcome. Consistent with these results, we found that 73% of ClinVar missense variants observed in mice, pigs, goats, and cattle were annotated as benign or likely benign, while this was 90% (P less than 2×10-8) for primate variants (as shown in Figure 50C) and 35% for the entire ClinVar database.

[0413] To establish that evolutionary distance, rather than domestication relics, is the primary driver of divergent selection coefficients across a wide range of evolutionary distances (as shown in Fig. 50D, Supplementary Table 4, and Supplementary Data File 2), we reran our analyses using fixed substitutions between pairs of close relatives that replace intraspecific polymorphisms. We found that depletion of human missense variants consistent with the state of intraspecific fixed substitutions increased with evolutionary branch length, with wild species not being significantly different from those of domesticated species. Consistent with earlier work in flies and yeast, the number of state-consistent fixed missense substitutions was found to be lower than expected by chance compared with diverging lineages.

[0414] Deep learning network for variant pathogenicity classification

[0415] The disclosed technology provides a deep learning network for variant pathogenicity classification. Based on the clinical application significance of variant classification, a large number of attempts have been made to use supervised machine learning to solve this problem, however, these efforts have been hampered by the lack of a real data set of suitable size that includes variants clearly marked as benign and pathogenic for training.

[0416] Existing databases of human-curated variants do not represent the entire genome, with 50% of the variants in the ClinVar database coming from just 200 genes (1% of human protein-coding genes). Furthermore, system learning determines that the supporting evidence for most human annotations is questionable, underestimating the difficulty of interpreting rare variants observed in only one patient. Despite the increasing rigor of human interpretation, a large number of classification guidelines have been developed around consensus practices, which risk reinforcing existing trends. To reduce human interpretation bias, recent classifiers are trained on common human polymorphisms or fixed human-chimpanzee surrogates, which also serve as their input: the prediction scores of earlier classifiers trained on human-curated databases. In the absence of independent, unbiased real-world datasets, objective benchmarks of the performance of these different approaches are elusive.

[0417] Variants from six non-human primates (chimpanzees, bonobos, gorillas, orangutans, rhesus monkeys, and marmosets) provide more than 300,000 unique missense variants that do not overlap with common human variants, largely representing common variants with benign outcomes screened by purifying selection, greatly expanding the training data set available for machine learning methods. In general, after excluding variants with uncertain meanings and variants with conflicting annotations, each primate species provides more variants than the entire ClinVar database (42,000, as of November 2017). In addition, this content is free of human interpretation bias.

[0418] By using a panel that included common human variants (AF greater than 0.1%) and primate variants (Supplementary Table 5 (e.g. Fig.58 We trained a new deep residual network, PrimateAI, for inputting other species (such as Figure 2 and 3 The amino acid sequences and orthologous sequences flanking the variant of interest are aligned in Figure 1. Unlike existing classifiers that use human-engineered features, our deep learning network learns accurate features directly from the original sequence. To incorporate information about protein structure, we trained separate networks for predicting secondary structure and solvent accessibility of the sequence alone, and then included these as subnetworks in the full model (see Figure 1). Figure 5 and 6 ). Given the small number of human proteins that have been successfully crystallized, inferring the structure from the original sequence has the advantage of avoiding biases due to incomplete protein structure and functional domain annotation. Including the protein structure, the total depth of the network is 36 convolutional layers, including approximately 400,000 trainable parameters.

[0419] In order to train the classifier only using variants with benign signatures, we formulate the prediction problem as to whether a given observed mutation is common in the population. Several factors influence the probability of observing a variant at high allele frequency, where we are only interested in its toxicity. Other factors include mutation rate, technical error such as sequencing range, and factors that influence neutral genetic drift such as gene conversion.

[0420] We matched each variant in the benign training set to a missing missense variant from 123,136 exomes from the ExAC database, controlled for various confounding factors, and trained a deep learning network to distinguish benign variants from matched controls (e.g. Fig.24 As the number of unlabeled variants far exceeds the size of the labeled benign training dataset, we trained eight networks simultaneously, each using a different set of unlabeled variants matched to the benign training dataset, to obtain a consensus prediction.

[0421] Using only the primary amino acid sequence as its input, the deep learning network accurately assigns high pathogenicity scores to useful protein functional domain residues, e.g. Fig. 20Shown is SCN2A for the voltage-gated sodium channel: an important disease gene for epilepsy, autism, and intellectual disability. The structure of SCN2A includes four homologous repeats, each of which includes six transmembrane helices (S1-6). In response to membrane depolarization, the positively charged S4 transmembrane helix moves toward the extracellular side of the membrane, allowing the S5 / S6 pore-forming domain to open via the S4-S5 connector. Mutations in the S4 and S4-5 connectors and S5 domains that are clinically associated with early-onset epileptic encephalopathy were predicted by the network to have the highest pathogenicity scores in the gene, and these mutations were also removed from healthy populations. We also found that the network identifies important amino acid positions within the domain and assigns the highest pathogenicity scores to mutations at these positions, such as DNA-containing residues of transcription factors and catalytic residues of enzymes (such as Fig.25A , 25B , 25C and 26).

[0422] To better understand how deep learning networks can gain insight into protein structure and function from primary sequences, we visualized the trainable parameters from the first three layers of the network. In these three layers, we observed that the network learned correlations between different amino acid weights, which were close to existing measures of amino acid distance such as the Grantham score (e.g. Fig. 27 The outputs of these initial layers become the inputs of subsequent layers, allowing deep networks to build progressively higher-order data displays.

[0423] Using the 10,000 common primate variants held out for training, we compared network performance with existing classification algorithms. Since 50% of all recently occurring missense singletons are filtered out by purifying selection at common allele frequencies (e.g. Fig.49A We determined the 50% score for each classifier on a set of 10,000 randomly selected variants that were matched to 10,000 common primate variants by mutation rate and sequencing coverage, and evaluated the accuracy of each classifier at this threshold (as shown in Fig.21D , 28A and Supplementary Data File 4). Our learned network outperforms other classifiers (91% accuracy) in assigning benign outcomes to 10,000 conserved common primate variants (the best model below has 80% accuracy).

[0424] Only with human variant data (e.g. Fig.21DCompared to the accuracy of networks trained with 100% serotonin (Figure ​(Figure4),), almost half of the improvement over existing methods comes from the use of deep learning networks and half from the augmentation of the primate variant training dataset. To test the classification of variants of uncertain significance in a clinical setting, we evaluated the ability of deep learning networks to distinguish de novo mutations in patients with neurodevelopmental disorders compared with healthy controls. Neurodevelopmental disorders are one of the largest categories of rare genetic diseases in terms of prevalence, and recent triple sequencing studies have implicated a central role for de novo missense and protein-truncating mutations.

[0425] We classified each identified de novo missense variant in 4293 affected individuals from the Deciphering Developmental Disorders (DDD) cohort against that in 2517 unaffected clan members from the Sample Database for Children (SSC) cohort and assessed the difference in prediction scores between the two assignments using a Wilcoxon rank-sum test ( Fig.21E , 29A In this task, the deep trained network significantly outperformed other classifiers (P less than 10-28, as shown in Figure 29B). Fig.21F and 28B The performance of the different classifiers on the primate variant-preserving dataset was correlated with the DDD case vs. control case dataset (correlation coefficient ρ = 0.57, P < 0.01), indicating that despite using completely different sources and methods (e.g. Fig. 30A The two datasets for evaluating pathogenicity were in good agreement.

[0426] Next, we aimed to assess the accuracy of the deep learning network in classifying benign and pathogenic variants in the same gene. Given that the DDD cohort consists primarily of index cases with affected children without affected first-degree relatives, it is important that the accuracy of the classifier is not inflated, favoring pathogenicity in genes with a de novo dominant inheritance pattern. We restricted the analysis to the 605 genes with nominally significant disease associations in the DDD study, calculating only protein-truncating variants (P < 0.05). In these genes, de novo missense variants were enriched 3 / 1 compared to expectations, (e.g., Fig.22A ), indicating that 67% of the mutations were pathogenic.

[0427] In the same gene set (such as Fig. 22B As shown in Figure 2, (P less than 10-15), the deep learning network can distinguish between pathogenic and benign de novo variants, which is much stronger than other methods (such as Fig. 22C and 28C At a binary critical value greater than 0.803 (as shown in Fig.22D and 30B), 65% of the de novo missense mutations in the cases were classified as pathogenic by the deep learning network, compared with 14% in the controls, which is consistent with a classification accuracy of 88% (as shown in Fig.22E and 30C Given the incomplete penetrance frequencies and heterogeneity of expressivity of neurodevelopmental disorders, this number may underestimate the accuracy of our classifier due to the inclusion of partially penetrant pathogenic variants in the control group.

[0428] Novel candidate gene discovery

[0429] Applying a threshold of ≥0.803 to stratify pathogenic missense variants increased enrichment of de novo missense variants in DDD patients by 1.5-2.2-fold, which is similar to protein-truncating variants (2.5-fold), discarding variants that were enriched beyond expectation by less than one-third of the total. This substantially increased statistical power, allowing the discovery of 14 additional candidate genes in intellectual disability. Previously, statistical power did not reach the gene-wide significance threshold in the original DDD study (Table 1).

[0430] Comparison with human expert selections

[0431] We tested the performance of various classifiers on the most recent human-curated variants from the ClinVar database and found no significant correlation between classifier performance on the ClinVar database and on the retained primate variant dataset or the DDD case-versus-control dataset (P = 0.12 and 0.34, respectively). Fig.31A and 31B We assume that existing classifiers are biased toward human expert selection, and that these human heuristics are often directionally correct but not optimal. In one example, the Grantham mean difference between pathogenic and benign variants in ClinVar is twice the difference between DDD cases and controls in 605 disease-associated genes (e.g. Figure 2 In contrast, human curation appears to underutilize protein structure, and in particular, the importance of surface-exposed residues where they can interact with other molecules. We observed that both ClinVar pathogenic and DDD de novo variants were associated with predicted solvent-exposed residues, but the difference in solvent accessibility between benign and pathogenic ClinVar variants was only half that of DDD cases versus controls. These findings suggest that confirmation bias is beneficial for human experts to more directly interpret factors such as Grantham and conservation. Machine learning classifiers trained on human curated databases are expected to reinforce these trends.

[0432] The results suggest that systematic primate population sequencing is an effective strategy for classifying the millions of human variants of uncertain significance that currently limit clinical genomic interpretation. The accuracy of the deep learning network for both retained common primate variants and clinical variants increased as the number of benign variants used to train the network increased (e.g. Fig.23A Furthermore, training the variants from each of the six non-human primate species individually helped improve the network's performance, whereas training the variants from more distant mammals had a negative impact on network performance (e.g. Fig. 23B and 23C These results support the assertion that common primate variants for penetrant Mendelian diseases are mostly benign in humans, whereas this is not the case for variants in more distant species.

[0433] Although the number of nonhuman primate genomes examined in this study is small relative to the number of sequenced human genomes and exomes, it is noteworthy that these additional primates provide a disproportionate amount of information on common benign variants. ExAC simulations suggest that only a few hundred individuals (e.g. Fig.56 After the rapid stagnation of common human variants (less than 0.1% allele frequency) in the primate population, sequencing millions of healthy populations has provided primarily additional rare variants. Unlike common variants, which are known to be largely clinically benign based on allele frequency, rare variants in healthy populations can cause recessive or dominant diseases with incomplete penetrance. Because each primate species carries a different set of common variants, sequencing dozens of primates from each species is an effective strategy for systematically cataloging benign missense variants in primate lineages. Indeed, the 134 individuals from six nonhuman primate species examined in this study provided nearly four times as many common missense variants as the 123,136 individuals in the ExAC study (Supplementary Table 5 (e.g., Fig.58 This is also true for sequencing studies of primate species involving hundreds of individuals, even when the number of unrelated individuals in wildlife sanctuaries and zoos is relatively small. Minimizing disturbance to wild populations is important from the standpoint of conservation and human treatment of nonhuman primates.

[0434] The genetic diversity of the human population is much lower than that of nonhuman primate species, about half the number of single nucleotide variants carried by chimpanzees, gorillas, and gibbons, and one-third the number of variants carried by orangutans. Although the level of genetic diversity of most nonhuman primate species is unknown, based on the large number of nonhuman primate species available, we infer that most likely benign human missense positions are covered by common variants in at least one primate species, allowing pathogenic variants to be identified by systematic elimination (e.g., Fig.23D). Even if only a small set of these species were sequenced, increasing the size of the training data allowed for more accurate prediction of missense outcomes via machine learning. Finally, while our findings focus on missense variants, this strategy can also be applied to inferring outcomes of noncoding variants, especially in conserved regulatory regions, where the human and primate genomes are sufficiently aligned to unambiguously determine if the variants are identical.

[0435] Of the 504 known non-human primate species, approximately 60% are threatened with extinction due to illegal poaching and widespread habitat loss. Population decline and potential species extinction will cause irreparable losses to genetic diversity, and it is urgent to encourage species conservation efforts worldwide, which will benefit rare and irreplaceable species and humans themselves.

[0436] Data generation and comparison

[0437] Coordinates for this disclosure relate to human genome build UCSC hg19 / GRCh37, including coordinates of other species variants mapped to hg19 using multiple sequence alignments. Canonical transcripts of protein-coding DNA sequences and multiple sequence alignments of 99 vertebrate genomes and branch lengths were downloaded from the UCSC genome browser.

[0438] We obtained human exome polymorphism data from the Exome Assembly Consortium (ExAC) / Genome Assembly Database (gnomAD Exome) v2.0. We obtained primate variant data from the Great Ape Genome Sequencing Project, which includes full genome sequencing data and genotypes for 24 chimpanzees, 13 bonobos, 27 gorillas, and 10 orangutans. We also included variants from 35 chimpanzees from separate studies of chimpanzees and bonobos, but due to different variant calling methods, we excluded variants from population analyses and used only these variants to train the deep learning model. In addition, 16 rhesus macaque individuals and 9 marmoset individuals were used to test the original genome project variants of these species, and individual-level information was not available. We obtained variant data for rhesus macaques, marmosets, pigs, cattle, goats, mice, chickens, and zebrafish from dbSNP. Since genotype information for individuals used for population analysis was not available, dbSNP also included additional orangutan variants that were used only for training the deep learning model. To avoid balancing selection effects, we also excluded variants in the extended major histocompatibility complex region (chr6:28,477,797-33,448,354) used for population analysis.

[0439] We used a multispecies comparison of 99 vertebrates to confirm one-to-one mapping of orthologs to human protein-coding regions and to prevent mapping of pseudogenes. Variant status was consistent if it occurred in the reference / variant orientation. To ensure that variants had the same predicted protein-coding consequence in both humans and other species, the other two nucleotides in the codon were required to be consistent between species for missense variants and synonymous mutation variants. The polymorphisms for each species included in the analysis are listed in the Supplementary Data file, and detailed metrics are shown in Supplementary Table 1.

[0440] For each of the four allele frequency categories (e.g. Fig.49A We used variation in intragenic regions to estimate the expected number of synonymous and missense variants in each of the 96 possible trinucleotide contexts and corrected for mutation rates (as shown in Fig.51 and Supplementary Tables 7 and 8 (e.g. Fig. 9 We also analyzed the consensus CpG and non-CpG dinucleotide variants separately and confirmed that the missense / synonymous mutation rate was flat across the two allele spectrums, indicating that our analysis is applicable to CpG and non-CpG variants despite their very different mutation rates (e.g. Fig.52A , 52B , 52C and 52D).

[0441] Attrition of state-consistent human missense variants that are polymorphic in other species

[0442] To evaluate whether variants in other species are tolerated at the common allele frequency in humans (greater than 0.1%), we identified state-consistent human variants with variability in other species. For each variant, we assigned it to one of four categories based on its allele frequency in the human population (singleton, greater than singleton-0.01% singleton, 0.01-0.1% and greater than 0.1% singleton) and assessed the reduction in missense / synonymous mutation rate (MSR) between rare (less than 0.1%) and common (greater than 0.1%) variants. At the common allele frequency (greater than 0.1%), the loss of state-consistent variants suggests that this portion of variants from other species with considerable toxicity is being screened by natural selection at the common allele frequency in humans.

[0443]

[0444] Calculate missense / synonymous mutation rates and attrition rates for each species, as Fig.50B and Supplementary Table 2. In addition, lower than the common chimpanzee variants (such as Fig.49B ), a single chimpanzee variant (e.g. Fig.49C ) and mammalian variants (e.g. Fig.50A), we performed a chi-square test (x2) on the 2*2 contingency table to test whether the difference in missense / synonymous mutation rates between rare variants and common variants was significant.

[0445] Since only a limited number of individuals from the Great Ape Genome Project were sequenced, we used the human allele frequency spectrum from ExAC to evaluate the subset of sampled variants that were either rare (less than 0.1%) or common (greater than 0.1%) in the general chimpanzee population. We sampled a set of 24 individuals based on ExAC allele frequencies and identified missense variants observed once or multiple times in this set. 99.8% of variants observed multiple times were likely to be common (greater than 0.1%) variants in the general population, while 69% of variants observed once were likely to be common variants in the general population. To verify that the observed depletion of missense variants in more distant mammals was not caused by the confounding effects of better conserved genes, we repeated the above analysis, limiting only genes with greater than 50% average nucleotide identity in multiple sequence alignments of 11 primates and 50 mammals compared to humans.

[0446] Removing 7% of human protein-coding genes from the analysis did not substantially affect the results. Second, to confirm that the results were not affected by variant calling events or domestication relics (because most species selected from dbSNP are domesticated species), we repeated the analysis using fixed surrogates from closely related species pairs instead of intraspecific polymorphisms (as shown in Figure 50D, Supplementary Table 4, and Supplementary Data File 2).

[0447] ClinVar analysis of polymorphisms in humans, primates, mammals, and other vertebrates

[0448] To test the clinical impact of variants consistent with status in other species, we downloaded the ClinVar database, excluding variants with conflicting pathogenic annotations or variants identified only as of uncertain significance. Below, Supplementary Table 9 shows the filtering steps, totaling 24,853 missense variants in the pathogenic category and 17,775 missense variants in the benign category.

[0449] We counted pathogenic and benign ClinVar variants consistent with variant status in humans, nonhuman primates, mammals, and other vertebrates. For humans, we simulated a set of 30 individuals sampled from ExAC allele frequencies. Supplementary Table 10 shows the number of benign and pathogenic variants for each species.

[0450] Generation of benign and unlabeled variants for model training

[0451] We constructed a large benign training dataset of common benign missense variants from humans and non-human primates for machine learning. The dataset includes common human variants (greater than 0.1% allele frequency, 83,546 variants) and variants from chimpanzees, bonobos, gorillas, orangutans, rhesus monkeys, and marmosets (301,690 unique primate variants). Supplementary Table 5 shows the number of benign training variants provided by each source.

[0452] We trained a deep learning network to distinguish between a set of labeled benign variants and a set of unlabeled variants matched with trinucleotide context controls, sequencing coverage, and alignability between species and humans. To obtain an unlabeled training dataset, we started with all possible missense variants in normalized coding regions. We excluded variants observed in 123,136 exomes from ExAC and variants in start and stop codons. A total of 68,258,623 unlabeled missense variants were generated. When selecting matched unlabeled variants for primate variants, variants were filtered to correct for regions with poor sequencing coverage and regions where the human and primate genomes did not align one-to-one.

[0453] By training eight models using the same set of labeled benign variants and eight groups of randomly sampled unlabeled variants and obtaining the average of their predictions, we obtained consensus predictions. We also selected a series of 10,000 primate variants randomly sampled for verification and testing, which were retained by training (as shown in Supplementary Table Data Document 3). For each data set, we sampled 10,000 non-labeled variants of trinucleotide context matches, and when comparing different classification algorithms (as shown in Supplementary Data Document 4), we used these variants to normalize the threshold of each classifier. In other embodiments, fewer models or additional models can be used in the set, ranging from 2 to 500.

[0454] We evaluated the classification accuracy of two versions of the deep learning network, one trained only on common human variants and the other trained on a complete benign labeled dataset that includes both common human and primate variants.

[0455] Deep Learning Network Architecture

[0456] For each variant, the pathogenicity prediction network takes as input a 51-length amino acid sequence centered on the variant of interest and substitutes the secondary structure and solvent accessibility network with the missense variant at the central position (see Figure 2). Figure 2 and 3Three 51-length position frequency matrices were generated from the multiple sequence alignment of 99 vertebrates, one for 11 primates, one for 50 mammals excluding primates, and one for 38 vertebrates excluding primates and mammals.

[0457] The secondary structure deep learning network predicts the three-state secondary structure for each amino acid position: alpha helix (H), beta sheet (B), and loop (C) (as shown in Supplementary Table 11). The solvent accessibility network describes the three-state solvent accessibility for each amino acid position: buried (B), intermediate (I), and exposed (E) (as shown in Supplementary Table 12). The networks all take the flanking amino acid sequence as their input and are trained using known non-redundant crystal structure labels from the protein database (as shown in Supplementary Table 13). For the input of the pre-trained three-state secondary structure and three-state solvent accessibility networks, we used a single-length position frequency matrix of length 51 and depth 20 generated from the multiple sequence alignment of all 99 vertebrates. After pre-training the network on known crystal structures in the protein database, the last two layers of the secondary structure and solvent models were removed and the network output was directly connected to the pathogenicity model input. The best test accuracy of the three-state secondary structure prediction model was 79.86% (Supplementary Table 14). Comparison of neural network predictions using DSSP-annotated structural landmarks for nearly 4000 human proteins with crystal structures (protein-defined secondary structure) showed no substantial differences compared to using predicted structural landmarks alone (shown in Supplementary Table 15).

[0458] The deep learning network for pathogenicity prediction (primateAI) and the deep learning network for predicting secondary outcomes and solvent accessibility both use residual block architectures. Figure 3 and Supplementary Table 16 ( Figure 4A , 4B and 4C) describe the detailed architecture of primateAI. Figure 6 and Supplementary Tables 11 and 12 ( Fig. 7A , 7B )( Fig. 8A , 8B ) describes a detailed architecture for predicting secondary structure and solvent accessibility networks.

[0459] Classifier performance on a hold-out test set of 10,000 primate variants.

[0460] We benchmark the deep learning network against 20 other previously published classifiers using the 10,000 conserved primate variants in the test dataset, for which we obtain prediction scores from the database dbNSFP. Fig.28AThe performance of each classifier on the 10,000 held-out primate variant test set is also shown. Since different classifiers have very different score distributions, we used 10,000 randomly selected unlabeled variants matched to the test set by trinucleotide context to identify a 50% threshold for each classifier. We benchmarked each classifier on a subset of variants in the 10,000 held-out primate variant test set that were classified as benign at the 50% threshold for that classifier to ensure fair comparison between methods.

[0461] For each classifier, Fig.28A and Supplementary Table 17 (e.g. Fig.34 As shown in Figure 2, the predicted benign variants are retained in the primate test using a threshold of 50%. It also shows that the performance of primateAI is strong in terms of the number of species compared at the variant position. As long as there is sufficient mammalian conservation information, primateAI usually performs well, which is correct for most protein coding sequences (such as Fig.57 ).

[0462] Analysis of de novo variants in DDD studies

[0463] We obtained the published novel variants from the DDD study and the novel variants in the healthy member control group from the SSC autism study. The DDD study provides a confidence level for the novel variants, and due to the variant calling error, we excluded the DDD data set variants with a threshold value less than 0.1 as potential false positives. In one embodiment, there are a total of 3512 missense novel variants from DDD-affected individuals and 1208 missense novel variants from the healthy control group. The canonical transcript annotation used by UCSC for 99-vertebrate multiple sequence alignment is slightly different from the transcript annotation used by DDD, making the total number of missense variants slightly different. We evaluated the ability of these classification methods to distinguish novel missense variants in DDD-affected individuals from novel missense variants in the unaffected member control group from autism studies. For each classifier, we reported the P value of the Wilcoxon rank sum test of the difference between the prediction scores of the two assignments (such as Supplementary Table 17 (such as Fig.34 )).

[0464] To measure the accuracy of various classifiers in distinguishing benign from pathogenic variants in the same disease genes, we repeated the analysis on a subset of 605 genes enriched for de novo protein-truncating variants in the DDD group (P < 0.05, Poisson's exact test) (Supplementary Table 18). Among these 605 genes, based on the 3 / 1 enrichment above expectation for de novo missense variants, we estimated that two-thirds of the de novo variants in the DDD dataset were pathogenic and one-third were benign. We assumed minimal incomplete penetrance and assumed that de novo missense variants in healthy controls were benign. For each classifier, we identified a threshold that yielded the same number of benign or pathogenic predictions as the empirical proportion observed in these data sets and used this threshold as a binary cutoff to assess the accuracy of each classifier in distinguishing de novo variants in cases versus controls. To construct receiver operating characteristic curves, we classified de novo DDD variants as pathogenic as true positive calls and classified de novo pathogenic variants in healthy controls as false positive calls. Because the DDD dataset includes one-third of benign de novo variants, the area under the curve (AUC) for a theoretical perfect classifier is less than 1. Therefore, a classifier that perfectly classifies benign and pathogenic variants would classify 67% of de novo variants in DDD patients as true positives, 33% of de novo variants in DDD patients as false negatives, and 100% of de novo variants in controls as true negatives, yielding a maximum possible AUC of 0.837 (e.g. Fig.29A and 29B And Supplementary Table 19 Fig.35 shown).

[0465] Novel candidate gene discovery

[0466] We tested for enrichment of de novo variants in genes by comparing the number of observed de novo mutations to the number expected under a null mutation model. We repeated the enrichment analysis from the DDD study and reported genes with novel significance genome-wide when counting only de novo missense variants with a primateAI score greater than 0.803. We adjusted the genome-wide expectation of de novo damaging missense variants by the fraction of missense variants that met a primateAI threshold greater than 0.803 (approximately one fifth of all possible missense variants genome-wide). According to the DDD study, four tests were required for each gene, one for protein truncation enrichment and one for protein-altering de novo mutation enrichment, both of which were used only in the larger meta-analysis of the DDD cohort and the neurodevelopmental triple sequencing cohort. Protein-altering de novo mutation enrichment was combined with a cluster test for missense de novo mutations in coding sequences by the Fisher method (Supplementary Tables 20, 21). The P value for each gene was taken from the minimum of the four tests, and significance was determined genome-wide for P values ​​less than 6.757 × 10-7. (α=0.05, four tests using 18,500 genes).

[0467] ClinVar classification accuracy

[0468] Since most existing classifiers are trained directly or indirectly with ClinVar content, such as using prediction scores from classifiers trained with ClinVar, we restricted our analysis of the ClinVar dataset to only ClinVar variants added since 2017. Due to the large overlap between recent ClinVar variants and other databases, we filtered out variants at common allele frequencies (greater than 0.1%) in ExAC, or variants in HGMD (Human Gene Mutation Database), LOVD (Leiden Open Variant Database), or Uniprot (Universal Protein Resource). After excluding variants annotated with uncertain significance and variants with conflicting annotations, 177 missense variants annotated as benign and 969 missense variants annotated as pathogenic remained. We scored ClinVar variants using deep learning networks and other classification methods. For each classifier, we identified a threshold that produced the same number of benign or pathogenic predictions as the empirical proportion observed in these data sets and used this threshold as a binary cutoff to assess the accuracy of each classifier (e.g. Fig.31A and 31B shown).

[0469] The impact of increasing the amount of training data and using different sources of training data

[0470] To evaluate the effect of the amount of training data on the performance of deep learning networks, we randomly sampled a subset of variants from 385,236 primate variants and common human variants from a labeled benign training set, keeping the basic deep learning network architecture the same. To show that variants from each individual primate species are beneficial to classification accuracy, while variants from each individual mammalian species reduce classification accuracy, according to one embodiment, we trained the deep learning network using a training dataset that included 83,546 human variants plus a certain number of randomly selected variants from each species, keeping the underlying network architecture the same. The constant number of variants added to the training set (23,380) is the total number of available variants in the species with the fewest missense variants, such as bonobos. We repeated the training process 5 times to obtain the median performance of each classifier.

[0471] Saturation of all possible human missense mutations with increasing numbers of sequenced primate populations

[0472] By simulating variants based on the common missense variants observed in humans in ExAC (with allele frequencies greater than 0.1%), we investigated the predicted saturation of all 70 million possible human missense mutations with common variants from 504 existing primate species. For each primate species, we simulated four common missense variants observed in humans (83,500 missense variants with allele frequencies greater than 0.1%), since humans have approximately half the number of variants per individual as other primate species, and filtered out 50% of human missense variants by purifying selection at allele frequencies greater than 0.1% (e.g., Fig.49A ).

[0473] To simulate the fraction of common missense variants in humans (allele frequency greater than 0.1%) found in an increasing number of human cohorts examined (e.g. Fig.56 ), we sampled genotypes based on ExAC allele frequencies and report on a subset of common variants that were observed at least once in these simulation groups.

[0474] In one embodiment, for practical application of the primateAI score, a control group (e.g. Fig.21D ), based on the enrichment of new variant cases, in variants with dominant inheritance pattern, the preferred threshold is greater than 0.8, which may be pathogenic, less than 0.6, which may be benign, and between 0.6-0.8, which may be intermediate, while in genes with recessive inheritance pattern, the preferred threshold is greater than 0.7, which may be pathogenic, and less than 0.5, which may be benign.

[0475] Figure 2 Represents an example architecture of a deep residual network for pathogenicity prediction, referred to in this paper as primateAI. Figure 2In the above, 1D refers to a one-dimensional convolutional layer. The predicted pathogenicity ranges from 0 (benign) to 1 (pathogenic). The network takes as input the human amino acid (AA) reference sequence and variable sequence (51AAs) centered on the variant, the position weight matrix (PWM) conservation spectrum calculated from 99 vertebrate species, and the secondary structure and solvent accessibility prediction deep learning network to predict the three-state protein secondary structure (H for helix, B for beta sheet, C for coiled) and the three-state solvent accessibility (B for buried, I for intermediate, and E for exposed).

[0476] Figure 3 Schematic diagram of primateAI, which is a deep learning network architecture for pathogenicity classification. The model input includes 51 amino acids (AA) of the flanking sequence of the reference sequence and the sequence replaced by the variant, the conservation represented by three 51-AA-length position weight matrices from primate, mammal and vertebrate alignments, and the input of the pre-trained secondary structure network and solvent accessibility network (also 51AA in length).

[0477] Figure 4A , 4B and 4C are Supplementary Table 16, which show the details of the example model architecture of the pathogenicity prediction deep learning model primateAI. Shape specifies the shape of the output tensor of each layer of the model, and activation is the activation given to the neurons of this layer. The model input refers to the position-specific frequency matrix of the flanking amino acid sequence around the variant (51AA long, 20 deep), the one-hot encoded human reference sequence and variable sequence (51AA long, 20 deep), and the secondary structure and solvent accessibility model (51AA long, 40 deep).

[0478] The examples described use one-dimensional convolutions. In other embodiments, the model may use different types of convolutions, such as 2D convolutions, 3D convolutions, dilated or atrous convolutions, transposed convolutions, separable convolutions, and depthwise convolutions. Certain layers also use ReLU activation functions that greatly speed up the convergence of stochastic gradient descent relative to saturating nonlinearities such as S-curves or hyperbolic tangents. Other examples of activation functions used in the disclosed technology include parametric ReLUs, leaky ReLUs, and exponential linear units (ELUs).

[0479] Some layers also use batch normalization (Ioffe and Szegedy 2015). In batch normalization, during training, the distribution of each layer in a convolutional neural network (CNN) is changed so that the distribution is different from one layer to another. This reduces the convergence rate of the optimization algorithm. Batch normalization is a technique to address this problem. The input of a batch normalization layer is annotated with x and its output is annotated with z. Batch normalization applies the following transformation of x:

[0480]

[0481] Batch Normalization performs mean-variance normalization of the input using μ and σ and linearly scales and shifts the unique input using γ and β. The normalization parameters μ and σ for the current layer of the training set are calculated using the exponential moving average method. In other words, these parameters are non-trainable parameters. In contrast, gamma and beta are trainable parameters. During training, the values ​​of μ and σ are used in the forward method during inference.

[0482] Figure 5 and 6 A deep training network for predicting secondary structure and solvent accessibility of proteins was demonstrated. The model input was a conserved position weighted matrix generated by RaptorX software (for training protein database sequences) or 99 vertebrate alignments (for training and perturbation of human protein sequences). The output of the second to last layer, 51 AAs in length, was the input to the deep learning network for pathogenicity classification.

[0483] Fig. 7A and 7B To supplement Table 11, an example model architecture details for a 3-state secondary structure prediction deep training (DL) model is shown. Shape specifies the shape of the output tensor of each layer of the model, and activation is the activation given to the neurons of this layer. The model input is a position-specific frequency matrix (55AA long, 20 deep) for the amino acid sequence flanking the variant.

[0484] Fig. 8A and 8B To supplement Table 12, an example model architecture details for a 3-state solvent accessibility prediction deep learning model is shown. Shape specifies the shape of the output tensor of each layer of the model, and activation is the activation given to the neurons of that layer. The model input is a position-specific frequency matrix (51AA long, 20 deep) for the amino acid sequence flanking the variant.

[0485] Fig. 20 Figure 3. Predicted pathogenicity scores for each amino acid position in the SCN2A gene, with annotations for key functional domains. The bottleneck primateAI score for missense substitutions at each amino acid position is plotted along the gene.

[0486] Fig.21D Figure 2. Comparison of classifier predictions for benignity on a test set of 10,000 common primate variants held out for training. The y-axis represents the percentage of primate variants correctly classified as benign after normalizing the threshold of each classifier to its 50% score on a set of 100,000 random variants with matching mutation rates.

[0487] Fig.21EThe distribution of primateAI prediction scores for de novo missense variants in decipherable developmental disorders (DDD) patients compared to unaffected members with corresponding Wilcoxon rank sum values ​​P is illustrated.

[0488] Fig.21F Depicted is a comparison of classifiers when isolating de novo missense variants in DDD cases compared to controls. Wilcoxon rank sum test P values ​​for each classifier are shown.

[0489] Fig.22A , 22B , 22C, 22D, and 22E illustrate the accuracy of the classifier among 605 DDD genes, with P less than 0.05. Fig.22A Figure 5 shows an above-expected enrichment of de novo missense mutations in affected individuals from the DDD group with 605 associated genes that were significant (P less than 0.05) for de novo protein-truncating variants. Fig. 22B Figure 3 Distribution of primateAI prediction scores for de novo missense variants in DDD patients versus unaffected members of 605 relevant genes, with corresponding Wilcoxon rank sum P values.

[0490] Fig. 22C Figure 2 shows a comparison of classifiers for separation of de novo missense variants in DDD individuals versus controls in 605 genes. The y-axis shows the Wilcoxon rank sum test value P for each classifier.

[0491] Fig.22D A comparison of various classifiers is depicted using receiver operating characteristic curves, with each classifier represented by the AUC.

[0492] Fig.22E The classification accuracy and AUC of each classifier are explained. Using thresholds, based on Fig.22A The enrichment shown shows that the classifier predicts the same number of pathogenic and benign variants as expected, and the classification accuracy shown is the average of the true positive rate and true negative rate. The dashed line shows the maximum achievable AUC for a perfect classifier, considering that 33% of the DDD de novo missense variants represent background.

[0493] Fig.23A , 23B , 23C and 23D show the effect of the data used to train the classifier accuracy. The deep learning network trained with an increased number of common primate and human variants was based on the full dataset (385,236 variants in total). Fig.23A , the classification performance of each network is based on the accuracy of 10,000 preserved primate variants and de novo variants in DDD cases versus controls.

[0494] According to one embodiment, Fig. 23B and 23C The performance of a network trained on a dataset consisting of 83,546 common human variants and 22,380 variants from a single primate or mammalian species is shown in Table 1. Fig. 23B ) and novel missense variants in DDD cases versus controls (as shown Fig.23C ), with results shown for each network trained with different sources of common variation.

[0495] Fig.23D Depicts expected saturation of all possible human benign missense positions for state-consistent common variants (>0.1%) across 504 extant primate species. The y-axis represents the fraction of human missense variants observed in at least one species, with CpG missense variants in green and all missense variants in blue. To simulate common variants in each primate species, we sampled the set of all possible single-nucleotide alternatives by substitution, which matched the observed trinucleotide context assignments for common human variants (>0.1% allele frequency) in ExAC.

[0496] Fig.24 Correction for the effect of sequencing extent in identifying common primate variants is demonstrated. The probability of observing a given variant in a nonhuman primate species is negatively correlated with the sequencing depth at that position in the ExAC / gnomAD exome datasets. In contrast, smaller gnomAD read depths do not affect the probability of observing a common human variant at that position (>0.1% allele frequency) because the large number of sequenced human exomes ensures that most common variants are identified. When selecting matching variants for each primate variant used to train the network, the effect of sequencing depth adjusts the probability of selecting variants, matching on trinucleotide context to control for mutation rate and gene conversion.

[0497] Fig.25A , 25B , 25C and 26 describe the protein motifs recognized by the neural network of the present disclosure. Fig.25A , 25B and 25C illustrate the neural network recognition of protein domains, showing the evaluation of primateAI scores for variants at each amino acid position in three different protein domains. Fig.25A COL1A2 collagen strands with glycine in the repeated GXX motif are highlighted. Clinically recognized mutations in collagen genes are mainly due to missense mutations of glycine in the GXX repeats, as these mutations interfere with the normal integration of collagen and have a strong dominant inhibitory effect. Fig.25B In the figure, the active site of the IDS sulfatase is highlighted, including the post-translational modification of the active site cysteine ​​to formylglycine. Fig.25C In Figure 1, the bHLHzip domain of the MYC transcription factor is shown. The base domain contacts the DNA through positively charged arginine and lysine residues (highlighted), which interact with the negatively charged sugar-phosphate backbone. The leucine zipper domain includes leucine residues separated by 7 amino acids (highlighted) that are essential for dimerization.

[0498] Fig.26 It is a line graph showing the interference effect of each position in and around the variant on the variant prediction deep learning score. We systematically zeroed the input of the adjacent amino acids around the variant (positions -25 to +25) and measured the change in pathogenicity predicted by the neural network. The figure shows the average change in the predicted pathogenicity score of 5000 randomly selected variants at each adjacent amino acid position.

[0499] Fig. 27 The association model of the weighted simulated BLOSUM62 and Grantham score models is described. The weighted correlation model of the first three layers of the secondary structure deep learning network shows the correlation between amino acids similar to the BLOSUM62 and Grantham score matrices. The left heat map shows the correlation of the parameter weights of the first convolutional layer, which is immediately followed by the two initial upsampling layers of the secondary structure deep learning network between amino acids encoded by the one-hot representation. The middle heat map shows the BLOSUM62 scores between amino acid pairs. The right heat map shows the Grantham distance between amino acids. The Pearson correlation between the deep learning weights and the BLOSUM62 score is 0.63 (P = 3.55 × 10–9). The correlation between the deep learning weights and the Grantham score is –0.59 (P = 4.36 × 10–8). The correlation between the BLOSUM62 and Grantham score is –0.72 (P = 8.09 × 10–13).

[0500] Fig.28A , 28B and 28C represent the performance evaluation of the deep learning network primateAI and other classifiers. Fig.28ADepicted is the accuracy of the deep learning network primateAI in predicting benign outcomes on a test set of 10,000 primate variants held out for training relative to other classifiers, including SIFT, PolyPhen-2, CADD, REVEL, M-CAP, LRT, MutationTaster, MutationAssessor, FATHMM, PROVEAN, VEST3, MetaSVM, MetaLR, MutPred, DANN, FATHMM-MKL_coding, Eigen, GenoCanyon, integrated_fitCons, and GERP. The y-axis represents the percentage of primate variants classified as benign based on normalizing the threshold of each classifier to its 50% score, with mutation rate and gene conversion controlled for trinucleotide context using a set of 10,000 randomly selected variants matched to primate variants.

[0501] Fig.28B Depicted is a comparison of the performance of the primateAI network in separating de novo missense variants in DDD cases vs controls, with a list of 20 existing methods above. The y-axis represents the P value of the Wilcoxon rank sum test for each classifier.

[0502] Fig.28C The performance of the primateAI network in separating de novo missense variants in DDD cases versus unaffected controls in 605 disease-associated genes is compared with 20 existing methods. The y-axis represents the P value of the Wilcoxon rank sum test for each classifier.

[0503] Fig.29A and 29B The distribution of prediction scores of the four classifiers is illustrated. Shown are histograms of the prediction scores of the four classifiers: SIFT, PolyPhen-2, CADD, and REVEL, for de novo missense variants in DDD cases versus unaffected controls, with corresponding Wilcoxon rank sum P values.

[0504] Fig. 30A , 30B and 30C compared the accuracy of the primateAI network and other classifiers in separating pathogenic and benign variants in 605 disease-related genes. Fig. 30A The middle scatter plots represent the performance of each classifier on DDD cases versus controls (y-axis) and benign prediction accuracy on the holdout primate dataset (x-axis). Fig. 30BDifferent classifiers for separating de novo missense variants in cases versus controls in 605 genes were compared, represented by receiver operating characteristic (ROC) curves, with the area under the curve (AUC) representing each classifier. Fig. 30C Represents the classification accuracy and AUC of the primateAI network, Fig.28A \28Bji 28C lists 20 classifiers. Using thresholds, based on Fig.22A The enrichment shown, the classifier can predict the same number of pathogenic and benign variants as expected, and the classification accuracy shown is the average of the true positive rate and true negative rate. Assuming that 67% of the de novo missense variants in DDD cases are pathogenic variants and 33% are diamond variants, while 100% of the de novo missense variants in controls are benign, the dashed line represents the maximum achievable AUC for a perfect classifier.

[0505] Fig.31A and 31B The correlation between the performance of the classifier on ClinVar variants selected by human experts and the performance of the classifier on empirical datasets is demonstrated. Fig.31A The classification accuracy of each of the 20 other classifiers on the 10,000 retained primate variants (x-axis) ClinVar variants and the primateAI network trained on human only or human plus primate data is shown. The Spearman correlation coefficient rho and the correlation value P are shown. To limit the evaluation to data not used to train the classifier, we only used ClinVar variants added between January and October 2017 and excluded common human variants from ExAC / gnomAD (greater than 0.1% allele frequency). Using a threshold, the classifier predicts the same number of pathogenic and benign variants observed by ClinVar, and the ClinVar classification accuracy is the average of the true positive rate and the true negative rate.

[0506] Scatter Fig.31B Classifier preparation showing ClinVar variants (y-axis) and the full dataset (x-axis) for DDD cases versus controls for each of the 20 classifiers and the primateAI network trained using either human-only or human plus primate data.

[0507] Fig.32 To supplement Table 14, the performance of the 3-state secondary structure and 3-state solvent accessibility prediction models on annotated samples from the protein database is shown, using 3637 unrelated protein sequences for training, 400 for validation, and 500 for testing. Only proteins with less than 25% sequence similarity in the protein database were selected. We report the accuracy of the deep learning network as the performance metric because the three classes are severely unbalanced in terms of secondary structure or solvent accessibility.

[0508] Fig.33 Supplementary Table 15 shows the performance comparison of deep learning networks using secondary structure markers annotated from the DSSP database, compared to deep learning networks using predicted secondary structure markers.

[0509] Fig.34 To supplement Table 17, the accuracy values ​​of the 10,000 retained primate variants evaluated and the P values ​​of the novel variants in DDD cases versus controls for each of the 20 classifiers are shown. The primateAI model with only human data is a deep learning network using a labeled benign training dataset that includes only common human variants (83.5K variants in the 0.1% of the major hospital population), while the primateAI model with human and primate data is a deep learning network trained on the full set of 385K labeled benign variants, including common human variants and primate variants.

[0510] Fig.35 Supplementary Table 19 shows a comparison of the performance of different classifiers for de novo variants in the DDD case-control dataset, limited to 605 disease-associated genes. To normalize the different methods, for each classifier we identified a threshold at which the classifier would predict the same number of pathogenic and benign variants based on enrichment in DDD and controls. The classification accuracy shown is the average of the true positive and true negative error rates at that threshold.

[0511] Fig.49A , 49B , 49C, 49D, and 49E describe the missense / synonymous mutation rates across the human allele spectrum. Fig.49A Represents missense and synonymous variants observed in 123,136 human individuals in the ExAC / gnomAD database divided into four categories by allele frequency. The shaded grey stripes represent the number of synonymous variants in each category, and the dark green stripes represent missense variants. The height of each stripe is reduced to the number of synonymous variants in each allele frequency category, and after adjusting for mutation rate, the number of missense / synonymous mutations and the missense / synonymous mutation rate are displayed. Fig.49B and 49C Describes variants common in chimpanzees (e.g. Fig.49B ) and chimpanzee single variants (e.g. Fig.49C ) Allelic spectrum of human missense and synonymous variants with consistent IBS status (IBS). The red box indicates the depletion of chimpanzee missense variants at common allele frequency (greater than 0.1%) relative to the rare human allele frequency (less than 0.1%) and also indicates the chi-squared (χ2) test value P.

[0512] Fig.49DRelative to ClinVar variants in a panel of 30 human samples sampled from ExAC / gnomAD allele frequencies (middle row) and variants observed in primates (bottom row), Fig.49E The number of benign and pathogenic missense variants in the entire ClinVar database (top row) is shown. Confirmed conflicting benign and pathogenic variants and variants annotated only as undetermined meaning are excluded.

[0513] Fig.50A , 50B , 50C and 50D indicate purifying selection for missense variants consistent with the status in other species. Fig.50A Depicts the allele spectrum of human missense and synonymous variants that are consistent with variant status in four non-primate mammalian species (mouse, pig, goat, and cattle). Red boxes indicate attrition of missense variants at common human allele frequency (greater than 0.1%) and also indicate the chi-square (χ2) test value P.

[0514] Fig.50B A scatter plot showing the loss of missense substitutions observed in other species at common human allele frequencies (greater than 0.1%) compared to the evolutionary distance from the human species, expressed in units of average branch length (average number of substitutions per nucleotide position). Next to the species name, the total branch length between each species and humans is indicated. Loss values ​​for individual and common variants are shown for species with variant frequencies, with the exception of gorillas, which include related individuals.

[0515] Figure 50C illustrates the number of benign and pathogenic missense variants sampled from a set of ExAC / gnomAD allele frequencies (top row) relative to variants observed in primates (middle row) and variants observed in mice, pigs, goats, and cattle (bottom row). Determined conflicting benign and pathogenic variants and variants annotated only as undetermined sense were excluded.

[0516] Figure 50D is a scatter plot showing the depletion of fixed missense substitutions (expressed in units of average branch length) observed in pairs of closely related species at common human allele frequency (greater than 0.1%) compared to evolutionary distance from the human species.

[0517] Fig.51Represents the missense: synonymous mutation rate expected for the entire human allele spectrum in the absence of purifying selection. The shaded gray stripes represent the number of synonymous mutation variants and the dark green stripes represent the number of missense variants. The dotted line represents the baseline for synonymous mutation variant formation. The missense: synonymous mutation rate for each allele frequency category is indicated. According to one embodiment, the expected number of missense and synonymous mutations in each allele category is calculated by selecting intragenic variants from the ExAC / gnomAD dataset including 123,136 exomes and using these variants to estimate the partial variants expected to belong to each of the four allele categories based on the trinucleotide context of the variants used to control mutation rate and GC bias in genetic variation.

[0518] Fig.52A , 52B , 52C, and 52D describe the missense:synonymous mutation rates for CpG and non-CpG variants. Fig.52A and 52B represents the CpG variants in the entire human allele spectrum (e.g. Fig.52A ) and non-CpG variants (e.g. Fig.52A ) with all variants from the ExAC / gnomAD exome. Fig.52C and 52D represents the CpG variants in the entire human allele spectrum (e.g. Fig.52C ) and non-CpG variants (e.g. Fig.52D ) that is restricted only to human variants consistent with the common polymorphic state in chimpanzees.

[0519] Fig.53 , 54 and 55 illustrate missense:synonymous mutation rates for human variants consistent with six primate states. For ExAC / gnomAD variants consistent with variant states, the pattern of missense:synonymous mutation rates for human allele spectrum appears in chimpanzees, bonobos, gorillas, orangutans, rhesus monkeys, and marmosets.

[0520] Fig.56 Saturation of new common missense variants discovered by increasing the size of the investigated human cohort was simulated. In the simulations, the genotypes of each sample were sampled according to the gnomAD allele frequencies. The fraction of discovered gnomAD common variants was averaged over 100 simulations for each sample size between 10 and 10,000 variants.

[0521] Fig.57Represents the accuracy of primateAI in different conservation profiles in the genome. The x-axis represents the percent alignment of 51AA around the sequence with 99 vertebrate alignments. The y-axis represents the classification performance of primateAI accuracy for variants within each conserved region based on a test dataset of 10,000 conserved primate variants.

[0522] Fig.58 To supplement Table 5, the impact of a benign training dataset on labeling of common human variants and variants present in non-human primates is shown.

[0523] Fig.59 The effect of allele frequency on the expected missense:synonymous mutation rate is shown in Supplementary Table 8. The expected number of synonymous and missense variants was calculated based on the allele frequency spectrum of variants in regions at least 20-30n away from the gene using the trinucleotide context for mutation rate and gene conversion bias.

[0524] Fig.60 To supplement Table 9, ClinVar analysis is shown. According to one embodiment, missense variants with conflicting annotations are filtered out from variants downloaded from the ClinVar database constructed in October 2017 and variants with uncertain significance are excluded, leaving 17,775 benign variants and 24,853 pathogenic variants.

[0525] According to one embodiment, Fig.61 The number of missense variants from other species found in ClinVar is shown to supplement Table 10. Variants are required to be consistent with the corresponding human variant status and have identical nucleotides at two other positions in the reading frame to ensure the same coding structure.

[0526] Fig.62 TABLE 1 shows one implementation of 14 additional candidate genes found in intellectual disability that did not previously meet genome-wide significance thresholds in the original DDD study.

[0527] Fig.63 Table 2 shows one embodiment of the mean difference in Grantham scores between pathogenic variants and benign variants in ClinVar, which is twice the difference in de novo variants between DDD cases and controls in 605 disease-associated genes.

[0528] Data Generation

[0529] All coordinates used in this article refer to the human genome construct UCSC hg19 / GRCh37, including coordinates of variants in other species, and were mapped to hg19 using multiple sequence alignments using the methods described in this section. Human protein-coding DNA sequences and multiple sequence alignments to 99 vertebrate genomes were downloaded from the UCSC genome browser for use in the hg19 build (http: / / hgdownload.soe.ucsc.edu / goldenPath / hg19 / multiz100way / alignments / knownCanonical.exonNuc.fa.gz). For genes with multiple canonical gene annotations, the longest coding transcript was selected in this article.

[0530] We downloaded human exome polymorphism data from the Exome Aggregation Database (ExAC) / Genome Aggregation Database (gnomAD) v2.0, which collects whole-genome sequencing data (WES) of 123,136 individuals from eight subpopulations worldwide (123: / / gnomad.broadinstitute.org / ). We excluded variants that were annotated in the ExAC VCF files, did not pass the default quality control filter, or fell outside the canonical coding region. To avoid the effects of balancing selection, we also excluded variants from the extended MHC region (chr6: 28,477,797-33,448,354) analyzed in primates. The Great Ape Genome Sequencing Project provided whole-genome sequencing data and genotypes for 24 chimpanzees, 13 bonobos, 27 gorillas, and 10 orangutans (including 5 from the Sumatran subspecies and 5 from the Bornean subspecies, which we collapsed for downstream analysis). For chimpanzee and bonobo studies, we provided genome sequences of 35 additional chimpanzees. However, because variants in these additional chimpanzees were not identified using the same methods as the Great Ape Genome Sequencing Project, we excluded them from the allele spectrum analysis and only used them for training the deep learning model. Variants from these primate diversity studies have been mapped to the human reference sequence (hg19). In addition, for marmosets and rhesus macaques, 16 rhesus macaque individuals and 9 marmoset individuals were used to determine the variation in the original genome sequences of these species, but individual-level information is not yet available.

[0531] The Great Ape Genome Sequencing Project4 provided whole-genome sequencing data and genotypes for 24 chimpanzees, 13 bonobos, 27 gorillas, and 10 orangutans (including 5 from the Sumatran subspecies and 5 from the Bornean subspecies, which we collapsed for downstream analyses). For studies on chimpanzees and bonobos, we provided genome sequences for an additional 35 chimpanzees. However, because variants in these additional chimpanzees were not identified using the same methods as the Great Ape Genome Sequencing Project, we excluded them from the allele spectrum analysis and only used them for training the deep learning model. Variants from these primate diversity studies have been mapped to the human reference sequence (hg19). In addition, for marmosets and rhesus macaques, 16 rhesus macaque individuals and 9 marmoset individuals were used to determine variation in the original genome sequences of these species, but individual-level information is not yet available.

[0532] For comparison with other primates and mammals, we also downloaded SNPs from other species from dbSNP, including those of rhesus macaques, marmosets, pigs, cattle, goats, mice, chickens, and zebrafish. dbSNP also includes other orangutan variants, which we only used for training the deep learning model because individual genotype information is not available for allele spectrum analysis. We dropped other species, such as dogs, cats, or sheep, because dbSNP only provides a limited number of variants for these species.

[0533] To map variants to humans, we used a multi-species alignment of 99 vertebrates to ensure 1:1 mapping with human protein coding region orthologs. Mapping variants using orthologous multi-species alignments is critical to eliminate artifacts caused by pseudogenes or retrotransposition sequences, which occur when mapping SNPs directly between species using tools that map many-to-one, such as liftOver. If the genome structure of the species in dbSNP did not match the genome structure of the species in the 99 vertebrate multiple sequence alignment, we used liftOver to update the variant to the genome structure used in the multiple sequence alignment. We accepted variants as state-consistent if they appeared in a reference / alternative manner, e.g., if it was G in the human reference sequence and A in the alternative allele, the variant was considered to be in the same state in other species with A in the reference sequence and G in the alternative allele. To ensure that variants had the same predicted protein coding outcome in humans and other species, we required that the other two nucleotides in the codon were identical in both missense and synonymous variants between species. The polymorphisms for each species included in the analysis are listed in Supplementary Data File 1 and detailed metrics are shown in Supplementary Table 1 .

[0534] To ensure that the variants in each dbSNP submission batch were of high quality and aligned to human variants, we calculated the missense / synonymous ratio for each batch and confirmed that this was below the expected ratio of 2.2:1; most species have ratios below 1:1, especially zebrafish and mouse, which have very large effective populations. We excluded two batches of SNPs from cows with unusually high missense / synonymous ratios (snpBatch_1000_BULL_GENOMES_1059190.gz with a ratio of 1.391 and snpBatch_COFACTOR_GENOMICS_1059634.gz with a ratio of 2.568) from further analysis. The average missense / synonymous ratio for the remaining cow batches was 0.8:1.

[0535] Corrected allele frequency for missense / synonymous ratio, mutation rate, genetic drift and GC-based genes The impact of conversion

[0536] In addition to the effects of purifying selection, the observed attrition of missense variants in humans at high allele frequencies may also be influenced by factors unrelated to natural selection. The probability of a neutral mutation arising at a particular allele frequency in a population is a function of mutation rate, gene conversion, and genetic drift, and these factors could potentially introduce biases in missense / synonymous mutation rates across the allele spectrum even in the absence of selective forces.

[0537] To calculate the expected missense / synonymous mutation ratios in each allele frequency category in the absence of protein-coding selection, we selected variants within intronic regions 31–50 bp upstream and 21–50 bp downstream of each exon. These regions were chosen to be sufficiently distant to avoid the influence of extended splice motifs. Because these regions are near the edge of the exome capture sequence of the ExAC / gnomAD exomes, we removed all chrX regions to ensure clear identification of variants and excluded regions with an average read depth of <30. Each variant and its immediately upstream and downstream nucleotides belonged to one of 64 trinucleotide contexts. If we mutated the intervening nucleotide to three other bases, a total of 64 × 3 = 192 trinucleotide configurations were possible. Since trinucleotide configurations and their reverse complements are equivalent, 96 trinucleotide contexts were valid. We observed a very strong effect of trinucleotide context on mutation rate and a smaller effect on GC-biased gene conversion, making trinucleotide context valid for constructing these variants.

[0538] In these intronic regions, we selected each variant from 126,136 ExAC / gnomAD exomes and classified them into 4 × 192 categories based on four categories of allele frequency (singleton, more than singleton, >0.1%) and 192 trinucleotide backgrounds. We normalized the number of variants observed in 4 × 192 categories (allele frequency × trinucleotide background) by dividing the total number of possible variants by this trinucleotide background (obtained by replacing each nucleotide in the intronic sequence in three different ways). Thus, for each of the 192 trinucleotide backgrounds, we have obtained the expected proportion of variants falling into the four allele frequency categories in the absence of protein-coding selection. This implicitly models the effects of mutation rate, GC-biased gene conversion, and genetic drift due to differences in trinucleotide backgrounds (Supplementary Table 7).

[0539] To obtain the expected missense / synonymous mutation rate in each allele frequency category, we calculated the total number of synonymous and missense mutations obtainable by single-nucleotide substitutions in the human genome and assigned them to each of the 192 trinucleotide contexts. For each case, we used a 4 × 192 table to calculate the number of each variant expected to fall into the four allele frequency categories. Finally, we summed the number of synonymous and missense variants in the 192 trinucleotide contexts to obtain the expected total number of synonymous and missense variants in each of the four allele frequency categories ( Fig.51 and Supplementary Table 8 ( Fig.59 )).

[0540] The expected missense / synonymous ratio is nearly constant across the allele spectrum and close to the 2.23:1 ratio expected for de novo variants in the absence of natural selection, with the exception of singleton variants, for which the expected missense / synonymous ratio is 2.46:1. This suggests that variants with the singleton allele frequency category in ExAC / gnomAD have an expected missense / synonymous ratio that is 10% higher than de novo mutations in the default case due to factors unrelated to protein-coding selective pressure (mutation rate, gene conversion, genetic drift). To address this issue, we reduced the missense / synonymous ratio of singletons by 10% in our allele frequency analysis ( Fig.49A , 49B, 49C, 49D and 49E, Fig.50A , 50B, 50C, and 50D). This small adjustment reduces the expected missense mutation loss of common human variants present in primates and other mammals by about 3.8% (e.g. Fig.49A , 49B, 49C, 49D and 49E and Fig.50A , 50B, 50C, and 50D). The higher missense / synonymous mutation rate for the singleton variants is due to the higher allele frequency of transition mutations (which are more likely to produce synonymous changes) due to their higher mutation rate than transversion mutations (which are more likely to produce missense changes).

[0541] Furthermore, this explains the missense / synonymous ratio of 2.33:1 observed for singleton variants in ExAC / gnomAD, which exceeds the expected ratio for de novo mutations, which is 2.23:1. After accounting for the effect of allele spectrum on the missense / synonymous ratio, this actually reflects a 5.3% reduction in singleton variants compared with expectation, which may be due to selection against pathogenic missense variants with a de novo dominant inheritance pattern. Indeed, when we consider only haploinsufficient genes with high loss-of-function cases (pLI>0.9), the missense / synonymous ratio of ExAC / gnomAD singleton variants is 2.04:1, indicating an approximate loss of haploinsufficient genes of ∼17%. This result is consistent with previous estimates that 20% of missense variants are equivalent to loss-of-function mutations, assuming some degree of incomplete penetrance.

[0542] Because their mutation rates vary widely, we also specifically investigated missense / synonymous mutation rates for CpG and non-CpG variants across the human allele spectrum ( Fig.52A , 52B, 52C, and 52D). We demonstrated that human variants identical to the common polymorphic state in chimpanzees have nearly invariant missense / synonymous mutation rates across the allele spectrum for both CpG and non-CpG mutations.

[0543] Attrition of human missense variants that are state-consistent with polymorphisms in other species

[0544] To assess whether variants from other species are tolerated at common allele frequencies in humans (>0.1%), we identified human variants that were consistent with the status of the variant in other species. For each variant, we assigned it to one of four categories based on its allele frequency in the human population (singleton, more than singleton, >0.1%) and estimated the reduction in missense / synonymous ratios (MSRs) between rare (<0.1%) and common (>0.1%) variants. The depletion of state-consistent missense variants at common human allele frequency (>0.1%) indicates the proportion of deleterious variants in other species that could have been filtered out by natural selection.

[0545]

[0546] The missense / synonymous mutation rate and attrition percentage were calculated for each species and displayed in Fig.50B and Supplementary Table 2. In addition, for chimpanzee common variants ( Fig.49A ), chimpanzee singleton variant ( Fig.49C ) and mammalian variants ( Fig.50A), we performed chi-square (χ2) homogeneity tests on 2 × 2 contingency tables to test whether the difference in missense / synonymous mutation rates between rare and common variants was significant.

[0547] Because sequencing was performed on only a small number of individuals in the Great Apes Diversity Project, we used the human allele frequency spectrum in EXAC / gnomAD to estimate the proportion of sampled variants that are rare (<0.1%) or common (>0.1%) in chimpanzee populations. We sampled a cohort of 24 individuals based on ExAC / gnomAD allele frequencies and identified missense variants that were observed once or more in that cohort. Variants observed more than once had a 99.8% probability of being in the common population (>0.1%), whereas variants observed only once had only a 69% probability. Fig.49B and 49C In this study, we show that because some chimpanzee singleton variants are rare deleterious mutations, we observed attrition of singleton chimpanzee variants at high human allele frequencies, but not at common chimpanzee variants. In a cohort of 24 individuals, about half of the chimpanzee variants were observed only once, and about half were observed multiple times.

[0548] To confirm that the observed depletion of missense variants in more ancient mammals was not due to confounding effects of more conserved genes, and thus to more accurately align, we repeated the above analysis but restricted it to genes with an average nucleotide identity >50% in a multiple sequence alignment of 11 primates and 50 mammals compared to humans (see Supplementary Table 3). Excluding approximately 7% of human protein-coding genes from this analysis had little effect on the results.

[0549] Fixed substitutions among primates, mammals, and distantly related vertebrates

[0550] To ensure that our results using dbSNP variants were not affected by variant data issues or domestication artifacts (because most of the species selected from dbSNP are domesticated), we also repeated the analysis of intraspecific polymorphisms using fixed substitutions from closely related species pairs. We downloaded a phylogenetic tree of 100 vertebrate species from the UCSC genome browser (http: / / hgdownload.soe.ucsc.edu / goldenPath / hg19 / multiz100way / hg19.100way.commonNames.nh), and their phylogenetic distances were expressed as branch lengths (average number of nucleotide substitutions per position). We selected closely related species pairs (branch length < 0.25) for further analysis. To identify fixed substitutions between closely related species pairs, we downloaded coding regions of multiple sequence alignments of 99 vertebrate genomes with humans and coding regions of alignments of 19 mammalian (16 primate) genomes with humans from the UCSC genome browser. Because some primates, such as bonobos, are not present in the 99 vertebrate alignments (http: / / hgdownload.soe.ucsc.edu / goldenPath / hg38 / multiz20way / alignments / knownCanonical.exonNuc.fa.gz), a multi-species alignment of 19 mammals was necessary. As shown in Figure 50D and Supplementary Table 4, we obtained a total of 15 closely related species pairs, including 5 primate species pairs.

[0551] We performed multiple sequence alignments of 19 mammalian or 99 vertebrate genomes to the human genome within canonical coding regions and listed the substituted nucleotides acquired between each selected pair of vertebrates in Supplementary Data File 2. These substituted nucleotides were mapped to the human genome, requiring that the other two nucleotides in the codon did not change between humans and other species and that the variants were accepted as reference or otherwise. Using human variants with fixed substitution state concordance for the relevant species pair, we calculated missense / synonymous mutation rate ratios for variants in rare (<0.1%) and common (>0.1%) allele frequency categories, thereby obtaining the proportion of nucleotides substituted under negative selection, as shown in Supplementary Table 4.

[0552] ClinVar analysis of polymorphism data for humans, primates, mammals, and other vertebrates

[0553] To examine the clinical impact of variants consistent with status in other species, we downloaded the variant summary release of the ClinVar database (ftp: / / ftp.ncbi.nlm.nih.gov / pub / clinvar / clinvar_20171029.vcf.gz released on November 2, 2017)12. This database contains 324,698 variants on the hg19 genome structure, of which 122,884 variants are missense single-nucleotide variants that map to our list of protein-coding genes (Supplementary Table 9). Most variants in the ClinVar database do not have missense consequences and were therefore excluded. Next, we filtered variants with interpretations that contradicted pathogenicity and retained only variants with benign, likely benign, pathogenic, and likely pathogenic annotations. We merged “variants” with “benign” or “likely benign” annotations into one category and merged “variants” with “pathogenic” or “likely pathogenic annotations”. Following the filtering steps shown in Supplementary Table 9, a total of 24,853 variants were in the pathogen category and 17,775 variants were in the benign category; the rest were excluded because they were variants of unknown significance or with conflicting annotations.

[0554] To obtain a baseline of ClinVar missense variants in the human population, we investigated ClinVar missense variants in 30 individuals sampled from the ExAC / gnomAD allele frequency cohort. The size of this cohort was chosen to roughly reflect the number of individuals sequenced in the Primate Diversity Project study. We report on 100 such simulated variants in the 30 human cohorts ( Fig.49E ). Because common human variants with benign consequences have been systematically annotated in ClinVar by curators, we excluded variants with an allele frequency greater than 1% to avoid this bias.

[0555] We analyzed ClinVar variants with consistent status in primates, mammals, and other vertebrates. The number of benign and pathogenic variants in each species is shown in Supplementary Table 10. The number of ClinVar variants present in humans, primates, and older mammals is shown in Fig.49E and 50B The chi-square (χ2) test for homogeneity showed differences in the ratios of benign to pathogenic variants.

[0556] Generation of benign variants for model training

[0557] Variants that are common in the human population are largely unaffected, except for rare cases of founder effects or balancing selection, making them suitable as benign training datasets for machine learning that are not subject to human interpretation biases. We used allele frequency data from 123,136 exomes in the ExAC / gnomAD database (version v2.0), excluding variants that did not pass the filter, which gave us 83,546 missense variants with an overall population allele frequency >= 0.1% of canonical protein-coding transcripts.

[0558] Based on our previous results that variants in primates are largely benign in humans, we created a benign training dataset for machine learning that included common human variants (>0.1% allele frequency), variants in chimpanzees, bonobos, gorillas, and orangutans from the Great Ape Diversity Project and other primate sequencing projects, and rhesus macaques, orangutans, and marmosets from dbSNP. A total of 301,690 unique primate variants were added to the benign training set according to one implementation. The number of benign training variants from different sources is shown in Supplementary Table 5.

[0559] It is important to note that although the majority of primate variants are common in their respective populations, a minority are rare variants. Because the number of sequenced individuals for nonhuman primates is limited, we expect the set of variants identified to generally represent common variation. Indeed, we found that the missense / synonymous ratio of variants in each primate species was less than half of the 2.23:1 ratio expected for de novo mutations, suggesting that these are the majority of common variants that have passed the filter. Furthermore, for the chimpanzee cohort, we estimate that approximately 84% of the identified variants are present in their respective populations at a common allele frequency (>0.1%). Since approximately 50% of the newly arising missense variants are acquired by amplifying variants at common human allele frequency (>0.1%) ( Fig.49A ) were filtered by purification selection, so this number is similar to The rare variants in the human genome correspond to 8.8% of the observed missense variants in humans, which are state-consistent with the observed primate variants ( Fig.49D ).

[0560] Applying the estimate that approximately 20% of human missense mutations are equivalent to loss of function, primate variants are expected to contain 3.2% fully pathogenic mutations, 91.2% benign mutations (tolerable allele frequency >0.1%), and 5.6% intermediate mutations that do not completely abrogate gene function but are sufficient to be filtered out at the common allele frequency (>0.1%). Despite the known deficiencies of this training dataset, the classification accuracy of the deep learning network was much better when trained on a benign training dataset containing both common human variants and primate variants compared to only common human variants. Therefore, at current classification accuracy, the amount of training data available is a more stringent limitation. As a large number of individuals within each primate species have been sequenced, it is possible to prepare training datasets containing a higher proportion of common primate variants, thereby reducing the contamination of the training dataset with pathogenic variants and further improving classification performance.

[0561] Generate unlabeled variants to supplement benign training datasets

[0562] All possible missense variants were generated by substituting three other nucleotides at each base position in the canonical coding region. We excluded variants observed in 123,136 exomes from ExAC / gnomAD, as well as variants in start or stop codons. A total of 68,258,623 unlabeled variants were generated. We assigned each unlabeled variant to one of 96 different trinucleotide context categories. We trained a deep learning network using a semisupervised approach by sampling variants in the unlabeled dataset that matched variants in the benign dataset (according to the trinucleotide context) from this unlabeled dataset and training a classifier to distinguish between benign and unlabeled training instances.

[0563] Filtering of unmarked variants

[0564] By showing examples of benign and unlabeled variants along with flanking amino acid sequences, deep learning networks can help learn protein regions that are highly intolerant of mutations. However, the absence of common variants in regions of protein sequence could be due to a high degree of purifying selection or could be due to technical artifacts that prevent variants in that region from being identified. To correct for the latter, we removed variants from both the benign and unlabeled datasets in regions where the average coverage of the ExAC / gnomAD datasets was <1. Similarly, when unlabeled variants in the benign dataset matched primate variants during training, we excluded unlabeled variants from regions where primates had no orthologous comparable sequences to humans in multiple sequence alignments.

[0565] Primate variants are retained for validation and testing, as well as de novo variants in affected and unaffected individuals

[0566] To validate and test the deep learning network, we randomly sampled two sets of 10,000 primate variants for validation and testing, but we did not train them. The remaining primate variants were used together with common human variants (>0.1% allele frequency) as a benign dataset for training the deep learning network. In addition, we sampled two sets of 10,000 unlabeled variants that matched the retained primate variants for validation and test sets.

[0567] We used 10,000 held-out primate variants in the validation set and a matching set of 10,000 unlabeled variants, and by measuring the ability of the network to distinguish between the two sets, we were able to monitor the performance of the deep learning network. Once the network performance reached saturation, we were able to determine the stopping point of training and avoid overfitting.

[0568] We used the 10,000 primate variants held out in the test dataset to benchmark the deep learning network and 20 other classifiers. Because the score distributions of different classifiers vary greatly, we used these unlabeled variants to identify the 50th percentile threshold for each classifier. We benchmarked the proportion of variants in the 10,000 held-out primate variant test set that were classified as benign at the 50th percentile threshold of that classifier to ensure reasonable comparisons between methods.

[0569] To evaluate the performance of deep learning networks using de novo variants in patients with neurodevelopmental disorders and healthy controls in a clinical setting, we downloaded de novo variants from the Demystifying Developmental Disorders (DDD) study and healthy sibling controls from the Simons Simplex Collection (SSC) autism study. The DDD study provides confidence scores for de novo variants, and we excluded variants with a threshold < 0.1 from the DDD dataset as potential false positives due to variant misidentification. We had a total of 3,512 missense de novo variants from DDD-infected individuals and 1,208 missense de novo variants from healthy controls.

[0570] To better simulate the real-world clinical situation of distinguishing benign from pathogenic variants of unknown significance in candidate disease gene sets, we restricted our analysis to de novo variants in the 605 genes associated with the disease in the DDD study (p < 0.05), calculated solely based on protein truncating variants (Supplementary Table 18). We assessed gene-specific enrichment of protein-truncating de novo variants by calculating statistical significance under the assumption-free condition of the expected number of de novo mutations given the gene-specific mutation rate and the number of chromosomes considered. We selected the 605 genes with a nominal P value < 0.05. We calculated the 605 genes ( Fig.22A) as the ratio of observed to expected de novo mutation counts and the difference between observed de novo mutations and expected de novo mutations. Of these 605 genes, we observed 380 de novo missense mutations in DDD-infected individuals ( Fig.22A ). For each classifier, including our own, a small fraction of variants were not predicted, typically because they did not map to the same transcript model used by the classifier. Therefore, for our deep learning network, we performed downstream analyses using 362 de novo missense mutations from DDD-infected individuals and 65 de novo missense mutations from healthy controls, as Fig.22A , as shown in 22B, 22C, 22D and 22E.

[0571] Saturation of all possible human missense mutations as the number of primate populations sequenced increases

[0572] We investigated saturation of all ~70M human missense variants by common variants in 504 primate species. For each primate, we simulated four times the number of common missense variants observed in humans (~83,500 missense variants with >0.1% allele frequency), because humans have about half as many variants per individual as other primates and because approximately 50% of human missense variants have been filtered out by purifying selection at >0.1% allele frequency ( Fig.49A ). We assigned simulated variants based on the observed distribution of common human missense variants across 96 trinucleotide contexts. For example, if 2% of common human missense variants are from CCG>CTG trinucleotide contexts, then we would require that 2% of simulated variants be randomly sampled CCG>CTG mutations. Using trinucleotide contexts controls for the effects of mutation rate, genetic drift, and gene conversion bias.

[0573] Assume that we identify all common variants (>0.1% allele frequency) in each primate species, Fig.23D The curve in Figure 2 shows the cumulative saturation of the approximately 70M possible human missense mutations caused by common variants in 504 primate species. Fig.49A As can be seen in Figure 2, approximately 50% of human missense mutations are sufficiently deleterious in humans and other primates to prevent them from rising to common allele frequency (>0.1%). Fig.23D The curve in represents the proportion of non-deleterious human missense mutations that become saturated with common primate variants as the number of primate species increases. We show that most non-deleterious human missense mutations will be saturated across the 504 primate species, and that non-deleterious CpG mutations will be saturated across the less diverse species due to their higher mutation rates.

[0574] To simulate the proportion of common missense variants (>0.1% allele frequency) found in humans as the number of human cohorts investigated increases ( Fig.36 ), we sampled genotypes based on gnomAD allele frequencies. For sample sizes ranging from 100 to 100K, the proportion of gnomAD common missense variants found in the 100 simulated samples was averaged.

[0575] Prediction of secondary structure and solvent accessibility

[0576] The deep learning network for pathogenicity prediction contains 36 total convolutional layers, including 19 convolutional layers for the secondary structure and solvent accessibility prediction network, and 17 convolutional layers for the main pathogenicity prediction network, which takes the results of the secondary structure and solvent accessibility network as input. Because the crystal structures of most human proteins are unknown, we trained two models to enable the network to learn protein structure from primary sequence. Figure 6 Both models use the same network architecture and inputs as shown. The input to the secondary structure and solvent accessibility networks is a 51-length × 20 amino acid position frequency matrix encoding conservation information from a multiple sequence alignment of human with 99 other vertebrates.

[0577] The secondary structure network was trained to predict three states of secondary structure: alpha helix (H), beta sheet (B), and coil (C). The solvent accessibility network was trained to predict three states of solvent accessibility: hidden (B), intermediate (I), and exposed (E). Both networks took only primary sequence as input and were trained using labels of known crystal structures from the Protein DataBank. The model predicts one state for each amino acid residue.

[0578] Data preparation for prediction of secondary structure and solvent accessibility

[0579] We used unrelated crystal structures from the Protein Data Bank to train the model. Amino acid sequences with more than 25% sequence similarity were removed. A total of 6,367 protein sequences were used for training, 400 protein sequences were used for validation, and 500 protein sequences were used for testing (Supplementary Table 13). The RaptorX dataset is available from the RaptorX website ( http: / / raptorx.uchicago.edu / download / ) to obtain the data used for training, which includes amino acid sequence, secondary structure and solvent accessibility labels.

[0580] Most of the solved crystal structures are non-human proteins, so in order to pre-train secondary structure and solvent models, we used the RaptorX suite (based on PSI-BLAST) to obtain relevant sequences, because human-based multiple sequence alignments are generally not available. We used the CNFsearch1.66_release tool of RaptorX to generate a multiple sequence alignment for the protein and calculated the amino acid at each position from the 99 closest alignments to form a position frequency matrix. For example, the specific command to retrieve the multiple sequence alignment of the 1u7lA.fasta protein using RaptorX is as follows: %. / buildFeature -i 1u7lA.fasta -c 10 -o . / TGT / 1u7lA.tgt %. / CNFsearch -a 30 -q 1u7lA

[0581] For each amino acid position in the dataset, we extracted a window from the position frequency matrix corresponding to the 51 amino acids flanking it and used it to predict a secondary structure or solvent accessibility signature for the amino acid in the center of the 51-amino acid sequence. The secondary structure and relative solvent accessibility signatures were obtained directly from the known 3D crystal structure of the protein using DSSP software and did not require prediction from the primary sequence. To include the secondary structure and solvent accessibility networks as part of the pathogenicity prediction network, we computed the position frequency matrix from a human-based 99 vertebrate multiple sequence alignment. Although the conservation matrices generated from the two methods were generally similar, we back-propagated them through the secondary structure model and the solvent accessibility model during the training process for pathogenicity prediction, allowing for fine-tuning of parameter weights.

[0582] Model architecture and training

[0583] We trained two separate deep convolutional neural network models to predict secondary structure and relative solvent accessibility of proteins. Both models had the same architecture and input data, but different prediction states. We performed a detailed hyperparameter search to optimize the models for best performance. Our deep learning network for pathogenicity prediction and our deep learning network for predicting secondary structure and solvent accessibility both adopted the residual block architecture that has been widely adopted for image classification. The residual block consists of repeated units of convolutions interspersed with skip connections that allow information from earlier layers to skip the residual block. In each residual block, the input layer is first batch normalized and then activated with a rectified linear unit (ReLU). After activation, it passes through a 1D convolutional layer. The intermediate layer output from the 1D convolutional layer is again batch normalized and ReLU activated before another 1D convolutional layer is performed. At the end of the second 1D convolution, we add its output with the original input to the residual block, which acts as a skip connection by allowing the original input information to bypass the residual block. In such an architecture, which its authors call a deep residual learning network, the input remains in its original state and the residual connections are not affected by the nonlinear activations of the model, allowing for efficient training of deeper networks. Figure 6 and Supplementary Table 11 ( Fig. 7A and 7B ) and 12( Fig. 8A and 8B ) provides a detailed architecture.

[0584] After the residual block, the softmax layer calculates the probability of three states for each amino acid, where the largest softmax probability determines the state of the amino acid. The model is trained on the entire protein sequence using the ADAM optimizer with a cumulative categorical cross entropy loss function. After pre-training the network on secondary structure and solvent accessibility, we no longer use the output of the network directly as the input of the pathogenicity prediction network, but put it before the softmax layer so that more information can pass through the pathogenicity prediction network.

[0585] The best test accuracy of the three-state secondary structure prediction model was 79.86% (Supplementary Table 14), which is similar to the state-of-the-art accuracy predicted by the DeepCNF model30. The best test accuracy of the three-state solvent accessibility prediction model was 60.31% (Supplementary Table 14), which is similar to the current best accuracy predicted by RaptorX on a similar training dataset. We also compared the predictions of the neural network using DSSP-annotated structural labels of approximately 4000 human proteins with crystal structures with the standard PrimateAI model using only predicted structural labels. We did not improve the accuracy of pathogenicity predictions when using labels with DSSP annotations (Supplementary Table 15).

[0586] Input features of deep learning models for pathogenicity prediction

[0587] After filtering, the training dataset for the pathogenicity prediction network contained 385,236 variants labeled as benign and 68,258,623 unlabeled variants. For each variant, we generated the following input features. The first input feature of each variant is its 51-length flanking amino acid sequence, i.e., 25 amino acids on each side of the variant obtained from the reference sequence of hg19, to provide sequence context for the deep learning model of the variant. The length of this flanking reference sequence is 51 amino acids in total. Through empirical observation, we found that using amino acids to represent protein sequences is more effective than using nucleotides to represent protein coding sequences.

[0588] The second feature is a human flanking amino acid sequence of length 51 with the alternative amino acid replaced by the variant at the central position. The alternative flanking sequence is identical to the reference flanking sequence in the first feature, except that the central position of the sequence contains the alternative amino acid instead of the reference amino acid. Both the reference and alternative human amino acid sequences are converted to one-bit efficient encoding vectors of length 51×20, where each amino acid is represented by a vector of 19 amino acids (value 0) and a single amino acid (value 1).

[0589] Three position frequency matrices (PFMs) were generated from multiple sequence alignments of variants from 99 vertebrates, one for 11 primates, one for 50 mammals (excluding primates), and one for 38 vertebrates (excluding primates and mammals). Each PFM was of size Lx20, where L is the length of the flanking sequence around the variant (in our case, L represents 51 amino acids).

[0590] For the input of the pre-trained three-state secondary structure and three-state solvent accessibility networks, we used a single PFM matrix generated from a multiple sequence alignment of 99 vertebrates, also of length 51 and depth 20. After pre-training the network on known crystal structures from ProteinDataBank, the last two layers of the secondary structure and solvent models (the entire maxpool layer and the output layer) were deleted, and the 51×40-shaped output of the previous layer was used as the input of the pathogenicity prediction network. We allowed backpropagation through the structured layers of the network to fine-tune the parameters.

[0591] Semi-supervised learning

[0592] Because semi-supervised learning algorithms use both labeled and unlabeled examples during training, they can produce classifiers that perform better than fully supervised learning algorithms that have only a small amount of labeled data available for training. The principle of semi-supervised learning is that the inherent knowledge in unlabeled data can be used to enhance the predictive power of supervised models that use only labeled examples, providing a potential advantage for semi-supervised learning. Supervised classifiers whose model parameters are learned from a small amount of labeled data can be shifted to a more realistic distribution (more similar to the distribution of the test data) using unlabeled data.

[0593] Another challenge that is prevalent in bioinformatics is the problem of imbalanced data. Imbalanced data occurs when the class to be predicted is underrepresented in the data, either because instances belonging to that class are rare (notable cases) or difficult to obtain. Ironically, those few classes are often the most important to learn from, as they may be associated with special cases.

[0594] One algorithm for dealing with imbalanced data distributions is one based on classifier ensembles. A limited amount of labeled data naturally leads to weaker classifiers, but an ensemble of weaker classifiers will often exceed the performance of any individual component classifier. Furthermore, an ensemble of classifiers will often improve the predictive accuracy obtained from a single classifier by a factor that justifies the effort and cost associated with learning multiple models. Intuitively, aggregating several classifiers should lead to better overfitting control because when the high variability of the individual classifiers is averaged out, the overfitting of the classifiers is also averaged out.

[0595] Due to the lack of datasets of sufficient size to reliably label pathogenic variants, we employed a semisupervised learning strategy. Although the ClinVar database has more than 300,000 entries, after removing variants with significant uncertainty, only about 42,000 missense variants remained with no conflicting pathogenic interpretations.

[0596] The systematic review also found that these entries often do not have sufficient clinical evidence to support their annotated pathogenicity. Moreover, most variants in human curated databases tend to be located in very few genes, making them mismatched with variants in benign training datasets, which are determined genome-wide using human common variants or chimpanzee-human fixed surrogates. Given the differences in the determined datasets, training supervised learning models with human curated variants as the pathogenic set and genome-wide common variants as the benign set could introduce significant bias.

[0597] We trained a deep learning network to distinguish between a set of labeled benign and unlabeled variants that were carefully matched to eliminate bias. According to one implementation, the 385,236 labeled benign variants included common human variants (>0.1% allele frequency) in the ExAC / gnomAD database and variants from six non-human primates.

[0598] We sampled a set of unlabeled variants that needed to be matched to benign variants on a trinucleotide background (to control for mutation rate, genetic drift, and gene conversion) and adjusted for the effects of comparability and sequence coverage on variant confirmation. Since the number of unlabeled variants greatly exceeded the number of labeled benign variants, we averaged the training set of eight models using the same labeled benign variants and eight models using randomly sampled unlabeled variants to obtain consistent predictions.

[0599] The motivation for choosing semi-supervised learning is that human curated variant databases are unreliable and contaminated, especially lacking reliable pathogenic variants. We obtained a set of reliable benign variants from gnomAD and primate variants. For pathogenic variants, we adopted an iterative balanced sampling method to sample pathogenic variants from a set of unknown variants (VUS variants, no annotated clinical significance).

[0600] To reduce sampling bias, we trained an ensemble of eight models that used the same set of benign training variants and eight different sets of pathogenic variants. Initially, we randomly sampled unknown variants to represent pathogenic variants. Next, the model ensemble was used to iteratively score a set of unknown variants that had not participated in the previous training cycle. The highest-scoring pathogenic variants were then obtained to replace 5% of the random unknown variants from the previous cycle. Note that we retained 25% more of the highest-scoring pathogenic variants than needed to be retained, so we could sample eight sets of pathogenic variants with different scores to replace the unknown variants, which increased the randomness of the eight models. A new pathogenic training set was then formed and a new training cycle was performed. This process was repeated until all of the initially randomly sampled unknown variants were replaced by high-confidence pathogenic variants predicted by the ensemble model. Fig.42 The iterative equalization sampling process is illustrated.

[0601] Balancing benign and unknown training sets

[0602] The sampling scheme of unknown variants matched with benign variants reduces the bias in our model training. When unknown variants are randomly sampled, deep learning models often extract biased information and give simple solutions. For example, if the amino acid substitution K->M occurs more frequently in unknown variants than in benign variants, deep learning models tend to classify the K->M substitution as pathogenic. Therefore, it is important to balance the distribution of amino acid substitutions between the two training sets.

[0603] Higher mutation types such as CpG transitions have a huge representation bias in common benign variants. Orthologous variants from other primates also follow human mutation rates, which means enrichment of highly expressed mutation types in the overall benign training set. If the sampling process of unknown variants is not well controlled and balanced, deep learning models will tend to classify CpG transitions as benign compared to less represented classes such as translocations or non-CpG transitions.

[0604] To prevent the deep learning model from converging to a trivial non-biological solution, we consider a trinucleotide background that balances benign variants and unknown variants. The trinucleotide is formed by the base before the variant, the reference base of the variant, and the base after the variant. And the reference base of the variant can be changed to other three nucleotides. There are 64x3 trinucleotide backgrounds in total.

[0605] Iterative Equalization Sampling

[0606] Cycle 1

[0607] We sampled the unknown variants to match the exact number of benign variants per trinucleotide context. In other words, in the first cycle, we reflected the benign and pathogenic training sets from the three-nucleotide context of the variants. The fact behind this sampling approach is that variants with the same mutation rate between the benign and unknown groups are represented in the same way, which prevents the model from converging to a trivial solution based on mutation rate.

[0608] Cycle 2-20

[0609] For cycle 2, we applied the cycle 1 trained model to score a set of unknown variants not implicated in cycle 1 and replaced 5% of the unknown variants with the most predicted pathogenic variants. This set of variants was entirely generated by the model, so we did not perform any balancing for trinucleotide background in this set. The remaining 95% of unknown variants required for training were sampled as 95% of the background counts for each trinucleotide in the benign variants.

[0610] Objectively, since cycle 1 uses a fully matched training set, the most predicted pathogenic variants generated do not have any mutation rate bias. Therefore, there is no need to consider any bias in this group. The remaining 95% of the data is still controlled by the trinucleotide background mutation rate to prevent the model from converging to a trivial solution.

[0611] For each cycle, the percentage of unknown variants replaced increased by 5%. For cycle 3, we replaced 5% of unknown variants with the variant with the highest pathogenicity predicted in the cycle 3 model. Cumulatively, the proportion of pathogenic variants increased to 10% and unknown variants reflected in the trinucleotide background decreased to 90%. The sampling process was similar for the remaining cycles.

[0612] Cycle 21

[0613] For cycle 21 (the last cycle), the entire pathogenicity training set consists entirely of the variants predicted by the deep learning model to be the most pathogenic. Since we have explicitly controlled the mutation rate bias at each cycle, the pathogenic variants can be used as training data, which is very reliable and not affected by the mutation rate bias. Therefore, the last cycle of training produces the final deep learning model for pathogenicity prediction.

[0614] Matching labeled benign training sets and unlabeled training sets

[0615] Balanced sampling of unlabeled variants is critical to eliminate biases that are unrelated to the deleteriousness of the variants. When confounding effects are not properly controlled, deep learning can easily inadvertently introduce biases that affect the distinction between classes. Variants common to humans tend to be enriched for variants from high-mutation types, such as those on CpG promoter regions. Similarly, primate polymorphisms follow the mutation rate of humans, meaning that there is an enrichment of variants from high-mutation types throughout the benign training set. If the sampling process of unlabeled variants is not well controlled and balanced, deep learning networks tend to rely on mutation rate biases to classify variants, and therefore they are more likely to classify CpG transitions as benign than less represented types such as translocations or non-CpG transitions. We sampled exactly the same number of unlabeled variants as labeled benign variants in each of the 96 trinucleotide contexts (discussed earlier).

[0616] When matching unlabeled variants to primate variants in the labeled benign dataset, we did not allow selection of variants from human genomic regions that had not undergone multiple sequence alignment across primate species, since no variant in that primate species could be identified at that locus.

[0617] In the 96 trinucleotide contexts, we corrected for primate variants by sequencing coverage. Variants that are common in the human population are frequently observed due to the large number of sequenced human populations, such that they are well defined even in regions with low sequencing coverage. This is not the case for primate variants, as only a small number of individuals have been sequenced. We divided the genome into 10 regions based on sequencing coverage of the ExAC / gnomAD exomes. For each region, we measured the proportion of primate variants in the labeled benign dataset relative to the unlabeled dataset. We used linear regression to calculate the probability of a primate variant being in the labeled benign dataset based solely on sequencing coverage ( Fig.24 ). When selecting unlabeled variants to match primate variants in the labeled benign dataset, we weighted the probability of sampling the variant according to the sequencing coverage at that site using a regression coefficient.

[0618] Generation of benign and unknown variants

[0619] Common variants in the human population

[0620] Recent studies have shown that common variants in human populations are generally benign. According to one embodiment, gnomAD provides 90,958 non-synonymous SNPs with a minor allele frequency (MAF) >= 0.1% within canonical coding regions. Those variants that pass the filter are retained, excluding indels. Variants occurring in start codons or stop codons and variants with protein truncation are deleted. According to one embodiment, by examining subpopulations, the total number of missense variants with MAF>= 0.1% in each subpopulation is increased to 245,360. These variants constitute part of the benign variant training set.

[0621] Common polymorphisms in great apes

[0622] Since coding regions are known to be highly conserved, it is straightforward to hypothesize whether a polymorphism segregates at high frequency in ape populations, which may also have a minor impact on the human fit. Polymorphism data for bonobos, chimpanzees, gorillas, and orangutans from the Great Ape Genome Project and other studies were merged with SNPs from rhesus macaques and marmosets in dbSNP.

[0623] Generation of unknown variants

[0624] All possible variants can be generated from each base position of the canonical coding region by replacing the nucleotide at each base position with three other nucleotides. When a new codon is formed, it causes a potential change in the amino acid at that position. Synonymous mutations will be filtered out.

[0625] Variants observed in the gnomAD dataset were removed. Variants occurring in the start codon or stop codon and variants forming the stop codon were removed. For SNPs with multiple gene annotations, the canonical gene annotation was selected to represent the annotation of the SNP. According to one embodiment, a total of 68,258,623 unknown variants were generated.

[0626] Filtering variants

[0627] Some regions of the human genome are difficult to map and interpret. Including these regions can have confounding effects on the training and test datasets. For example, regions under high selective pressure tend to have a limited number of polymorphisms. However, regions that are difficult to sequence also have fewer polymorphisms. To avoid confounding our models, we removed variants from genes that were not sequenced by gnomAD.

[0628] Typically, benign variants are found in well-sequenced regions that tend to be conserved across multiple species. Unknown variants are randomly sampled across the genome, including some regions with poor coverage. This leads to a bias in the certainty between the benign and unknown sets. To reduce the bias, we filtered variants in gnomAD with an interpretation depth of less than 10. We also filtered all variants with more than 10% missing data in the flanking sequence alignments across all mammalian species.

[0629] Validation and test data

[0630] To validate and test the pathogenicity models, according to one embodiment, we randomly sampled two sets of 10,000 benign variants each for validation and testing from a large set of benign variants. The remainder of the benign variants were used to train the deep learning model. These variants were specifically sampled from orthologous primate variants to ensure fair comparisons between methods, as those methods are trained on common human variants. According to one embodiment, we also randomly sampled two sets of 10,000 unknown variants each for validation and testing. We ensured that the number of unknown variants in the 192 trinucleotide background matched the number of benign variants used for validation and testing, respectively.

[0631] We evaluated the performance of multiple methods in a clinical setting using de novo variants from children with autism or degenerative developmental disorder (DDD) and their unaffected siblings. According to one embodiment, there were a total of 3821 missense de novo variants from DDD cases and 2736 missense de novo variants from autism cases. According to one embodiment, there were 1231 unaffected siblings with missense de novo variants.

[0632] Deep Learning Network Architecture

[0633] The pathogenicity prediction network receives five direct inputs and two indirect inputs through the secondary structure and solvent accessibility networks. The five direct inputs are amino acid sequences of length 51 × depth 20 (encoding 20 different amino acids), including the human reference amino acid sequence without variants (1a), the human alternative amino acid sequence with variants substituted in (1b), the PFM of the multiple sequence alignment of primate species (1c), the PFM of the multiple sequence alignment of mammalian species (1d), and the PFM of the multiple sequence alignment of more ancient vertebrate species (1e). The secondary structure and solvent accessibility networks receive PFMs from the multiple sequence alignment (1f) and (1g), respectively, as input and output them as input to the main pathogenicity prediction network as its input. The secondary structure and solvent accessibility networks have been pre-trained on known protein crystal structures from Protein DataBank and allowed to back-propagate during pathogenicity model training.

[0634] The five direct input channels are passed through an upsampling convolutional layer of 40 kernels with linear activation. The human reference amino acid sequence (1a) is merged with the PFM from the multiple sequence alignment of primates, mammals, and vertebrates (Merge 1a). Similarly, the human alternative amino acid sequence (1b) is merged with the PFM from the multiple sequence alignment of primates, mammals, and vertebrates (Merge 1b). This creates two parallel tracks, one for the reference sequence and the other with the alternative sequence, in which the variant has been replaced.

[0635] The merged feature maps of the reference and alternate channels (Merge 1a and Merge 1b) are passed through a series of six residual blocks (layers 2a to 7a, Merge 2a and layers 2b to 7b, Merge 2b). The outputs of the residual blocks (Merge 2a and Merge 2b) are concatenated together to form a feature map of size (51, 80) (Merge 3a, Merge 3b) that fully mixes the data in the reference and alternate channels. Next, the data has two parallel paths through the network, either through a series of six residual blocks, where each residual block contains two convolutional layers as defined in Section 2.1 (Merge 3 to 9, layers 9 to 46, excluding layers 21, 34), or through skip connections, which connect the outputs of every two residual blocks after 1D convolution (layers 21, 37, 47). Finally, the merged activations (Merge 10) are fed to another residual block (layers 48 to 53, Merge 11). The activations from pooling 11 are assigned to a 1D convolution (layer 54) with filter size 1 and sigmoid activation, and then passed through a global max pooling layer, which selects a value that represents the network's prediction of the pathogenicity of the variant. A schematic diagram of the model is shown in Figure 3 and Supplementary Table 16 ( Figure 4A, 4B and 4C).

[0636] Model Overview

[0637] We developed a semi-supervised deep convolutional neural network (CNN) model to predict the pathogenicity of variants. The input features of the model include protein sequences and conservation maps flanking the variants, as well as the loss of missense variants in specific gene regions. We also predicted the changes in secondary structure and solvent accessibility caused by the variants through a deep learning model and integrated them into our pathogenicity prediction model. To train the model, we generated benign variants from common variants in human subpopulations and orthologous variants from primates. However, we still lack a reliable source of pathogenic variants. We initially trained the model with benign and unknown variants, and then used a semi-supervised iterative balanced sampling (IBS) algorithm to gradually replace the unknown variants with a set of pathogenic variants with high confidence. Finally, we demonstrated that our model outperforms existing methods in distinguishing between de novo variants that cause human developmental disorders and benign diseases.

[0638] The use of residual blocks

[0639] Fig.17 The residual block is shown. Both our deep learning model for pathogenicity prediction and the deep learning model for predicting secondary structure and solvent accessibility adopt the definition of the residual block, which was first described in this paper. The structure of the residual block is shown in the figure below. The input layer is first batch normalized and then normalized with the non-linear activation "ReLU". After the activation, it passes through a 1D convolution layer. The intermediate layer output from the 1D convolution layer is batch normalized and ReLU activated again, and then another 1D convolution layer is performed. At the end of the second 1D convolution, we merge its output with the original input. In such an architecture, the input remains in its original state, and the residual connection is not affected by the non-linear activation of the model.

[0640] Atrous convolutions allow for large receptive fields with few trainable parameters. Atrous convolutions are a type of convolution in which the kernel is applied to an area larger than its length by skipping input values ​​by a certain step size (also called the atrous convolution rate or dilation factor). Atrous convolutions increase the spacing between convolution filter / kernel elements so that larger spacings of neighboring input entries (e.g., nucleotides, amino acids) are considered when performing the convolution operation. This allows for inclusion of long-range contextual correlations in the input. Atrous convolutions retain part of the convolution computation for reuse when processing neighboring nucleotides.

[0641] Novelty of our model

[0642] Our approach differs from existing methods for predicting variant pathogenicity in three aspects. First, our method employs a novel semi-supervised deep convolutional neural network architecture. Second, reliable benign variants are obtained from human common variants from gnomAD and primate variants, and the high-confidence pathogenicity training set is generated by iterative balanced sampling and training to avoid cyclic training and testing of the model using the same human-curated variant database. Third, deep learning models of secondary structure and solvent accessibility are integrated into the architecture of our pathogenicity model. The information obtained from the structure and solvent models is not limited to label predictions of specific amino acid residues. Moreover, the readout layer is removed from the structure and solvent models, and the pre-trained model is merged with the pathogenicity model. When training the pathogenicity model, the structure and solvent pre-trained layers are also back-propagated to minimize the error. This allows the pre-trained structure and solvent models to be dedicated to the pathogenicity prediction problem.

[0643] Training secondary structure and solvent accessibility models

[0644] Data preparation

[0645] We trained a deep convolutional neural network to predict the three-state secondary structure and three-state solvent accessibility of proteins. Protein annotations from PDB were used to train the model. According to one embodiment, sequences with more than 25% similarity to the sequence profile were removed. In total, according to one embodiment, 6,293 protein sequences were used for training, 392 protein sequences were used for validation, and 499 protein sequences were tested.

[0646] Position-specific scoring matrices (PSSM) conservation profiles of proteins were generated by searching UniRef90 by running PSI-BLASTE with an E-value threshold of 0.001 and 3 iterations. Any unknown amino acid and its secondary structure were set as blank. We also performed PSI-BLAST with similar parameter settings for all human genes to collect their PSSM conservation profiles. These matrices were used to integrate the structural models into pathogenicity prediction. The amino acids of the protein sequence were then converted into one-bit efficient encoding vectors. The protein sequence and PSSM matrix were reshaped into an Lx20 matrix, where L is the length of the protein. The three predicted labels for secondary structure include helix (H), beta sheet (B), and coil (C). The three labels for solvent accessibility include buried (B), intermediate (I), and exposed (E). One label corresponds to one amino acid residue. The labels are encoded as one-bit efficient encoding vectors of dimension 3.

[0647] Model architecture and training

[0648] We trained two end-to-end deep convolutional neural network models to predict the three-state secondary structure and three-state solvent accessibility of proteins, respectively. The two models have similar configurations, including two input channels, one for protein sequence and the other for protein conservation profile. The dimension of each input channel is Lx20, where L represents the length of the protein.

[0649] Each input channel passes through a 1D convolutional layer (Layer 1a and 1b) with 40 kernels and linear activations. This layer is used to upsample the input dimension from 20 to 40. Note that all other layers throughout the model use 40 kernels. The activations of the two layers (1a and 1b) are merged together by summing each of the 40 dimensions (i.e., merge mode = "sum"). The output of the merge node passes through a layer of 1D convolution (Layer 2) followed by a linear activation.

[0650] The activations of layer 2 pass through a series of 9 residual blocks (layers 3 to 11) as defined above. The activations of layer 3 are fed back to layer 4, the activations of layer 4 are fed back to layer 5, and so on. There are also skip connections that directly add the outputs of every 3rd residual block (layers 5, 8, and 11). The combined activations are then fed back into 2 1D convolutions with ReLU activations (layers 12 and 13). The activations of layer 13 are given to a softmax readout layer. This softmax computes the probabilities of the three class outputs given an input.

[0651] For the best secondary structure model, the 1D convolution was used with a dilation ratio of 1. For the solvent accessibility model, the last 3 residual blocks (layers 9, 10, and 11) were used with a dilation ratio of 2 to increase the kernel coverage. The secondary structure of proteins strongly depends on the interactions between immediately adjacent amino acids. Therefore, models with higher kernel coverage slightly improve the performance. On the other hand, solvent accessibility is affected by large-scale interactions between amino acids. As a result, models with high kernel coverage using dilated convolutions are more than 2% more accurate than models with short coverage.

[0652] The following table provides detailed information about the activation values ​​and parameters of each layer of the three-state secondary structure prediction model according to one embodiment:

[0653]

[0654]

[0655] According to one embodiment, the following table shows the details of the solvent accessibility

[0656]

[0657]

[0658]

[0659] The secondary structure category of a specific amino acid residue is determined by the maximum predicted softmax probability. The ADAM optimizer is used to train the model of the entire protein sequence with the cumulative categorical cross entropy loss function to optimize the back propagation algorithm.

[0660] The best test accuracy of the three-state secondary structure prediction model was 80.32%, which was similar to the best accuracy predicted by the DeepCNF model on a similar training dataset.

[0661] The best test accuracy of the three-state solvent accessibility prediction model is 64.83%, which is similar to the current best accuracy predicted by RaptorX on a similar training dataset.

[0662] We integrated the pre-trained three-state secondary structure and solvent accessibility prediction models into our pathogenicity prediction model as described below.

[0663] Training models to predict pathogenicity of variants

[0664] Input features of pathogenicity prediction models

[0665] As mentioned above, for the pathogenicity prediction problem, there is a training set of benign variants and a training set of unknown variants to train the pathogenicity model. For each variant, we prepared the following input features to feed into our model.

[0666] The first input feature of each variant is its flanking amino acid sequence, i.e., 25 amino acids on each side of the variant obtained from the reference sequence of hg19, to provide a deep learning model of the sequence context of the variant. In total, the length of the flanking reference sequence is 51 amino acids.

[0667] The second feature is the alternative amino acid that forms the variant. Instead of providing a reference alternative amino acid pair directly, we provide an alternative flanking sequence for the model. The alternative flanking sequence is identical to the reference flanking sequence in the first feature, except that the middle position of the sequence contains an alternative amino acid instead of the reference amino acid.

[0668] Both sequences were then converted into one-bit efficient coded vectors of length 51x20, where each amino acid is represented by a vector of 20 0s or 1s.

[0669] Three position weight matrices (PWMs) were then generated from the multiple sequence alignment (MSA) of the 99 vertebrate variants, one PWM for the 12 primates, one PWM for the 47 mammals (excluding primates), and one PWM for the 40 vertebrates (excluding primates and mammals). Each PWM has a size of L x 20, where L is the length of the flanking sequence around the variant (in this case, L represents 51 amino acids). It contains the counts of amino acids seen in each class of species.

[0670] We also generated a PSSM matrix of the 51 amino acid variant-flanking sequences from Psi blast. This was used to integrate the three-state secondary structure and solvent accessibility prediction models for pathogenicity prediction.

[0671] We trained pathogenicity and solvent accessibility models using reference sequences (input 1), alternate sequences (input 2), PWM matrices for primates (input 3), mammals (input 4), vertebrates (input 5), and three-state secondary structure information.

[0672] Deep learning model training

[0673] Fig.19 A block diagram is provided to provide an overview of the deep learning model workflow. The pathogenicity training model includes five direct inputs and four indirect inputs. The five direct input features include the reference sequence (1a), alternative sequence (1b), primate conserved sequence (1c), mammalian conserved sequence (1d), and vertebrate conserved sequence (1e). The indirect inputs include secondary structure based on the reference sequence (1f), secondary structure based on the alternative sequence (1g), solvent accessibility based on the reference sequence (1h), and solvent accessibility based on the alternative sequence (1i).

[0674] For indirect inputs 1f and 1g, we loaded the pre-trained layers of the secondary structure prediction model, excluding the softmax layer. For input 1f, the pre-trained layers are based on the human reference sequence of the variant and the PSSM of the variant generated by PSI-BLAST. Similarly, for input 1g, the pre-trained layers of the secondary structure prediction model take the human alternative sequence as well as the PSSM matrix as input. Inputs 1h and 1i correspond to similar pre-trained channels containing solvent accessibility information for the reference and alternative sequences of the variant, respectively.

[0675] The five direct input channels are passed through an upsampling convolutional layer with 40 kernels with linear activation values. Layers 1a, 1c, and 1h are merged with the values ​​of the 40 feature dimensions added together to produce layer 2a. In other words, the feature map of the reference sequence is merged with the three types of conservative feature maps. Similarly, layers 1b, 1c, 1d, and 1e are merged with the values ​​of the 40 feature dimensions added together to produce layer 2b, that is, the features of the alternative sequence are merged with the three types of conservative features.

[0676] Layers 2a and 2b are batch normalized using ReLU activation values, and each layer passes through a 1D convolutional layer with a filter size of 40 (3a and 3b). The outputs of layers 3a and 3b are merged with 1f, 1g, 1h, and 1i, where the feature maps are associated with each other. In other words, the feature maps of the reference sequence with conserved profiles and the alternative sequences with conserved profiles are merged with the secondary structure feature maps of the reference and alternative sequences and the solvent accessibility feature maps of the reference and alternative sequences (layer 4).

[0677] The output of layer 4 passes through six residual blocks (layers 5, 6, 7, 8, 9, 10). The last three residual blocks of the 1D convolution have a dilation rate of 2 to provide higher coverage for the kernel. The output of layer 10 passes through a 1D convolution with a filter size of 1 and a sigmoid activation (layer 11). The output of layer 11 passes through a global maxpool that selects a single value for the variant. This value represents the pathogenicity of the variant. Details of one embodiment of the pathogenicity prediction model are shown in the following table.

[0678]

[0679]

[0680]

[0681]

[0682] gather

[0683] In one embodiment, for each cycle of our method, we run eight different models trained on the same benign dataset and eight different unknown datasets, and average the predictions of the evaluation dataset across the eight models. When there are multiple randomly sampled sets of unknown variants in the model, sampling bias can be reduced and controlled.

[0684] Additionally, employing an ensemble approach improves the performance of the models on our evaluation datasets. CADD uses an ensemble of 10 models and takes the average score of all 10 models to score the variants. Here we attempt to use a similar ensemble approach. We benchmark the results using one ensemble and then increase the number of ensembles to evaluate the performance gain. Note that each ensemble has 8 models, which are trained on the same benign dataset and 8 different unknown datasets. The seed value of the random number generator is different for different ensembles, so that the sets of random variants are drawn differently from each other.

[0685] The following table shows detailed results according to one embodiment.

[0686]

[0687] Compared to one ensemble, 5 and 10 ensembles produced more significant p-values ​​when evaluated using the DDD dataset. However, increasing the number of ensembles did not further improve performance, indicating saturation of the ensembles. Ensembles reduce sampling bias with a large number of unknown variants. However, we also needed to match 192 trinucleotide contexts between benign and pathogenic classes, which significantly limited our sampling space and led to rapid saturation. We conclude that the ensemble approach significantly improved the performance of the model and further enriched our understanding of the model.

[0688] Early stopping for training pathogenicity models

[0689] Defining stopping criteria for model training is challenging due to the lack of reliable annotated pathogenic variant samples. To avoid using pathogenic variants in model evaluation, in one embodiment, we used 10,000 benign validation variants from orthologous primates and 10,000 trinucleotide contexts matching unknown variants. After each stage of training the model, we evaluated the benign validation variants and unknown validation variants. We used the Wilcoxon rank sum test to evaluate the difference in the probability distributions of the two validation variant sets.

[0690] As the model's ability to distinguish between benign variants and a set of position-unknown variants improves, the p-values ​​of the test become more significant. We stop training if the model's...

Claims

1. A computer-implemented method for identifying pathogenic single nucleotide polymorphisms (SNPs), the method comprising: Using an ensemble of convolutional neural networks (CNNs) to process unlabeled SNPs in a cyclic manner, wherein the cycle includes a first stage and a second stage; The unlabeled SNPs include SNP variants that were generated in combination, not labeled as benign, and have not been added to the elite predicted pathogenic SNPs; The processing in the first phase of the cycle includes: Selecting separate pseudo-pathogenic training examples from the unlabeled SNPs for each CNN in the set to be used in conjunction with the elite predicted pathogenic SNPs, if any; and Using the individual pseudo-pathogenic training examples and elite predicted pathogenicity SNPs, if any, a set of CNNs trained with balanced general SNP benign training examples, thereby updating the parameters of the CNNs in the set; Processing in the second phase of the cycle includes: Applying an ensemble of CNNs with updated parameters to generate pathogenicity scores for unlabeled SNPs not involved in the first stage; Use pathogenicity scores to identify outlier SNPs; and Accumulation of abnormal SNPs among elite predicted pathogenic SNPs; Executing the first and second phases of processing in the first cycle includes: Starting with no elite predicted pathogenic SNPs; and End with the outlier SNPs identified and accumulated in the second-stage processing as elite predicted pathogenic SNPs; The first and second phases of the treatment are repeated in the second and subsequent cycles, including: Starting with elite predicted pathogenic SNPs identified as abnormal in a previous cycle of Phase II processing; and End with the accumulation of additionally identified outlier SNPs with the elite predicted pathogenic SNPs; and These cycles are repeated until a final cycle, after which the number of elite predicted pathogenic SNPs accumulated in the second stage of processing is within a predetermined range of the number of SNPs in the generic benign training examples.

2. The computer-implemented method of claim 1 , wherein: Each of the CNN sets is a deep convolutional neural network, abbreviated as DCNN.

3. The computer-implemented method of claim 2, wherein: The set includes 4 to 16 DCNNs.

4. The computer-implemented method of claim 2, wherein: The set includes 8 DCNNs.

5. The computer-implemented method of any one of claims 2 to 4, further comprising: An ensemble of DCNNs are trained at each epoch in the loop.

6. The computer-implemented method of any one of claims 2 to 4, further comprising: The pathogenicity scores were used to select outlier SNPs by summing the scores from the ensemble of DCNNs.

7. The computer-implemented method of any one of claims 2 to 4, further comprising: The pathogenicity scores were used to select outlier SNPs by taking the maximum average for each SNP scored by the ensemble of DCNNs.

8. The computer-implemented method of any one of claims 1 to 4, wherein: Sampling without replacement during the current cycle results in separate training sets that do not intersect with the pseudo-pathogenic SNPs during the current cycle.

9. The computer-implemented method of any one of claims 1 to 4, wherein: The final cycle is the predetermined number of cycles.

10. The computer-implemented method of claim 9, wherein: The predetermined number of cycles is 20.

11. The computer-implemented method of any one of claims 1 to 4, wherein: The parameters are at least convolutional filter weights and a learning rate.

12. The computer-implemented method of any one of claims 1 to 4, further comprising selecting one CNN in the set of CNNs to use as a pathogenicity classifier.

13. The computer-implemented method of claim 12, further comprising selecting the one CNN in response to the one CNN being different from other CNNs in the set of CNNs on a validation sample.

14. The computer-implemented method of any one of claims 1 to 4, wherein: The universal elite set completed in the cycle has at least 400,000 elite predicted pathogenic SNPs.

15. The computer-implemented method of any one of claims 1 to 4, further comprising, in each cycle, matching trinucleotide contexts between generic benign training examples and separate pseudo pathogenic training examples to prevent mutation rate bias in elite predicted pathogenic SNPs.

16. The computer-implemented method of any one of claims 1 to 4, wherein: In each successive cycle, the sampling of individual pseudo-pathogenic training examples from unlabeled SNPs was reduced by 5%.

17. A computer-implemented system for identifying pathogenic single nucleotide polymorphisms (SNPs), comprising one or more processors coupled to a memory, wherein the memory is loaded with computer instructions for the computer-implemented method for identifying pathogenic single nucleotide polymorphisms (SNPs) according to any one of claims 1-16.

18. A non-transitory computer-readable storage medium having printed thereon computer program instructions for a computer-implemented method for identifying pathogenic single nucleotide polymorphisms (SNPs) according to any one of claims 1-16.

Images