Semi-supervised learning for training an ensemble of deep convolutional neural networks

Abstract

To disclose a technique relating to constructing a convolutional neural network-based classifier for variant classification.SOLUTION: A disclosed technique relates to training a convolutional neural network-based classifier regarding training data by using a backpropagation-based gradient update technique that progressively collates an output of a convolutional neural network-based classifier with a corresponding ground truth level. The convolutional neural network-based classifier includes groups of residual blocks. Each group of residual blocks is parameterized by the number of convolution filters in the residual blocks, a convolution window size of the residual blocks, and an expansion convolution rate of the residual blocks, the convolution window size varies between the groups of the residual blocks, and the expansion convolution rate varies between the groups of the residual blocks. Training data includes a benign training example and a pathogenic training example of a translated sequence pair generated from a benign variant and a pathogenic variant.SELECTED DRAWING: Figure 19

Description

[Technical Field]

[0001] appendix The appendix includes a bibliography of potentially relevant references listed in papers written by the inventors whose subject matter is addressed in U.S. provisional applications to which this application claims priority / benefit. These references may be made available to counsel upon request or may be available via Global Dossier. The papers are the first listed references.

[0002] Priority Application This application is a continuation of U.S. Provisional Patent Application No. 62 / 573,144, filed October 16, 2017, by Hong Gao, Kai-How Farh, Laksshman Sundaram, and Jeremy Francis McRae, entitled "Training a Deep Pathogenicity Classifier Using Large-Scale Benign Training Data" (Attorney Docket No. ILLM 1000-1 / IP-1611-PRV), and U.S. Provisional Patent Application No. 62 / 573,149, filed October 16, 2017, by Kai-How Farh, Laksshman Sundaram, Samskruthi Reddy Padigepati, and Jeremy Francis McRae, entitled "Pathogenicity Classifier Based On Deep Convolutional Neural Networks (CNNS)" (Attorney Docket No. ILLM This application claims priority to or benefit of U.S. Provisional Patent Application No. 62 / 573,153, entitled "Deep Semi-Supervised Learning that Generates Large-Scale Pathogenic Training Data," filed October 16, 2017, by Hong Gao, Kai-How Farh, Laksshman Sundaram, and Jeremy Francis McRae (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 November 7, 2017, by Hong Gao, Kai-How Farh, and Laksshman Sundaram (Attorney Docket No. ILLM 1000-4 / IP-1618-PRV). These provisional applications are incorporated herein by reference for all purposes.

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

[0004] International Patent Application No. PCT / US18 / 55840, entitled "DEEP LEARNING-BASED TECHNIQUES FOR TRAINING DEEP CONVOLUTIONAL NEURAL NETWORKS," filed October 15, 2018, by Hong Gao, Kai-How Farh, Laksshman Sundaram, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-8 / IP-1611-PCT), later published as PCT Application No. WO (TBA).

[0005] International Patent Application No. PCT / US18 / 55878, entitled "DEEP CONVOLUTIONAL NEURAL NETWORKS FOR VARIANT CLASSIFICATION," filed October 15, 2018, by Laksshman Sundaram, Kai-How Farh, Hong Gao, Samskruthi Reddy Padigepati, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-9 / IP-1612-PCT), later published as PCT Application No. WO(TBA).

[0006] and concurrently filed U.S. non-provisional patent application entitled "DEEP LEARNING-BASED TECHNIQUES FOR TRAINING DEEP CONVOLUTIONAL NEURAL NETWORKS" by Hong Gao, Kai-How Farh, Laksshman Sundaram, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-5 / IP-1611-US).

[0007] Co-filed U.S. non-provisional patent application entitled "DEEP CONVOLUTIONAL NEURAL NETWORKS FOR VARIANT CLASSIFICATION" by Laksshman Sundaram, Kai-How Farh, Hong Gao, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-6 / IP-1612-US).

[0008] and concurrently filed U.S. non-provisional patent application entitled "SEMI-SUPERVISED LEARNING FOR TRAINING AN ENSEMBLE OF DEEP CONVOLUTIONAL NEURAL NETWORKS" by Laksshman Sundaram, Kai-How Farh, Hong Gao, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-7 / IP-1613-US).

[0009] Document 1 - AVDOord, 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

[0010] Document 2 - SOArik, 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

[0011] Document 3 - F. Yu and V. Koltun, "MULTI-SCALE CONTEXT AGGREGATION BY DILATED CONVOLUTIONS", arXiv:1511.07122, 2016

[0012] Document 4 - K. He, X. Zhang, S. Ren, and J. Sun, "DEEP RESIDUAL LEARNING FOR IMAGE RECOGNITION," arXiv:1512.03385, 2015

[0013] Document 5 - R.K. Srivastava, K. Greff, and J. Schmidhuber, "HIGHWAY NETWORKS," arXiv:1505.00387, 2015

[0014] Document 6 - G. Huang, Z. Liu, L. van der Maaten, and KQWeinberger, "DENSELY CONNECTED CONVOLUTIONAL NETWORKS", arXiv:1608.06993, 2017

[0015] Document 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

[0016] Document 8 - S. Ioffe, and C. Szegedy, "BATCH NORMALIZATION: ACCELERATING DEEP NETWORK TRAINING BY REDUCING INTERNAL COVARIATE SHIFT," arXiv:1502.03167, 2015

[0017] Document 9 - J.M.W. Olterink, T.Leiner, M.A.Viergever, and I.Isgum, "DILATED CONVOLUTIONAL NEURAL NETWORKS FOR CARDIOVASCULAR MR SEGMENTATION IN CONGENITAL HEART DISEASE," arXiv:1704.03669, 2017

[0018] Document 10 - LCPiqueras, "AUTOREGRESSIVE MODEL BASED ON A DEEP CONVOLUTIONAL NEURAL NETWORK FOR AUDIO GENERATION", Tampere University of Technology, 2016

[0019] Document 11 - J.Wu, "Introduction to Convolutional Neural Networks", Nanjing University, 2017

[0020] Document 12 - I.J.Goodfellow, D.Warde-Farley, M.Mirza, A.Courville, and Y.Bengio, "CONVOLUTIONAL NETWORKS," Deep Learning, MIT Press, 2016

[0021] 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

[0022] Document 1 describes a deep convolutional neural network architecture that uses a group of residual blocks with convolution filters having the same convolution window size, a batch normalization layer, a rectified linear unit (abbreviated as ReLU) layer, a dimensionality transformation layer, a dilated convolution layer with an exponentially increasing atrous convolution rate, a skip connection, and a softmax classification layer to accept an input sequence and generate an output sequence that scores 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 parameters of the neural network components described in Document 1. For example, unlike Document 1, the dilated convolution rate in the disclosed technology increases 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.

[0023] Document 2 provides details of the deep convolutional neural network architecture described in Document 1.

[0024] Document 3 describes dilated convolution used by the disclosed technology. In this specification, dilated convolution is also referred to as "dilated convolution." Dilated / dilated convolution enables a large receptive field with a small number of trainable parameters. Dilated / dilated convolution is a convolution in which the kernel is applied over an area longer than its length by skipping input values ​​using a step, also called the dilated convolution rate or dilation factor. Dilated / dilated convolution adds separation between elements of the convolution filter / kernel so that a longer interval of adjacent input entries (e.g., nucleotides, amino acids) is considered when the convolution operation is performed. This allows for the incorporation of long-range contextual dependencies in the input. Dilated convolution saves partial convolution calculation results for reuse as adjacent nucleotides are processed.

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

[0026] Document 5 describes skip connections used by the disclosed technology. In this specification, skip connections are also referred to as "highway networks."

[0027] Document 6 describes the densely connected convolutional network architecture used by the disclosed technology.

[0028] Document 7 describes a dimensional transformation convolution layer and a module-based processing pipeline used by the disclosed technology. An example of a dimensional transformation convolution is a 1x1 convolution.

[0029] Document 8 describes the batch normalization layer used by the disclosed technique.

[0030] Document 9 also describes the dilation / expansion convolution used by the disclosed technique.

[0031] 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 dilated / extended convolutions.

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

[0033] Document 12 provides details of various convolution operations that may be used by the disclosed technology.

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

[0035] Incorporation by reference of tables submitted electronically with this application The following table files in ASCII text format are submitted herewith and incorporated by reference. The file names, creation dates, and sizes are as follows:

[0036] SupplementaryTable1.txt October 2, 2018 13KB

[0037] SupplementaryTable2.txt October 2, 2018 13KB

[0038] SupplementaryTable3.txt October 2, 2018 11KB

[0039] SupplementaryTable4.txt October 2, 2018 13KB

[0040] SupplementaryTable6.txt October 2, 2018 12KB

[0041] SupplementaryTable7.txt October 2, 2018 44KB

[0042] SupplementaryTable13.txt October 2, 2018 119KB

[0043] SupplementaryTable18.txt October 2, 2018 35KB

[0044] SupplementaryTable20.txt October 2, 2018 1027KB

[0045] SupplementaryTable20Summary.txt October 2, 2018 9KB

[0046] SupplementaryTable21.txt October 2, 2018 24KB

[0047] SupplementaryTable21.txt October 2, 2018 24KB

[0048] SupplementaryTable18.txt October 4, 2018 35KB

[0049] DataFileS1.txt, October 4, 2018, 138MB

[0050] DataFileS2.txt, October 4, 2018, 980MB

[0051] DataFileS3.txt October 4, 2018 1.01MB

[0052] DataFileS4.txt October 4, 2018 834KB

[0053] Pathogenicity_prediction_model.txt October 4, 2018 8.24KB

[0054] Supplementary Table 1: Details of the variants from each species used in the analysis. This table contains intermediate results in the pipeline for each of these data sources. Note that this table is provided in SupplementaryTable1.txt.

[0055] Supplementary Table 2: Depletion rates of missense variants present in other species at common human allele frequencies. The depletion rates were calculated based on the missense:synonymous ratio in common variants (>0.1%) compared to rare variants (<0.1%), using variants that were identical-by-state between humans and other species. Note that this table is provided in SupplementaryTable2.txt.

[0056] Supplementary Table 3: Depletion rates of missense variants present in other species at common human allele frequencies, constrained to only genes with an average nucleotide conservation rate of greater than 50% between humans and other mammals. The depletion rates were calculated based on the missense:synonymous ratio of common variants (>0.1%) compared to rare variants (<0.1%), using variants that were identical between humans and other species. Note that this table is provided in SupplementaryTable3.txt.

[0057] Supplementary Table 4: Depletion rates of missense variants present as fixed substitutions in related species pairs at common human allele frequencies. This depletion rate was calculated based on the missense:synonymous ratio in common variants (>0.1%) compared to rare variants (<0.1%), using variants that were identical between humans and related species pairs. Note that this table is provided in SupplementaryTable4.txt.

[0058] Supplementary Table 6: Domain-specific annotation of the SCN2A gene. Wilcoxon rank sum p-values ​​indicate the difference in PrimateAI scores in specific domains compared to the entire protein. The domains highlighted in bold cover approximately 7% of the protein but have the majority of ClinVar pathogenicity annotations. This corresponds well with the average PrimateAI score for those domains, which are the top three pathogenic domains according to the PrimateAI model. Note that this table is provided in SupplementaryTable6.txt.

[0059] Supplementary Table 7: Raw counts used in calculating the effect of allele frequency on the expected missense:synonymous ratio. Expected counts of synonymous and missense variants were calculated based on variants in intronic regions, using trinucleotide context to account for mutation rate and gene conversion. Note that this table is provided in SupplementaryTables.xlsx.

[0060] Supplementary Table 13: List of protein names from the Protein DataBank (PDB) used to train deep learning models for three-state secondary structure and three-state solvent accessibility prediction. The label column indicates whether the protein is used in the training / validation / validation stages of model training. Note that this table is provided in SupplementaryTable13.txt.

[0061] Supplementary Table 18: List of 605 genes that were nominally significant for association with disease in the DDD study, calculated from protein-truncating mutations only (p<0.05). Note that this table is provided in SupplementaryTable18.txt.

[0062] Supplementary Table 20: Gene-by-gene de novo mutation (DNM) enrichment test results for all genes with at least one observed DNM. P values ​​are given when including all DNMs and after removing missense DNMs with a PrimateAI score less than 0.803. FDR-corrected P values ​​are also given. Counts of observed protein-truncating (PTV) DNMs and missense DNMs from the DDD cohort alone and from the full meta-analysis cohort are also included. Similar counts of observed and predicted missense DNMs are also included, first when including all missense DNMs and second after removing all missense DNMs with a PrimateAI score less than 0.803. Note that this table is provided in SupplementaryTable20.txt and SupplementaryTable20Summary.txt.

[0063] Supplementary Table 21: Results of testing the enrichment of de novo mutations in genes with FDR<0.1. Counts of observed protein-truncating (PTV) de novo mutations and counts of other protein-transforming de novo mutations are included, once for all missense de novo mutations and once for damaging missense mutations only. P-values ​​are given when including all missense sites compared to P-values ​​after excluding low-scoring missense sites. Note that this table is provided in SupplementaryTable21.txt.

[0064] DataFileS1: List of all variants present in other species. The "ClinVar Significance" column contains the available, non-conflicting ClinVar annotations. Note that this table is provided in DataFileS1.txt.

[0065] DataFileS2: List of all fixed substitutions from related species pairs. Note that this table is provided in DataFileS2.txt.

[0066] DataFileS3: List of withheld benign test variants with primates and IBS. A benign test variant is an uncommon human variant with one or more primate species and IBS. Note that this table is provided in DataFileS3.txt.

[0067] DataFileS4: List of unlabeled variants in primates and IBS that match the retained benign test variants. Unlabeled variants match the benign test variants for mutation rate, coverage bias, and alignability to primate species. Note that this table is provided in DataFileS4.txt.

[0068] Pathogenicity_prediction_model: Code in the Python programming language that enables the disclosed techniques according to one implementation. Note that this code file is provided in Pathogenicity_prediction_model.txt.

[0069] Field of the Disclosed Technology The disclosed technology relates to artificial intelligence type computers and digital data processing systems and corresponding data processing methods and products for the emulation of intelligence (i.e., knowledge-based systems, inference systems, knowledge acquisition systems), including systems for reasoning with uncertainty (e.g., fuzzy logic systems), adaptive systems, machine learning systems, and artificial neural networks. Specifically, the disclosed technology relates to using deep learning-based techniques to train deep convolutional neural networks. [Background technology]

[0070] The subject matter discussed in this section should not be deemed to be prior art merely as a result of its mention in this section. Similarly, the problems mentioned in this section or associated with the subject matter provided as background should not be deemed to have been previously recognized in the prior art. The subject matter in this section merely represents different approaches, which may themselves correspond to implementations of the claimed technology.

[0071] [Machine Learning] In machine learning, input variables are used to predict an output variable. The input variables are often called features and are expressed as X=(X1,X2,...,X k ), and for each X i are the features. The output variables are often called the response or dependent variables, and the variable Y i The relationship between Y and the corresponding X can be written in the following general form: Y=f(x)+∈

[0072] In the above equation, f is the feature (X1, X2, ..., X k ), where ∈ is a random error term. This error term is independent of X and has a mean of 0.

[0073] In practice, feature X is available without Y or without knowing the exact relationship between X and Y. The error term has zero mean, so the goal is to estimate f.

[0074]

number

[0075] In the above formula,

[0076]

number

[0077] is an estimate of ∈, which is often considered a black box,

[0078]

number

[0079] This means that only the relationship between the input and output of the equation is known, but the question of why this works remains unanswered.

[0080] function

[0081]

number

[0082] is discovered using learning. Supervised learning and unsupervised learning are two methods used in machine learning for this task. In supervised learning, labeled data is used for training. By indicating the inputs and the corresponding outputs (=labels), the function

[0083]

number

[0084] is optimized to approximate the output. In unsupervised learning, the goal is to find hidden structure from unlabeled data. The algorithm has no measure of accuracy on the input data, which distinguishes it from supervised learning.

[0085] [Neural Networks] FIG. 1A shows one implementation of a fully connected neural network with multiple layers. 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 the hidden layer, and two neurons in the output layer. The hidden layer has an activation function f(·), and the output layer has an activation function g(·). These connections are weighted by numerical weights (e.g., w ) that are adjusted during the training process so that a properly trained network will respond correctly when presented with an image to recognize. 11 , w 21 , w 12 , w 31 , w 22 , w 32 , v 11 , v 22 ). The input layer processes the raw input, and the hidden layer processes the output from the input layer based on the weights of the connections between the input and hidden layers. The output layer takes the output from the hidden layer and processes it based on the weights of the connections between the hidden and output layers. The network contains multiple layers of feature detection neurons, each with a number of neurons that correspond to different combinations of inputs from the previous layer. These layers are constructed so that the first layer detects a set of basic patterns in the input image data, the second layer detects patterns of patterns, and the third layer detects patterns of those patterns.

[0086] An overview of the application of deep learning in genetics can be found in the following publications: · T. Ching et al., Opportunities And Obstacles For Deep Learning In Biology And Medicine, www.biorxiv.org:142760, 2017 · Angermueller C, Parnamaa T, Parts L, Stegle O, Deep Learning For Computational Biology. Mol Syst Biol. 2016;12:878 · Park Y, Kellis M, 2015 Deep Learning For Regulatory Genomics. Nat. Biotechnol. 33, 825-826, (doi:10.1038 / nbt.3313) Min S, Lee B, and Yoon S, Deep Learning in Bioinformatics. Brief. Bioinform. bbw068 (2016) · Leung MK, Delong A, Alipanahi B, et al., Machine Learning In Genomic Medicine: A Review of Computational Problems and Data Sets, 2016 · Libbrecht MW, Noble WS, Machine Learning Applications In Genetics and Genomics. Nature Reviews Genetics 2015;16(6):321-32 [Prior art documents] [Patent documents]

[0087] [Patent Document 1] International Patent Application No. PCT / US18 / 55840 [Patent Document 2] International Patent Application No. PCT / US18 / 55878 [Patent Document 3] U.S. Patent Application No. (Undefined) (Attorney Docket No. ILLM 1000-5 / IP-1611-US) [Patent Document 4] U.S. Patent Application No. (Undefined) (Attorney Docket No. ILLM 1000-6 / IP-1612-US) [Patent Document 5] U.S. Patent Application No. (Undefined) (Attorney Docket No. ILLM 1000-7 / IP-1613-US) [Patent Document 6] International Patent Application Publication No. WO07010252 [Patent Document 7] International Patent Application No. PCTGB2007 / 003798 [Patent Document 8] U.S. Patent Application Publication No. 2009 / 0088327 [Patent Document 9] U.S. Patent Application Publication No. 2016 / 0085910 [Patent Document 10] U.S. Patent Application Publication No. 2013 / 0296175 [Patent Document 11] International Patent Application Publication No. WO 04 / 018497 [Patent Document 12] U.S. Patent No. 7,057,026 [Patent Document 13] International Patent Application Publication No. WO 91 / 06678 [Patent Document 14] International Patent Application Publication No. WO 07 / 123744 [Patent Document 15] U.S. Patent No. 7,329,492 [Patent Document 16] U.S. Patent No. 7,211,414 [Patent Document 17] U.S. Patent No. 7,315,019 [Patent Document 18] U.S. Patent No. 7,405,281 [Patent Document 19] U.S. Patent Application Publication No. 2008 / 0108082 [Patent Document 20] U.S. Patent No. 5,641,658 [Patent Document 21] U.S. Patent Application Publication No. 2002 / 0055100 [Patent Document 22] U.S. Patent No. 7,115,400 [Patent Document 23] U.S. Patent Application Publication No. 2004 / 0096853 [Patent Document 24] U.S. Patent Application Publication No. 2004 / 0002090 [Patent Document 25] U.S. Patent Application Publication No. 2007 / 0128624 [License 26] U.S. Patent and Trademark Publication No. 2008 / 0009420 [License 27] U.S. Patent and Trademark Publication No. 2007 / 0099208A1 [License 28] U.S. Patent and Trademark Publication No. 2007 / 0166705A1 [License 29] U.S. Patent and Trademark Publication No. 2008 / 0280773A1 [License 30] U.S. Patent Application No. 13 / 018255 [Patent Document 31] International Patent Application Publication No. WO 2014 / 142831 [Non-licensed literature]

[0088] [Non-licensed Document 1] AVDOord, S.Dieleman, H.Zen, K.Simonyan, O.Vinyals, A.Graves, N.Kalchbrenner, A.Senior, K.Kavukcuoglu, "WAVENET: A GENERATIVE MODEL FOR RAW AUDIO", arXiv:1609.03499, 2016 [Non-licensed Document 2] SOArik, M.Chrzanowski, A.Coates, G.Diamos, A.Gibiansky, Y.Kang, X.Li, J.Miller, A.Ng, J.Raiman, S.Sengupta, およびM.Shoeybi, "DEEP VOICE: REAL-TIME NEURAL TEXT-TO-SPEECH", arXiv:1702.07825, 2017 [Non-licensed Document 3] F.YuおよびV.Koltun, "MULTI-SCALE CONTEXT AGGREGATION BY DILATED CONVOLUTIONS", arXiv:1511.07122, 2016 [Non-licensed Document 4] K. He, X. Zhang, S. Ren, and J. Sun, “DEEP RESIDUAL LEARNING FOR IMAGE RECOGNITION”, arXiv:1512.03385, 2015 [Non-Patent Document 5] R.K. Srivastava, K. Greff, and J. Schmidhuber, "HIGHWAY NETWORKS," arXiv:1505.00387, 2015 [Non-patent document 6] G. Huang, Z. Liu, L. van der Maaten, and KQWeinberger, “DENSELY CONNECTED CONVOLUTIONAL NETWORKS”, arXiv:1608.06993, 2017 [Non-Patent Document 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 [Non-patent document 8] S. Ioffe, and C. Szegedy, “BATCH NORMALIZATION: ACCELERATING DEEP NETWORK TRAINING BY REDUCING INTERNAL COVARIATE SHIFT,” arXiv:1502.03167, 2015 [Non-Patent Document 9] JMWolterink, T. Leiner, MAViergever, and I. Isgum, “DILATED CONVOLUTIONAL NEURAL NETWORKS FOR CARDIOVASCULAR MR SEGMENTATION IN CONGENITAL HEART DISEASE”, arXiv:1704.03669, 2017 [Non-Patent Document 10] LCPiqueras, “AUTOREGRESSIVE MODEL BASED ON A DEEP CONVOLUTIONAL NEURAL NETWORK FOR AUDIO GENERATION”, Tampere University of Technology, 2016 [Non-Patent Document 11] J. Wu, "Introduction to Convolutional Neural Networks", Nanjing University, 2017 [Non-Patent Document 12] I. J. Goodfellow, D. Warde-Farley, M. Mirza, A. Courville, and Y. Bengio, "CONVOLUTIONAL NETWORKS," Deep Learning, MIT Press, 2016 [Non-Patent 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 [Non-Patent Document 14] T. Ching et al., Opportunities And Obstacles For Deep Learning In Biology And Medicine, www.biorxiv.org:142760, 2017 [Non-Patent Document 15] Angermueller C, Parnamaa T, Parts L, Stegle O, Deep Learning For Computational Biology. Mol Syst Biol. 2016;12:878 [Non-Patent Document 16] Park Y, Kellis M, 2015 Deep Learning For Regulatory Genomics. Nat. Biotechnol. 33, 825-826, (doi:10.1038 / nbt.3313) [Non-Patent Document 17] Min S, Lee B, and Yoon S, Deep Learning In Bioinformatics. Brief. Bioinform. bbw068 (2016) [Non-Patent Document 18] Leung MK, Delong A, Alipanahi B, et al. Machine Learning In Genomic Medicine: A Review of Computational Problems and Data Sets, 2016 [Non-Patent Document 19] Libbrecht MW, Noble WS, Machine Learning Applications In Genetics and Genomics. Nature Reviews Genetics 2015;16(6):321-32 [Non-Patent Document 20] K. He, X. Zhang, S. Ren, and J. Sun, “DEEP RESIDUAL LEARNING FOR IMAGE RECOGNITION”, arXiv:1512.03385, 2015 [Non-Patent Document 21] Bentley et al., Nature 456:53-59 (2008) [Non-Patent Document 22] Lizardi et al., Nat. Genet. 19:225-232 (1998) [Non-Patent Document 23] Dunn, Tamsen and Berry, Gwenn and Emig-Agius, Dorothea and Jiang, Yu and Iyer, Anita and Udar, Nitin and Stromberg, Michael, 2017, Pisces: An Accurate and Versatile Single Sample Somatic and Germline Variant Caller, 595-595, 10.1145 / 3107411.3108203 [Non-Patent Document 24] T Saunders, Christopher and Wong, Wendy and Swamy, Sajani and Becq, Jennifer and J Murray, Lisa and Cheetham, Keira, 2012, Strelka: Accurate somatic small-variant calling from sequenced tumor-normal sample pairs, Bioinformatics (Oxford, England), 28, 1811-7, 10.1093 / bioinformatics / bts271 [Non-Patent Document 25] Kim, S, Scheffler, K, Halpern, AL, Bekritsky, MA, Noh, E, Kallberg M, Chen, X, Beyter, D, Krusche, P, and Saunders, CT, 2017, Strelka2: Fast and accurate variant calling for clinical sequencing applications. [Non-Patent Document 26] As described in Stromberg, Michael and Roy, Rajat and Lajugie, Julien and Jiang, Yu and Li, Haochen and Margulies, Elliott, 2017, Nirvana: Clinical Grade Variant Annotator, 596-596, 10.1145 / 3107411.3108204, Illumina Inc. Summary of the Invention [Means for solving the problem]

[0089] In the drawings, like reference characters generally refer to like parts throughout the various views. Also, the drawings are not necessarily to scale, emphasis instead generally being placed upon 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: [Brief explanation of the drawings]

[0090] [Figure 1A] FIG. 1 illustrates an implementation of a feedforward neural network with multiple layers. [Figure 1B] FIG. 1 is a diagram of one implementation of the operation of a convolutional neural network. [Figure 1C] FIG. 1 is a block diagram for training a convolutional neural network according to one implementation of the disclosed technology. [Figure 1D] FIG. 10 is a diagram of an implementation of a subsampling layer (average / max pooling) in accordance with an implementation of the disclosed technology. [Figure 1E] FIG. 1 illustrates an implementation of a ReLU nonlinear layer in accordance with an implementation of the disclosed technology. [Figure 1F] FIG. 1 illustrates one implementation of two layers of convolution. [Figure 1G] FIG. 10 illustrates residual connections that reinject previous information downstream via the addition of feature maps. [Figure 1H]FIG. 1 illustrates one implementation of residual blocks and skip connections. [Figure 1I] FIG. 1 illustrates a batch normalization forward pass. [Figure 1J] FIG. 1 illustrates batch normalization transformation during testing. [Figure 1K] FIG. 1 illustrates a batch normalization backward pass. [Figure 1L] FIG. 1 illustrates the use of batch normalization layers after and before convolutional or densely connected layers. [Figure 1M] FIG. 1 illustrates one implementation of 1D convolution. [Figure 1N] FIG. 1 illustrates how Global Average Pooling (GAP) works. [Figure 1O] FIG. 1 illustrates dilated convolution. [Figure 1P] FIG. 1 illustrates an implementation of stacked dilated convolution. [Figure 1Q] FIG. 1 illustrates an exemplary computing environment in which the disclosed technology can operate. [Figure 2] FIG. 1 illustrates an exemplary architecture of a deep residual network for pathogenicity prediction, referred to herein as "PrimateAI." [Figure 3] Schematic of PrimateAI, a deep learning network architecture for pathogenicity classification. [Figure 4A] Supplementary Table 16 details an exemplary model architecture for the pathogenicity prediction deep learning model PrimateAI. [Figure 4B] Supplementary Table 16 details an exemplary model architecture for the pathogenicity prediction deep learning model PrimateAI. [Figure 4C] Supplementary Table 16 details an exemplary model architecture for the pathogenicity prediction deep learning model PrimateAI. [Figure 5] FIG. 1 illustrates the deep learning network architecture used to predict protein secondary structure and solvent accessibility. [Figure 6] FIG. 1 illustrates the deep learning network architecture used to predict protein secondary structure and solvent accessibility. [Figure 7A] Supplementary Table 11 shows details of an exemplary model architecture for a three-state secondary structure prediction deep learning (DL) model. [Figure 7B] Supplementary Table 11 shows details of an exemplary model architecture for a three-state secondary structure prediction deep learning (DL) model. [Figure 8A] Supplementary Table 12 shows details of an exemplary model architecture for a three-state solvent accessibility prediction deep learning model. [Figure 8B] Supplementary Table 12 shows details of an exemplary model architecture for a three-state solvent accessibility prediction deep learning model. [Figure 9] FIG. 1 illustrates one implementation of generating reference and alternative protein sequences from benign and pathogenic variants. [Figure 10] FIG. 1 illustrates one implementation of aligning a reference protein sequence with an alternative protein sequence. [Figure 11] A diagram showing one implementation of generating a position-specific frequency matrix (abbreviated as PFM), also called a position-specific weight matrix (abbreviated as PWM) or a position-specific scoring matrix (abbreviated as PSSM). [Figure 12] FIG. 1 shows the processing of secondary structure and solvent accessibility sub-networks. [Figure 13] FIG. 1 shows the processing of secondary structure and solvent accessibility sub-networks. [Figure 14] FIG. 1 shows the processing of secondary structure and solvent accessibility sub-networks. [Figure 15] FIG. 1 shows the processing of secondary structure and solvent accessibility sub-networks. [Figure 16]

[0023] Figure 1 illustrates the operation of the variant pathogenicity classifier. As used herein, the term variant also refers to single nucleotide polymorphisms (abbreviated as SNPs), and generally refers to single nucleotide variants (abbreviated as SNVs). [Figure 17] FIG. 2 illustrates a residual block. [Figure 18] FIG. 1 shows the neural network architecture of the secondary structure and solvent accessibility subnetworks. [Figure 19] FIG. 1 shows the neural network architecture of the variant pathogenicity classifier. [Figure 20] FIG. 1 shows predicted pathogenicity scores at each amino acid location in the SCN2A gene, annotated for important functional domains. [Figure 21D] FIG. 1 shows a comparison of classifiers in predicting benign outcomes on a test set of 10,000 common primate variants that were withheld from training. [Figure 21E] Distribution of PrimateAI prediction scores for de novo missense variants occurring in patients with Deciphering Developmental Disorders (DDD) compared with unaffected siblings, along with corresponding Wilcoxon rank sum P values. [Figure 21F] Comparison of classifiers in separating de novo missense variants in DDD cases vs. controls. Wilcoxon rank sum test P values ​​are shown for each classifier. [Figure 22A] FIG. 1 shows the enrichment of de novo missense mutations in affected individuals from the DDD cohort that was greater than expected within 605 associated genes that was significant for de novo protein-truncating mutations (P<0.05). [Figure 22B] Distribution of PrimateAI prediction scores for de novo missense variants occurring in DDD patients vs. unaffected siblings within 605 relevant genes, along with corresponding Wilcoxon rank sum P-values. [Figure 22C] FIG. 1 shows a comparison of different classifiers in separating de novo missense variants in cases vs. controls within 605 genes. [Figure 22D] FIG. 1 shows a comparison of different classifiers plotted on receiver operating characteristic curves, with the area under the curve (AUC) shown for each classifier. [Figure 22E] FIG. 1 shows the classification accuracy and area under the curve (AUC) for each classifier. [Figure 23A] FIG. 1 illustrates the impact on classification accuracy of data used for training. [Figure 23B] FIG. 1 illustrates the impact on classification accuracy of data used for training. [Figure 23C] FIG. 1 illustrates the impact on classification accuracy of data used for training. [Figure 23D] FIG. 1 illustrates the impact on classification accuracy of data used for training. [Figure 24] FIG. 1 shows correcting for the effect of sequencing coverage on the ascertainment of common primate variants. [Figure 25A] FIG. 1 illustrates protein motif recognition by the disclosed neural network. [Figure 25B] FIG. 1 illustrates protein motif recognition by the disclosed neural network. [Figure 25C] FIG. 1 illustrates protein motif recognition by the disclosed neural network. [Figure 26] 1 illustrates protein motif recognition by the disclosed neural network, including line plots showing the effect of perturbing each location in and around the variant on the predicted deep learning score for the variant. [Figure 27] FIG. 1 shows that the correlation pattern of weights follows the BLOSUM62 and Grantham scoring matrices. [Figure 28A] Figure 1 shows performance evaluation of deep learning networks PrimateAI and other classifiers. [Figure 28B]Figure 1 shows performance evaluation of deep learning networks PrimateAI and other classifiers. [Figure 28C] Figure 1 shows performance evaluation of deep learning networks PrimateAI and other classifiers. [Figure 29A] FIG. 1 shows the distribution of prediction scores for four classifiers. [Figure 29B] FIG. 1 shows the distribution of prediction scores for four classifiers. [Figure 30A] Figure 1 compares the accuracy of the PrimateAI network and other classifiers in separating pathogenic and benign variants in 605 disease-associated genes. [Figure 30B] Figure 1 compares the accuracy of the PrimateAI network and other classifiers in separating pathogenic and benign variants in 605 disease-associated genes. [Figure 30C] Figure 1 compares the accuracy of the PrimateAI network and other classifiers in separating pathogenic and benign variants in 605 disease-associated genes. [Figure 31A] FIG. 1 shows the correlation between classifier performance on expert-curated ClinVar variants and performance on an empirical dataset. [Figure 31B] FIG. 1 shows the correlation between classifier performance on expert-curated ClinVar variants and performance on an empirical dataset. [Figure 32] Supplementary Table 14 shows the performance of the three-state secondary structure prediction model and the three-state solvent accessibility prediction model on annotated samples from the Protein Databank. [Figure 33] Supplementary Table 15 shows the performance comparison of deep learning networks using annotated secondary structure labels of human proteins from the DSSP database. [Figure 34]Supplementary Table 17 shows accuracy values ​​for the 10,000 retained primate variants and p-values ​​for de novo variants in DDD cases vs. controls for each of the 20 classifiers evaluated. [Figure 35] Supplementary Table 19 shows a comparison of the performance of different classifiers on de novo variants in the DDD case vs. control datasets, constrained to 605 disease-associated genes. [Figure 36] FIG. 1 illustrates a computing environment for the disclosed semi-supervised learner. [Figure 37] FIG. 1 illustrates the various cycles of the disclosed semi-supervised learning. [Figure 38] FIG. 1 illustrates the various cycles of the disclosed semi-supervised learning. [Figure 39] FIG. 1 illustrates the various cycles of the disclosed semi-supervised learning. [Figure 40] FIG. 1 illustrates the various cycles of the disclosed semi-supervised learning. [Figure 41] FIG. 1 illustrates the various cycles of the disclosed semi-supervised learning. [Figure 42] FIG. 1 illustrates an iterative balanced sampling process. [Figure 43] FIG. 1 illustrates one implementation of a computing environment used to generate a benign dataset. [Figure 44] FIG. 1 shows one implementation of generating benign human missense SNPs. [Figure 45] 1 shows one implementation of a human orthologous missense SNP. A missense SNP in a non-human species with a canonical codon and an alternative codon that matches in humans. [Figure 46] FIG. 1 illustrates one implementation of classifying SNPs in non-human primate species (eg, chimpanzees) with matching reference codons in humans as benign. [Figure 47] FIG. 1 illustrates one implementation of calculating enrichment scores and comparing them. [Figure 48] FIG. 1 illustrates one implementation of a benign SNP dataset. [Figure 49A] FIG. 1 shows missense:synonymous ratios across the human allele frequency spectrum. [Figure 49B] FIG. 1 shows missense:synonymous ratios across the human allele frequency spectrum. [Figure 49C] FIG. 1 shows missense:synonymous ratios across the human allele frequency spectrum. [Figure 49D] FIG. 1 shows missense:synonymous ratios across the human allele frequency spectrum. [Figure 49E] FIG. 1 shows missense:synonymous ratios across the human allele frequency spectrum. [Figure 50A] FIG. 1 shows the purifying selection against missense variants in the same state as in other species. [Figure 50B] FIG. 1 shows the purifying selection against missense variants in the same state as in other species. [Figure 50C] FIG. 1 shows the purifying selection against missense variants in the same state as in other species. [Figure 50D] FIG. 1 shows the purifying selection against missense variants in the same state as in other species. [Figure 51] FIG. 1 shows the expected missense:synonymous ratio across the human allele frequency spectrum in the absence of purifying selection. [Figure 52A] FIG. 1 shows missense:synonymous ratios for CpG and non-CpG variants. [Figure 52B] FIG. 1 shows missense:synonymous ratios for CpG and non-CpG variants. [Figure 52C] FIG. 1 shows missense:synonymous ratios for CpG and non-CpG variants. [Figure 52D] FIG. 1 shows missense:synonymous ratios for CpG and non-CpG variants. [Figure 53] FIG. 1 shows the missense:synonymous ratios of human variants in the same state as six primates. [Figure 54] FIG. 1 shows the missense:synonymous ratios of human variants in the same state as six primates. [Figure 55] FIG. 1 shows the missense:synonymous ratios of human variants in the same state as six primates. [Figure 56] Simulations showing the saturation of new common missense variants discovered by increasing the size of the human cohort studied. [Figure 57] FIG. 1 shows the accuracy of PrimateAI across different conservation profiles in genomes. [Figure 58] Supplementary Table 5 shows contributions to the labeled benign training dataset from common human variants and variants present in non-human primates. [Figure 59] Supplementary Table 8 shows the effect of allele frequency on the expected missense:synonymous ratio. [Figure 60] Supplementary Table 9 shows the ClinVar analysis. [Figure 61] 10 is Supplementary Table 10 showing the number of missense variants from other species found in ClinVar, according to one implementation. [Figure 62] Table 1 shows one implementation of the discovery of 14 additional gene candidates in intellectual disability. [Figure 63] Table 2 shows one implementation of the mean difference in Grantham score between pathogenic and benign variants in ClinVar. [Figure 64] FIG. 1 illustrates one implementation of gene-by-gene enrichment analysis. [Figure 65] FIG. 1 illustrates one implementation of genome-wide enrichment analysis. [Figure 66] FIG. 1 is a simplified block diagram of a computer system that can be used to implement the disclosed techniques. DETAILED DESCRIPTION OF THE INVENTION

[0091] The following discussion is presented to enable any person skilled in the art to make and use the disclosed technology and is provided in the context of a particular application and its requirements. Various modifications to the disclosed implementations will be readily 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 is to be accorded the widest scope consistent with the principles and features disclosed herein.

[0092] [introduction] [Convolutional Neural Networks] 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 the input feature space, while convolutional layers learn local patterns. In the case of images, patterns are found in small 2D windows of the input. This important 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.

[0093] Regarding the first property, after learning a pattern in the bottom right corner of a photo, the convolutional layer can recognize it anywhere, for example, in the top left corner. A densely connected network would have to re-learn the pattern if it appeared in a new location. This makes convolutional neural networks data-efficient because they require fewer training samples to learn a representation with generalization capabilities.

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

[0095] Convolutional neural networks learn highly nonlinear mappings by interconnecting layers of artificial neurons arranged in many different layers with activation functions that make the layers interdependent. Convolutional neural networks contain one or more convolutional layers interspersed with one or more subsampling and nonlinear layers, which are usually followed by one or more fully connected layers. Each element in 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 identical weights. These local, shared weights reduce the network's complexity, allowing convolutional neural networks to avoid the complexity of data reconstruction during the feature extraction and regression or classification processes when multidimensional input data enters the network.

[0096] The convolution is performed over a 3D tensor called a feature map, which has two spatial axes (height and width) and a depth axis (also called the channel axis). In an RGB image, the depth axis has dimension 3 because the image has three color channels: red, green, and blue. In a black-and-white photograph, the depth is 1 (level of gray). The convolution operation extracts patches from the input feature map, applies the same transformation to all of these patches, and generates an output feature map. This output feature map is still a 3D tensor, with width and height. Its depth can be arbitrary because the output depth is a parameter of the layer, and different channels in the depth axis no longer represent specific colors as in the RGB input, but rather filters. Filters encode specific aspects of the input data; at a high level, a single filter can encode the concept of, for example, the presence of a face in the input.

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

[0098] Convolution is defined by two key parameters: (1) the size of the patch extracted from the input, which is usually 1x1, 3x3, or 5x5, and (2) the depth of the output feature map, the number of filters computed by the convolution. Often, these start at a depth of 32, continue to a depth of 64, and end at a depth of 128 or 256.

[0099] Convolution works by sliding these windows of size 3x3 or 5x5 across the 3D input feature map, stopping at each location and extracting a 3D patch (shape (window_height, window_width, input_depth)) of surrounding features. Each such 3D patch is then converted (via a tensor product with the same learned weight matrix, called the convolution kernel) into a 1D vector of shape (output_depth). All of these vectors are then spatially reassembled into a 3D output map of shape (height, width, output_depth). Each spatial location in the output feature map corresponds to the same location 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, in a 3x3 window, the vector output [i,j,:] comes from the 3D patch input [i-1:i+1,j-1:J+1,:]. The complete process is detailed in Figure 1B.

[0100] A convolutional neural network (CNNN) comprises a convolutional layer that performs convolution operations between input values ​​and a convolutional filter (a matrix of weights), which is learned through numerous gradient update iterations during training. If (m,n) is the filter size and W is the matrix of weights, the convolutional layer convolves 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 by which the convolutional filter slides across the input is called the stride, and the filter area (m×n) is called the receptive field. The same convolutional filter is applied across different locations of the input, which reduces the number of weights to be learned. This also enables position-invariant learning, i.e., if an important pattern is present in the input, the convolutional filter learns the important pattern regardless of where it is located in the sequence.

[0101] [Training a convolutional neural network] 1C shows a block diagram of training a convolutional neural network according to one implementation of the disclosed technology. The convolutional neural network is tuned or trained so that input data leads to a particular output estimate. The convolutional neural network is tuned using backpropagation based on a comparison of the output estimate to ground truth until the output estimate asymptotically matches or approaches the ground truth.

[0102] Convolutional neural networks are trained by adjusting the weights between neurons based on the difference between the ground truth and the actual output, which can be expressed mathematically as follows:

[0103]

number

[0104] where δ = (ground truth) - (actual output)

[0105] In one implementation, the training rules are defined as follows: w nm ←w nm +α(t m -φ m )α n

[0106] In the above equation, the arrows indicate the update of the value, and t m is the target value of neuron m, and φ m is the calculated current output of neuron m, and α n is the input n and α is the learning rate.

[0107] An intermediate step in training involves generating feature vectors from input data using convolutional layers. Starting at the output, gradients are calculated for the weights in each layer. This is called a backward pass, or going back. Weights in the network are updated using a combination of negative gradients and previous weights.

[0108] In one implementation, the convolutional neural network uses a stochastic gradient update algorithm (such as ADAM) that performs backpropagation of errors via gradient descent. An example of a sigmoid function-based backpropagation algorithm is written as follows:

[0109]

number

[0110] In the sigmoid function above, h is the weighted sum calculated by the neuron. The sigmoid function has the following derivative:

[0111]

number

[0112] The algorithm involves computing the activations of all neurons in the network and producing an output for the forward pass. The activation of neuron m in the hidden layer can be written as:

[0113]

number

[0114] This is done for all hidden layers to obtain activations written as

[0115]

number

[0116] Then the error and correction weights are calculated for each layer. The error at the output is calculated as follows: δ ok =(t k -φ k )φ k (1-φ k )

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

[0118]

number

[0119] The output layer weights are updated as follows: v mk ←v mk +αδ ok φ m

[0120] The hidden layer weights are updated using a learning rate α as follows: v nm ←w nm +αδ hm a n

[0121] In one implementation, the convolutional neural network uses gradient descent optimization to calculate the error across all layers. In such optimization, given an input feature vector x and a predicted output

[0122]

number

[0123] When the goal is y,

[0124]

number

[0125] The loss function is defined as the cost of predicting l for

[0126]

number

[0127] The expected output is

[0128]

number

[0129] is transformed from the input feature vector x using a function f, which is parameterized by the weights of the convolutional neural network, i.e.

[0130]

number

[0131] The loss function is

[0132]

number

[0133] , 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:

[0134]

number

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

[0136] where α is the learning rate, and the loss is calculated as the average over a set of n data pairs. The calculation terminates when the learning rate α becomes small enough for linear convergence. In other implementations, for computational efficiency, the gradient is calculated using only selected data pairs that are fed into the Nesterov accelerated gradient and adaptive gradient methods.

[0137] In one implementation, the convolutional neural network uses stochastic gradient descent (SGD) to compute the cost function. SGD calculates the gradient with respect to the weights in the loss function for one randomized data pair z, as described below: t It is approximated by calculating only from v t+1 =μv-α∇wQ(z t ,w t ) w t+1 =w t +v t+1

[0138] where α is the learning rate, μ is the momentum, and t is the current weight state before the update. The convergence speed of SGD is approximately O(1 / t) when the learning rate α decreases sufficiently fast and when it decreases sufficiently slow. In other implementations, convolutional neural networks use different loss functions, such as Euclidean loss and softmax loss. In further implementations, the Adam stochastic optimizer is used by convolutional neural networks.

[0139] [Convolutional layer] The convolutional layer of a convolutional neural network functions as a feature extractor. It acts as an adaptive feature extractor that can learn and decompose input data into hierarchical features. In one implementation, a convolutional layer takes two images as input and produces a third image as output. In such an implementation, 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 to the input image to produce the output image. Thus, for an input vector f of length n and a kernel g of length m, the convolution of f and g, f*g, is defined as follows:

[0140]

number

[0141] A convolution operation involves sliding a kernel across the input image. For each location in the kernel, overlapping values ​​from the kernel and input image are multiplied and the results are added. The sum of these products is the value of the output image at the point in the input image where the kernel is centered. The different outputs obtained from multiple kernels are called feature maps.

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

[0143] Subsampling Layer 1D is an implementation of a subsampling layer according to one implementation of the disclosed technology. The subsampling layer reduces the resolution of features extracted by the convolutional layer to make the extracted features or feature maps robust to noise and distortion. In one implementation, the subsampling layer utilizes two types of pooling operations: average pooling and max pooling. The pooling operations divide the input into non-overlapping two-dimensional spaces. In average pooling, the average of four values ​​in a region is calculated. In max pooling, the maximum of the four values ​​is selected.

[0144] In one implementation, a subsampling layer includes a pooling operation on a set of neurons in the previous layer by mapping its output to only one of the inputs in max pooling, and to the average of the inputs in average pooling. In max pooling, the output of a pooling neuron is φ0=max(φ1,φ2,...,φ N ) is the maximum value present in the input, as described by

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

[0146] In average pooling, the output of the pooling neuron is

[0147]

number

[0148] is the average value of the input values ​​present with the set of input neurons, as described by

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

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

[0151] Nonlinear Layer FIG. 1E shows an implementation of a nonlinear layer according to one implementation of the disclosed technology. The nonlinear layer uses different nonlinear trigger functions to signal distinct identities of likely features on each hidden layer. The nonlinear layer uses various unique functions to implement nonlinear triggering, including rectified linear unit (ReLU), hyperbolic tangent, absolute hyperbolic tangent, sigmoid, and continuous trigger (nonlinear) functions. In one implementation, the ReLU activation implements the function y=max(x,0), keeping the input and output sizes of the layer the same. The advantage of using ReLU is that convolutional neural networks are trained more times and faster. ReLU is a discontinuous, nonsaturating activation function that is linear with respect to the input when the input is greater than or equal to 0 and 0 otherwise. Mathematically, the ReLU activation function is written as follows: φ(h)=max(h,0)

[0152]

number

[0153] In another implementation, the convolutional neural network: φ(h)=(a+bh) c We use a power unit activation function, which is a continuous, non-saturating function described by:

[0154] where a, b, and c are parameters that control the shift, scale, and power, respectively. The power activation function can produce activations that are asymmetric in x and y when c is odd, and y-symmetric when c is even. In some implementations, this unit produces non-normalized linear activations.

[0155] In yet another implementation, the convolutional neural network is configured to compute the logistic function

[0156]

number

[0157] We use a sigmoid unit activation function, which is a continuous saturating function described by:

[0158] where β = 1. The sigmoid unit activation function does not produce negative activations and is asymmetric only about the y-axis.

[0159] [Convolution example] FIG. 1F shows an implementation of two convolution layers. In FIG. 1F, an input with a dimension of 2048 is convolved. In Convolution 1, the input is convolved by a convolution layer consisting of two channels of 16 kernels of size 3×3. The resulting 16 feature maps are then normalized by a ReLU activation function with ReLU1 and then pooled in Pool 1 by average pooling using a 16-channel pooling layer with a kernel of size 3×3. In Convolution 2, the output of Pool 1 is then convolved by another convolution layer consisting of 16 channels of 30 kernels with a dimension of 3×3. This is followed by average pooling in Pool 2 with yet another ReLU2 and a kernel size of 2×2. The convolution layers use variable strides and padding, for example, 0, 1, 2, and 3. The resulting feature vector is 512-dimensional according to one implementation.

[0160] In other implementations, the convolutional neural network uses a different number 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 a larger number of neurons per layer, for example, 1, 2, or 3 fully connected layers with 100 to 200 neurons per layer. In another implementation, the convolutional neural network is a deep network with more layers and fewer neurons per layer, for example, 5, 6, or 8 fully connected layers with 30 to 50 neurons per layer.

[0161] [Forward Pass] The output of the neuron at row x, column y in the lth convolutional layer and kth feature map for f convolutional cores in the feature map is determined by the following formula:

[0162]

number

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

[0164]

number

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

[0166]

number

[0167] [Backpropagation] The output deviation of the kth neuron in the output layer is determined by the following formula:

[0168]

number

[0169] The input deviation of the kth neuron in the output layer is determined by the following formula:

[0170]

number

[0171] The weight and bias variance of the kth neuron in the output layer is determined by the following formula:

[0172]

number

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

[0174]

number

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

[0176]

number

[0177] The weight and bias variances for row x, column y in the mth feature map of the previous layer, which receives input from k neurons in the hidden layer, are determined by the following formula:

[0178]

number

[0179] The output bias for row x, column y in the mth feature map of subsample layer S is determined by the following formula:

[0180]

number

[0181] The input bias for row x, column y in the mth feature map of subsample layer S is determined by the following formula:

[0182]

number

[0183] The weight and bias variance in row x, column y in the mth feature map of subsample layer S and convolutional layer C is determined by the following formula:

[0184]

number

[0185] The output bias for row x, column y in the kth feature map of convolutional layer C is determined by the following formula:

[0186]

number

[0187] The input bias for row x, column y in the kth feature map of convolutional layer C is determined by the following formula:

[0188]

number

[0189] The weight and bias variance in row r, column c of the mth convolutional core of the kth feature map of the lth convolutional layer C:

[0190]

number

[0191] Residual Connection FIG. 1G illustrates residual connections that reinject previous information downstream via feature map addition. Residual connections comprise reinjecting previous representations into the downstream flow of data by adding past output tensors to later output tensors, which helps prevent information loss along the data processing flow. 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 comprise making the outputs of earlier layers available as inputs to later layers, effectively creating a shortcut in the sequential network. Earlier outputs are summed with later activations rather than concatenated to them, which assumes both activations are the same size. If they are different sizes, a linear transformation can be used to reshape the earlier activations into the target shape. 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 by reference for all purposes as if fully set forth herein.

[0192] Residual learning and skip connections Figure 1H shows one implementation of the residual block and skip connections. The main idea of ​​residual learning is that the residual mapping is much easier to learn than the original mapping. The residual network stacks several residual units to mitigate the degradation of training accuracy. The residual block utilizes an extra skip connection to eliminate gradient vanishing in deep neural networks. At the beginning of the residual block, the data flow is separated into two streams: the first stream carries the block's unchanged input, while the second stream applies weights and nonlinearities. At the end of the block, the two streams are combined using element-wise summation. The main advantage of such a construction is that it allows gradients to flow more easily through the network. Additional information about residual blocks and skip connections can be found in A. V. Dord, 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.

[0193] By benefiting from residual networks, deep convolutional neural networks (CNNs) can be easily trained to achieve improved accuracy in image classification and object detection. A convolutional feedforward network connects the output of the lth layer as the input to the (l+1)th layer, which leads to the following layer transitions: x l =H l (x l-1 ) The residual block is the identity function: x l =H l (x l-1 )+x l-1 The advantage of the residual block is that gradients can flow directly from later layers to earlier layers through the identity function. However, the identity function and H lThe outputs of are combined by addition, which can hinder the information flow in the network.

[0194] [Dilated Fold] Figure 1O shows a dilated convolution. Dilated convolution, sometimes called dilated convolution, literally means "with holes." Its name comes from the French algorithme a trous, which computes the fast dyadic wavelet transform. In these types of convolutional layers, the inputs corresponding to the receptive field of the filter are not adjacent points. This is shown in Figure 1O. The distance between the inputs depends on the dilation factor.

[0195] [WaveNet] WaveNet is a deep neural network for generating raw audio waveforms. WaveNet is distinguished from other convolutional networks because it can capture a relatively large "visual field" at low cost. Furthermore, it is possible to add local and global signal conditioning, which allows WaveNet to be used as a text-to-speech (TTS) engine with multiple voices, where TTS imposes local, voice-specific, and global conditioning.

[0196] The main building block of WaveNet is the causal dilated convolution. As an extension of causal dilated convolution, WaveNet allows for the stacking of these convolutions, as shown in Figure 1P. To obtain the same receptive field using dilated convolutions in this diagram, another dilated layer is required. Stacking is a repetition of dilated convolutions, connecting the outputs of dilated convolution layers to a single output. This allows WaveNet to obtain a large "visual" field for one output node at a relatively low computational cost. For comparison, to obtain a visual field of 512 inputs, a fully convolutional network (FCN) requires 511 layers. For a dilated convolutional network, eight layers are required. Stacked dilated convolutions require only seven layers for a two-layer stack or six layers for a four-layer stack. To get an idea of ​​the difference in computational power required to cover the same visual field, the following table shows the number of weights required in a network under the assumption of one filter per layer and a filter width of 2. It is further assumed that the network uses 8-bit binary encoding.

[0197] [Table 1]

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

[0199] [Batch Normalization] Batch normalization is a method for accelerating deep network training by making data normalization an integral part of the network architecture. Batch normalization can adaptively normalize data even as the mean and variance change over time during training. Batch normalization works by internally maintaining an exponential moving average of the mean and variance for each batch of data seen during training. The main impact of batch normalization is that, much like residual connections, it aids gradient propagation, thereby enabling deep networks. Some ultra-deep networks can only be trained if they contain multiple batch normalization layers. Additional information about batch normalization can be found in S. Ioffe and C. Szegedy, "BATCH NORMALIZATION: ACCELERATING DEEP NETWORK TRAINING BY REDUCING INTERNAL COVARIATE SHIFT," arXiv:1502.03167, 2015, incorporated herein by reference for all purposes as if fully set forth herein.

[0200] Batch normalization can be viewed as yet another layer that can be inserted into a model architecture, like a fully connected or convolutional layer. Batch normalization layers are typically used after convolutional or densely connected layers. Batch normalization layers can also be used before convolutional or densely connected layers. Both implementations can be used with the disclosed technology and are shown in FIG. 1L. Batch normalization layers take an axis argument, which specifies the feature axis to be normalized. This argument defaults to 1, which is the last axis in the input tensor. This is the correct value when using Dense, Conv1D, RNN, and Conv2D layers with data_format set to "channels_last." However, in the niche use case of Conv2D layers with data_format set to "channels_first," the feature axis is axis 1, and the axis argument in batch normalization can be set to 1.

[0201] Batch normalization provides a definition for feeding forward inputs and computing gradients with respect to parameters and its own inputs through 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 an activation function. In a convolutional layer, different elements of the same feature map at different positions, i.e., activations, are normalized in the same way to follow the nature of convolution. Therefore, all activations in a mini-batch are normalized across all positions, not just activation-by-activation.

[0202] Internal covariate shifts are a major reason why deep architectures are notoriously slow to train: this is due to the fact that deep networks not only have to learn new representations at each layer, but also account for changes in their distributions.

[0203] Covariate shift in general is a known problem in the domain of deep learning and frequently occurs in real-world problems. A common covariate shift problem is the difference in distributions between training and test sets, which can lead to suboptimal generalization performance. This problem is usually addressed by a standardization or whitening preprocessing step. However, the whitening operation in particular is computationally intensive and therefore impractical in online situations, especially when covariate shift occurs across different layers.

[0204] Internal covariate shift is the phenomenon where the distribution of network activations changes across multiple layers due to changes in network parameters during training. Ideally, each layer should be transformed into a space where each layer has the same distribution but the functional relationships remain the same. To avoid expensive computation of covariance matrices to decorrelate and whiten the data at each layer and step, we normalize the distribution of each input feature at each layer across each mini-batch to have a mean of 0 and a standard deviation of 1.

[0205] Forward Pass During the forward pass, the mean and variance of the mini-batch 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 scale and shift parameters. Batch normalization forward pass f BN is illustrated in Figure 1I.

[0206] In Figure 1I, μ B is the batch mean,

[0207]

number

[0208] is the batch variance. The learned scale and shift parameters are denoted as γ and β, respectively. For simplicity, the batch normalization procedure is described here on a per-activation basis, omitting the corresponding indices.

[0209] Because normalization is a differentiable transformation, errors are propagated to these learned parameters, and thus it is possible to restore the network's recall ability by learning an identity transformation. Conversely, by learning scale and shift parameters identical to the corresponding batch statistics, where that was the optimal operation to perform, the batch normalization transformation has no effect on the network. At testing time, the batch mean and variance are replaced by the statistics of their respective populations, since the inputs do not depend on other samples from the mini-batch. An alternative approach is to keep averaging the batch statistics during training and use them to calculate the network output at testing time. At testing time, the batch normalization transformation can be represented as shown in Figure 1J. In Figure 1J, μ D and

[0210]

number

[0211] are not batch statistics, but represent the population mean and variance, respectively.

[0212] Backward Pass Since normalization is a differentiable operation, the backward pass can be calculated as illustrated in FIG. 1K.

[0213] [1D convolution] 1D convolution extracts local 1D patches or subsequences from a sequence, as shown in Figure 1M. 1D convolution obtains each output time step from a temporal patch in the input sequence. 1D convolution layers recognize local patterns in a sequence. Because the same input transformation is performed for each patch, a pattern learned at one location in the input sequence can be recognized later at a different location, making 1D convolution layer transforms invariant to temporal transformations. For example, a 1D convolution layer processing a sequence of bases using a convolution window of size 5 should be able to learn base sequences of length 5 or less and recognize base motifs in any context in the input sequence. Therefore, base-level 1D convolution can learn about base morphology.

[0214] [Global Average Pooling] Figure 1N shows how global average pooling (GAP) works. Global average pooling can be used to replace a fully connected (FC) layer for classification by spatially averaging features in the final layer for scoring. This reduces the training load and bypasses the overfitting problem. Global average pooling applies a structural transformation before the model, which is equivalent to a linear transformation with predefined weights. Global average pooling reduces the number of parameters and eliminates the fully connected layer. The fully connected layer is usually the layer with the most parameters and connections, and global average pooling provides a much lower-cost approach to achieving similar results. The main idea of ​​global average pooling is to generate the average value from the feature maps of each final layer as a reliability coefficient for scoring and feed it directly into the softmax layer.

[0215] Global average pooling has three advantages: (1) there are no redundant parameters in the global average pooling layer, so overfitting in the global average pooling layer is avoided; (2) the output of global average pooling is the average of the entire feature map, so global average pooling is more robust to spatial transformations; and (3) due to the large number of parameters in the fully connected layer, which usually accounts for more than 50% of all parameters in the entire network, replacing them with a global average pooling layer can significantly reduce the model size, which makes global average pooling very useful in model compression.

[0216] Global average pooling makes sense because stronger features in the final layer are expected to have higher mean values. In some implementations, global average pooling can be used as a proxy for classification scores. The feature map under global average pooling can be interpreted as a confidence map, which can enforce a correspondence between the feature map and the category. Global average pooling can be particularly effective when the features in the final layer are sufficiently abstracted for direct classification. However, global average pooling itself is not sufficient when multi-level features are to be combined into groups such as submodels, which is best performed by adding a simple fully connected layer or other classifier after global average pooling.

[0217] [Deep Learning in Genomics] Genetic variation can help explain many diseases. Each human has a unique genetic code, and within a group of individuals, there are many genetic variants. The majority of harmful genetic variants have been depleted from the genome by natural selection. It is important to identify which genetic variants are likely to be pathogenic or harmful. This will help researchers focus on genetic variants that are likely to be pathogenic, accelerating the diagnosis and treatment of many diseases.

[0218] Modeling the properties and functional impact (e.g., pathogenicity) of variants is an important but challenging task in the field of genomics. Despite rapid advances in functional genome sequencing technologies, interpreting the functional consequences of variants poses significant challenges due to the complexity of cell-type-specific transcriptional regulatory systems.

[0219] Advances in biochemical technology over the past few decades have led to next-generation sequencing (NGS) platforms that rapidly generate genomic data at a much lower cost than ever before. Such overwhelming amounts of sequenced DNA remain challenging to annotate. Supervised machine learning algorithms typically perform well when large amounts of labeled data are available. In bioinformatics and many other data-rich training methods, the process of labeling instances is expensive. However, unlabeled instances are inexpensive and readily available. In scenarios where the amount of labeled data is relatively small and the amount of unlabeled data is significantly larger, semi-supervised learning offers a cost-effective alternative to manual labeling.

[0220] Opportunities arise to use semi-supervised algorithms to build deep learning-based pathogenicity classifiers that accurately predict the pathogenicity of variants, resulting in a database of pathogenic variants free from human diagnostic bias.

[0221] With respect to pathogenicity classifiers, deep neural networks are a type of artificial neural network that uses multiple nonlinear complex transformation layers to continuously model high-level features. Deep neural networks provide feedback via backpropagation, which conveys the difference between observed and predicted outputs, to adjust parameters. Deep neural networks have evolved with the availability of large training datasets, parallel and distributed computing capabilities, and sophisticated training algorithms. Deep neural networks have driven significant advances in many areas, such as computer vision, speech recognition, and natural language processing.

[0222] 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 due to their architecture, which includes convolutional layers, nonlinear layers, and pooling layers. Recurrent neural networks are designed to exploit continuous information in input data using cyclic connections between building blocks such as perceptrons, long-term memory units, and gated recurrent units. In addition, many other emerging deep neural networks, such as deep spatiotemporal neural networks, multidimensional recurrent neural networks, and convolutional autoencoders, have been proposed for limited contexts.

[0223] The goal of training a deep neural network is to optimize the weight parameters in each layer, gradually synthesizing simpler features into more complex ones so that the best hierarchical representation can be learned from the data. A single cycle of the optimization process is organized as follows: First, given a training dataset, a forward pass calculates the output in each layer in turn and propagates a function signal forward through the network. At the final output layer, an objective loss function measures the error between the inferred output and a given label. To minimize the training error, a backward pass uses the chain rule to backpropagate the error signal and calculates the gradients for all weights throughout 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 examples. Several optimization algorithms are derived from stochastic gradient descent. For example, the Adagrad and Adam training algorithms perform stochastic gradient descent while adaptively modifying the learning rate based on the update frequency and the moments of the gradient for each parameter, respectively.

[0224] Another core element in training deep neural networks is regularization, which refers to a strategy intended to achieve good generalization performance by avoiding overfitting. For example, weight decay adds a penalty term to the objective loss function so that weight parameters converge to smaller absolute values. Dropout randomly removes hidden units from a neural network during training, which can be viewed as an ensemble of possible subnetworks. To improve the dropout capability, a new activation function, maxout, and a variant of dropout for recurrent neural networks, called rnnDrop, have been proposed. Furthermore, batch normalization provides a new regularization method through normalization of scalar features for each activation within a mini-batch and learning their respective means and variances as parameters.

[0225] Given the multidimensional and high-dimensional nature of sequencing data, deep neural networks hold great promise for bioinformatics research due to their broad applicability and high predictive power. Convolutional neural networks have been adapted 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, repetitive local patterns in DNA that are predicted to have significant biological functions. A distinctive feature of convolutional neural networks is the use of convolutional filters. Unlike traditional classification methods based on carefully designed, human-crafted features, convolutional filters perform adaptive feature learning, a process similar to mapping raw input data into useful representations of knowledge. In this sense, convolutional filters function as a series of motif scanners, because a set of such filters can recognize relevant patterns in the input and update them during the training procedure. Recurrent neural networks can capture long-range dependencies in continuous data of variable length, such as protein or DNA sequences.

[0226] Therefore, powerful computational models for predicting variant pathogenicity could have enormous benefits for both basic science and translational research.

[0227] Common polymorphisms represent natural experiments whose health has been tested by natural selection over many generations. By comparing allele frequency distributions for missense and synonymous substitutions in humans, we found that the presence of a missense variant at high allele frequency in non-human primate species predicts with high confidence that the variant is also under natural selection in human populations. In contrast, common variants in more distantly related species are subject to negative selection as evolutionary distance increases.

[0228] We utilize common mutations from six non-human primate species to train a semi-supervised deep learning network that accurately classifies clinical de novo missense mutations using sequence alone. With over 500 known species, the primate lineage contains enough common mutations to systematically model the effects of most human variants of unknown significance.

[0229] The human reference genome harbors 70 million potentially protein-altering missense substitutions, the majority of which are rare mutations whose impact on human health has not been characterized. These variants of unknown significance pose challenges to genome interpretation for clinical applications and are an obstacle to the long-term adoption of sequencing for population-wide screening and personalized medicine.

[0230] Cataloging common mutations across diverse human populations is a useful strategy for identifying clinically benign mutations, but the common mutations available from modern humans are limited by bottleneck events in our species' distant past. Humans and chimpanzees share 99% sequence identity, suggesting that natural selection acting on chimpanzee variants may model the effects of identical variants in humans. Because the average coalescent time for natural polymorphisms in human populations is a fraction of the species' divergence time, naturally occurring chimpanzee mutations largely span a mutation space that does not overlap with human mutations, except for rare cases of haplotypes maintained by balancing selection.

[0231] The recent availability of aggregated exon data from 60,706 individuals made it possible to test this hypothesis by comparing the allele frequency spectrum for missense and synonymous mutations. Singleton variants in ExAC closely match the expected 2.2:1 missense:synonymous ratio predicted by de novo mutations after adjusting for mutation rate using trinucleotide context, but at higher allele frequencies, the number of observed missense variants decreases due to the removal of deleterious variants by natural selection. The pattern of missense:synonymous ratios across the allele frequency spectrum indicates that the majority of missense variants with population frequencies below 0.1% are mildly deleterious, i.e., not sufficiently pathogenic to warrant their immediate removal from the population, nor sufficiently neutral to be tolerated at high allele frequencies, consistent with previous observations on more limited population data. These findings support the widespread empirical practice in clinics of eliminating variants with allele frequencies higher than 0.1%–1% as likely benign for penetrant genetic disorders, with a few well-documented exceptions caused by balancing selection and founder effects.

[0232] When we repeated this analysis for a subset of human variants identical to common chimpanzee variants (observed more than once in sequencing chimpanzee populations), we found that the missense:synonymous ratio remained roughly constant across the allele frequency spectrum. The high allele frequencies of these variants in chimpanzee populations indicate that these variants have already been subjected to natural selection in chimpanzees, and their neutral effects on human population fitness provide remarkable evidence that selection pressures against missense variants are highly congruent in the two species. The lower missense:synonymous ratio observed in chimpanzees is consistent with a larger effective population size in the ancestral chimpanzee population, which allowed for efficient elimination of mildly deleterious variants.

[0233] In contrast, rare chimpanzee variants (observed only once in sequencing chimpanzee populations) show a modest decline in the missense:synonymous ratio at higher allele frequencies. Simulating a cohort of the same size from human mutation data, we estimate that only 64% of variants observed once in a cohort of this size have an allele frequency higher than 0.1% in the overall population, compared with 99.8% of variants seen multiple times in the cohort having an allele frequency higher than 0.1% in the overall population, indicating that not all rare chimpanzee variants are subject to selection. Overall, we estimate that 16% of confirmed chimpanzee missense variants have an allele frequency lower than 0.1% in the overall population and are subject to negative selection at higher allele frequencies.

[0234] We next characterized human variants that are identical to those observed in other nonhuman primate species (bonobos, gorillas, orangutans, rhesus macaques, and marmosets). Similar to chimpanzees, we observed that missense:synonymous ratios were roughly equivalent across the allele frequency spectrum, with a slight depletion of missense mutations at high allele frequencies, likely due to the inclusion of a small number of rare variants (approximately 5-15%). These results suggest that selective pressures for missense variants are largely consistent within the primate lineage, at least up to New World monkeys, estimated to have diverged from the ancestral human lineage approximately 35 million years ago.

[0235] Human missense variants with identical status to variants in other primates are strongly enriched for benign outcomes in ClinVar. After excluding variants with unknown or conflicting annotations, we observe that human variants with primate orthologs are annotated as benign or likely benign in ClinVar approximately 95% of the time, compared with 45% for missense mutations overall. The small proportion of ClinVar variants classified as pathogenic from nonhuman primates is comparable to the proportion of pathogenic ClinVar variants that would be observed by identifying rare variants from similarly sized cohorts of healthy humans. A significant proportion of these variants annotated as pathogenic or likely pathogenic before the advent of large allele frequency databases may be evaluated differently today.

[0236] The field of human genetics has long relied on model organisms to infer the clinical impact of human mutations, but the long evolutionary distance to most genetically tractable animal models raises concerns about the generalizability of these findings to humans. To investigate the alignment of natural selection against missense variants in humans and more distantly related species, we extended our analysis beyond the primate lineage to include largely common mutations from four additional mammalian species (mouse, pig, goat, and cow) and two more distantly related vertebrates (chicken and zebrafish). In contrast to previous primate analyses, we observed a significant depletion of missense mutations at common allele frequencies compared to rare allele frequencies, particularly at greater evolutionary distances, indicating that a significant proportion of common missense mutations in more distantly related species will be subject to negative selection in human populations. Nevertheless, observation of missense variants in more distantly related vertebrates increases the probability of a benign outcome, because the proportion of common missense variants depleted by natural selection is much lower than the approximately 50% depletion rate for the reference human missense variant. Consistent with these results, we found that human missense variants observed in mice, dogs, pigs, and cattle have an approximately 85% probability of being annotated as benign or likely benign in ClinVar, compared with 95% for primate variants and 45% for the entire ClinVar database.

[0237] The existence of closely related species pairs across various evolutionary distances also provides an opportunity to assess the functional consequences of fixed missense substitutions in human populations. Within closely related species pairs (branch lengths <0.1) on the mammalian tree, fixed missense mutations are observed to be depleted at common allele frequencies relative to rare allele frequencies, indicating that a significant proportion of cross-species fixed substitutions are non-neutral in humans, even within the primate lineage. Comparison of the degree of missense depletion indicates that cross-species fixed substitutions are far less neutral than intraspecific polymorphisms. Interestingly, cross-species mutations between closely related mammals are significantly less pathogenic in ClinVar (83% more likely to be annotated as benign or likely benign) than intraspecific common polymorphisms, suggesting that these changes do not abolish protein function but rather result in modulation of protein function that confers species-specific adaptive benefits.

[0238] The large number of potential variants of unknown significance and the critical importance of accurate variant classification for clinical applications have led to numerous attempts to solve the problem using machine learning, but these efforts have been severely limited by the insufficient amount of common human variants and the questionable quality of annotations in curated databases. Mutations from six nonhuman primate species contributed over 300,000 unique missense variants that do not overlap with common human mutations and have mostly benign consequences, greatly expanding the size of training datasets available for machine learning methods.

[0239] Unlike previous models that utilize numerous human-crafted features and meta-classifiers, we apply a simple deep learning residual network that takes as input only the amino acid sequence flanking the variant of interest and orthologous sequence alignments in other species. To provide the network with information about protein structure, we train two separate networks to learn secondary structure and solvent accessibility from sequence alone, and incorporate these as subnetworks into a larger deep learning network to predict their impact on protein structure. Using sequence as a starting point avoids potential biases in protein structure and functional domain annotations that may be incompletely validated or inconsistently applied.

[0240] We use semi-supervised learning to overcome the problem of the training set containing only benign-labeled variants by first training an ensemble of networks to separate likely benign primate variants from random unknown variants that are matched for mutation rate and sequencing coverage. This ensemble of networks is used to score the full set of unknown variants and to influence the selection of unknown variants to seed the next iteration of the classifier by biasing them toward unknown variants with a predicted outcome of being more pathogenic, taking gradual steps at each iteration to prevent the model from prematurely converging to a suboptimal outcome.

[0241] Common primate variants also provide a clean evaluation dataset, completely independent of previously used training data, for evaluating existing methods, which are becoming difficult to objectively evaluate due to the proliferation of meta-classifiers. We evaluated the performance of our model alongside four other popular classification algorithms (Sift, Polyphen2, CADD, and M-CAP) using 10,000 provided common primate variants. Because approximately 50% of all human missense variants are eliminated by natural selection at common allele frequencies, we calculated the 50th percentile score for each classifier on a set of randomly selected missense variants matched by mutation rate to the 10,000 provided common primate variants and used that threshold to evaluate the submitted common primate variants. The accuracy of our deep learning model was significantly better than other classifiers on this independent evaluation dataset, whether using a deep learning network trained only on common human variants or on both common human and primate variants.

[0242] A recent trio-sequencing study catalogued thousands of de novo mutations in patients with neurodevelopmental disorders and their healthy siblings, allowing for evaluation of the strength of various classification algorithms in separating de novo missense mutations in cases vs. controls. For each of the four classification algorithms, we scored each de novo missense variant in cases vs. controls and reported p-values ​​from Wilcoxon rank-sum tests of the difference between the two distributions. In this clinical scenario, we performed a deep learning method trained on primate variants (p ∼10). -33 ) is compared with other classifiers (p approx. 10 -13 From 10 -19 ) performed much better than previously reported expectations for this cohort. Based on the approximately 1.3-fold enrichment of de novo missense variants, and on previous estimates that approximately 20% of missense variants produce loss-of-function effects, a perfect classifier would have a p of approximately 10 -40This indicates that our classifier still has room for improvement.

[0243] The accuracy of the deep learning classifier scales with the size of the training dataset, and mutation data from each of the six primate species independently contribute to improving the classifier's accuracy. The large number and diversity of nonhuman primate species, along with evidence that selection pressures on protein-altering variants are largely consistent within primate lineages, suggests systematic primate population sequencing as an effective strategy for classifying the millions of human variants of unknown significance that currently constrain clinical genome interpretation. Of the 504 known nonhuman primate species, approximately 60% are endangered due to hunting and habitat loss, motivating urgent global conservation efforts that will benefit both these unique and irreplaceable species and ourselves.

[0244] Although whole genome data are not available in as aggregated a form as exon data, limiting our ability to detect the effects of natural selection in deep intronic regions also allowed us to calculate the observed and expected counts of cryptic splicing variants far from exon regions. Overall, we observed 60% deletions in cryptic splicing variants located more than 50 nt from exon-intron boundaries. The signal attenuation is likely due to a combination of smaller sample sizes in whole genome data compared to exons and the greater difficulty in predicting the effects of deep intronic variants.

[0245] [term] All literature and similar materials cited in this application, including but not limited to patents, patent applications, articles, books, papers, and web pages, are expressly incorporated by reference in their entirety, regardless of the format of such literature and similar materials. In the event that one or more of the cited literature and similar materials, including but not limited to defined terms, term usage, techniques described, etc., differs from or conflicts with this application, this application controls.

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

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

[0248] This application uses the terms "protein" and "translated sequence" interchangeably.

[0249] This application uses the terms "codon" and "base triplet" interchangeably.

[0250] This application uses the terms "amino acid" and "translation unit" interchangeably.

[0251] 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.

[0252] The term "chromosome" refers to the genetic information carrier of a living cell, which is derived from chromatin strands comprising DNA and protein components (especially histones). The conventional internationally recognized individual human genome chromosome numbering system is utilized herein.

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

[0254] The term "sample" as used herein typically refers to a sample derived from a bodily fluid, cell, tissue, organ, or organism that contains a nucleic acid or mixture of nucleic acids containing at least one nucleic acid sequence to be sequenced and / or phased. Such samples include, but are not limited to, saliva / oral fluid, amniotic fluid, blood, blood fragments, fine needle biopsy samples (e.g., direct biopsy, fine needle biopsy, etc.), urine, peritoneal fluid, pleural fluid, tissue explants, organ cultures, and any other tissue or cell specimen, or a portion thereof or a derivative thereof, or a material isolated therefrom. While samples are often taken from human subjects (e.g., patients), samples can be taken from any organism that has chromosomes, including, but not limited to, dogs, cats, horses, goats, sheep, cows, pigs, etc. Samples may be used directly as obtained from a biological source or may be used after pretreatment to modify the sample's characteristics. For example, such pretreatment may include preparing plasma from blood, diluting mucus, etc. Pretreatment methods may also involve, but are not limited to, filtration, precipitation, dilution, distillation, mixing, centrifugation, freezing, lyophilization, concentration, amplification, nucleic acid fragmentation, inactivation of interfering factors, addition of reagents, lysis, and the like.

[0255] The term "sequence" includes or refers to a chain of nucleotides linked together. The nucleotides may be based on DNA or RNA. It should be understood that a sequence may include multiple subsequences. For example, a single sequence (e.g., a PCR amplicon) may have 350 nucleotides. A sample read may include multiple subsequences within these 350 nucleotides. For example, a sample read may include first and second flanking subsequences, e.g., having 20-50 nucleotides. The first and second flanking subsequences may be located on either side of a repetitive segment with a corresponding subsequence (e.g., 40-100 nucleotides). Each of the flanking subsequences may include (or include a portion of) a primer subsequence (e.g., 10-30 nucleotides). For ease of reading, the term "subsequence" is referred to as "sequence," but it should be understood that the two sequences are not necessarily separate from each other on a common strand. To distinguish between the various sequences described herein, the sequences may be given different labels (e.g., target sequence, primer sequence, flanking sequence, reference sequence, etc.) Other terms, such as "allele," may be given different labels to distinguish between similar entities.

[0256] 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 genome rearrangements and repeat segments, as well as gene fusions and novel transcripts. Paired-end sequencing methodologies are described in International Patent Application Publication No. WO07010252, International Patent Application No. PCTGB2007 / 003798, and US Patent Application Publication No. 2009 / 0088327, each of which is incorporated herein by reference. In one example, a series of operations can be performed as follows: (a) generate a cluster of nucleic acids; (b) linearize the nucleic acids; (c) hybridize a first sequencing primer and perform repeated cycles of extension, scanning, and deblocking as described above; (d) "reverse" the target nucleic acid on the flow cell surface by synthesizing a complementary copy; and (e) linearize the resynthesized strand. (f) Hybridizing a second sequencing primer and performing repeated cycles of extension, scanning, and deblocking as described above, where the reverse operation can be performed to introduce reagents as described above for a single cycle of bridge amplification.

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

[0258] The term "read" refers to a collection of sequence data describing a fragment of a nucleotide sample or standard. The term "read" can refer to a sample read and / or a reference read. Typically, but not necessarily, a read represents a short sequence of consecutive base pairs in a sample or standard. A read can be represented in letters (ATCG) by the base pair sequence of the sample or standard fragment. A read can be stored in a memory device and processed accordingly to determine whether the read matches a reference sequence or meets other criteria. A read can be obtained directly from a sequencing device or indirectly from stored sequence information about the sample. In some cases, a read is a DNA sequence of sufficient length (e.g., at least about 25 bp) that can be used to identify a larger sequence or region, for example, that can be aligned and specifically assigned to a chromosome or genomic region or gene.

[0259] Next-generation sequencing methods include, for example, sequencing by synthesis (Illumina), pyrosequencing (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 can vary from about 30 bp to more than 10,000 bp. For example, Illumina sequencing methods using SOLiD sequencers generate nucleic acid reads of about 50 bp. In another example, Ion Torrent sequencing can generate nucleic acid reads of up to 400 bp, while 454 pyrosequencing can generate nucleic acid reads of about 700 bp. In yet another example, single-molecule real-time sequencing methods can generate reads of 10,000 bp to 15,000 bp. Thus, in some implementations, the nucleic acid sequence reads have a length of 30-100 bp, 50-200 bp, or 50-400 bp.

[0260] The terms "sample read," "sample sequence," or "sample fragment" refer to sequence data of a genomic sequence of interest from a sample. For example, a sample read comprises sequence data from a PCR amplicon having a forward primer sequence and a reverse primer sequence. The sequence data can be obtained from any sequence selection method. A sample read can be, for example, from a sequencing-by-synthesis (SBS) reaction, a sequencing-by-ligation reaction, or any other suitable sequencing method for which it is desired to determine the length and / or identity of repetitive elements. A sample read can be a consensus (e.g., average or weighted) sequence derived from multiple sample reads. In some implementations, providing a reference sequence comprises identifying a locus of interest based on a primer sequence of a PCR amplicon.

[0261] The term "raw fragment" refers to sequence data for a portion of a target genome sequence that at least partially overlaps with a specified location or secondary location of interest within a sample read or sample fragment. Non-limiting examples of raw fragments include duplex stitched fragments, simplex stitched fragments, duplex un-stitched fragments, and simplex un-stitched fragments. The term "raw" is used to indicate that a raw fragment contains sequence data that has some relationship to the sequence data in the sample read, regardless of whether it exhibits a supporting variant that corresponds to and proves or confirms the authenticity of a potential variant in the sample read. The term "raw fragment" does not indicate that a fragment necessarily contains a supporting variant that validates a variant call in the sample read. For example, when a variant calling application determines that a sample read exhibits a first variant, the variant calling application may determine that one or more raw fragments lack a corresponding type of "supporting" variant that would be expected to be present if that variant were present in the sample read.

[0262] The terms "mapping," "aligned," "alignment," or "aligning" refer to the process of comparing a read or tag to a reference sequence, thereby determining whether the reference sequence contains the read sequence. If the reference sequence contains the read, the read may be mapped to the reference sequence, or in some implementations, may be mapped to a specific location within the reference sequence. In some cases, alignment simply indicates whether the read is a member of a specific reference sequence (i.e., whether the read is present or absent within the reference sequence). For example, alignment of a read to the reference sequence for human chromosome 13 indicates whether the read is present within the reference sequence for chromosome 13. A tool that provides this information may be called a set membership tester. In some cases, alignment additionally indicates the location within the reference sequence to which the read or tag maps. For example, if the reference sequence is the entire human genome sequence, alignment may indicate that the read is present on chromosome 13 and may further indicate that the read is located on a specific strand and / or site of chromosome 13.

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

[0264] The term "variant" refers to a nucleic acid sequence that differs from a reference nucleic acid. Typical nucleic acid sequence variants include, but are not limited to, single nucleotide polymorphisms (SNPs), short deletion and insertion polymorphisms (indels), copy number variations (CNVs), microsatellite markers, or short tandem repeats and structural variations. Somatic variant calling attempts to identify variants present at low frequencies in DNA samples. Somatic variant calling is of interest in the context of cancer treatment. Cancer is caused by the accumulation of DNA mutations. DNA samples from tumors are generally heterogeneous, containing some normal cells, some cells at early stages of cancer progression (with a small number of mutations), and some cells at later stages (with a large number of mutations). Due to this heterogeneity, somatic mutations often appear at low frequencies when sequencing tumors (e.g., from FFPE samples). For example, an SNV may be found in only 10% of reads containing a given base. A variant to be classified as somatic or germline by a variant classifier is also referred to herein as a "variant under test."

[0265] The term "noise" refers to erroneous variant calls that result from one or more errors in the sequencing process and / or variant calling application.

[0266] The term "variant frequency" refers to the relative frequency of an allele (genetic variant) at a particular locus in a population, expressed as a proportion or percentage. For example, this proportion or percentage can be the proportion of all chromosomes in the population that have that allele. As an example, a sample variant frequency refers to the relative frequency of an allele / variant at a particular locus / location along a subject's genome sequence across a "population" corresponding to the number of reads and / or samples obtained for the subject's genome sequence from an individual. As another example, a reference variant frequency refers to the relative frequency of an allele / variant at a particular locus / location along one or more reference genome sequences, where the "population" corresponds to the number of reads and / or samples obtained for one or more reference genome sequences from a population of normal individuals.

[0267] The term "variant array frequency (VAF)" refers to the percentage of sequenced reads at a target location that are observed to match a variant, divided by the total coverage. VAF is a measure of the proportion of sequenced reads that carry a variant.

[0268] The terms "location," "designated location," and "locus" refer to the position or coordinate of one or more nucleotides within a sequence of nucleotides. The terms "location," "designated location," and "locus" also refer to the position or coordinate of one or more base pairs within a sequence of nucleotides.

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

[0270] The term "threshold" herein refers to a numeric or non-numeric value used as a cutoff to characterize a sample, nucleic acid, or portion thereof (e.g., a read). The threshold may vary based on empirical analysis. The threshold may be compared to a measured or calculated value to determine whether a source suggesting such a value should be classified in a particular manner. The threshold may be identified empirically or analytically. The choice of threshold depends on the level of confidence a user desires to have in making a classification. The threshold may be chosen for a specific purpose (e.g., to balance sensitivity and selectivity). As used herein, the term "threshold" refers to a point at which the course of an analysis may be altered and / or an action may be initiated. A threshold need not be a predetermined number. Instead, the threshold may be, for example, a function based on multiple factors. The threshold may be adaptive to the situation. Moreover, a threshold may indicate an upper limit, a lower limit, or a range between limiting values.

[0271] In some implementations, a measure or score based on the sequencing data may be compared to a threshold. As used herein, the terms "measure" or "score" may include a value or result determined from the sequencing data or a function based on a value or result determined from the sequencing data. Like a threshold, a measure or score may be adaptive to the situation. For example, the measure or score may be a normalized value. As an example of a score or measure, one or more implementations may use a count score when analyzing data. The count score may be based on the number of sample reads. The sample reads may have undergone one or more filtering stages so that the sample reads have at least one common characteristic or quality. For example, each of the sample reads used to determine the count score may be aligned with a reference sequence or assigned as a potential allele. The number of sample reads with the common characteristic may be counted to determine a read count. The count score may be based on the read count. In some implementations, the count score may be a value equal to the read count. In other implementations, the count score may be based on the read count and other information. For example, the count score may be based on the read count for a particular allele of a locus and the total number of reads for the locus. In some implementations, the count score may be based on the read count for the locus and previously obtained data. In some implementations, the count score may be a normalized score between multiple predetermined values. The count score may also be a function of the read counts from other loci in the sample, or the read counts from other samples run simultaneously with the subject sample. For example, the count score may be a function of the read count for a particular allele and the read counts for other loci in the sample and / or from other samples. As an example, the read counts from other loci and / or from other samples may be used to normalize the count score for a particular allele.

[0272] The term "coverage" or "fragment coverage" refers to a count or other measure of the number of sample reads to the same fragment of a sequence. The read count may represent a count of the number of reads that cover the corresponding fragment. Alternatively, coverage may be determined by multiplying the read count by a specified factor based on historical knowledge, sample knowledge, locus knowledge, etc.

[0273] The term "read depth" (a number conventionally followed by "x") refers to the number of sequenced reads with overlapping alignments at the target location. It is often expressed as an average or as a percentage above a cutoff across a set of intervals (such as exons, genes, or panels). For example, a clinical report may state that the panel average coverage is 1.105x, with 98% of the target bases covered >100x.

[0274] The term "base call quality score" or "Q score" refers to a PHRED-scaled probability ranging from 0 to 20 that is inversely proportional to the probability that a single sequenced base is correct. For example, a T base call with a Q of 20 is considered likely to be correct with a confidence P-value of 0.01. Any base call with a Q < 20 should be considered low quality, and any identified variant where a significant proportion of sequenced reads supporting the variant are of low quality should be considered a likely false positive.

[0275] The term "variant read" or "number of variant reads" refers to the number of sequenced reads that support the presence of a variant.

[0276] [Sequencing process] The implementations described herein may be applicable to analyzing nucleic acid sequences to identify sequence variations. The implementations may be used to analyze potential variants / alleles of a genetic location / locus and determine the genotype of the locus, or in other words, to provide a genotype call for the locus. For example, nucleic acid sequences may be analyzed according to the methods and systems described in U.S. Patent Application Publication Nos. 2016 / 0085910 and 2013 / 0296175, the entire subject matter of which is expressly incorporated herein by reference.

[0277] In one implementation, the sequencing process includes receiving a sample containing or suspected of containing nucleic acids, such as DNA. The sample can be from a known or unknown source, such as an animal (e.g., a human), a plant, bacteria, or a fungus. The sample can be collected directly from the source. For example, blood or saliva can be collected directly from an individual. Alternatively, the sample may not be collected directly from the source. One or more processors then instruct the system to prepare the sample for sequencing. This preparation can include removing extraneous material and / or isolating some material (e.g., DNA). The biological sample can be prepared to contain features for a particular assay. For example, the biological sample can be prepared for sequencing-by-synthesis (SBS). In some implementations, preparing can include amplifying several regions of a genome. For example, preparing can include amplifying predetermined loci known to contain STRs and / or SNRs. The loci can be amplified using predetermined primer sequences.

[0278] The one or more processors then instruct the system to sequence the sample. Sequencing can be performed through a variety of known sequencing protocols. In certain implementations, sequencing includes SBS. In SBS, multiple fluorescently labeled nucleotides are used to sequence multiple clusters (potentially millions of clusters) of amplified DNA present on a surface of an optical substrate (e.g., a surface that at least partially defines a channel in a flow cell). The flow cell can contain a nucleic acid sample for sequencing, where the flow cell is placed in a suitable flow cell holder.

[0279] Nucleic acids can be prepared with known primer sequences flanking an unknown target sequence. To initiate the first SBS sequencing cycle, one or more differentially labeled nucleotides and a DNA polymerase can be flowed into or through a flow cell by a fluidics subsystem. Either a single type of nucleotide can be added at a time, or the nucleotides used in the sequencing procedure can be specifically designed with invertible ends, allowing each cycle of the sequencing reaction to occur simultaneously in the presence of several types of labeled nucleotides (e.g., A, C, T, G). The nucleotides can contain detectable label moieties, such as fluorescent dyes. When four nucleotides are mixed together, the polymerase can select the correct base to incorporate, extending each sequence by one base. Unincorporated nucleotides can be washed away by flowing a wash solution through the flow cell. One or more lasers can excite the nucleic acids to induce fluorescence. The fluorescence emitted from the nucleic acids is based on the fluorescent dye of the incorporated base, and different fluorescent dyes can emit light at different wavelengths. A deblocking reagent can be added to the flow cell to remove the invertible terminal group from the extended and detected DNA strand. The deblocking reagent can then be washed away by passing a wash solution through the flow cell. The flow cell is then ready for another cycle of sequencing, starting with the introduction of a labeled nucleotide as described above. The fluidic and detection operations can be repeated several times to complete the sequencing run.Exemplary sequencing methods are described, for example, in Bentley et al., Nature 456:53-59 (2008), International Patent Application Publication No. WO 04 / 018497, U.S. Patent No. 7,057,026, International Patent Application Publication No. WO 91 / 06678, International Patent Application 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.

[0280] In some implementations, nucleic acids can be attached to a surface and amplified before or during sequencing.For example, amplification can be performed using bridge amplification 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), as described, for example, 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.

[0281] One exemplary SBS protocol utilizes modified nucleotides with removable 3' blocks, as described, for example, in International Patent Application 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, repeated cycles of SBS reagents can be introduced into a flow cell to which target nucleic acids are attached, for example, as a result of a bridge amplification protocol. Nucleic acid clusters can be converted into single-stranded forms using a linearization solution. The linearization solution can include, for example, a restriction endonuclease capable of cleaving one strand of each cluster. Other methods of cleavage can be used as alternatives to restriction enzymes or nicking enzymes, including, inter alia, chemical cleavage (e.g., cleavage of diol bonds using periodate), cleavage of abasic sites by cleavage with endonucleases (e.g., "USER", part number M5505S, as supplied by NEB, Ipswich, Massachusetts, USA) upon exposure to heat or alkali, cleavage of ribonucleotides that would otherwise be incorporated into amplification products consisting of deoxyribonucleotides, photochemical cleavage, or cleavage of peptide linkers. After the linearization procedure, a sequencing primer can be introduced into the flow cell under conditions for hybridization of the sequencing primer to the target nucleic acid to be sequenced.

[0282] The flow cell can then be contacted with an SBS extension reagent having a modified nucleotide with a removable 3' block and a fluorescent label under conditions to extend the primers hybridized to each target nucleic acid by adding a single nucleotide. Only a single nucleotide is added to each primer because, once the modified nucleotide is incorporated into the growing polynucleotide chain complementary to the region of the template being sequenced, the polymerase cannot add additional nucleotides because there is no free 3'-OH group available to direct further sequence extension. The SBS extension reagent can be removed and replaced with a scanning reagent containing a component that protects the sample under excitation by radiation. Exemplary components of the scanning reagent are described in U.S. Patent Application Publication No. 2008 / 0280773A1 and U.S. Patent Application No. 13 / 018255, each of which is incorporated herein by reference. The extended nucleic acid can then be detected by fluorescence in the presence of the scanning reagent. Once fluorescence is detected, the 3' block can be removed using a deblocking reagent appropriate for the blocking group used. Exemplary deblocking reagents useful for each blocking group are described in International Patent Application Publication No. WO004018497, U.S. Patent Application Publication No. 2007 / 0166705A1, and U.S. Patent No. 7,057,026, each of which is incorporated herein by reference. The deblocking reagent can be washed away, leaving the target nucleic acid hybridized to the extended primer with a 3'OH group, which then becomes available for the addition of additional nucleotides. Thus, cycles of adding an extension reagent, a scanning reagent, and a deblocking reagent, with optional washing between one or more of the operations, can be repeated until the desired sequence is obtained. The above cycles can be performed using a single extension reagent introduction operation per cycle, known to correspond to a specific base, when each modified nucleotide is labeled with a different label. The different labels facilitate the differentiation of the nucleotides added during each incorporation operation.Alternatively, each cycle may involve a separate operation of extending reagent introduction followed by a separate operation of scanning reagent introduction and detection, in which case two or more of the nucleotides may have the same label and can be distinguished based on a known order of introduction.

[0283] Although the sequencing operations are discussed above with respect to specific SBS protocols, it will be understood that other protocols for sequencing and any of a variety of other molecular analytical methods may be performed as desired.

[0284] One or more processors of the system then receive the sequencing data for subsequent analysis. The sequencing data can be formatted in various ways, such as a .BAM file. The sequencing data can include, for example, several sample reads. The sequencing data can include multiple sample reads with corresponding sample sequences of nucleotides. Although only one sample read is discussed, it should be understood that the sequencing data can include, for example, hundreds, thousands, hundreds of thousands, or millions of sample reads. Different sample reads can have different numbers of nucleotides. For example, the sample reads can range from 10 nucleotides to about 500 or more nucleotides. The sample reads can span the entire genome of the source. In one example, the sample reads are directed to a predetermined locus, such as a locus with a suspected STR or a suspected SNP.

[0285] Each sample read may include a sequence of nucleotides, which may be referred to as a sample sequence, sample fragment, or target sequence. The sample sequence may include, for example, a primer sequence, a flanking sequence, and a target sequence. The number of nucleotides in a sample sequence may include 30, 40, 50, 60, 70, 80, 90, 100, or more. In some implementations, one or more of the sample reads (or sample sequences) include at least 150 nucleotides, 200 nucleotides, 300 nucleotides, 400 nucleotides, 500 nucleotides, or more. In some implementations, a sample read may include more than 1,000 nucleotides, more than 2,000 nucleotides, or more. A sample read (or sample sequence) may include a primer sequence at one or both ends.

[0286] Next, one or more processors analyze the sequencing data to obtain potential variant calls and sample variant frequencies for the sample variant calls. This operation may also be referred to as a variant call application or variant caller. Thus, the variant caller identifies or detects variants, and the variant classifier classifies the detected variants as somatic or germline. Alternative variant callers may be utilized in accordance with implementations herein, where different variant callers may be used based on the type of sequencing operation being performed, the characteristics of the sample of interest, etc. One non-limiting example of a variant calling application is the Pisces™ application by Illumina Inc. (San Diego, CA), hosted at https: / / github.com / Illumina / Pisces and described in the article Dunn, Tamsen and Berry, Gwenn and Emig-Agius, Dorothea and Jiang, Yu and Iyer, Anita and Udar, Nitin and Stromberg, Michael, (2017), Pisces: An Accurate and Versatile Single Sample Somatic and Germline Variant Caller, 595-595, 10.1145 / 3107411.3108203, the full subject matter of which is incorporated herein by reference in its entirety.

[0287] Such a variant calling application may comprise four sequentially executed modules.

[0288] (1) Pisces Read Stitcher: Reduces noise by stitching paired reads (read 1 and read 2 of the same molecule) in a BAM into a consensus read. The output is a stitched BAM.

[0289] (2) Pisces Variant Caller: Calls small SNVs, insertions, and deletions. Pisces includes a variant-collapsing algorithm into coalescent variants that are decomposed by read boundaries, a basic filtering algorithm, and a simple Poisson-based variant confidence scoring algorithm. The output is a VCF.

[0290] (3) Pisces Variant Quality Recalibrator (VQR): If variant calls predominantly follow patterns associated with heat damage or FFPE deamination, the VQR step lowers the variant Q score of suspicious variant calls. The output is an adjusted VCF.

[0291] (4) Pisces Variant Phaser (Scylla): Uses a read-backed greedy clustering method to assemble small numbers of variants from clonal subpopulations into complex alleles, allowing for more accurate determination of functional consequences by downstream tools. The output is an adjusted VCF.

[0292] Additionally or alternatively, this operation may utilize the variant calling application Strelka™ application by Illumina Inc., hosted at https: / / github.com / Illumina / strelka and described in the article T Saunders, Christopher and Wong, Wendy and Swamy, Sajani and Becq, Jennifer and J Murray, Lisa and Cheetham, Keira, (2012), Strelka: Accurate somatic small-variant calling from sequenced tumor-normal sample pairs, Bioinformatics (Oxford, England), 28, 1811-7, 10.1093 / bioinformatics / bts271, the full subject matter of which is expressly incorporated herein by reference in its entirety. Further, in addition, or instead, this operation may utilize the variant calling application Strelka2™ application by Illumina Inc., hosted at https: / / github.com / Illumina / strelka and described in the editorial Kim, S, Scheffler, K, Halpern, AL, Bekritsky, MA, Noh, E, Kallberg, M, Chen, X, Beyter, D, Krusche, P, and Saunders, CT, (2017), Strelka2: Fast and accurate variant calling for clinical sequencing applications, the full subject matter of which is expressly incorporated herein by reference in its entirety.Furthermore, in addition, or instead, this operation may utilize variant annotation / calling tools such as the Nirvana™ application by Illumina Inc., hosted at https: / / github.com / Illumina / Nirvana / wiki and described in the article Stromberg, Michael and Roy, Rajat and Lajugie, Julien and Jiang, Yu and Li, Haochen and Margulies, Elliott, (2017), Nirvana: Clinical Grade Variant Annotator, 596-596, 10.1145 / 3107411.3108204, the full subject matter of which is expressly incorporated herein by reference in its entirety.

[0293] Such variant annotation / calling tools can apply different algorithmic techniques, such as those disclosed in Nirvana.

[0294] a. Identifying all overlapping transcripts using an interval array. For functional annotation, we can identify all transcripts that overlap with a variant and use an interval tree. However, because the set of intervals can be static, it was possible to further optimize the interval tree 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, k is really small compared to n for most variants, so the effective runtime on an interval tree is O(k lg n). We improved this to O(lg n+k) by creating an interval array where all intervals are stored in a sorted array, so we only need to find the first overlapping interval and then count up through the remaining (k-1).

[0295] b. CNV / SV (Yu): Annotations for copy number and structural variations may be provided. Similar to annotations for small variants, transcripts overlapping with SVs and previously reported structural variations may also be annotated in online databases. Unlike small variants, not all overlapping transcripts need to be annotated because too many transcripts overlap with large SVs. Instead, all overlapping transcripts belonging to partially overlapping genes may be annotated. Specifically, affected introns, exons, and consequences caused by structural variations may be reported for these transcripts. While an option to allow output of all overlapping transcripts is possible, basic information about these transcripts, such as the gene name and a flag indicating whether they are canonical or partial overlaps with the transcript, may be reported. For each SV / CNV, it is also of interest to know whether these variants have been studied and their frequency in different populations. Therefore, we reported overlapping SVs in external databases such as 1000 genomes, DGV, and ClinGen. To avoid using an arbitrary cutoff for determining which SVs overlap, instead, all overlapping transcripts may be used and the mutual overlap rate calculated, i.e., the overlapping length divided by the shorter of the lengths of these two SVs.

[0296] c. Report Supplementary Annotations. There are two types of supplementary annotations: small variants and structural variants (SVs). SVs may be modeled as intervals, and the interval array discussed above can be used to identify overlapping SVs. Small variants are modeled as points and matched by location and (optionally) allele. They are therefore searched using a binary search-like algorithm. Because the supplementary annotation database can be very large, a much smaller index is created to map chromosome locations to file locations where the supplementary annotations reside. The index is a sorted array of objects (consisting of chromosome locations and file locations) that can be binary searched using the location. To keep the index size small, multiple locations (up to some maximum count) are compressed into a single object that stores only the value of the first location and the difference for subsequent locations. Using a binary search, the runtime is O(lg n), where n is the number of items in the database.

[0297] d.VEP cache file

[0298] e. Transcript Database: The transcript cache (cache) and supplementary database (SAdb) files are serialized dumps of data objects such as 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 the pre-computed interval array and perform lookups on it. Because the cache is loaded into memory and lookups are very fast (as explained above), finding duplicate transcripts is very fast in Nirvana (estimated to be less than 1% of the total runtime?).

[0299] f. Supplementary Database: The data sources for SAdb are listed under Supplementary Materials. For small variants, SAdb is produced by a k-way merge of all data sources, such that each object in the database (identified by reference name and location) holds all associated supplementary annotations. Issues encountered while parsing data source files are documented in detail on the Nirvana homepage. To limit memory usage, only the SA index is loaded into memory. This index allows for fast lookup of file locations for supplementary annotations. However, adding supplementary annotations has been identified as Nirvana's largest bottleneck (approximately 30% of the total runtime), since the data must be fetched from disk.

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

[0301] Such variant annotation tools may involve preprocessing. For example, Nirvana included numerous annotations from external data sources, such as ExAC, EVS, the 1000 Genomes project, dbSNP, ClinVar, Cosmic, DGV, and ClinGen. To fully utilize these databases, information from them must be sanitized. We implemented different strategies to address various inconsistencies present from various data sources. For example, if there are multiple dbSNP entries for the same location and alternative alleles, we add all IDs to a comma-separated list of IDs. If there are multiple entries with different CAF values ​​for the same allele, we use the first CAF value. For conflicting ExAC and EVS entries, we consider the number of sample counts, and the entry with the higher sample count is used. For the 1000 Genomes Project, we removed the allele frequencies of conflicting alleles. Another issue is inaccurate information. We primarily extracted allele frequency information from the 1000 Genome Projects. However, we noticed that for GRCh38, we did not exclude samples with genotypes for which the allele frequency reported in the information field was not available, leading to lower frequencies for variants not available for all samples. To ensure the accuracy of our annotations, we calculated true allele frequencies using all individual-level genotypes. As is known, the same variant may have different representations based on different alignments. To ensure accurate reporting of information about identified variants, we must preprocess variants from different resources to ensure consistent representation. For all external data sources, we trimmed alleles to remove duplicated nucleotides in both the reference and alternative alleles. For ClinVar, we directly parsed the xml file and performed a 5' alignment for all variants. This approach is often used in vcf files. Different databases may contain the same set of information.To avoid unnecessary duplication, we removed some redundant information, for example, we removed variants in the DGV that had the 1000 genome projects as their data source, because we had already reported these variants in the 1000 genomes with more detailed information.

[0302] According to at least some implementations, the variant calling application provides calls for low-frequency variants, germline calling, etc. As a non-limiting example, the variant calling application may be performed on tumor-only samples and / or tumor-normal paired samples. The variant calling application may search for single nucleotide variations (SNVs), multi-nucleotide variations (MNVs), indels, etc. The variant calling application filters discrepancies due to sequencing or sample preparation errors while identifying variants. For each variant, the variant caller identifies the reference sequence, the location of the variant, and the potential variant sequence (e.g., an A-to-C SNV or an AG-to-A deletion). The variant calling application identifies the sample sequence (or sample fragment), the reference sequence / fragment, and a variant call as an indication that the variant is present. The variant calling application may identify raw fragments and specify the raw fragment designation, a count of the number of raw fragments validating the potential variant call, the location within the raw fragment where the supporting variant occurred, and other relevant information. Non-limiting examples of raw fragments include duplex stitched fragments, simplex stitched fragments, duplex un-stitched fragments, and simplex un-stitched fragments.

[0303] The variant calling application may output calls in various formats, such as a .VCF file or a .GVCF file. By way of example only, the variant calling application may be included in a MiSeqReporter pipeline (e.g., when implemented on a MiSeq® sequencer instrument). Optionally, this application may be implemented using various workflows. An analysis may include a single protocol or a combination of protocols that analyze sample reads in a specified manner to obtain desired information.

[0304] The one or more processors then perform validation operations in relation to the potential variant calls. The validation operations may be based on quality scores and / or a tiered testing hierarchy, as described below. When the validation operations confirm or substantiate the potential variant call, the validation operations pass the variant call information (from the variant calling application) to the sample report generator. Alternatively, when the validation operations invalidate or disqualify the potential variant call, the validation operations pass a corresponding indication (e.g., a negative indicator, a no call indicator, an invalid call indicator) to the sample report generator. The validation operations may also pass a confidence score regarding the degree of confidence that the variant call or invalid call assignment is correct.

[0305] The one or more processors then generate and store a sample report. The sample report may include, for example, information about multiple loci for the sample. For example, for each locus of a predetermined set of loci, the sample report may at least one of provide a genotype call, indicate that a genotype call cannot be made, provide a confidence score for the certainty of the genotype call, or indicate potential problems with the assay for one or more loci. The sample report may also indicate the gender of the individual who provided the sample and / or indicate that the sample includes multiple sources. As used herein, a "sample report" may include digital data (e.g., a data file) of a locus or a predetermined set of loci and / or a printed report of a locus or set of loci. Thus, generating or providing may include creating a data file and / or printing the sample report or displaying the sample report.

[0306] The sample report may indicate that a variant call was determined but not validated. When a variant call is determined to be invalid, the sample report may indicate additional information regarding the basis for the determination that the variant call cannot be validated. For example, the additional information in the report may include a description of the raw fragments and the degree to which the raw fragments support or refute the variant call (e.g., counts). Additionally, or alternatively, the additional information in the report may include a quality score obtained according to implementations described herein.

[0307] [Variant Call Application] Implementations disclosed herein include analyzing sequencing data to identify potential variant calls. Variant calling can be performed on stored data from previously performed sequencing operations. Additionally, or alternatively, variant calling can be performed in real time while a sequencing operation is being performed. Each sample read is assigned a corresponding locus. A sample read can be assigned to a corresponding locus based on the sequence of nucleotides in the sample read, or in other words, the order of nucleotides (e.g., A, C, G, T) within the sample read. Based on this analysis, the sample read can be designated as containing a potential variant / allele of a particular locus. The sample read can be collected (or aggregated or pooled) with other sample reads designated as containing a potential variant / allele of a locus. The assignment operation may also be referred to as a calling operation, in which a sample read is identified as potentially associated with a particular genetic location / locus. The sample read can be analyzed to locate one or more identifying sequences of nucleotides (e.g., primer sequences) that distinguish the sample read from other sample reads. More specifically, the discriminating sequence can identify a sample read from other sample reads as being associated with a particular locus.

[0308] The assignment operation may include analyzing a series of n nucleotides of the identification sequence to determine whether the series of n nucleotides of the identification sequence substantially matches one or more of the selected sequences. In certain implementations, the assignment operation may include analyzing the first n nucleotides of the sample sequence to determine whether the first n nucleotides of the sample sequence substantially matches one or more of the selected sequences. The number n may have various values, and this value may be programmed into the protocol or entered by a user. For example, the number n may be defined as the number of nucleotides of the shortest selected sequence in the database. The number n may be a predetermined number. The predetermined number may 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, in other implementations, fewer or more nucleotides may be used. The number n may also be selected by an individual, such as a user of the system. The number n may be based on one or more conditions. For example, the number n can be defined as the number of nucleotides in the shortest primer sequence in the database or a specified number, whichever is smaller. In some implementations, a minimum value for n, such as 15, can be used so that any primer sequence that is less than 15 nucleotides can be designated as an exception.

[0309] In some cases, a series of n nucleotides of the discriminating sequence may not exactly match the nucleotides of the selected sequence. Nevertheless, the discriminating sequence may be substantially identical to the selected sequence if it is nearly identical to the selected sequence. For example, if a series of n nucleotides (e.g., the first n nucleotides) of the discriminating sequence matches the selected sequence without a specified number of mismatches (e.g., 3) and / or a specified number of shifts (e.g., 2), the sample read may be called for the locus. Rules may be established so that each mismatch or shift may be counted as a difference between the sample read and the primer sequence. If the number of differences is less than a specified number, the sample read may be called for the corresponding locus (i.e., assigned to the corresponding locus). In some implementations, a matching score may be determined based on the number of differences between the discriminating sequence of the sample read and the selected sequence associated with the locus. If the matching score passes a specified matching threshold, the locus corresponding to the selected sequence may be designated as a potential locus for the sample read. In some implementations, a subsequent analysis may be performed to determine whether the sample read is called for the locus.

[0310] If the sample read substantially matches (i.e., closely matches, or, as described above, nearly matches) one of the selected sequences in the database, the sample read is assigned or assigned to the locus correlated with the selected sequence. This is sometimes referred to as locus calling or preliminary locus calling, in which the sample read is called for a locus correlated with the selected sequence. However, as discussed above, the sample read may be called for more than one locus. In such implementations, further analysis may be performed to call or assign the sample read to only one of the potential loci. In some implementations, the sample read compared to the database of reference sequences is the first read from paired-end sequencing. When paired-end sequencing is performed, a second read (representing a raw fragment) correlated to the sample read is obtained. After assignment, subsequent analysis performed on the assigned read may be based on the type of locus called for the assigned read.

[0311] Next, the sample reads are analyzed to identify potential variant calls. Among other things, the results of the analysis identify potential variant calls, sample variant frequencies, reference sequences, and the locations where the variants occur within the target genome sequence. For example, if a locus is known to contain a SNP, the assigned reads called for the locus may undergo analysis to identify the SNPs in the assigned reads. If a locus is known to contain a polymorphic repetitive DNA element, the assigned reads may be analyzed to identify or characterize the polymorphic repetitive DNA elements within the sample read. In some implementations, if the assigned read substantially matches both an STR locus and an SNP locus, a warning or flag may be assigned to the sample read. The sample read may be designated as both an STR locus and an SNP locus. This analysis may include aligning the assigned reads according to an assignment protocol to determine the sequence and / or length of the assigned reads. Alignment protocols may include methods described in International Patent Application No. PCT / US2013 / 030867, filed March 15, 2013 (Publication No. WO 2014 / 142831), which is incorporated herein by reference in its entirety.

[0312] The one or more processors then analyze the raw fragments to determine whether a supporting variant is present at the corresponding location within the raw fragment. Various types of raw fragments may be identified. For example, the variant caller may identify a type of raw fragment that exhibits a variant that validates the original variant call. For example, the type of raw fragment may represent a duplex stitched fragment, a simplex stitched fragment, a duplex un-stitched fragment, or a simplex un-stitched fragment. Optionally, instead of or in addition to the foregoing examples, other raw fragments may be identified. In connection with identifying each type of raw fragment, the variant caller also identifies the location within the raw fragment where the supporting variant occurred, as well as a count of the number of raw fragments that exhibited the supporting variant. For example, the variant caller may output an indication that 10 reads of the raw fragment were identified as representing a duplex stitched fragment with a supporting variant at a specific location X. The variant caller may also output an indication that 5 reads of the raw fragment were identified as representing a simplex un-stitched fragment with a supporting variant at a specific location Y. The variant caller may also output the number of raw fragments that did not contain supporting variants that correspond to the reference sequence and thus provide evidence to validate the potential variant call in the genomic sequence of interest.

[0313] Next, a count of raw fragments is maintained, including the supporting variants and the locations where the supporting variants occur. Additionally or alternatively, a count of raw fragments that did not contain a supporting variant at the target location (relative to the location of the potential variant call in the sample read or sample fragment) can be maintained. Additionally or alternatively, a count of raw fragments that correspond to the reference sequence and do not confirm or validate the potential variant call can be maintained. The determined information, including the count and type of raw fragments that support the potential variant call, the location of the supporting variant in the raw fragment, the count of raw fragments that do not support the potential variant call, etc., is output to the variant call validation application.

[0314] When a potential variant call is identified, the process outputs an indication of the potential variant call, the variant sequence, the variant location, and the associated reference sequence. Because errors can cause the calling process to identify an incorrect variant, the variant call is designated to represent a "potential" variant. According to implementations herein, the potential variant call is analyzed to reduce or remove incorrect variants or false positives. Additionally or alternatively, the process analyzes one or more raw fragments associated with the sample read and outputs corresponding variant calls associated with the raw fragments.

[0315] [Generating a good training set] Although millions of human genomes and exons have been sequenced, their clinical application remains limited by the difficulty of distinguishing disease-causing mutations from benign genetic variations. Here, we demonstrate that common missense variants in other primate species are largely clinically benign in humans, allowing pathogenic variants to be systematically identified by a process of elimination. Using hundreds of thousands of common variants from population sequencing of six nonhuman primate species, we trained a deep neural network that identified pathogenic variants in patients with rare diseases with 88% accuracy and enabled the discovery of 14 new candidate genes in intellectual disability with genome-wide significance. Cataloging common variants from additional primate species will improve interpretation of millions of variants of uncertain significance and further advance the clinical application of human genome sequencing.

[0316] The clinical potential of diagnostic sequencing is limited by the difficulty of interpreting rare genetic variants in human populations and inferring their impact on disease risk. Clinically significant genetic variants tend to be extremely rare in the population due to their adverse health effects, and for the majority, their impact on human health has not been determined. The large number and rarity of these variants with uncertain clinical significance poses a formidable barrier to the adoption of sequencing for personalized medicine and population-wide health screening.

[0317] Because most penetrant Mendelian diseases have very low population prevalence, the observation of a variant at high frequency in a population is strong evidence in favor of a benign outcome. While assessing common variation across diverse human populations is a useful strategy for cataloging benign variants, the total number of common variants in modern humans is limited by bottleneck events in the recent history of our species, which resulted in the loss of much of our ancestral diversity. Studies of modern human populations have focused on effective population sizes (N) of less than 10,000 individuals within the past 15,000–65,000 years. e ), and the small pool of common polymorphisms stems from the limited mutational capacity in a population of this size: of the 70 million potential protein-altering missense substitutions in the reference genome, roughly only 1 in 1000 has an overall population allele frequency above 0.1%.

[0318] Outside of modern human populations, chimpanzees constitute the next closest extant species, sharing 99.4% amino acid sequence identity. The close homology of protein-coding sequences in humans and chimpanzees suggests that purifying selection acting on protein-coding variants in chimpanzees may also model the health consequences of mutations in the same state in humans.

[0319] Neutral polymorphisms contribute to the ancestral lineage of humans (approximately 4N e Because the average time that a chimpanzee persists (in generations) is a fraction of the species divergence time (approximately 6 million years ago), naturally occurring chimpanzee mutations span a mutation space that largely does not overlap except by chance, except for rare instances of haplotypes maintained by balancing selection. If identical polymorphisms affect fitness similarly in the two species, the presence of a variant at high allele frequency in the chimpanzee population should indicate a benign outcome in humans, expanding the catalog of known variants whose benign outcome has been established by purifying selection.

[0320] [Results - Common variants in other primates are mostly benign in humans] The recent availability of aggregated exon data, including 123,136 individuals collected in the Exome Aggregation Consortium (ExAC) and Genome Aggregation Database (gnomAD), allows us to measure the influence of natural selection on missense and synonymous mutations across the allele frequency spectrum. Rare singleton variants observed only once in a cohort closely match the expected missense / synonymous ratio of 2.2 / 1 predicted by de novo mutations after adjusting for the effect of trinucleotide context on mutation rate (Figures 49A, 51, and 52A, 52B, 52C, and 52D). However, at higher allele frequencies, the number of observed missense variants decreases due to the purging of deleterious mutations by natural selection. The gradual decline in the missense / synonymous ratio with increasing allele frequency is consistent with the fact that a significant portion of missense variants with population frequencies <0.1% have mildly deleterious consequences despite being observed in healthy individuals. These findings support the widespread empirical practice in clinics of eliminating variants with allele frequencies higher than 0.1% to approximately 1% as likely benign for penetrant genetic disorders, with a few well-documented exceptions caused by balancing selection and founder effects.

[0321] We identified common chimpanzee variants sampled more than once in a cohort of 24 unrelated individuals. We estimated that 99.8% of these variants are common in the general chimpanzee population (allele frequency (AF) >0.1%), indicating that these variants have already been subjected to purifying selection. We examined the human allele frequency spectrum for corresponding cognate human variants, excluding extended major histocompatibility complex regions as known regions of balancing selection, along with variants lacking one-to-one mapping in multiple sequence alignments (Figure 49B). For human variants cognate to common chimpanzee variants, the missense / synonymous ratio was largely constant across the human allele frequency spectrum (P > 0.5 by chi-square (χ2) test), consistent with the absence of negative selection against common chimpanzee variants in the human population and the matching selection coefficients for missense variants in the two species. The low missense / synonymous ratio observed in human variants identical to common chimpanzee variants is due to the larger effective population size (N e Approximately 73,000), which allows for more efficient removal of mildly deleterious mutations.

[0322] In contrast, for singleton chimpanzee variants (sampled only once in the cohort), a significant decrease in the missense / synonymous ratio was observed in the common allele frequency (P < 5.8 × 10 -6 , Figure 49C), which indicates that 24% of singleton chimpanzee missense variants would be filtered out by purifying selection in human populations at allele frequencies higher than 0.1%. This depletion indicates that the majority of chimpanzee singleton variants are rare deleterious mutations whose deleterious effects on health prevent them from reaching common allele frequencies in either species. We estimate that only 69% of singleton variants are common in the general chimpanzee population (AF > 0.1%).

[0323] We then identified human variants identical to mutations observed in at least one of six nonhuman primate species. The mutations in each of the six species were identified from the Great Ape Genome Project (chimpanzee, bonobo, gorilla, and orangutan) or submitted to the Single Nucleotide Polymorphism Database (dbSNP) from the Primate Genome Project (rhesus macaque, marmoset), and represent largely common variants based on the limited number of sequenced individuals and the low missense:synonymous ratio observed for each species (Supplementary Table 1). Similar to chimpanzees, we found that the missense / synonymous ratios of variants from six nonhuman primate species were roughly equivalent across the human allele frequency spectrum, except for a mild depletion of missense mutations at common allele frequencies (Figure 49D, Figure 53, Figure 54, and Figure 55, Supplementary Data File 1) expected due to the inclusion of a small proportion of rare variants (<0.1% allele frequency in chimpanzees and approximately 16% at lower allele frequencies due to fewer sequenced individuals in the other species). These results suggest that selection coefficients for identical missense variants are consistent within the primate lineage, at least up to New World monkeys, which are estimated to have diverged from the ancestral human lineage approximately 35 million years ago.

[0324] We found that human missense variants identical to observed primate variants were strongly enriched for benign outcomes in the ClinVar database. After excluding variants of uncertain significance and variants with conflicting annotations, ClinVar variants present in at least one nonhuman primate species were annotated as benign or likely benign in an average of 90% of cases, compared with 35% for ClinVar missense variants overall (P < 10 -40, Figure 49E). The pathogenicity of the ClinVar annotations for primate variants is slightly higher than that observed from sampling a similar-sized cohort of healthy humans, excluding human variants with allele frequencies higher than 1% to reduce curation bias (approximately 95% benign or likely benign results, P = 0.07).

[0325] The field of human genetics has long relied on model organisms to infer the clinical impact of human mutations, but the long evolutionary distance to most genetically tractable animal models raises concerns about the generalizability of model organism findings to humans. We extended our analysis beyond the primate lineage to include largely common mutations from four additional mammalian species (mouse, pig, goat, and cow) and two more distantly related vertebrates (chicken and zebrafish). We selected species with sufficient genome-wide mutation identification in dbSNP and confirmed these as largely common variants based on a missense / synonymous ratio much less than 2.2 / 1. In contrast to our primate analysis, human missense mutations identical to mutations in more distantly related species are significantly depleted in common allele frequency (Figure 50A), and the degree of this depletion increases at longer evolutionary distances (Figure 50B and Supplementary Tables 2 and 3).

[0326] Missense mutations that are deleterious in humans but tolerate high allele frequencies in more distantly related species indicate that the coefficients of selection for identical missense mutations differ considerably between humans and more distantly related species. Nevertheless, the presence of missense variants in more distantly related mammals increases the probability of a benign outcome, because the proportion of missense variants depleted by natural selection at common allele frequencies is lower than the approximately 50% depletion rate observed for human missense variants overall (Figure 49A). Consistent with these results, ClinVar missense variants observed in mice, pigs, goats, and cattle have a 73% probability of being annotated with a benign or likely benign outcome, compared with 90% for primate variants (P<2×10). -8 , Figure 50C), and found that for the entire ClinVar database it was 35%.

[0327] To confirm that evolutionary distance, rather than domestication artifacts, was the primary cause of differences in selection coefficients, we repeated the analysis using fixed substitutions between pairs of closely related species instead of intraspecific polymorphisms across a wide range of evolutionary distances (Figure 50D, Supplementary Table 4, and Supplementary Data File 2). We found that the depletion rate of human missense variants identical to interspecific fixed substitutions increased with evolutionary branch length, with no discernible difference for wild species compared with domesticated species. This is consistent with previous work in flies and yeast, which found that the number of identical fixed missense substitutions was lower than expected by chance in divergent lineages.

[0328] [Deep learning network for variant pathogenicity classification] The disclosed technology provides deep learning networks for pathogenicity classification of variants. The importance of variant classification to clinical applications has prompted many attempts to use supervised machine learning to address the problem, but these efforts have been hampered by the lack of an adequately sized truth dataset containing confidently labeled benign and pathogenic variants for training.

[0329] Existing databases of expert-curated variants are not genome-wide representative, with approximately 50% of the variants in the ClinVar database coming from only 200 genes (approximately 1% of human protein-coding genes). Furthermore, systematic studies have identified questionable evidence for many expert annotations, highlighting the difficulty of interpreting rare variants that may only be observed in a single patient. While expert interpretations are becoming increasingly rigorous, classification guidelines are largely developed around agreed-upon conventions, risking reinforcing existing trends. To reduce bias in human interpretation, classifiers have recently been trained on common human polymorphisms or fixed human-chimpanzee substitutions, but these also use as input the prediction scores of previous classifiers trained on human-curated databases. Objective benchmarking of the performance of these various methods is difficult to achieve without an independent, unbiased, truth dataset.

[0330] Variation from six nonhuman primate species (chimpanzee, bonobo, gorilla, orangutan, rhesus macaque, and marmoset) contributes over 300,000 unique missense variants that do not overlap with common human variants, representing mostly benign common variants that have been filtered out by purifying selection, greatly expanding the training dataset available for machine learning methods. On average, each primate species contributes more variants than the entire ClinVar database (approximately 42,000 missense variants as of November 2017, after excluding variants of uncertain significance and variants with conflicting annotations). Additionally, this content is free from the bias of human interpretation.

[0331] Using a dataset with common human variants (AF > 0.1%) and primate mutations (Supplementary Table 5 (Figure 5-8)), we trained a novel deep residual network, PrimateAI. PrimateAI takes as input the amino acid sequence flanking the variant of interest and orthologous sequence alignments in other species (Figures 2 and 3). Unlike existing classifiers that utilize human-curated features, our deep learning network learns to extract features directly from the original sequence. To incorporate information about protein structure, we trained separate networks to predict secondary structure and solvent accessibility from sequence alone and then included these as subnetworks in the full model (Figures 5 and 6). Given the small number of human proteins that have been successfully crystallized, inferring structure from the original sequence has the advantage of avoiding bias due to incomplete protein structure and functional domain annotations. The total depth of the network was 36 layers of convolution with approximately 400,000 trainable parameters, along with the included protein structures.

[0332] To train a classifier using only variants with benign labels, We formulated the prediction problem as whether a given mutation is likely to be observed as a common variant in a population. Several factors influence the probability of observing variants at high allele frequencies, and we are only interested in the deleteriousness of those variants. Other factors include mutation rate, technical artifacts such as sequencing coverage, and factors affecting neutral genetic drift, such as gene conversion.

[0333] We considered each of these indistinguishable factors, matched each variant in the benign training set with missense mutations absent from 123,136 exons from the ExAC database, and trained a deep learning network to distinguish between benign variants and matched controls (Figure 24). Because the number of unlabeled variants greatly exceeds the size of the labeled benign training dataset, we trained eight networks in parallel, each using a different set of unlabeled variants matched against the benign training dataset, to obtain consensus predictions.

[0334] Using only the original amino acid sequence as input, deep learning networks accurately assign high pathogenicity scores to residues in informative protein functional domains, as shown for the voltage-gated sodium channel SCN2A, a major disease gene in epilepsy, autism, and intellectual disability (Figure 20). The structure of SCN2A comprises four homologous repeats, each containing six transmembrane helices (S1–S6). Membrane depolarization causes the positively charged S4 transmembrane helix to move toward the extracellular side of the membrane, opening the S5 / S6 pore-forming domain via the S4-S5 linker. Mutations in the S4, S4-S5 linker, and S5 domain, which are clinically associated with early signs of epileptic encephalopathy, are predicted by the network as having the highest pathogenicity scores in the gene and are depleted relative to variants in healthy populations (Supplementary Table 6). We also found that the network recognizes the locations of critical amino acids within domains and assigns the highest pathogenicity scores to mutations in these locations, such as DNA-contacting residues in transcription factors and catalytic residues in enzymes (Figures 25A, 25B, 25C, and 26).

[0335] To better understand how well the deep learning network derives insights into protein structure and function from the original sequence, we visualized the trainable parameters from the first three layers of the network. Within these layers, we observed that the network learns correlations between the weights of different amino acids that approximate existing measures of amino acid distance, such as the Grantham score (Figure 27). The outputs of these early layers become the inputs for later layers, allowing the deep learning network to build progressively higher-order representations of the data.

[0336] We compared the performance of our network with existing classification algorithms using 10,000 common primate variants withheld from training. Because approximately 50% of all newly emerged human missense variants are filtered out by purifying selection at common allele frequencies (Figure 49A), we calculated the 50th percentile score for each classifier on a set of 10,000 randomly selected variants matched by mutation rate and sequencing coverage to the 10,000 common primate variants and assessed the accuracy of each classifier at that threshold (Figure 21D, Figure 28A, and Supplementary Data File 4). Our deep learning network (91% accuracy) outperformed other classifiers (the next best model was 80% accurate) in assigning benign outcomes to the 10,000 reserved common primate variants.

[0337] When compared to the accuracy of a network trained using only human mutation data (Figure 21D), roughly half of the improvement over existing methods comes from using a deep learning network, and half comes from augmenting the training dataset with primate mutations. To test the classification of variants of uncertain significance in a clinical scenario, we evaluated the deep learning network's ability to distinguish de novo mutations occurring in patients with neurodevelopmental disorders versus healthy controls. By prevalence, neurodevelopmental disorders constitute one of the largest categories of rare genetic diseases, and recent trio-sequencing studies suggest a central role for de novo missense and protein-truncating mutations.

[0338] We classified de novo missense variants in 4293 affected individuals from the Deciphering Developmental Disorders (DDD) cohort versus 2517 unaffected siblings from the Simon's Simplex Collection (SSC) cohort, each confidently called, and assessed the difference in prediction scores between the two distributions using a Wilcoxon rank-sum test (Figure 21E and Figures 29A and 29B). The deep learning network clearly outperformed the other classifiers for this task (P<10). -28 , Figure 21F and Figure 28B). Furthermore, the performance of various classifiers on the reserved primate variant dataset and the DDD case vs. control dataset was correlated (Spearman ρ = 0.57, P < 0.01), showing good agreement between the two datasets for assessing pathogenicity, despite using completely different sources and methods (Figure 30A).

[0339] Next, we attempt to estimate the accuracy of the deep learning network for classifying benign versus pathogenic mutations within the same gene. Given that the majority of DDD populations comprise an index case of an affected child without affected first-degree relatives, it is important to show that the classifier is not inflating accuracy by overestimating the pathogenicity of genes with a de novo dominant mode of inheritance. We constrained the analysis to 605 genes that were nominally significant for disease association in DDD studies, calculated from protein-truncating mutations only (P < 0.05). Within these genes, de novo missense mutations were enriched 3 / 1 compared to expectation (Figure 22A), indicating that approximately 67% are pathogenic.

[0340] The deep learning network was able to distinguish between pathogenic and benign de novo variants within the same set of genes (P < 10 -15, Figure 22B), significantly outperforming other methods (Figure 22C and Figure 28C). At a binary cutoff of 0.803 or greater (Figure 22D and Figure 30B), the deep learning network classified 65% of de novo missense mutations in the case group as pathogenic, compared with 14% of de novo missense mutations in the control group, corresponding to a classification accuracy of 88% (Figure 22E and Figure 30C). Given the frequent incomplete penetrance and variable expressivity in neurodevelopmental disorders, this figure likely underestimates the accuracy of our classifier due to the inclusion of partially penetrant pathogenic variants in the control group.

[0341] [Discovery of new gene candidates] Applying a threshold of 0.803 or higher to stratify pathogenic missense variants increases the enrichment of de novo missense variants in DDD patients from 1.5-fold to 2.2-fold, closer to protein-truncating mutations (2.5-fold), while discarding less than one-third of the total number of variants enriched beyond expectation. This significantly increases statistical power and allows the discovery of 14 additional candidate genes in intellectual disability that previously did not reach the genome-wide significance threshold in the original DDD study (Table 1).

[0342] [Compare with expert selections] We investigated the performance of various classifiers on recent expert-curated variants from the ClinVar database and found that classifier performance on the ClinVar dataset was not strongly correlated with either the withheld primate variant dataset or the DDD case vs. control dataset (P = 0.12 and P = 0.34, respectively) (Figure 31A and Figure 31B). We hypothesize that existing classifiers are biased by expert curation, and human heuristics, while tending in the right direction, can be suboptimal. One example is the difference in mean Grantham scores between pathogenic and benign variants in ClinVar, which is twice the difference between de novo variants in DDD cases vs. controls within 605 disease-associated genes (Table 2). In comparison, expert curation appears to under-exploit the importance of protein structure, particularly surface-exposed residues that may be able to interact with other molecules. We observed that both ClinVar pathogenic variants and DDD de novo mutations were associated with predicted solvent-exposed residues, but the difference in solvent accessibility between benign and pathogenic ClinVar variants was only half of the difference seen for DDD cases versus controls. These findings suggest an ascertainment bias favoring factors that are easier for experts to interpret, such as Grantham scores and conservation rates. Machine learning classifiers trained on human-curated databases are expected to reinforce these trends.

[0343] Our results suggest that systematic sequencing of primate populations is an effective strategy for classifying millions of human variants of uncertain significance that currently constrain clinical genome interpretation. The accuracy of our deep learning network for both reserved common primate variants and clinical variants increases with the number of benign variants used to train the network (Figure 23A). Furthermore, training on variants from each of six nonhuman primate species independently contributes to improved network performance, whereas training on variants from more distantly related mammals negatively impacts network performance (Figures 23B and 23C). These results support the assertion that while common primate variants are largely benign in humans with regard to penetrant Mendelian diseases, the same cannot be said for mutations in more distantly related species.

[0344] Although the number of nonhuman primate genomes surveyed in this study is small compared to the number of human genomes and exons that have been sequenced, it is important to note that these additional primates contribute a disproportionate amount of information about common benign mutations. Simulations using ExAC indicate that while the discovery of common human variants (>0.1% allele frequency) quickly plateaus after only a few hundred individuals (Figure 56), sequencing additional healthy populations of up to several million individuals will primarily contribute additional rare variants. Unlike common variants known to be largely clinically benign based on allele frequency, rare variants in healthy populations can cause recessive genetic disorders or dominant genetic disorders with incomplete penetrance. Because each primate species has a different pool of common variants, sequencing dozens of individuals of each species is an effective strategy for systematically cataloging benign missense mutations in the primate lineage. Indeed, the 134 individuals from the six nonhuman primate species surveyed in this study contributed nearly four times as many common missense mutations as the 123,136 humans from the ExAC study (Supplementary Table 5). Sequencing studies of primate populations involving hundreds of individuals may be feasible with relatively few unrelated individuals living in wild sanctuaries and zoos, thereby minimizing disturbance to wild populations, which is important from the perspective of conservation and ethical treatment of nonhuman primates.

[0345] Modern human populations are much less genetically diverse than most nonhuman primate species, with roughly half the number of single-nucleotide variants per individual compared to chimpanzees, gorillas, and gibbons, and one-third the number of variants per individual compared to orangutans. While the level of genetic diversity for most nonhuman primate species is unknown, the large number of extant nonhuman primate species allows us to estimate that the majority of potential benign human missense sites are likely covered by common variants in at least one primate species, allowing pathogenic variants to be systematically identified by a process of elimination (Figure 23D). Even with only a subset of these species sequenced, increasing the training data size will enable more accurate prediction of missense consequences using machine learning. Finally, although our findings focus on missense mutations, this strategy may also be applicable to inferring the consequences of noncoding mutations, particularly in conserved control regions where there is sufficient alignment between the human and primate genomes to unambiguously determine whether a variant is identical.

[0346] Of the 504 known non-human primate species, approximately 60% are threatened with extinction due to poaching and large-scale habitat loss. Population declines and possible extinction in these species represent an irreplaceable loss of genetic diversity, motivating urgent global conservation efforts that would benefit both these unique and irreplaceable species and ourselves.

[0347] [Data generation and alignment] Coordinates in the application refer to the human genome build UCSC hg19 / GRCh37, which includes coordinates for variants in other species that map to hg19 using multiple sequence alignments. Canonical transcript to protein-coding DNA sequences and multiple sequence alignments and branch lengths of 99 vertebrate genomes were downloaded from the UCSC Genome Browser.

[0348] Human exon polymorphism data were obtained from the Exome Aggregation Consortium (ExAC) / Genome Aggregation Database (gnomAD exomes) v2.0. Primate variation data were obtained from the Great Ape Genome Sequencing Project, which includes whole genome sequencing data and genotypes for 24 chimpanzees, 13 bonobos, 27 gorillas, and 10 orangutans. Variations from 35 chimpanzees from a separate chimpanzee and bonobo study were also included, but due to differences in variant calling methods, these were excluded from population analyses and used only for training deep learning models. In addition, 16 rhesus macaque individuals and 9 marmoset individuals were used to evaluate variation in the original genome projects for these species, but individual-level information was not available. Variation data for rhesus macaques, marmosets, pigs, cows, goats, mice, chickens, and zebrafish were obtained from dbSNP. dbSNP also contains additional orangutan variants, which we used only for training our deep learning model because individual genotype information was not available for the population analysis. To avoid the influence of balancing selection, we also excluded variants from within the extended major histocompatibility complex region (chr6:28,477,797-33,448,354) for the population analysis.

[0349] A multi-species alignment of 99 vertebrate species was used to ensure orthologous one-to-one mapping to human protein-coding regions and to prevent mapping to pseudogenes. Variants were accepted as identical if they occurred in either the canonical or alternate orientation. To ensure that variants have the same predicted protein-coding outcome in both humans and other species, we required that the other two nucleotides in the codon be identical between species for both missense and synonymous variants. Polymorphisms from each species included in the analysis are listed in Supplementary Data File 1, and detailed scales are provided in Supplementary Table 1.

[0350] For each of the four allele frequency categories (Figure 49A), we used mutations in intronic regions to estimate the expected number of synonymous and missense variants in each of the 96 potential trinucleotide contexts and to correct for mutation rates (Figure 51 and Supplementary Tables 7 and 8 (Figure 59)). We also analyzed identical CpG dinucleotide and non-CpG dinucleotide variants separately and verified that the missense / synonymous ratio was flat across the allele frequency spectrum for both classes. This indicates that our analysis is applicable to both CpG and non-CpG variants, despite the large differences in their mutation rates (Figures 52A, 52B, 52C, and 52D).

[0351] [Depletion rates of human missense variants identical to polymorphisms in other species] To assess whether variants present in other species are tolerant at common allele frequencies (>0.1%) in humans, we identified human variants that were identical to mutations 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%, 0.01%–0.1%, or >0.1%) and estimated the missense-to-synonymous ratio (MSR) decline between rare (<0.1%) and common (>0.1%) variants. The depletion rate of missense variants identical to the common allele frequency (>0.1%) in humans indicates the proportion of variants from other species that are sufficiently deleterious to be eliminated by natural selection at common allele frequencies in humans.

[0352]

number

[0353] The missense / synonymous ratio and depletion rate were calculated for each species and are shown in Figure 50B and Supplementary Table 2. In addition, a chi-squared test of homology (χ ) was performed on the 2 × 2 contingency tables to test whether the difference in the missense / synonymous ratio between rare and common variants was significant for chimpanzee common variants (Figure 49B), chimpanzee singleton variants (Figure 49C), and mammalian variants (Figure 50A). 2 ) was carried out.

[0354] Because sequencing was performed on only a limited number of individuals from the Great Ape Genome Project, we used the human allele frequency spectrum from ExAC to estimate the proportion of sampled variants that were rare (<0.1%) or common (>0.1%) in the general chimpanzee population. We sampled a cohort of 24 humans based on ExAC allele frequencies and identified missense variants that were observed once or more than once in this cohort. Variants observed more than once had a 99.8% probability of being common in the general population (>0.1%), whereas variants observed only once in the cohort had a 69% probability of being common in the general population. To verify that the observed depletion of missense variants in more distantly related mammals was not due to the indistinguishable influence of well-conserved and therefore more precisely aligned genes, we repeated the above analysis, restricting it to only genes with an average nucleotide identity of >50% in multiple sequence alignments of 11 primates and 50 mammals compared to humans (see Supplementary Table 3).

[0355] This removed approximately 7% of human protein-coding genes from the analysis without significantly affecting the results. Additionally, to ensure that the results were not affected by issues with variant calling or domestication artifacts (as the majority of species selected from dbSNP are domesticated), we repeated the analysis using fixed substitutions from pairs of closely related species instead of interspecific polymorphisms (Figure 50D, Supplementary Table 4, and Supplementary Data File 2).

[0356] [ClinVar analysis of polymorphism data for humans, primates, mammals, and other vertebrates] To investigate the clinical impact of variants with identical status to other species, we downloaded the ClinVar database, excluding variants that had conflicting pathogenic annotations or were labeled only as variants of uncertain significance. After the filtering step shown in Supplementary Table 9, there are a total of 24,853 missense variants in the pathogenic category and 17,775 missense variants in the benign category.

[0357] We counted the number of pathogenic and benign ClinVar variants that were identical to mutations in humans, nonhuman primates, mammals, and other vertebrates. For humans, we simulated a cohort of 30 humans sampled from ExAC allele frequencies. The number of benign and pathogenic variants for each species is shown in Supplementary Table 10.

[0358] [Generating benign and unlabeled variants for model training] For machine learning, we constructed a benign training dataset of mostly common benign missense variants from humans and nonhuman primates. This dataset includes common human variants (>0.1% allele frequency, 83,546 variants) as well as variants from chimpanzees, bonobos, gorillas, orangutans, rhesus macaques, and marmosets (301,690 unique primate variants). The number of benign training variants contributed by each source is shown in Supplementary Table 5.

[0359] A deep learning network was trained to distinguish between a set of collated labeled benign variants and an unlabeled set of variants, taking into account trinucleotide context, sequencing coverage, and their likelihood of alignment between species and humans. To obtain the unlabeled training dataset, we started with all potential missense variants in canonical coding regions. We excluded variants observed in 123,136 exons from ExAC, as well as variants in start or stop codons. In total, 68,258,623 unlabeled missense variants were generated. This was then filtered to correct for regions with poor sequencing coverage and regions that did not have a one-to-one alignment between the human and primate genomes when selecting collated unlabeled variants for primate variants.

[0360] Eight models using the same set of labeled benign variants and eight randomly sampled sets of unlabeled variants were trained and averaged to obtain consensus predictions. For validation and testing, two randomly sampled sets of 10,000 primate variants were also excluded and withheld from training (Supplementary Data File 3). For each of these sets, 10,000 unlabeled variants matched by trinucleotide context were sampled and used to normalize the thresholds of each classifier when comparing different classification algorithms (Supplementary Data File 4). In other implementations, fewer or additional models, ranging from 2 to 500, may be used in the ensemble.

[0361] We evaluated the classification accuracy of two versions of the deep learning network: one trained using only common human variants, and the other trained using a complete benign labelled dataset containing both common human and primate variants.

[0362] [Deep learning network architecture] For each variant, the pathogenicity prediction network takes as input the 51-amino acid sequence centered on the variant of interest and the outputs of the secondary structure and solvent accessibility networks (Figures 2 and 3), along with the missense variant substituted at the center location. Three 51-length location frequency matrices are generated from the multiple sequence alignments of 99 vertebrates, including one location frequency matrix for 11 primates, one location frequency matrix for 50 mammals excluding primates, and one location frequency matrix for 38 vertebrates excluding primates and mammals.

[0363] The secondary structure deep learning network predicts three states of secondary structure at each amino acid location: α-helix (H), β-sheet (B), and coil (C) (Supplementary Table 11). The solvent accessibility network predicts three states of solvent accessibility at each amino acid location: buried (B), intermediate (I), and exposed (E) (Supplementary Table 12). Both networks took only flanking amino acid sequences as input and were trained using labels from known, non-redundant crystal structures in the Protein DataBank (Supplementary Table 13). A single length-location frequency matrix, also 51 in length and 20 in depth, generated from multiple sequence alignments for all 99 vertebrate species was used as input to the pretrained three-state secondary structure network and three-state solvent accessibility network. After pretraining the network on known crystal structures from the Protein DataBank, the final two layers for the secondary structure and solvent models were removed, and the network output was directly connected to the input of the pathogenicity model. The best test accuracy achieved for the three-state secondary structure prediction model was 79.86% (Supplementary Table 14). Comparing neural network predictions using DSSP (Define Secondary Structure of Proteins) and annotated structural labels for approximately 4000 human proteins that had crystal structures, there was no significant difference when using only predicted structural labels (Supplementary Table 15).

[0364] Both our deep learning network for pathogenicity prediction (PrimateAI) and our deep learning network for predicting secondary structure and solvent accessibility employed a residue-block architecture. The detailed architecture of PrimateAI is described in Figure 3 and Supplementary Table 16 (Figure 4A, Figure 4B, and Figure 4C). The detailed architecture of the network for predicting secondary structure and solvent accessibility is described in Figure 6 and Supplementary Tables 11 (Figure 7A and Figure 7B) and 12 (Figure 8A and Figure 8B).

[0365] [Benchmarking classifier performance on a held-out test set of 10,000 primate variants] We used 10,000 reserved primate variants in the test dataset to benchmark our deep learning network and 20 other previously published classifiers, whose predicted scores were obtained from the database dbNSFP. The performance of each classifier on the 10,000 reserved primate variants test set is also given in Figure 28A. Because different classifiers had highly variable score distributions, we used 10,000 randomly selected unlabeled variants matched to the test set by trinucleotide context to identify the 50th percentile threshold for each classifier. To ensure a fair comparison between methods, we benchmarked each classifier on the proportion of variants in the 10,000 reserved primate variants test set that were classified as benign at the 50th percentile threshold for that classifier.

[0366] For each of the classifiers, the proportion of retained primate test variants predicted as benign using the 50th percentile threshold is also shown in Figure 28A and Supplementary Table 17 (Figure 34). The performance of PrimateAI is robust with respect to the number of aligned species at the variant location, and it also shows good overall performance as long as sufficient conservation information from mammals is available, which is the case for most protein-coding sequences (Figure 57).

[0367] [Analysis of de novo variants from the DDD study] We retrieved published de novo variants from the DDD study and de novo variants from the healthy sibling control group in the SSC autism study. The DDD study provided a confidence level for de novo variants, and we excluded variants with a threshold below 0.1 from the DDD dataset as potential false positives due to variant calling errors. In one implementation, there were a total of 3512 missense de novo variants from DDD-affected individuals and 1208 missense de novo variants from the healthy control group. The canonical transcript annotation used by UCSC for multiple sequence alignments of 99 vertebrate species differed slightly from the transcript annotation used by DDD, resulting in small differences in the total number of missense variants. We evaluated the ability of this classification method to distinguish between de novo missense variants in DDD-affected individuals and those in unaffected sibling controls from the autism study. For each classifier, we reported the P-value from the Wilcoxon rank sum test of the difference between the predicted scores for the two distributions (Supplementary Table 17 (Figure 34)).

[0368] To measure the accuracy of various classifiers in distinguishing between benign and pathogenic variants within the same disease gene, we repeated the analysis on a subset of 605 genes enriched for de novo protein-truncating mutations in the DDD cohort (P < 0.05, Poisson exact test) (Supplementary Table 18). Based on a greater-than-expected enrichment of de novo missense mutations within these 605 genes, we estimated that two-thirds of de novo variants in the DDD dataset were pathogenic and one-third were benign. We assumed minimal incomplete penetrance and that de novo missense mutations in healthy controls were benign. For each classifier, we identified a threshold that produced the same number of benign or pathogenic predictions as the empirical proportion observed in these datasets, and used this threshold as a binary cutoff to estimate each classifier's accuracy in distinguishing between de novo mutations in cases and controls. To construct receiver operating characteristic curves, we treated the classification of de novo DDD variants as pathogenic as a true-positive call and the classification of de novo variants as pathogenic in healthy controls as a false-positive call. Because the DDD dataset contains one-third benign de novo variants, the area under the curve (AUC) for a theoretically perfect classifier is less than 1. Thus, a classifier that perfectly separates 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 positives, and 100% of de novo variants in controls as true negatives, yielding the largest possible AUC of 0.837 (Figure 29A and Figure 29B and Supplementary Table 19 (Figure 35)).

[0369] [Discovery of new gene candidates] Gene enrichment for de novo mutations was tested by comparing the number of observed de novo mutations to the number expected under a null mutation model. We repeated the enrichment analysis performed in the DDD study and reported genes that were newly genome-wide significant when counting only de novo missense mutations with a PrimateAI score above 0.803. We adjusted the genome-wide expectation for de novo damaging missense mutations by the proportion of missense variants meeting a PrimateAI threshold above 0.803 (roughly one-fifth of all potential missense mutations genome-wide). For each DDD study, each gene required four tests: one for protein-cleaving enrichment and one for enrichment for protein-altering de novo mutations; both were tested for the DDD cohort alone and for the larger meta-analysis of the neurodevelopmental trio sequencing cohort. Enrichment of protein-altering de novo mutations was combined with a test for clustering of missense de novo mutations within coding sequences by Fisher's method (Supplementary Tables 20 and 21). The P value for each gene was taken from the minimum of the four tests, and genome-wide significance was determined to be P < 6.757 × 10 -7 was determined as (α = 0.05, 18,500 genes using four tests).

[0370] [Accuracy of ClinVar classification] Because the majority of existing classifiers are trained either directly or indirectly on ClinVar content, such as by using predicted scores from classifiers trained on ClinVar, we limited our analysis of the ClinVar dataset to only use ClinVar variants added since 2017. Because there was considerable overlap between recent ClinVar variants and other databases, we further filtered to remove variants found at common allele frequencies (>0.1%) in ExAC or present in the Human Gene Mutation Database (HGMD), Leiden Open Variation Database (LOVD), or Universal Protein Resource (Uniprot). After removing variants annotated only as of uncertain significance and variants with conflicting annotations, 177 missense variants with benign annotations and 969 missense variants with pathogenic annotations remained. These ClinVar variants were scored using both deep learning networks and other classification methods. For each classifier, a threshold that produced the same number of benign and pathogenic predictions was identified as the empirical proportion observed in these databases, and this threshold was used as a binary cutoff to estimate the accuracy of each classifier (Figures 31A and 31B).

[0371] The impact of increasing training data size and using different sources of training data. To evaluate the impact of training data size on the performance of deep learning networks, a subset of variants was randomly sampled from a training set of 385,236 primate and common human variants labeled as benign, while the underlying deep learning network architecture remained the same. To demonstrate that variants from each individual primate species contribute to classification accuracy, while variants from each individual mammalian species contribute to lower classification accuracy, according to one implementation, a deep learning network was trained using a training dataset comprising 83,546 human variants and a randomly selected fixed number of variants for each species, again keeping the underlying network architecture the same. The fixed number of variants (23,380) added to the training set was the total number of variants available in the species with the smallest number of missense variants, i.e., bonobos. The training procedure was repeated five times to obtain the median performance for each classifier.

[0372] Saturation of all potential human missense mutations with increasing numbers of primate populations sequenced. We investigated the expected saturation of all approximately 70 million potential human missense mutations by common variants present in 504 extant primate species by simulating variants based on the trinucleotide context of common human missense variants (>0.1% allele frequency) observed in ExAC. For each primate species, we simulated four times the number of common missense variants observed in humans (approximately 83,500 missense variants with allele frequencies above 0.1%). This is because humans have roughly half the number of variants per individual compared to other primate species, and purifying selection removes approximately 50% of human missense variants at allele frequencies above 0.1% (Figure 49A).

[0373] To model the proportion of common human missense variants (>0.1% allele frequency) discovered as the size of the human cohort studied increases (Figure 56), we sampled genotypes according to ExAC allele frequency and reported the proportion of common variants observed at least once in these simulated cohorts.

[0374] In one implementation, for practical application of the PrimateAI score, for genes with dominant inheritance mode, based on the enrichment of de novo variants in cases compared to controls (Figure 21D), a threshold of >0.8 is preferred for likely pathogenic classification, <0.6 for likely benign classification, and 0.6-0.8 is intermediate; for genes with recessive inheritance mode, a threshold of >0.7 is preferred for likely pathogenic classification and <0.5 for likely benign classification.

[0375] Figure 2 shows an exemplary architecture of a deep residual network for pathogenicity prediction, referred to herein as "PrimateAI." In Figure 2, 1D refers to a one-dimensional convolutional layer. Predicted pathogenicity is on a scale from 0 (benign) to 1 (pathogenic). The network takes as input alternative sequences (51 AA) centered around human amino acid (AA) standards and variants, position-specific weight matrix (PWM) conservation profiles calculated from 99 vertebrate species, and the output of a secondary structure and solvent accessibility prediction deep learning network. This deep learning network predicts three states of protein secondary structure (helix-H, β-sheet-B, and coil-C) and three states of solvent accessibility (buried-B, middle-I, and exposed-E).

[0376] Figure 3 shows a schematic of PrimateAI, a deep learning network architecture for pathogenicity classification. Inputs to the model include 51 amino acids (AA) of the flanking sequences for both the reference sequence and the sequence with the substituted variant, conservation rates represented by three 51-AA-long position-specific weight matrices from primate, mammal, and vertebrate alignments, and the outputs of pre-trained secondary structure and solvent accessibility networks (also 51 AA long).

[0377] Figures 4A, 4B, and 4C show details of an exemplary model architecture for the pathogenicity prediction deep learning model PrimateAI, as shown in Supplementary Table 16. The shape specifies the shape of the output tensor at each layer of the model, and the activation is the activation given to the neurons in the layer. The inputs to the model are a position-specific frequency matrix for the flanking amino acid sequences around the variant (51 AA long, 20 AA deep), one-hot encoded human reference and alternative sequences (51 AA long, 20 AA deep), and outputs from the secondary structure and solvent accessibility models (51 AA long, 40 AA deep).

[0378] The example shown uses 1D convolution. In other implementations, the model can use different types of convolution, such as 2D convolution, 3D convolution, dilated or expanded convolution, transposed convolution, separable convolution, and depthwise separable convolution. Some layers also use a ReLU activation function, which greatly accelerates the convergence of stochastic gradient descent compared to saturating nonlinearities such as sigmoid or hyperbolic tangent. Other examples of activation functions that can be used by the disclosed technology include parametric ReLU, leaky ReLU, and exponential linear unit (ELU).

[0379] Some layers also use batch normalization (Ioffe and Szegedy, 2015). Regarding batch normalization, the distributions of each layer in a convolutional neural network (CNN) change during training and vary from layer to layer. This slows down the convergence speed of optimization algorithms. Batch normalization is a technique to overcome this problem. If we denote the input of a batch normalization layer by x and its output using z, batch normalization applies the following transformation on x:

[0380]

number

[0381] Batch normalization applies mean-variance normalization to the input x using μ and σ, and linearly scales and shifts it using γ and β. The normalization parameters μ and σ are calculated for the current layer over the training set using a method called exponential moving average. In other words, they are not trainable parameters. In contrast, γ and β are trainable parameters. The values ​​of μ and σ calculated during training are used in the forward pass during inference.

[0382] Figures 5 and 6 show the deep learning network architecture used to predict protein secondary structure and solvent accessibility. The input to the model is a position-specific weight matrix using conservation ratios generated by RaptorX software (for training on Protein Data Bank sequences) or an alignment of 99 vertebrate species (for training and inference on human protein sequences). The output of the second to last layer, which is 51 AA in length, becomes the input to the deep learning network for pathogenicity classification.

[0383] Figure 7A and Figure 7B show details of an exemplary model architecture for a three-state secondary structure prediction deep learning (DL) model. The shape specifies the shape of the output tensor at each layer of the model, and the activation is the activation given to the neurons in the layer. The input to the model was a position-specific frequency matrix (length 51 AA, depth 20) for the flanking amino acid sequences around the variant.

[0384] Figure 8A and Figure 8B show details of an exemplary model architecture for a three-state solvent accessibility prediction deep learning model. The shape specifies the shape of the output tensor at each layer of the model, and the activation is the activation given to the neurons in the layer. The input to the model was a position-specific frequency matrix (length 51 AA, depth 20) for the flanking amino acid sequences around the variant.

[0385] Figure 20 shows the predicted pathogenicity score at each amino acid position in the SCN2A gene, annotated for important functional domains. Plotted along the gene is the average PrimateAI score for missense substitutions at each amino acid position.

[0386] Figure 21D shows a comparison of classifiers in predicting benign outcomes for a test set of 10,000 common primate variants that were withheld from training. The y-axis represents the percentage of primate variants correctly classified as benign after normalizing each classifier's threshold to the 50th percentile score for a set of 10,000 random variants matched for mutation rate.

[0387] Figure 21E shows the distribution of PrimateAI prediction scores for de novo missense variants occurring in Deciphering Developmental Disorders (DDD) patients compared to unaffected siblings, along with the corresponding Wilcoxon rank sum P-values.

[0388] Figure 21F shows a comparison of classifiers in separating de novo missense variants in DDD cases vs. controls. Wilcoxon rank sum test p-values ​​are shown for each classifier.

[0389] Figures 22A, 22B, 22C, 22D, and 22E show classification accuracy within the 605 DDD genes with P<0.05. Figure 22A shows greater-than-expected enrichment of de novo missense mutations in affected individuals from the DDD cohort within the 605 relevant genes that was significant for de novo protein-truncating mutations (P<0.05). Figure 22B shows the distribution of PrimateAI prediction scores for de novo missense variants occurring in DDD patients vs. unaffected siblings within the 605 relevant genes, along with the corresponding Wilcoxon rank sum P-values.

[0390] Figure 22C shows a comparison of various classifiers in separating de novo missense variants in cases versus controls within 605 genes. The y-axis shows the Wilcoxon rank-sum test P-value for each classifier.

[0391] FIG. 22D shows a comparison of various classifiers plotted on receiver operating characteristic curves, with the AUC shown for each classifier.

[0392] Figure 22E shows the classification accuracy and AUC for each classifier. The classification accuracy shown is the average true positive and true negative error rate using a threshold such that the classifier predicts the same number of pathogenic and benign variants as would be predicted based on the enrichment shown in Figure 22A. To account for the fact that 33% of DDD de novo missense variants are background, the maximum achievable AUC for a perfect classifier is shown as a dotted line.

[0393] Figures 23A, 23B, 23C, and 23D show the impact of the data used for training on classification accuracy. Deep learning networks are trained using increasing numbers of common primate and human variants, up to the full dataset (385,236 variants). In Figure 23A, the classification performance of each network is benchmarked for accuracy against 10,000 reserved primate variants and de novo variants in DDD cases versus controls.

[0394] Figures 23B and 23C show the performance of networks trained using a dataset with 83,546 common human variants and 23,380 variants from a single primate or mammalian species, according to one implementation. Results are shown for each network trained with a different source of common mutation, benchmarked against 10,000 reserved primate variants (Figure 23B) and de novo missense variants in DDD cases vs. controls (Figure 23C).

[0395] Figure 23D shows the expected saturation of all potential benign human missense positions by common variants (>0.1%) of the same state in 504 extant primate species. The y-axis shows the proportion of human missense variants observed in at least one primate species, with CpG missense variants shown in green and all missense variants shown in blue. To simulate common variants in each primate species, we sampled from the set of all potential single-base substitutions with replacements and matched the trinucleotide context distribution observed for common human variants (>0.1% allele frequency) in ExAC.

[0396] Figure 24 shows correcting for the effect of sequencing coverage on the confirmation of common primate variants. The probability of observing a given variant in non-human primate species is inversely correlated with the sequencing depth at that position in the ExAC / gnomAD exon dataset. In contrast, lower gnomAD read depth did not affect the probability of observing a common human variant (>0.1% allele frequency) at that position, because the large number of human exons sequenced makes the confirmation of common mutations almost guaranteed. When selecting a match for each primate variant to train the network, the probability of selecting a variant was adjusted for the effect of sequencing depth, in addition to matching against trinucleotide context to account for mutation rate and gene conversion.

[0397] Figures 25A, 25B, 25C, and 26 show the recognition of protein motifs by the disclosed neural network. With reference to Figures 25A, 25B, and 25C, the average PrimateAI scores for variants at each amino acid position within three different protein domains are shown to illustrate the neural network's recognition of protein domains. In Figure 25A, the collagen chain COL1A2, which contains a glycine within a repeating GXX motif, is highlighted. Clinically identified mutations in collagen genes are mostly missense mutations of glycines within the GXX repeats, as these interfere with the normal assembly of collagen and exert a strong dominant-negative effect. In Figure 25B, the active site of the IDS sulfatase enzyme is highlighted, which contains a cysteine ​​in the active site that is post-translationally modified to formylglycine. In Figure 25C, the bHLHzip domain of the MYC transcription factor is shown. The basic domain contacts DNA through positively charged arginine and lysine residues (highlighted) that interact with the negatively charged sugar phosphate backbone, and the leucine zipper domain contains leucine residues (highlighted) separated by seven amino acids that are critical for dimerization.

[0398] Figure 26 includes a line plot showing the effect of perturbing each position in and around a variant on the predicted deep learning score for the variant. Inputs at nearby amino acids (positions -25 to +25) around the variant were systematically set to 0, and the change in predicted pathogenicity of the variant by the neural network was measured. The plot shows the average change in predicted pathogenicity score for 5,000 randomly selected variants with perturbations at each nearby amino acid position.

[0399] Figure 27 shows that the weight correlation pattern mimics the BLOSUM62 score matrix and Grantham score matrix. The weight correlation pattern from the first three layers of the secondary structure deep learning network shows correlations between amino acids similar to the BLOSUM62 score matrix and Grantham score matrix. The left heatmap shows the correlation of parameter weights from the first convolutional layer, following the two initial upsampling layers of the secondary structure deep learning network, between amino acids encoded using one-hot representation. The middle heatmap shows the BLOSUM62 scores between pairs of amino acids. The right heatmap shows the Grantham distances between amino acids. The Pearson correlation between deep learning weights and BLOSUM62 scores 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 score and the Grantham score was -0.72 (P = 8.09 × 10 -13 )

[0400] Figures 28A, 28B, and 28C show performance evaluation of the deep learning network PrimateAI versus other classifiers. Figure 28A shows the accuracy of the deep learning network PrimateAI in predicting benign outcomes on a test set of 10,000 primate variants that were withheld from training, and a comparison with 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 normalized thresholds for each classifier to the 50th percentile score using a set of 10,000 randomly selected variants matched to primate variants on trinucleotide context to account for mutation rate and gene conversion.

[0401] Figure 28B shows a comparison of the performance of the PrimateAI network in separating de novo missense variants in DDD cases vs. controls, along with the 20 existing methods listed above. The y-axis shows the Wilcoxon rank-sum test P-value for each classifier.

[0402] Figure 28C shows a comparison of the performance of the PrimateAI network with the 20 methods listed above in separating de novo missense variants in DDD cases versus unaffected controls within 605 disease-associated genes. The y-axis shows the Wilcoxon rank-sum test P-value for each classifier.

[0403] Figures 29A and 29B show the distribution of prediction scores for four classifiers. Histograms of prediction scores for four classifiers, including SIFT, PolyPhen-2, CADD, and REVEL, for de novo missense variants occurring in DDD cases versus unaffected controls are shown, along with the corresponding Wilcoxon rank sum P-values.

[0404] Figures 30A, 30B, and 30C compare the accuracy of the PrimateAI network and other classifiers in classifying pathogenic and benign variants in 605 disease-associated genes. The scatter plot in Figure 30A shows the performance of each classifier for DDD cases vs. controls (y-axis) and the accuracy of benign prediction for the withheld primate dataset (x-axis). Figure 30B compares different classifiers in separating de novo missense variants in cases vs. controls within 605 genes, shown on receiver operating characteristic (ROC) curves, with the area under the curve (AUC) shown for each classifier. Figure 30C shows the classification accuracy and AUC for the PrimateAI network and 20 classifiers listed in Figures 28A, 28B, and 28C. The classification accuracy shown is the average true positive and true negative rate using a threshold such that the classifier predicts the same number of pathogenic and benign variants as expected based on the enrichment in Figure 22 A. The dotted line indicates the maximum achievable AUC for a perfect classifier, assuming 67% pathogenic and 33% benign de novo missense variants in the DDD case group and 100% benign de novo missense variants in the control group.

[0405] Figures 31A and 31B show the correlation between classifier performance on expert-curated ClinVar variants and performance on an empirical dataset. The scatterplot in Figure 31A shows the classification accuracy for 10,000 retained primate variants (x-axis) and the classification accuracy for ClinVar variants (y-axis) for each of 20 other classifiers and the PrimateAI network, trained on human-only data or human + primate data. The Spearman correlation coefficient rho and associated P-value are shown. To limit the evaluation to data that was not used to train the classifier, only ClinVar variants added between January 2017 and November 2017 were used, and common human variants (>0.1% allele frequency) from ExAC / gnomAD were excluded. The accuracy of ClinVar classification shown is the average true positive and true negative rates using thresholds such that the classifier predicts the same number of pathogenic and benign variants as observed in the ClinVar dataset.

[0406] The scatter plot in Figure 31B shows the full dataset of DDD cases vs. controls (x-axis) and classification accuracy for ClinVar variants (y-axis) for each of 20 other classifiers and the PrimateAI network, trained on human-only data or human+primate data.

[0407] Figure 32 is Supplementary Table 14 showing the performance of the three-state secondary structure prediction model and the three-state solvent accessibility prediction model on annotated samples from the Protein DataBank using 6367 unrelated protein sequences for training, 400 for validation, and 500 for testing. Only proteins with sequence identity less than 25% were selected from the Protein DataBank. Because the three classes are not significantly biased for either secondary structure or solvent accessibility, we report the accuracy of the deep learning network as a measure of performance.

[0408] Figure 33 is Supplementary Table 15 showing the performance comparison of deep learning networks using predicted secondary structure labels as well as annotated secondary structure labels of human proteins from the DSSP database when they were available.

[0409] Figure 34 shows the accuracy values ​​for the 10,000 retained primate variants and the p-values ​​for de novo variants in the DDD case vs. control groups for each of the 20 classifiers evaluated in Supplementary Table 17. The PrimateAI model using human data only is our deep learning network trained using a labeled benign training dataset consisting only of common human variants (83,500 variants that are >0.1% in the population), while the PrimateAI model using human + primate data is our deep learning network trained on the full set of 385,000 labeled benign variants, including both common human and primate variants.

[0410] Figure 35 is Supplementary Table 19, which shows a comparison of the performance of different classifiers for de novo variants in the DDD case vs. control dataset, constrained to 605 disease-associated genes. To normalize between different methods, we identified a threshold for each classifier such that the classifier predicted the same number of pathogenic and benign variants as would be expected based on enrichment in the DDD and control sets. The classification accuracy shown is the average of the true positive and true negative error rates at this threshold.

[0411] Figures 49A, 49B, 49C, 49D, and 49E show the missense / synonymous ratio across the human allele frequency spectrum. Figure 49A shows missense and synonymous variants observed in 123,136 humans from the ExAC / gnomAD database, divided into four categories by allele frequency. Gray-shaded bars represent the count of synonymous variants in each category, and dark green bars represent missense variants. The height of each bar is scaled to the number of synonymous variants in each allele frequency category, and the missense / synonymous counts and ratios are displayed after adjusting for mutation rate. Figures 49B and 49C show the allele frequency spectrum for human missense and synonymous variants with the same status (IBS) as common chimpanzee variants (Figure 49B) and chimpanzee singleton variants (Figure 49C). The depletion rate of chimpanzee missense variants with common human allele frequencies (>0.1%) compared with rare human allele frequencies (<0.1%) was correlated with the associated chi-square (χ 2 ) test are shown by red boxes along with the P value.

[0412] Figure 49D shows human variants observed in at least one non-human primate species. Figure 49E shows the counts of benign and pathogenic missense variants in the entire ClinVar database (top row), compared with the counts of ClinVar variants in a cohort of 30 humans sampled from ExAC / gnomAD allele frequencies (middle row), and the counts of variants observed in primates (bottom row). Inconsistent benign and pathogenicity determinations and variants annotated only as of uncertain significance were excluded.

[0413] Figures 50A, 50B, 50C, and 50D show purifying selection for missense variants that are consistent with other species. Figure 50A shows the allele frequency spectrum for human missense and synonymous variants that are consistent with variants present in four non-primate mammalian species (murine, porcine, caprine, and bovine). The depletion rate of missense variants with common human allele frequencies (>0.1%) is shown to be significantly higher than the associated chi-squared (χ 2 ) test are shown by red boxes along with the P value.

[0414] Figure 50B is a scatterplot showing the depletion rate of missense variants observed in other species at common human allele frequencies (>0.1%) compared with the evolutionary distance of the species from humans, expressed in units of branch length (average number of substitutions per nucleotide position). The total branch length between each species and humans is shown next to the species name. Depletion values ​​for singleton and common variants are shown for species for which variant frequencies were available, except for gorillas, which included related individuals.

[0415] Figure 50C shows the counts of benign and pathogenic missense variants in a cohort of 30 humans sampled from ExAC / gnomAD allele frequencies (top row), compared with the counts of variants observed in primates (middle row), and the counts of variants observed in mice, pigs, goats, and cattle (bottom row). Conflicting benign and pathogenic calls and variants annotated only as of uncertain significance were excluded.

[0416] Figure 50D is a scatter plot showing the depletion rate of fixed missense substitutions observed in pairs of closely related species at common human allele frequencies (>0.1%) compared with the evolutionary distance of the species from humans (expressed in units of average branch length).

[0417] Figure 51 shows the expected missense:synonymous ratio across the human allele frequency spectrum in the absence of purifying selection. Gray-shaded bars represent the number of synonymous variants, and dark green bars represent the number of missense variants. The dotted line indicates the baseline formed by synonymous variants. The missense:synonymous ratio is shown for each allele frequency category. According to one implementation, the expected missense and synonymous counts in each allele frequency category were calculated by taking intronic variants from the ExAC / gnomAD dataset, which comprises 123,136 exons, and using those intronic variants to estimate the proportion of variants expected to fall into each of the four allele frequency categories based on the trinucleotide context of the variants, which takes into account mutation rate and GC bias in gene conversion.

[0418] Figures 52A, 52B, 52C, and 52D show the missense:synonymous ratios for CpG variants and non-CpG variants. Figures 52A and 52B show the missense:synonymous ratios for CpG variants (Figure 52A) and non-CpG variants (Figure 52B) across the human allele frequency spectrum using all variants from the ExAC / gnomAD exons. Figures 52C and 52D show the missense:synonymous ratios for CpG variants (Figure 52C) and non-CpG variants (Figure 52D) across the human allele frequency spectrum, constrained to only human variants that are identical to common chimpanzee polymorphisms.

[0419] Figures 53, 54, and 55 show the missense:synonymous ratios of human variants that are identical to those present in six primate species. The patterns of missense:synonymous ratios across the human allele frequency spectrum are shown for ExAC / gnomAD variants that are identical to mutations present in chimpanzees, bonobos, gorillas, orangutans, rhesus monkeys, and marmosets.

[0420] Figure 56 is a simulation showing the saturation of new common missense variants discovered by increasing the size of the human cohort studied. In the simulation, the genotypes of each sample were sampled according to the gnomAD allele frequency. The proportion of common gnomAD variants discovered is averaged over 100 simulations for each sample size from 10 to 100,000.

[0421] Figure 57 shows the accuracy of PrimateAI across various conservation profiles in genomes. The x-axis represents the percentage alignability of 51 AA around a sequence with 99 vertebrate alignments. The y-axis represents the classification accuracy performance of PrimateAI for variants in each of the conservation bins, benchmarked on a test dataset of 10,000 retained primate variants.

[0422] Figure 58 is Supplementary Table 5 showing contributions to the labeled benign training dataset from common human variants and variants present in non-human primates.

[0423] Figure 59 shows the effect of allele frequency on the expected missense:synonymous ratio in Supplementary Table 8. Expected counts of synonymous and missense variants were calculated based on the allele frequency spectrum of variants in intronic regions at least 20–30 nt away from exon boundaries, using trinucleotide context to account for mutation rate and gene conversion bias.

[0424] Figure 60 is Supplementary Table 9 showing the ClinVar analysis. According to one implementation, variants downloaded from the November 2017 build of the ClinVar database were filtered to exclude missense variants with conflicting annotations and remove variants of uncertain significance, leaving 17,775 benign and 24,853 pathogenic variants.

[0425] Figure 61 is Supplementary Table 10 showing the number of missense variants from other species discovered in ClinVar, according to one implementation. Variants were required to be identical to the corresponding human variant and to have identical nucleotides in the other two positions in the reading frame to ensure the same coding outcome.

[0426] Figure 62 is Table 1 showing one implementation of the discovery of 14 additional gene candidates in intellectual disability that had not previously reached the threshold of genome-wide significance in the original DDD study.

[0427] Figure 63 is Table 2 showing one implementation of the difference in the mean Grantham score between pathogenic and benign variants in ClinVar, which is twice the difference in de novo variants in DDD cases vs. controls within 605 disease-associated genes.

[0428] [Data Generation] All coordinates used herein refer to the human genome build UCSC hg19 / GRCh37, including coordinates for variants in other species that were mapped to hg19 using multiple sequence alignments using the procedures described in this section. Protein-coding DNA sequences and multiple sequence alignments of 99 vertebrate genomes with humans were downloaded from the UCSC Genome Browser for 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.

[0429] Human exon polymorphism data were downloaded from the Exome Aggregation Consortium (ExAC) / genome Aggregation Database (gnomAD) v2.0, which collects whole-exome sequencing (WES) data for 123,136 individuals from eight subpopulations worldwide (http: / / gnomad.broadinstitute.org / ). Variants that did not pass the default quality control filters as annotated in the ExAC VCF file or that were located outside the canonical coding region were excluded. To avoid the effects of balancing selection, we also excluded variants from within the extended MHC region (chr6:28,477,797–33,448,354) for primate analysis. The Great Ape Genome Sequencing Project provides whole-genome sequencing data and genotypes for 24 chimpanzees, 13 bonobos, 27 gorillas, and 10 orangutans (including five from the Sumatran subspecies and five from the Borneo subspecies, which were collapsed for downstream analysis). Studies on chimpanzees and bonobos provide genome sequences for an additional 35 chimpanzees. However, because variants for these additional chimpanzees were not called using the same methods as the Great Ape Genome Sequencing Project, they were excluded from the allele frequency spectrum analysis and used only to train deep learning models. Variations from these primate diversity studies have already been mapped to the human reference (hg19). Additionally, for marmosets and rhesus macaques, 16 rhesus macaque individuals and nine marmoset individuals were used to assess variation in the original sequencing of these species' genomes, but individual-level information is not available.

[0430] The Great Ape Genome Sequencing Project 4 provides whole-genome sequencing data and genotypes for 24 chimpanzees, 13 bonobos, 27 gorillas, and 10 orangutans (including five from the Sumatran subspecies and five from the Borneo subspecies, which were collapsed for downstream analysis). Studies on chimpanzees and bonobos provide genome sequences for an additional 35 chimpanzees. However, because variants for these additional chimpanzees were not called using the same methods as the Great Ape Genome Sequencing Project, they were excluded from the allele frequency spectrum analysis and used only to train deep learning models. Variation from these primate diversity studies had already been mapped to the human reference (hg19). Additionally, for marmosets and rhesus macaques, 16 rhesus macaque individuals and nine marmoset individuals were used to assess variation in the original sequencing of these species' genomes, but individual-level information is not available.

[0431] For comparison with other primates and mammals, we also downloaded SNPs from dbSNP for other species, including rhesus macaques, marmosets, pigs, cows, goats, mice, chickens, and zebrafish. dbSNP also contains additional orangutan variants, which we used only to train our deep learning model because individual genotype information was not available for allele frequency spectrum analysis. Other species, such as dogs, cats, or sheep, were discarded because dbSNP provides only a limited number of variants for those species.

[0432] To map variants to humans, a multi-species alignment of 99 vertebrates was used to ensure orthologous 1:1 mapping to human protein-coding regions. Mapping variants using an orthologous multi-species alignment was essential to remove artifacts caused by pseudogenes or retrotransposed sequences. This artifact occurs when directly mapping SNPs between species using tools such as liftOver, which allows many-to-one mapping. When the genome build of a species in dbSNP did not match the genome build of the species in the multiple sequence alignment of 99 vertebrates, liftOver was used to update the variant to the genome build used in the multiple sequence alignment. If a variant occurred in the reference / alternate orientation, the variant was accepted as being in the same state. For example, if the human reference was G and the alternate allele was A, this was considered to be in the same state as a variant in another species where the reference was A and the alternate allele was G. To ensure that variants have the same predicted protein-coding outcome in both humans and other species, we required that the other two nucleotides in the codon be identical between species for both missense and synonymous variants. Polymorphisms from each species included in the analysis are listed in Supplementary Data File 1, and detailed measures are shown in Supplementary Table 1.

[0433] To ensure that the variants from each dbSNP submitter batch were of high quality and aligned correctly to humans, we calculated the missense:synonymous ratio for each batch and confirmed that this was less than the expected ratio of 2.2:1, and that most species, particularly zebrafish and mice, which are expected to have very large effective population sizes, had ratios less than 1:1. Two batches of SNPs from cattle were excluded from further analysis because they had 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). The average missense:synonymous ratio for the remaining cattle batches was 0.8:1.

[0434] [Correcting for the effects of allele frequency on missense:synonymous ratios, mutation rates, genetic drift, and GC-biased gene conversion] In addition to the action of purifying selection, the observed depletion rates of human missense variants at high allele frequencies may also be influenced by factors unrelated to natural selection. The probability of a natural variant appearing at a particular allele frequency in a population is a function of mutation rate, gene conversion, and genetic drift, and these factors can bias the missense:synonymous ratio across the allele frequency spectrum even in the absence of selective pressure.

[0435] To calculate the expected missense:synonymous ratio for each allele frequency category in the absence of protein-coding selection, we selected variants within the intronic regions 31–50 bp upstream and 21–50 bp downstream of each exon. These regions were chosen to be far enough away to avoid the influence of extended splice motifs. Because these regions are close to the exon capture sequence edges for the ExAC / gnomAD exons, we removed any chrX regions and removed regions with an average read depth of less than 30 to ensure unbiased ascertainment of variants. Each variant and its immediately upstream and downstream nucleotides are located within one of 64 trinucleotide contexts. If the middle nucleotide is mutated to three other bases, a total of 64 × 3 = 192 possible trinucleotide configurations are possible. Because trinucleotide configurations and their reverse complements are equivalent, there are effectively 96 trinucleotide contexts. We observed that trinucleotide context has a very strong effect on mutation rates and a smaller effect on GC-biased gene conversion, making trinucleotide context useful for modeling these variables.

[0436] Within these intronic regions, we extracted variants from 126,136 ExAC / gnomAD exons and separated them into 4 × 192 categories based on four allele frequency categories (singleton, greater than singleton ~0.01%, 0.01–0.1%, >0.1%) and 192 trinucleotide contexts. We normalized the number of variants observed in each of the 4 × 192 categories (allele frequency × trinucleotide context) by dividing the total number of potential variants by their trinucleotide context (obtained by substituting each nucleotide in the intron sequence in three different ways). For each of the 192 trinucleotide contexts, we thus obtained the expected proportion of variants falling into each of 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 context (Supplementary Table 7).

[0437] To obtain the expected missense:synonymous ratio for each allele frequency category, we counted the total number of potential synonymous and missense mutations in the human genome accessible by single-nucleotide substitution and assigned each of them to one of 192 trinucleotide contexts. For each context, we used a 4 × 192 table to calculate the expected number of variants falling into each of the four allele frequency categories. Finally, we summed the numbers of synonymous and missense variants across the 192 trinucleotide contexts to obtain the expected total number of synonymous and missense variants in each of the four allele frequency categories (Figure 51 and Supplementary Table 8 (Figure 59)).

[0438] With the exception of singleton variants, for which the expected missense:synonymous ratio was 2.46:1, the expected missense:synonymous ratio was nearly constant across the allele frequency spectrum, close to the 2.23:1 ratio expected for de novo variants in the absence of natural selection. This indicates that, due to the activity of factors independent of protein-coding selection pressure (mutation rate, gene conversion, and genetic drift), variants in the singleton allele frequency category in ExAC / gnomAD are expected to have a missense:synonymous ratio approximately 10% higher than de novo mutations by default. To correct for this, we adjusted the missense:synonymous ratio for singletons downward by 10% in the allele frequency analysis (Figures 49A, 49B, 49C, 49D, and 49E, and Figures 50A, 50B, 50C, and 50D). This small adjustment reduced the estimated missense depletion rate for common human variants present in primates and other mammals by approximately 3.8% (shown in Figures 49A, 49B, 49C, 49D, and 49E, and Figures 50A, 50B, 50C, and 50D). The higher missense:synonymous ratio for singleton variants is due to the higher mutation rate of transversion mutations (which are more likely to create missense changes) and the higher allele frequency of transition mutations (which are more likely to create synonymous changes).

[0439] Furthermore, this explains the observed missense:synonymous ratio of 2.33:1 for singleton variants in ExAC / gnomAD, which exceeds the 2.23:1 ratio expected for de novo mutations. After accounting for the effect of the allele frequency spectrum on the missense:synonymous ratio, this actually reflects a 5.3% depletion of singleton variants compared to expectation, likely due to selection against pathogenic missense variants with a de novo dominant mode of inheritance. Indeed, when considering only haploinsufficient genes with a high probability of loss of function (pLI > 0.9), the missense:synonymous ratio for ExAC / gnomAD singleton variants is 2.04:1, indicating a depletion rate of approximately 17% within haploinsufficient genes. 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.

[0440] We also specifically examined the missense:synonymous ratios for CpG and non-CpG variants across the human allele frequency spectrum due to the large differences in mutation rates (Figures 52A, 52B, 52C, and 52D). We found that for both CpG and non-CpG variants, human variants identical to common chimpanzee polymorphisms have nearly constant missense:synonymous ratios across the allele frequency spectrum.

[0441] [Depletion rates of human missense variants with identical status to polymorphisms in other species] To assess whether variants from other species are tolerated at common allele frequencies (>0.1%) in humans, we identified human variants that were identical to mutations 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%, 0.01%–0.1%, or >0.1%) and estimated the missense:synonymous ratio (MSR) decline between rare variants (<0.1%) and common variants (>0.1%). The depletion rate of identical missense variants at common human allele frequencies (>0.1%) indicates the proportion of variants from other species that are sufficiently deleterious to be eliminated by natural selection at common allele frequencies in humans.

[0442]

number

[0443] Missense:synonymous ratios and percentages of depletion were calculated for each species and are shown in Figure 50B and Supplementary Table 2. Additionally, chi-squared (χ2) tests of similarity were performed on 2 x 2 contingency tables to test whether the difference in missense:synonymous ratios between rare and common variants was significant for common chimpanzee variants (Figure 49A), chimpanzee singleton variants (Figure 49C), and mammalian variants (Figure 50A).

[0444] Because sequencing was performed only on a limited number of individuals from the Great Ape Diversity Project, we used the human allele frequency spectrum from ExAC / gnomAD to estimate the proportion of sampled variants that were rare (<0.1%) or common (>0.1%) in the general chimpanzee population. We sampled a cohort of 24 individuals based on ExAC / gnomAD allele frequencies and identified missense variants that were observed once or more than once in this cohort. Variants observed more than once had a 99.8% probability of being common in the general population (>0.1%), whereas variants observed only once in the cohort had a 69% probability of being common in the general population. Figures 49B and 49C show that depletion of singleton chimpanzee variants at high allele frequencies in humans is observed, but not common chimpanzee variants, as a result of a proportion of chimpanzee singleton variants being rare deleterious mutations. Approximately half of the chimpanzee variants observed in the cohort of 24 individuals were observed only once, and approximately half were observed more than once.

[0445] To confirm that the observed depletion rates for missense variants in more distantly related mammals were not due to the confounding effects of well-conserved, and therefore more precisely aligned, genes, we repeated the above analysis, constraining it to only genes with an average nucleotide identity of greater than 50% in multiple sequence alignments of 11 primates and 50 mammals compared to humans (see Supplementary Table 3). This removed approximately 7% of human protein-coding genes from the analysis without substantially affecting the results.

[0446] [Fixed substitutions among primates, mammals, and distantly related vertebrates] To ensure that our results using dbSNP variation were not affected by issues with variant data or artifacts due to domestication (since the majority of species selected from dbSNP are domesticated), we also repeated the analysis using fixed substitutions from pairs of closely related species instead of intraspecific polymorphisms. Phylogenetic trees of 100 vertebrate species, along with evolutionary distance measured by branch length (the average number of nucleotide substitutions per location), were downloaded from the UCSC Genome Browser (http: / / hgdownload.soe.ucsc.edu / goldenPath / hg19 / multiz100way / hg19.100way.commonNames.nh). Closely related species pairs (branch lengths <0.25) were selected for further analysis. To identify fixed substitutions between pairs of closely related species, we downloaded coding regions from the UCSC Genome Browser for multiple sequence alignments of 99 vertebrate genomes with humans and for alignments of 19 mammalian (16 primate) genomes with humans. An additional 19 mammalian multi-species alignment was necessary because some primate species, such as bonobos, were not present in the 99-vertebrate alignment (http: / / hgdownload.soe.ucsc.edu / goldenPath / hg38 / multiz20way / alignments / knownCanonical.exo nNuc.fa.gz). In total, we obtained 15 pairs of closely related species, including five primate pairs, as listed in Figure 50D and Supplementary Table 4.

[0447] We imported multiple sequence alignments of 19 mammalian genomes or 99 vertebrate genomes with humans within the canonical coding regions to obtain nucleotide substitutions between each selected pair of vertebrates, as listed in Supplementary Data File 2. These substitutions were mapped to the human genome, provided that the other two nucleotides in the codon remained unchanged between humans and the other species, and accepting variants in either the canonical or alternative orientation. Using human variants that were identical to the fixed substitutions from the related species pairs, we calculated the missense:synonymous ratio for variants in the rare (<0.1%) allele frequency category and the common (>0.1%) allele frequency category, yielding the proportion of substitutions fixed under negative selection, as shown in Supplementary Table 4.

[0448] [ClinVar analysis of polymorphism data for humans, primates, mammals, and other vertebrates] To investigate the clinical impact of variants with identical status in other species, we downloaded the release variant summary for the ClinVar database (published November 2, 2017, at ftp: / / ftp.ncbi.nlm.nih.gov / pub / clinvar / clinvar_20171029.vcf.gz). This database contains 324,698 variants on the hg19 genome build, of which 122,884 were missense single-nucleotide variants mapping to our list of protein-coding genes (Supplementary Table 9). The majority of variants in the ClinVar database did not result in missense confounds and were excluded. We then filtered variants with conflicting pathogenicity interpretations, leaving only variants with benign, likely benign, pathogenic, and likely pathogenic annotations. Variants with benign and likely benign annotations were combined into a single category, and variants with pathogenic or likely pathogenic annotations were also combined. After the filtering step shown in Supplementary Table 9, there were a total of 24,853 variants in the pathogenic category and 17,775 variants in the benign category; the remainder were excluded as variants with unknown significance or conflicting annotations.

[0449] To obtain a baseline for ClinVar missense variants in the human population, we examined ClinVar missense variants in a cohort of 30 individuals sampled from the ExAC / gnomAD allele frequency database. This cohort size was chosen to roughly reflect the number of individuals sequenced in the Primate Diversity Project study. We report the average number of pathogenic and benign variants in the 30-person cohort from 100 such simulations (Figure 49E). Because experts have systematically annotated common human variants in ClinVar with benign outcomes, we excluded variants with allele frequencies higher than 1% to avoid this curation bias.

[0450] ClinVar variants that were identical to mutations in primates, mammals, and other vertebrates were analyzed. The number of benign and pathogenic variants for each species is shown in Supplementary Table 10. A summary of the number of ClinVar variants present in humans, primates, and more distantly related mammals is shown in Figure 49E and Figure 50B, along with results from a chi-squared (χ2) test of homology for differences in the ratio of benign to pathogenic variants.

[0451] [Generating benign variants for model training] Common variants in human populations are largely neutral except for rare cases of founder effects or balancing selection, making them suitable as benign training datasets for machine learning that are not influenced by biases from human interpretation. Using allele frequency data from 123,136 exons from the ExAC / gnomAD database (release v2.0), excluding variants that did not pass the filter, 83,546 missense variants with an overall population allele frequency of 0.1% or higher remained within canonical protein-coding transcripts.

[0452] Based on our previous results showing that variants present in primates are mostly benign in humans, we created a benign training dataset for machine learning that includes common human variants (>0.1% allele frequency), variants from chimpanzees, bonobos, gorillas, and orangutans from the Great Ape Diversity Project and additional primate sequencing, and variants from rhesus macaques, orangutans, and marmosets from dbSNP. In total, 301,690 unique primate variants were added to the benign training set, according to one implementation. The number of benign training variants contributed by each source is shown in Supplementary Table 5.

[0453] It is important to note that the majority of primate variants are common in their respective populations, with a minority being rare variants. Because the number of sequenced individuals in nonhuman primate species was limited, the set of identified variants is expected to generally represent common mutations. Indeed, the missense:synonymous ratio for variants from each primate species was found to be less than half the expected 2.23:1 ratio for de novo mutations, indicating that the majority of these are common variants that have already been subjected to selection. Furthermore, for the chimpanzee cohort, approximately 84% of the identified variants were estimated to be present at common allele frequencies (>0.1%) in their respective populations. Because approximately 50% of newly discovered missense mutations are filtered by purifying selection at common human allele frequencies (>0.1%) (Figure 49A), this figure is consistent with approximately 16% rare variants accounting for the observed depletion rate of 8.8% of human missense variants, which is identical to the observed primate mutations (Figure 49D).

[0454] Applying the estimate that approximately 20% of human missense mutations are loss-of-function equivalent, we expect primate variants to contain 3.2% fully pathogenic mutations, 91.2% benign mutations (tolerating allele frequencies >0.1%), and 5.6% intermediate mutations that do not completely abolish gene function but are not sufficiently deleterious to be eliminated at common allele frequencies (>0.1%). Although this training dataset is known to be incomplete, the classification accuracy of deep learning networks trained on a benign training dataset containing both common human and primate variants was significantly better than that trained on common human variants alone. Thus, the amount of available training data appears to be a stronger limitation on current classification accuracy. As larger numbers of individuals are sequenced in each primate species, it will be possible to prepare training datasets containing a higher proportion of common primate variants, reducing contamination from pathogenic variants in the training dataset and further improving classification performance.

[0455] [Generating unlabeled variants to supplement the benign training dataset] All potential missense variants were generated from each base position in the canonical coding region by substituting the nucleotide at that position with three other nucleotides. Variants observed in 123,136 exons from ExAC / gnomAD and variants at start or stop codons were excluded. Overall, 68,258,623 unlabeled variants were generated. Each unlabeled variant was assigned to one of 96 different trinucleotide context categories. A deep learning network was trained using a semi-supervised approach by sampling variants from this unlabeled dataset that matched variants in the benign dataset by trinucleotide context and training a classifier to distinguish between benign and unlabeled training examples.

[0456] Additional filtering for unlabeled variants By presenting examples of benign and unlabeled variants along with the flanking amino acid sequences, deep learning networks learn regions of proteins that are highly intolerant of mutations. However, the absence of common variants in a region of a protein sequence may be due to strong purifying selection or to technical artifacts that prevent variants from being called in the region. To correct for the latter, we removed variants from both the benign and unlabeled datasets from regions where the ExAC / gnomAD dataset had an average coverage less than 1. Similarly, when matching unlabeled variants to primate variants in the benign dataset during training, we excluded unlabeled variants from regions where primates did not have orthologous alignable sequences with humans in multiple sequence alignments.

[0457] [Primate variants reserved for validation and testing, as well as de novo variants from affected and unaffected individuals] For validation and testing of the deep learning network, two sets of 10,000 primate variants were randomly sampled and held back for training. The remainder of the primate variants, along with common human variants (>0.1% allele frequency), were used as a benign dataset for training the deep learning network. In addition, two sets of 10,000 unlabeled variants matched against the held-back primate variants were also sampled for the validation and testing sets.

[0458] We used the 10,000 retained primate variants in the validation set and matched the 10,000 unlabeled variants to monitor the performance of the deep learning network during the training process by measuring the network's ability to distinguish between variants in the two sets, allowing us to determine a stopping point in training once the network's performance saturates, thus avoiding overfitting.

[0459] We used the 10,000 reserved primate variants in the test dataset to benchmark the deep learning network as well as 20 other classifiers. Because different classifiers had widely varying score distributions, we used these unlabeled variants to identify the 50th percentile threshold for each classifier. To ensure a fair comparison between methods, we benchmarked each classifier on the proportion of variants in the 10,000 reserved primate variants test set that were classified as benign at the 50th percentile threshold for that classifier.

[0460] To evaluate the performance of deep learning networks in a clinical setting using de novo variants in affected individuals with neurodevelopmental disorders and healthy controls, we downloaded de novo variants from the Deciphering Developmental Disorders (DDD) study and from healthy sibling controls in the Simons Simplex Collection (SSC) autism study. The DDD study provided a confidence level for de novo variants, and we excluded variants from the DDD dataset with a threshold below 0.1 as potential false positives due to variant calling errors. Overall, there were 3512 missense de novo variants from affected individuals with DDD and 1208 missense de novo variants from healthy controls.

[0461] To better model the real-world clinical scenario of distinguishing benign and pathogenic variants of uncertain significance among a pool of disease gene candidates, we limited our analysis to de novo variants within 605 genes associated with disease (p<0.05) in DDD studies, calculated exclusively from protein-truncating mutations (Supplementary Table 18). Gene-specific enrichment of protein-truncating de novo mutations was assessed by calculating statistical significance under the null hypothesis of the expected number of de novo mutations, given the gene-specific mutation rate and the number of chromosomes considered. We selected 605 genes with a nominal P value <0.05. The excess of synonymous and missense de novo mutations within the 605 genes (Figure 22A) was calculated as the ratio of observed to expected de novo mutation counts, as well as the difference between observed and expected de novo mutations. Within these 605 genes, we observed 380 de novo missense mutations from individuals affected by DDD (Figure 22A). For each of the classifiers, including our own, a small percentage of variants had no predictions, generally because they did not map to the same transcript model used by the classifier. Therefore, our deep learning network performed downstream analyses in Figures 22A, 22B, 22C, 22D, and 22E using 362 de novo missense mutations from individuals affected by DDD and 65 de novo missense mutations from healthy controls.

[0462] [Saturation of all potential human missense mutations with increasing numbers of sequenced primate populations] We investigated the expected saturation of all 70 million potential human missense mutations by common variants present in 504 extant primate species. For each primate species, we simulated four times the number of common missense variants observed in humans (approximately 83,500 missense variants with allele frequencies greater than 0.1%). This is because humans appear to have roughly half the number of variants per individual compared to other primate species, and approximately 50% of human missense variants are removed by purifying selection at allele frequencies greater than 0.1% (Figure 49A). We assigned simulated variants based on the observed distribution of common human missense variants in 96 trinucleotide contexts. For example, if 2% of common human missense variants were from the CCC>CTG trinucleotide context, we required that 2% of the simulated variants were randomly sampled CCG>CTG mutations. This has the effect of using trinucleotide context to take into account the effects of mutation rate, genetic drift, and gene conversion bias.

[0463] The curve in Figure 23D shows the cumulative saturation of approximately 70 million potential human missense mutations by common variants present in any of the 504 primate species, assuming that all common variants (>0.1% allele frequency) in each primate species are identified. Since Figure 49A shows that approximately 50% of human missense mutations are sufficiently deleterious to prevent them from reaching common allele frequencies (>0.1%) in both humans and other primates, the curve in Figure 23D indicates that the proportion of harmless human missense mutations becomes saturated by common primate mutations as the number of primate species increases. With 504 primate species, the majority of harmless human missense mutations become saturated, while harmless CpG mutations become saturated in far fewer species due to their higher mutation rates.

[0464] To model the proportion of common human missense variants (>0.1% allele frequency) discovered by increasing the size of the human cohort studied (Figure 36), genotypes were sampled according to gnomAD allele frequency. The proportion of common gnomAD missense variants discovered was averaged over 100 simulations for each sample size from 10 to 100,000.

[0465] [Secondary structure and solvent accessibility prediction] The deep learning network for pathogenicity prediction contains 36 convolutional layers in total, 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, two models were trained to enable the network to learn protein structure from primary sequence. Both models used the same network architecture and inputs, shown in Figure 6. The input to the secondary structure and solvent accessibility network is a position-specific frequency matrix of amino acids, 51 × 20 in length, that encodes conservation information from a multiple sequence alignment of humans with 99 other vertebrate species.

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

[0467] Data preparation for secondary structure and solvent accessibility prediction To train the model, unrelated crystal structures from the Protein Databank were used. Amino acid sequences with sequence affinity greater than 25% were removed. In total, 6367 protein sequences were used for training, 400 for validation, and 500 for validation (Supplementary Table 13). The data used for training, including amino acid sequences and secondary structure and solvent accessibility labels, are available on the RaptorX website: http: / / raptorx.uchicago.edu / download / .

[0468] Because most solved crystal structures are of non-human proteins, we used the RaptorX suite (based on PSI-BLAST) to retrieve related sequences to pretrain secondary structure and solvent models, since human-based multiple sequence alignments were not publicly available. We generated multiple sequence alignments for proteins using the CNFsearch1.66_release tool from RaptorX and formed position-specific frequency matrices by counting amino acids at each position from the 99 closest alignments. For example, the specific command using RaptorX to retrieve the multiple sequence alignment for the 1u71A.fasta protein was as follows: % . / buildFeature-i 1u71A.fasta-c 10-o . / TGT / 1u71A.tgt % . / CNFsearch-a 30-q 1u71A

[0469] For each amino acid location in the dataset, we captured a window from the position-specific frequency matrix corresponding to the 51 flanking amino acids and used it to predict labels for either secondary structure or solvent accessibility for the amino acid at the center of the 51-amino acid sequence. Labels for secondary structure and relative solvent accessibility were obtained directly from the known 3D crystal structure of the protein using DSSP software, eliminating the need for prediction from the primary sequence. To incorporate the secondary structure and solvent accessibility network as part of the pathogenicity prediction network, we calculated position-specific frequency matrices from a human-based multiple sequence alignment of 99 vertebrate species. While the conservation matrices generated from these two methods are generally similar, we enabled backpropagation through the secondary structure and solvent accessibility models during training for pathogenicity prediction to allow for precise tuning of parameter weights.

[0470] [Model Architecture and Training] We trained two separate deep convolutional neural network models to predict protein secondary structure and relative solvent accessibility. The two models have identical architectures and input data but differ in their prediction states. Extensive hyperparameter searches were conducted to optimize the models for best performance. Both our deep learning network for pathogenicity prediction and our deep learning network for predicting secondary structure and solvent accessibility employ a residual block architecture, widely adopted for its success in image classification. The residual block comprises recurrent convolutional units interspersed with skip connections, which allow information from earlier layers to skip the residual block. In each residual block, the input layer is first batch normalized, followed by an activation layer using rectified linear units (ReLU). The activation is then passed through a 1D convolutional layer. This intermediate output from the 1D convolutional layer is again batch normalized and ReLU activated, followed by another 1D convolutional layer. At the end of the second 1D convolution, its output is summed with the original input to form a residual block, which acts as a skip connection by allowing the original input information to bypass the residual block. In such architectures, which the authors term deep residual learning networks, the input is preserved in its original state and the residual connections are kept free of nonlinear activations from the model, allowing for the effective training of deeper networks. The detailed architecture is provided in Figure 6 and Supplementary Table 11 (Figures 7A and 7B) and Figure 12 (Figures 8A and 8B).

[0471] Following the residual block, a softmax layer calculates the three-state probability for each amino acid, and the largest softmax probability among them determines the state of the amino acid. The model is trained using an ADAM optimizer with an accumulated categorical cross entropy loss function for the entire protein sequence. Once the network is pre-trained on secondary structure and solvent accessibility, instead of directly taking the output of the network as input for the pathogenicity prediction network, the layer before the softmax layer is taken, thereby allowing more information to pass through to the pathogenicity prediction network.

[0472] The best-tested accuracy achieved for the three-state secondary structure prediction model was 79.86% (Supplementary Table 14), similar to the state-of-the-art accuracy predicted by the DeepCNF model. The best-tested accuracy for the three-state solvent accessibility prediction model was 60.31% (Supplementary Table 14), similar to the current best accuracy predicted by RaptorX on a similar training dataset. We also compared neural network predictions using DSSP-annotated structural labels for approximately 4,000 human proteins that had crystal structures with the standard PrimateAI model, which uses only predicted structural labels. No further improvement in pathogenicity prediction accuracy was obtained when using DSSP-annotated labels (Supplementary Table 15).

[0473] [Input features of deep learning models for pathogenicity prediction] The training dataset for the pathogenicity prediction network contained 385,236 benign-labeled variants and 68,258,623 unlabeled variants after filtering. For each variant, the following input features were generated: The first input feature for each variant is a 51-amino acid flanking sequence, i.e., 25 amino acids on each side of the variant derived from the hg19 reference sequence, to provide the sequence context of the variant to the deep learning model. Overall, this flanking reference sequence is 51 amino acids long. Through empirical observation, we found that representing protein sequences using amino acids was more effective than representing protein-coding sequences using nucleotides.

[0474] The second feature is a flanking human amino acid sequence of length 51 with the alternative amino acid substituted at the central location by the variant. The alternative flanking sequence is the same as the reference flanking sequence in the first feature, except that the central location 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 into one-hot coded vectors of length 51 x 20, where each amino acid is represented by a vector of 19 amino acids with value 0 and a single amino acid with value 1.

[0475] Three position-specific frequency matrices (PFMs) are generated from the multiple sequence alignment of 99 vertebrates against the variants, including one PFM for 11 primates, one PFM for 50 mammals excluding primates, and one PFM for 38 vertebrates excluding primates and mammals. Each PFM has dimensions L × 20, where L is the length of the flanking sequences around the variant (in our case, L represents 51 amino acids).

[0476] A single PFM matrix generated from a multiple sequence alignment for all 99 vertebrate species, also of length 51 and depth 20, was used as input to the pretrained three-state secondary structure and three-state solvent accessibility networks. After pretraining the network on known crystal structures from the Protein DataBank, the last two layers for the secondary structure and solvent models were removed (global max-pooling and output layers), and the 51 × 40 geometry of the output of the previous layers was used as input to the pathogenicity prediction network. Backpropagation through the structural layers of the network was allowed to fine-tune the parameters.

[0477] [Semi-supervised learning] Because semi-supervised learning algorithms use both labeled and unlabeled instances in the training process, they can produce classifiers that achieve better performance than fully supervised learning algorithms that have only a small amount of labeled data available for training. The principle behind semi-supervised learning is that inherent knowledge in unlabeled data can be exploited to enhance the predictive ability of a supervised model that uses only labeled instances, thereby providing a potential benefit of semi-supervised learning. Model parameters learned by a supervised classifier from a small amount of labeled data can be guided by the unlabeled data toward a more realistic distribution (which closely resembles the distribution of the test data).

[0478] Another pervasive challenge in bioinformatics is the problem of data imbalance. The data imbalance phenomenon occurs when one of the classes to be predicted is underrepresented in the data because instances belonging to that class are rare (notable cases) or difficult to obtain. Ironically, the minority classes are usually the most important for learning because they can be associated with special cases.

[0479] An algorithmic approach to address imbalanced data distribution is based on ensembles of classifiers. While a limited amount of labeled data naturally leads to weaker classifiers, an ensemble of weak classifiers tends to exceed the performance of any single constituent classifier. Moreover, ensembles typically improve the accuracy of predictions obtained from a single classifier by a factor that justifies the effort and cost associated with training multiple models. Intuitively, aggregating several classifiers leads to better control of overfitting, since averaging out the high variability of individual classifiers also averages out the overfitting of classifiers.

[0480] In the absence of a suitably sized dataset of highly confidently labeled pathogenic variants, we pursued a semi-supervised learning strategy. The ClinVar database has 300,000 entries, but after removing variants of uncertain significance, only approximately 42,000 missense variants remained for which pathogenic interpretations were consistent.

[0481] A systematic review also found that these entries did not have sufficient clinical evidence to support the annotated pathogenicity. Furthermore, the majority of variants in human-curated databases tend to be in a very small set of genes, which makes them inconsistent with common human variants or variants in genome-wide validated benign training datasets using fixed chimpanzee-human substitutions. Given how disparately validated the datasets are, training a supervised learning model using human-curated variants as the pathogenic set and genome-wide common variants as the benign set is likely to introduce significant bias.

[0482] A deep learning network was trained to distinguish between a set of labeled benign variants, carefully matched to remove bias, and an unlabeled set of variants. According to one implementation, the set of 385,236 labeled benign variants included common human variants (>0.1% allele frequency) from the ExAC / gnomAD database and variants from six non-human primate species.

[0483] We sampled a set of unlabeled variants, conditioning on matches to benign variants in trinucleotide context (to account for mutation rate, genetic drift, and gene conversion), and adjusted for the effects of alignability and sequence coverage on variant ascertainment. Because the number of unlabeled variants greatly exceeds the number of labeled benign variants, we trained eight models using the same set of labeled benign variants and eight randomly sampled sets of unlabeled variants and averaged their predictions to obtain a consensus prediction.

[0484] The motivation for choosing semi-supervised learning is that human-curated variant databases are unreliable and noisy, particularly lacking reliable pathogenic variants. We derived a set of reliable benign variants from common human and primate variants from gnomAD. For pathogenic variants, we employed an iterative balanced sampling approach to initially sample pathogenic variants from the set of unknown variants (VUS variants not annotated with clinical significance).

[0485] To reduce sampling bias, an ensemble of eight models was trained using the same set of benign training variants and eight different sets of pathogenic variants. First, unknown variants were randomly sampled to represent pathogenic variants. Next, the ensemble of models was iteratively used to score a set of unknown variants that were not involved in the previous training cycle. The top-scoring pathogenic variants were then taken to replace 5% of the random unknown variants in the previous cycle. Note that because 25% more top-scoring pathogenic variants than required were retained, eight different sets of scored pathogenic variants could be sampled to replace unknown variants, increasing 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 initial randomly sampled unknown variants were replaced by reliable pathogenic variants predicted by the ensemble model. Figure 42 illustrates the iterative balanced sampling process.

[0486] Balancing the benign and unknown training sets A sampling method for unknown variants that is consistent with benign variants is useful for reducing bias in our model training. When unknown variants are randomly sampled, deep learning models often extract biased information and present trivial 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 always classify the K→M substitution as pathogenic. Therefore, it is important to balance the distribution of amino acid sequences between the two training sets.

[0487] More mutable classes, such as CpG transitions, have a large representation bias in common benign variants. Orthologous variants from other primates also follow human mutation rates, suggesting enrichment of highly mutable classes in the overall benign training set. When the sampling procedure for unknown variants is less controlled and unbalanced, deep learning models tend to be more likely to classify CpG transitions as benign compared to less commonly represented classes, such as transversions or non-CpG transitions.

[0488] To prevent deep learning models from converging to trivial, nonbiological solutions, we consider balancing the trinucleotide contexts of benign and unknown variants. A trinucleotide is formed by the base before the variant, the reference base of the variant, and the base after the variant. The reference base of the variant can then be changed to three other nucleotides. In total, there are 64 × 3 trinucleotide contexts.

[0489] [Repeated Balanced Sampling] Cycle 1 The unknown variants were sampled to match the exact number of benign variants for each trinucleotide context. In other words, in the first cycle, the benign and pathogenic training sets were mirrored with respect to the trinucleotide context of the variants. The intuition behind such a sampling method is that there is an equal representation of variants with the same mutation rate between the benign and unknown sets. This prevents the model from converging to a trivial solution based on mutation rate.

[0490] Cycle 2 to Cycle 20 For cycle 2, the trained model from cycle 1 was applied to score a set of unknown variants not involved in cycle 1, replacing 5% of the unknown variants with the top predicted pathogenic variants. This set was generated purely by the model, and no balancing was applied to the trinucleotide context within this set. The remaining 95% of unknown variants required for training were sampled to represent 95% of the counts of each nucleotide context among benign variants.

[0491] Intuitively, because cycle 1 uses a perfectly matched training set, the top predicted pathogenic variants are generated without any mutation rate bias. Therefore, there is no need to consider any bias in this set. The remaining 95% of the data is still controlled for trinucleotide context mutation rates to prevent the model from converging to a trivial solution.

[0492] For each cycle, the percentage of replaced unknown variants increases by 5%. In cycle 3, 5% of unknown variants were replaced with the top predicted pathogenic variant from the cycle 3 model. Cumulatively, the proportion of pathogenic variants increases to 10%, and the proportion of unknown variants mirrored with trinucleotide context falls to 90%. The sampling process is similar for the remaining cycles.

[0493] Cycle 21 In the final cycle, cycle 21, the entire pathogenicity training set consists of the top pathogenic variants predicted purely from the deep learning model. Because we have explicitly considered mutation rate bias in each cycle, the pathogenic variants are reliable and unaffected by mutation rate bias for use as training data. Thus, the final cycle of training yields the final deep learning model for pathogenicity prediction.

[0494] [Matching labeled benign training sets to unlabeled training sets] Balanced sampling of unlabeled variants is crucial to removing biases unrelated to variant deleteriousness. Without proper control for confounding influences, deep learning can easily discriminate between classes through inadvertently introduced biases. Common human variants tend to be enriched for variants from classes with high mutagenic potential, such as variants on CpG islands. Similarly, primate polymorphisms also follow human mutation rates, suggesting enrichment of highly mutagenic variants across benign training sets. If the unlabeled variant sampling procedure is not well controlled and balanced, deep learning networks tend to rely on mutation rate bias to classify variants, and are therefore more likely to classify CpG transitions as benign compared to less common classes such as transversions 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 above).

[0495] When matching unlabeled variants against primate variants in the labeled benign dataset, we did not allow unlabeled variants to be selected from regions of the human genome where no primate species was aligned in the multiple sequence alignment, because it was not possible to call a variant in that primate species at that location.

[0496] Within each of the 96 trinucleotide contexts, we corrected for sequencing coverage for primate variants. Due to the large number of humans sequenced, common variants in the human population are observed frequently enough that they are well confirmed, even in areas with low sequencing coverage. The same was not true for primate variants, because only a small number of individuals were sequenced. We divided the genome into 10 bins based on sequencing coverage in the ExAC / gnomAD exons. For each bin, we measured the proportion of primate variants in the labeled benign dataset versus the unlabeled dataset. Based on sequencing coverage alone, we used linear regression to calculate the probability that a primate variant was a member of the labeled benign dataset (Figure 24). When selecting unlabeled variants to match against primate variants in the labeled benign dataset, we used a regression coefficient to weight the probability of sampling a variant based on sequencing coverage at that location.

[0497] [Generation of benign and unknown variants] Common variants in the human population Recent studies have demonstrated that common variants in human populations are generally benign. According to one implementation, gnomAD provides 90,958 nonsynonymous SNPs with a minor allele frequency (MAF) of 0.1% or greater within canonical coding regions. Variants that pass the filter are retained. Indels are filtered out. Variants occurring at start or stop codons, as well as protein-truncating variants, are removed. Upon probing the subpopulations, the total number of missense variants with a MAF of 0.1% or greater within each subpopulation increases to 245,360 according to one implementation. These variants form part of the training set of benign variants.

[0498] Common polymorphisms in great apes Because coding regions are known to be highly conserved, it is easy to hypothesize that if polymorphisms segregate at high frequencies in great ape populations, they may also have mild effects on fitness. Polymorphism data for bonobos, chimpanzees, gorillas, and orangutans from the Great Ape Genome Project and other studies were integrated with SNPs for rhesus macaques and marmosets from dbSNP.

[0499] Generation of unknown variants All potential variants are generated from each base position in the canonical coding region by substituting the nucleotide at that position with three other nucleotides. A new codon is formed, leading to a potential change in the amino acid at that position. Synonymous changes are filtered out.

[0500] Variants observed in the gnomAD dataset are removed. Variants occurring at the start or stop codon, as well as variants forming stop codons, are removed. For SNPs with multiple gene annotations, a canonical gene annotation is selected to represent the annotation of the SNP. In total, according to one implementation, 68,258,623 unknown variants are generated.

[0501] Additional filtering of variants Some regions of the human genome are known to be difficult to align reads to, and including these regions would cause confounding effects in the training and testing datasets. For example, regions that are subject to high selective pressure tend to have a limited number of polymorphisms. On the other hand, regions that are difficult to sequence also have fewer polymorphisms. To avoid such confounding input to our model, we removed variants from genes that were not sequenced by gnomAD.

[0502] Generally, benign variants are found in well-sequenced regions that tend to be conserved across multiple species. Meanwhile, unknown variants are randomly sampled across the genome, including some low-coverage regions. This leads to ascertainment bias between the benign and unknown sets. To reduce bias, we removed variants with a read depth of less than 10 in gnomAD. We also removed all variants with more than 10% missing data in flanking sequence alignments across all mammalian species.

[0503] Validation and verification data For the validation and testing of the pathogenicity model, according to one implementation, two sets of 10,000 benign variants were randomly sampled from the large pool of benign variants for validation and testing, respectively. The remaining benign variants were used to train the deep learning model. These variants were specifically sampled from orthologous primate variants to ensure fair comparison between methods, since some methods are trained on common human variants. According to one implementation, two sets of 10,000 unknown variants were randomly sampled separately for validation and testing. Ensure that the number of unknown variants in each of the 192 trinucleotide contexts matches the number of benign variants in the validation set and the testing set, respectively.

[0504] The performance of several methods was evaluated in a clinical setting using de novo variants in affected children with autism or developmental disorder (DDD) and their unaffected siblings. Overall, according to one implementation, there were 3821 missense de novo variants from the DDD case group and 2736 missense de novo variants from the autism case group. According to one implementation, there were 1231 missense de novo variants for unaffected siblings.

[0505] [Deep Learning Network Architecture] The pathogenicity prediction network receives five direct inputs and two indirect inputs via secondary structure and solvent accessibility networks. The five direct inputs are 51 amino acid sequences long and 20 deep (encoding 20 different amino acids) and include a reference human amino acid sequence without variants (1a), an alternative human amino acid sequence substituted with variants (1b), PFMs from a multiple sequence alignment of primate species (1c), a PFM from a multiple sequence alignment of mammalian species (1d), and a PFM from a multiple sequence alignment of more distantly related vertebrate species (1e). The secondary structure and solvent accessibility networks each receive PFMs from multiple sequence alignments (1f) and (1g) as inputs and provide outputs as inputs to the main pathogenicity prediction network. The secondary structure and solvent accessibility networks are pre-trained on known protein crystal structures from the Protein DataBank, enabling backpropagation during pathogenicity model training.

[0506] 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 PFMs from a multiple sequence alignment of primates, mammals, and vertebrates (merge 1a). Similarly, the human alternative amino acid sequence (1b) is merged with PFMs from a multiple sequence alignment of primates, mammals, and vertebrates (merge 1b). This creates two parallel paths, one for the reference sequence and one for the alternative sequence substituted with the variant.

[0507] The combined feature maps (Merge 1a and Merge 1b) for both the reference and alternate channels are passed through a series of six residual blocks (Layers 2a-7a, Merge 2a, and Layers 2b-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), which thoroughly blends the data from the reference and alternate channels. The data then has two paths through the network in parallel: either through a series of six residual blocks, each containing two convolutional layers as defined in Section 2.1 (Merge 3-9, Layers 21, Layers 9-46 except Layer 34), or via skip connections (Layers 21, 37, Layer 47) connecting the outputs of every other residual block after passing through a 1D convolution. Finally, the combined activation (Merge 10) is fed to another residual block (Layers 48-53, Merge 11). The activations from the merger are fed into a 1D convolution and sigmoid activation with filter size 1 (layer 54), then passed through a global max-pooling layer that selects a single value representing the network's prediction of variant pathogenicity. A schematic illustration of the model is shown in Figure 3 and Supplementary Table 16 (Figure 4A, Figure 4B, and Figure 4C).

[0508] [Model Overview] We developed a semi-supervised deep convolutional neural network (CNN) model to predict variant pathogenicity. The input features to the model included the protein sequence and conservation profile of flanking variants and the depletion rate of missense variants in specific gene regions. The deep learning model predicted the changes in secondary structure and solvent accessibility caused by the variants 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, a reliable source of pathogenic variants remains lacking. We first 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 predicted with high confidence. Finally, we demonstrated that our model outperforms existing methods in distinguishing de novo variants causing developmental disorders from benign variants in humans.

[0509] [Use residual blocks] Figure 17 shows the residual block. Both our deep learning model for pathogenicity prediction and our deep learning model for predicting secondary structure and solvent accessibility employ the definition of the residual block first presented in [1]. The structure of the residual block is shown in the following diagram. The input layer is first batch normalized, followed by a nonlinear activation (ReLU). The activation is then passed through a 1D convolutional layer. This intermediate output from the 1D convolutional layer is again batch normalized and ReLU activated, followed by another 1D convolutional layer. At the end of the second 1D convolution, the output is integrated with the original output. In such an architecture, the input is kept in its original state, and the residual connections are kept free of the nonlinear activation of the model.

[0510] Dilated / dilated convolutions allow for large receptive fields with a small number of trainable parameters. Dilated / dilated convolutions are convolutions in which the kernel is applied over an area greater than its length by skipping input values ​​using some step, also called the dilated convolution rate or dilation factor. Dilated / dilated convolutions add spacing between elements of the convolution filter / kernel so that a longer interval of neighboring input entries (e.g., nucleotides, amino acids) is considered when the convolution operation is performed. This allows for the incorporation of long-range contextual dependencies in the input. Dilated convolutions save partial convolution calculation results for reuse when neighboring nucleotides are processed.

[0511] [Novelty of our model] Our method differs from existing methods for predicting variant pathogenicity in three ways. First, our method employs a novel architecture of a semi-supervised deep convolutional neural network. Second, reliable benign variants are obtained from common human and primate variants from gnomAD, while a high-confidence pathogenic training set is generated through iterative balanced sampling and training to avoid cyclical training and testing of models using the same human-curated variant database. Third, deep learning models for secondary structure and solvent accessibility are integrated into our pathogenicity model architecture. Information obtained from the structure and solvent model is not limited to label prediction for specific amino acid residues. Rather, the readout layer is removed from the structure and solvent model, and the pre-trained model is integrated with the pathogenicity model. During pathogenicity model training, the pre-trained structure and solvent layers are also backpropagated to minimize error. This helps the pre-trained structure and solvent model focus on the pathogenicity prediction problem.

[0512] Training of secondary structure and solvent accessi...

Claims

1. 1. A system comprising: one or more processors; When executed by the one or more processors, the system further comprises: Towards a pathogenicity prediction neural network, a reference amino acid sequence encoded by the protein; an encoded alternative amino acid sequence for the protein comprising the variant amino acid; and At least a conservation profile generated from amino acids of related species aligned with a reference amino acid sequence corresponding to the encoded reference amino acid sequence. and the encoded alternative amino acid sequence, the encoded reference amino acid sequence, and Based on at least said stored profile, generating, as an output of the pathogenicity prediction neural network, a pathogenicity prediction indicating the likelihood that the variant amino acid is benign or pathogenic; A non-transitory storage medium containing computer instructions for causing the Including, the system.

2. The system of claim 1 , wherein the conservation profile comprises a position-specific frequency matrix generated from a multiple sequence alignment of the related species.

3. The system of claim 1 , wherein the conservation profile comprises a position-specific scoring matrix generated from a multiple sequence alignment of the related species.

4. At least the stored profile comprises: a primate conservation profile generated from amino acids of primate species aligned with the reference amino acid sequence; a mammalian conservation profile generated from amino acids of mammalian species aligned with the reference amino acid sequence; or A vertebrate conservation profile generated from the amino acids of vertebrate species aligned with the reference amino acid sequence The system of claim 1 , comprising one or more of:

5. 2. The system of claim 1, further comprising computer instructions that, when executed by the one or more processors, cause the system to input three-state secondary structure predictions for each amino acid in the encoded reference amino acid sequence into the pathogenicity prediction neural network.

6. 6. The system of claim 5, further comprising computer instructions that, when executed by the one or more processors, cause the system to generate the three-state secondary structure predictions based on the conservation profile using a secondary structure neural network.

7. 2. The system of claim 1, further comprising computer instructions that, when executed by the one or more processors, cause the system to input three-state solvent accessibility predictions for each amino acid in the encoded reference amino acid sequence into the pathogenicity prediction neural network.

8. 8. The system of claim 7, further comprising computer instructions that, when executed by the one or more processors, cause the system to generate the three-state solvent accessibility prediction based on the storage profile using a solvent accessibility neural network.

9. 2. The system of claim 1, wherein the variant amino acids of the encoded alternative amino acid sequence are caused by nucleotide substitutions of nucleotides in the reference amino acid sequence.

10. When executed by the one or more processors, the system further comprises: inputting the encoded reference amino acid sequence into a first convolutional layer of the pathogenicity prediction neural network; inputting the encoded alternative amino acid sequence into a second convolutional layer of the pathogenicity prediction neural network; inputting the stored profile into a third convolutional layer of the pathogenicity prediction neural network; The system of claim 1 , further comprising computer instructions for:

11. 1. A computer-implemented method, comprising: Towards a pathogenicity prediction neural network, a reference amino acid sequence encoded by the protein; an encoded alternative amino acid sequence for the protein comprising the variant amino acid; and At least a conservation profile generated from amino acids of related species aligned with a reference amino acid sequence corresponding to the encoded reference amino acid sequence. and the encoded alternative amino acid sequence, the encoded reference amino acid sequence, and Based on at least said stored profile, generating, as an output of the pathogenicity prediction neural network, a pathogenicity prediction indicating the likelihood that the variant amino acid is benign or pathogenic; 23. A computer-implemented method comprising:

12. The computer-implemented method of claim 11 , wherein the conservation profile comprises a position-specific frequency matrix generated from a multiple sequence alignment of the related species.

13. The computer-implemented method of claim 11 , wherein the conservation profile comprises a position-specific scoring matrix generated from a multiple sequence alignment of the related species.

14. At least the stored profile comprises: a primate conservation profile generated from amino acids of primate species aligned with the reference amino acid sequence; a mammalian conservation profile generated from amino acids of mammalian species aligned with the reference amino acid sequence; or A vertebrate conservation profile generated from the amino acids of vertebrate species aligned with the reference amino acid sequence The computer-implemented method of claim 11 , comprising one or more of:

15. 12. The computer-implemented method of claim 11, further comprising inputting into the pathogenicity prediction neural network a three-state secondary structure prediction for each amino acid in the encoded reference amino acid sequence.

16. 16. The computer-implemented method of claim 15, further comprising utilizing a secondary structure neural network to generate the three-state secondary structure prediction based on the conservation profile.

17. 12. The computer-implemented method of claim 11, further comprising inputting into the pathogenicity prediction neural network three-state solvent accessibility predictions for each amino acid in the encoded reference amino acid sequence.

18. 20. The computer-implemented method of claim 17, further comprising utilizing a solvent accessibility neural network to generate the three-state solvent accessibility prediction based on the conservation profile.

19. 20. The computer-implemented method of claim 17, wherein the variant amino acids of the encoded alternative amino acid sequence are caused by nucleotide substitutions of nucleotides in the reference amino acid sequence.

Images