Variant pathogenicity classifier based on recurrent neural network

Abstract

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

Description

[0001] This application is a divisional application of the application with application number 201880043817.3, application date October 15, 2018, and invention name “Deep Convolutional Neural Network Training Method Based on Deep Learning”.

[0002] appendix

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

[0004] Priority application

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

[0006] Merge

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

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

[0009] PCT Patent Application No. PCT / US2018 / _______ (Attorney Docket No. ILLM 1000-10 / IP-1613-PCT), entitled “Semi-Supervised Learning for Training a Set of Deep Convolutional Neural Networks,” filed concurrently by Laksshman Sundaram, Kai-How Farh, Hong Gao, and Jeremy Francis McRae on October 15, 2018, and subsequently published as PCT Publication No. WO______.

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

[0011] U.S. nonprovisional patent application titled “Deep Convolutional Neural Networks for Variant Classification” filed concurrently by Laksshman Sundaram, Kai-How Farh, Hong Gao, Samskruthi ReddyPadigepati, and Jeremy Francis McRae (Attorney Docket No. ILLM 1000-6 / IP-1612-US).

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

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

[0014] Reference 2 - S. Arik, 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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[0029] Document 4 describes the residual blocks and residual connections used by the disclosed technology.

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

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

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

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

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

[0035] Reference 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 / atrous convolutions.

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

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

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

[0039] Incorporation by Reference Form for Electronically Submitted Applications

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

[0041]

[0042]

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

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

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

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

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

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

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

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

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

[0052] Supplementary Table 21: Results of enrichment tests for de novo mutations in genes with an FDR < 0.1. This table includes counts of observed protein-truncating (PTV) de novo mutations and counts of other protein-altering de novo mutations, once for all missense de novo mutations and once for only damaging missense mutations. P values ​​are presented for inclusion of all missense sites versus exclusion of low-scoring missense sites. Note: This table is listed in Supplementary Table 21.txt.

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

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

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

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

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

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

[0059] The subject matter discussed in this section should not be admitted to be prior art simply by virtue of being mentioned in this section. Similarly, issues mentioned in this section or associated with the subject matter provided as background should not be assumed to have been previously recognized in the prior art. The subject matter in this section merely represents different approaches that themselves may also correspond to implementations of the claimed technology.

[0060] Machine Learning

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

[0062] Y=f(X)+∈

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

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

[0065]

[0066] In the above equation, is an estimate of ∈, which is usually considered as a black box, meaning that only The relationship between the input and output is known, but the question of why it works remains unanswered.

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

[0068] Neural Networks

[0069] Figure 1A One implementation of a fully connected neural network with multiple layers is shown. A neural network is a system of interconnected artificial neurons (e.g., a1, a2, a3) that exchange messages between each other. The neural network shown has three inputs, two neurons in 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(·). The connections have digital weights (e.g., w) that are tuned during the training process. 11 、w 21 、w 12 、w 31 、w 22 、w 32 、v 11 、v 22 ) so that a properly trained network responds correctly to recognition when fed an image. 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 includes multiple layers of feature detection neurons. Each layer has many neurons that respond to different combinations of inputs from the previous layer. The layers are constructed so that the first layer detects a set of raw patterns in the input image data, the second layer detects patterns of patterns, and the third layer detects patterns of these patterns.

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

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

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

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

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

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

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

[0077] In the drawings, like reference numerals generally refer to like parts throughout the different views. Moreover, the drawings are not necessarily drawn to scale, with emphasis 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, in which:

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

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

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

[0081] Figure 1Dis an implementation of a subsampling layer (averaging / max pooling) according to one implementation of the disclosed technology.

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

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

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

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

[0086] Figure 1I The batch normalization forward pass is shown.

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

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

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

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

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

[0092] Figure 1O Illustrated dilated convolution.

[0093] Figure 1P Shown is an implementation of stacked dilated convolutions.

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

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

[0096] Figure 3 A schematic diagram depicting the deep learning network structure of PrimateAI, the pathogenicity classification framework.

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

[0098] Figure 5 and 6 Describes a deep learning network for predicting protein secondary structure and solvent accessibility.

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

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

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

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

[0103] Figure 11 It is an implementation method for generating a position frequency matrix (abbreviated PFM), where PFM is also called a position weighting matrix (abbreviated PWM) or a position specific scoring matrix (abbreviated PSSM).

[0104] Figure 12 、 13 , 14 and 15 show the treatment of the secondary structure and solvent accessibility subnetworks.

[0105] Figure 16 Operation of the variant pathogenicity classifier. As used herein, the term variant is also referred to as single nucleotide polymorphism (abbreviated SNP), often also referred to as single nucleotide variant (abbreviated SNV).

[0106] Figure 17 Residual box plots are illustrated.

[0107] Figure 18 The neural network architecture of the secondary structure and solvent accessibility subnetworks was described.

[0108] Figure 19 The neural network architecture of the variant pathogenicity classifier is shown.

[0109] Figure 20 The predicted pathogenicity score at each amino acid position in the SCN2A gene is shown, with annotations of the major functional domains.

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

[0111] Figure 21B Figure 2. Distribution of PrimateAI prediction scores for de novo missense variants arising in patients with developmental disorders (DDD) compared with their unaffected siblings, along with corresponding Wilcoxon rank-sum test P values.

[0112] Figure 21C Depicted is a comparison of classifiers for separation of de novo missense variants in DDD cases compared to controls. Wilcoxon rank sum test P values ​​for each classifier are shown.

[0113] Figure 22A Figure 3: 605 genes associated with significant (P < 0.05) de novo protein-truncating variants, showing an unexpected enrichment of de novo missense mutations in affected individuals of the DDD cohort.

[0114] Figure 22B Figure 5. Distribution of PrimateAI prediction scores for de novo missense variants in individuals with DDD compared with their unaffected siblings in 605 relevant genes, along with corresponding Wilcoxon rank-sum test P values.

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

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

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

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

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

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

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

[0122] Figure 28A 、 28B and 28C show the performance evaluation of the deep learning network PrimateAI and other classifiers.

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

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

[0125] Figure 31A and 31B The correlation between the classifier's performance on human-selected ClinVar variants and its performance on an empirical implementation is shown.

[0126] Figure 32 Supplementary Table 14 shows the performance of the 3-state secondary structure and 3-state solvent accessibility prediction models on annotated samples of the Protein Data Bank.

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

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

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

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

[0131] Figure 37 、 38 , 39, 40 and 41 show the respective cycles of the disclosed semi-supervised learning.

[0132] Figure 42 is an illustration of the iterative balanced sampling process.

[0133] Figure 43 An implementation of a computing environment for generating a benign implementation is described.

[0134] Figure 44 One implementation for generating benign human missense SNPs is described.

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

[0136] Figure 46 One implementation is described for classifying SNPs from non-human primate species (e.g., chimpanzees) that have reference codons that match humans as benign.

[0137] Figure 47 An implementation of computing enrichment scores and comparing them is described.

[0138] Figure 48 An implementation of one implementation of a benign SNP is described.

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

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

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

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

[0143] Figure 53 、 54 and 55 illustrate missense:synonymous ratios for human variants that are identical to the six primate states.

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

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

[0146] Figure 58 Supplementary Table 5 shows the contribution of common human variants and variants present in non-human primates to labeling one implementation of benign training.

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

[0148] Figure 60 is Supplementary Table 9, showing the ClinVar analysis.

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

[0150] Figure 62 Table 1 shows one implementation method for discovering 14 additional candidate genes for intellectual disability.

[0151] Figure 63 Table 2 shows one implementation of the mean difference in Grantham scores between pathogenic and benign variants in ClinVar.

[0152] Figure 64 One implementation of the enrichment analysis for each gene is shown.

[0153] Figure 65 One implementation of gene set enrichment analysis is shown.

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

[0155] 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 specific application and its requirements. Various modifications to the disclosed implementations will be apparent to those skilled in the art, and the general principles defined herein can 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 should be given the widest scope consistent with the principles and features disclosed herein.

[0156] introduce

[0157] Convolutional Neural Networks

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

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

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

[0161] Convolutional neural networks learn highly nonlinear mappings by interconnecting layers of artificial neurons arranged in many different layers with activation functions that correlate these layers. They consist of one or more convolutional layers interspersed with one or more subsampling and nonlinear layers, typically followed by one or more fully connected layers. Each element of a convolutional neural network receives input from a set of features from the previous layer. Convolutional neural networks learn simultaneously because neurons in the same feature map have the same weights. These locally shared weights reduce the complexity of the network, allowing convolutional neural networks to avoid the complexity of feature extraction and data reconstruction during regression or classification when multidimensional input data enters the network.

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

[0163] 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): it computes 32 filters on its input. Each of these 32 output channels contains a grid of 26x26 values, which is a map of the filter's response on the input, indicating the response of that filter pattern at different locations in the input. This is what the term "feature map" means: each dimension in the depth axis is a feature (or filter), and the 2D tensor output [:,:,n] is a 2D spatial map of the response of this filter on the input.

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

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

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

[0167] Training a Convolutional Neural Network

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

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

[0170] △w i =x i δ

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

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

[0173]

[0174] In the above equation: the arrows indicate the updating of the values; t m is the target value of neuron m; is the current output of the calculated neuron m; a n is the input n; and α is the learning rate.

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

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

[0177]

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

[0179]

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

[0181]

[0182]

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

[0184]

[0185]

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

[0187]

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

[0189]

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

[0191]

[0192] The following learning rate α is used to update the weights of the hidden layer:

[0193] vnm←wnm+αδhman

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

[0195]

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

[0197] In the equation above, α is the learning rate. Furthermore, the loss is calculated as the average over a set of n data pairs. Upon linear convergence, the computation terminates when the learning rate α becomes sufficiently small. In other implementations, gradients are computed using only selected data pairs fed into the Nesterov accelerated gradient and adaptive gradient algorithms to inject computational efficiency.

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

[0199]

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

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

[0202] Convolutional layer

[0203] The convolutional layers of a convolutional neural network act as feature extractors. The convolutional layers act as adaptive feature extractors 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 on the input image, thereby producing 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 is f*g

[0204] is defined as:

[0205]

[0206] The convolution operation involves sliding a kernel over the input image. For each kernel location, the overlapping values ​​of the kernel and the 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 centered on the kernel. The different outputs produced by many kernels are called feature maps.

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

[0208] Subsampling layer

[0209] Figure 1D This is an implementation of a subsampling layer according to one embodiment of the disclosed technology. The subsampling layer reduces the resolution of features extracted by the convolutional layer, making the extracted features or feature maps more robust to noise and deformation. In one implementation, the subsampling layer employs two types of pooling operations: average pooling and maximum pooling. Pooling divides the input into non-overlapping two-dimensional regions. Average pooling calculates the average of four values ​​within the region. Maximum pooling selects the maximum of the four values.

[0210] In one implementation, a subsampling layer consists of a pooling operation on a group of neurons in the previous layer: mapping their output to just one input in max pooling, or mapping their output to the average of the inputs in average pooling. In max pooling, the output of the pooling neuron is the maximum value residing in the input, as described below:

[0211]

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

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

[0214]

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

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

[0217] Nonlinear layer

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

[0219]

[0220]

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

[0222]

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

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

[0225]

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

[0227] Convolution Example

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

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

[0230] Forward pass

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

[0232]

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

[0234]

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

[0236]

[0237] Backpropagation

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

[0239]

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

[0241]

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

[0243]

[0244]

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

[0246]

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

[0248]

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

[0250]

[0251]

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

[0253]

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

[0255]

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

[0257]

[0258]

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

[0260]

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

[0262]

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

[0264]

[0265]

[0266] Residual Connection

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

[0268] Residual learning and skip connections

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

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

[0271] Dilated convolution

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

[0273] WaveNet

[0274] Wavenet is a deep neural network for generating raw audio waveforms. Wavenet differs from other convolutional networks because it can achieve a relatively large "field of view" at a low cost. Furthermore, it can condition the signal both locally and globally, allowing Wavenet to be used as a text-to-speech (TTS) engine with multiple voices, where TTS applies both local conditioning and voice-specific global conditioning.

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

[0276]

[0277]

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

[0279] Batch Normalization

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

[0281] Batch Normalization can be thought of as another layer that can be inserted into the model architecture, just like a fully connected or convolutional layer. Batch Normalization layers are typically used after convolutional or densely connected layers. It can also be used before convolutional or densely connected layers. Both implementations can be used by the disclosed techniques and are Figure 1L . The Batch Normalization layer takes an axis argument that 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 the data format set to "channels_last". However, in the niche use case of a Conv2D layer with the data format set to "channels_first", the feature axis is axis 1; the axis argument in Batch Normalization can be set to 1.

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

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

[0284] Covariate shift is a known problem in deep learning and frequently occurs in real-world problems. A common covariate shift problem is a difference in the distribution of the training and test sets, which can lead to suboptimal generalization performance. This problem is often addressed through a normalization or whitening preprocessing step. However, whitening, in particular, is computationally expensive and therefore impractical in online settings, especially if covariate shift occurs throughout different layers.

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

[0286] Forward pass

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

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

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

[0290] Backward pass

[0291] Since normalization is a differentiable operation, it can be done as Figure 1K The backward pass is computed as described in .

[0292] 1D convolution

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

[0294] Global average pooling

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

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

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

[0298] Deep Learning in Genomics

[0299] Genetic variation can help explain many diseases. Every person has a unique genetic code, and many genetic variants exist within a population. Most undesirable genetic variants have been removed from the genome through natural selection. It is important to determine which genetic variants are likely pathogenic, or harmful. This will help researchers focus on genetic variants that may be causative and accelerate the diagnosis and treatment of many diseases.

[0300] Modeling the nature and functional effects of variants, such as pathogenicity, is an important but challenging task in genomics. Despite rapid advances in functional genomic sequencing technologies, interpreting the functional consequences of variants remains a significant challenge due to the complexity of cell-type-specific transcriptional regulatory systems.

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

[0302] There is an opportunity to use semi-supervised algorithms to build deep learning-based pathogenicity classifiers that can accurately predict the pathogenicity of variants. This may result in a database of pathogenic variants that is free of human confirmation bias.

[0303] Regarding pathogenicity classifiers, deep neural networks are a type of artificial neural network that successfully model high-level features using multiple layers of complex, nonlinear transformations. Deep neural networks provide feedback via backpropagation, which uses the difference between observed and predicted outputs to adjust parameters. Deep neural networks have evolved with the availability of large-scale training datasets, the power of parallel and distributed computing, and sophisticated training algorithms. They have already facilitated significant advances in numerous fields, including computer vision, speech recognition, and natural language processing.

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

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

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

[0307] Given the multidimensional and high-dimensional nature of sequencing data, deep neural networks hold great promise in bioinformatics research due to their broad applicability and high predictive power. Convolutional neural networks have been used to address sequence-based problems in genomics, such as motif discovery, pathogenic variant identification, and gene expression inference. Convolutional neural networks utilize a weight-sharing strategy that is particularly useful for studying DNA, as they can capture sequence motifs—short, recurring local patterns in DNA that are hypothesized to have significant biological functions. A hallmark of convolutional neural networks is the use of convolutional filters. Unlike traditional classification methods, which rely on carefully designed and handcrafted features, convolutional filters perform adaptive feature learning, akin to mapping raw input data into an informational representation of knowledge. In this sense, convolutional filters act as a series of motif scanners, as a set of such filters can identify relevant patterns in the input and self-update during training. Recurrent neural networks can capture long-range dependencies in sequence data of varying lengths, such as protein or DNA sequences.

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

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

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

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

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

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

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

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

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

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

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

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

[0320] The large number of potential variants of unknown significance and the critical importance of accurate variant classification for clinical applications have prompted multiple attempts to address the problem using machine learning, but these efforts have been largely limited by the insufficient number of common human variants and the questionable quality of annotation in curated databases. Variants from six nonhuman primates have contributed over 300,000 unique missense variants that do not overlap with common human variants and are largely benign, significantly expanding the size of the training datasets available for machine learning approaches.

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

[0322] We used semi-supervised learning to overcome the problem of the training set containing only variants with benign labels by initially training an ensemble of networks to distinguish likely benign primate variants from random unknown variants matched in mutation rate and sequencing coverage. This ensemble of networks was used to score the entire set of unknown variants and influence the selection of unknown variants for the next iteration of the seeding classifier by biasing them towards variants with greater prognostic value. Gradual steps were taken at each iteration to prevent the model from prematurely converging to suboptimal results.

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

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

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

[0326] While less aggregated whole-genome data are available as exome data, limiting our ability to detect the effects of natural selection in deep intronic regions, we were also able to calculate observed vs. expected counts for cryptic splicing variants located far from exonic regions. Overall, we observed a 60% depletion of cryptic splicing variants at distances >50 nt from exon-intron boundaries. This attenuated signal is likely due to a combination of smaller sample sizes with whole-genome data compared to exome data and the greater difficulty in predicting the effects of deep intronic variants.

[0327] the term

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

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

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

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

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

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

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

[0335] The term "chromosome" refers to the genetic carrier of a living cell, derived from strands of chromatin containing DNA and protein components (particularly histones). This document uses the commonly accepted, internationally recognized chromosome numbering system for the human genome.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[0352] 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, multiple loci, or an entire chromosome, depending on the number of recombination events (if any) that have occurred between a given set of loci.

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

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

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

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

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

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

[0359] Sequencing process

[0360] The embodiments described herein can be applied to analyzing nucleic acid sequences to identify sequence variations. Implementations can be used to analyze potential variants / alleles of a genetic location / locus and determine the genotype of the genetic locus, or in other words, provide genotype identification for the locus. For example, nucleic acid sequences can be analyzed according to the methods and systems described in U.S. Patent Application Publication No. 2016 / 0085910 and U.S. Patent Application Publication No. 2013 / 0296175, the entire subject matter of which is expressly incorporated herein by reference in its entirety.

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

[0362] Next, one or more processor guidance systems are ordered by the sample. Order-checking can be carried out by a variety of known sequencing schemes. In a specific implementation, order-checking includes SBS. In SBS, a plurality of fluorescently labeled nucleotides are used to order-check the DNA clusters (possibly millions of clusters) of a plurality of amplifications present on an optical substrate surface (e.g., at least partially limiting the surface of the channel in the flow cell). The flow cell can comprise a nucleic acid sample for order-checking, wherein the flow cell is placed in a suitable flow cell support.

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

[0364] In some embodiments, nucleic acid can be attached to a surface and amplified before or during sequencing. For example, bridge amplification can be used to amplify to form nucleic acid clusters on the 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 surfaces is rolling circle amplification (RCA), for example, as described in Lizardi et al., Nat. Genet. 19: 225-232 (1998), and U.S. Patent Application Publication No. 2007 / 0099208 A1, each of which is incorporated herein by reference in its entirety.

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

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

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

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

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

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

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

[0372] (1) Pisces Interpretation 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.

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

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

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

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

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

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

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

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

[0381] d.VEP cache file

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

[0383] f. Supplementary Database: The data sources for SAdb are listed under Supplementary Materials. A small variant of SAdb is produced by performing a k-way merge of all data sources so that each object in the database (identified by a reference name and location) has all associated supplementary annotations. The issues encountered when parsing the 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 lookups of file locations to retrieve supplementary annotations. However, since the data must be fetched from disk, adding supplementary annotations has been identified as Nirvana's largest bottleneck (configured to be 30% of the total runtime).

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

[0385] Such variant annotation tools can include preprocessing. For example, Nirvana includes a large number of annotations from external data sources such as ExAC, EVS, the 1000 Genomes Project, dbSNP, ClinVar, Cosmic, DGV, and ClinGen. To fully utilize these databases, we must clean the information from them. We implemented different strategies to deal with different conflicts from different data sources. For example, in the case of multiple dbSNP entries for the same position and alternating alleles, we merged all IDs into a comma-separated list of IDs; if there were multiple entries with different CAF values ​​for the same allele, the first CAF value was used. For conflicting ExAC and EVS entries, we took into account the number of sample counts and used the entry with the higher sample count. In the 1000 Genomes Project, we removed the allele frequencies of conflicting alleles. Another problem is incorrect information. We primarily extracted allele frequency information from the 1000 Genomes Project. However, we noted that for GRCh38, the allele frequencies reported in the info field do not exclude samples for which genotypes are unavailable, resulting in lower frequencies for variants that are not available for all samples. To ensure annotation accuracy, we used all individual-level genotypes to calculate true allele frequencies. It is well known that the same variant can have different representations based on different alignments. To ensure accurate reporting of information for identified variants, we preprocessed variants from different resources to achieve consistent representation. For all external data sources, we trimmed alleles to remove duplicate nucleotides in the reference and alternative alleles. For ClinVar, we directly parsed the XML files and performed a five-prime alignment on all variants, which is commonly used in VCF files. Different databases can contain the same set of information. To avoid unnecessary duplication, we removed some duplicate information. For example, we removed variants from the DGV with the 1000 Genomes Project data source because we already report these variants in 1000 Genomes and provide more detailed information.

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

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

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

[0389] Next, one or more processors generate and store sample reports.Sample report can include, for example, information about a plurality of genetic loci about a sample.For example, for each genetic locus of a predetermined group of genetic loci, sample report can at least provide one of the following: genotype identification is provided; Indication can not carry out genotype identification; The confidence score of the certainty of genotype identification is provided; Or point out the possible problems of the determination method of relevant one or more genetic loci.Sample report can also indicate the individual's sex and / or indicator sample that provides sample and include multiple sources.As used herein, " sample report " can include the digital data (for example, data file) of the predetermined set of genetic locus or genetic locus and / or the printed report of the set of genetic locus or genetic locus.Therefore, generating or providing can include creating data file and / or printing sample report, or displaying sample report.

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

[0391] Variant Identification Application

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

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

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

[0395] If sample is interpreted one of the selection sequences in the effective matching database (that is, fully matched or almost matched as described above), then sample is interpreted and assigned to or assigned to the genetic locus relevant to the selection sequence. This can be called locus identification or temporary locus identification, wherein interprets for the genetic locus identification sample relevant to the selection sequence. However, as mentioned above, it may be necessary to carry out sample interpretation more than one locus. In such implementation, further analysis can be performed to identify or distribute the sample interpretation only for one of potential genetic locus. In some implementations, the sample interpretation compared with the reference sequence database is the first interpretation from paired end sequencing. When performing paired end sequencing, the second interpretation (representing the original fragment) relevant to the sample interpretation will be obtained. After the distribution, the subsequent analysis performed using the interpretation distributed can be based on the type of the genetic locus identified for the interpretation distributed.

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

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

[0398] Next, a count of raw fragments that include supported variants is maintained, as well as the locations where the supported variants occur. Additionally or alternatively, a count of raw fragments that do not include supported variants at the location of interest (relative to the sample interpretation or the location of potential variant identification in the sample fragment) can be maintained. Additionally or alternatively, a count of raw fragments that correspond to the reference sequence and do not authenticate or confirm potential variant identification can be maintained. The determined information will be output to the variant identification verification application, including the count and type of raw fragments that support potential variant identification, the location of the supported variant in the raw fragment, the count of raw fragments that do not support potential variant identification, etc.

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

[0400] Benign training set generation

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

[0402] The clinical applicability of diagnostic sequencing is limited by the difficulty in interpreting rare genetic variants in human populations and inferring their impact on disease risk. Because rare genetic variants have deleterious health effects, clinically significant genetic variants are extremely rare in human populations, and for the vast majority of rare genetic variants, their impact on human health is uncertain. This large number of rare variants of uncertain clinical significance severely hinders the use of sequencing for personalized medicine and population-wide health screening.

[0403] Most penetrant Mendelian disorders have low population prevalence, so observing a variant at high frequency in a population is strong evidence for a favorable outcome. Assessing common variation across diverse human populations is an effective strategy for cataloging benign variants, but the total number of common variants in present-day humans is limited due to a bottleneck in our species' recent history, during which much of the original diversity was lost. Studies of present-day human populations have increased significantly from effective population sizes (Ne) of fewer than 10,000 individuals over the past 15,000–65,000 years. A small set of common polymorphisms accounts for the limited capacity for variation in populations of this size. Of the over 700 million potential protein-altering missense substitutions in the reference genome, only approximately one in a thousand represents an allele frequency greater than 0.1% of the total population.

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

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

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

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

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

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

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

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

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

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

[0414] To confirm that evolutionary distance, rather than domestication artifacts, is the primary driver of diverging selection coefficients, we reran our analyses using fixed substitutions between pairs of closely related individuals that replace intraspecific polymorphisms across a wide range of evolutionary distances (as shown in Figure 50D, Supplementary Table 4, and Supplementary Data File 2). We found that the attrition of human missense variants consistent with the state of fixed substitutions within species increased with branch length and was not significantly different in wild species compared to those in domesticated species. Consistent with earlier work in flies and yeast, the number of fixed missense substitutions with consistent state was lower than expected for diverging lineages.

[0415] Deep learning network for variant pathogenicity classification

[0416] The present disclosure provides a deep learning network for variant pathogenicity classification. Based on the clinical application significance of variant classification, a large number of attempts have been made to use supervised machine learning to solve this problem. However, these efforts have been hampered by the lack of a real dataset of suitable size that includes variants clearly labeled as benign and pathogenic for training.

[0417] Existing databases of human-curated variants do not represent the entire genome; 50% of the variants in the ClinVar database come from just 200 genes (1% of human protein-coding genes). Furthermore, system learning identifies most human annotations as having questionable supporting evidence, underestimating the difficulty of interpreting rare variants observed in only a single patient. While human interpretation is becoming increasingly rigorous, numerous classification guidelines have been developed around consensus practices, which risk reinforcing existing trends. To mitigate human interpretation bias, recent classifiers are trained on common human polymorphisms or fixed human-chimpanzee surrogates, which also serve as their input: the prediction scores of earlier classifiers trained on human-curated databases. In the absence of independent, unbiased real-world datasets, objective benchmarks of the performance of these different approaches are difficult to come by.

[0418] Variants from six non-human primates (chimpanzee, bonobo, gorilla, orangutan, rhesus macaque, and marmoset) provide over 300,000 unique missense variants that do not overlap with common human variants. These variants largely represent common variants that have been screened for benign outcomes through purifying selection, significantly expanding the training dataset available for machine learning methods. In general, after excluding variants of uncertain significance and variants with conflicting annotations, each primate species provides more variants than the entire ClinVar database (42,000 as of November 2017). Furthermore, this content is free of human interpretation bias.

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

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

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

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

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

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

[0425] Only with human variant data (e.g. Figure 21ACompared to the accuracy of a network trained on a dataset of primate variants (shown), nearly half of the improvement over existing methods stems from the use of a deep learning network, and half from the expansion of the primate variant training dataset. To test the classification of variants of uncertain significance in a clinical setting, we evaluated the ability of the deep learning network to distinguish de novo mutations in patients with neurodevelopmental disorders compared with healthy controls. Neurodevelopmental disorders are one of the largest categories of rare genetic diseases in terms of prevalence, and recent triple sequencing studies have implicated a central role for de novo missense and protein-truncating mutations.

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

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

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

[0429] Novel candidate gene discovery

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

[0431] Comparison with human expert selections

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

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

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

[0435] The genetic diversity of the human population is much lower than that of nonhuman primate species, with approximately half the number of single nucleotide variants carried by chimpanzees, gorillas, and gibbons, and one-third the number of variants carried by orangutans. Although the level of genetic diversity of most nonhuman primate species is unknown, based on the large number of nonhuman primate species available, we infer that most likely benign human missense positions are covered by common variants in at least one primate species, allowing the identification of pathogenic variants by systematic elimination (e.g., Figure 23DEven when sequencing only a small subset of these species, increasing the size of the training data allows for more accurate prediction of missense outcomes using machine learning. Finally, while our findings focus on missense variants, this strategy can also be applied to inferring the outcomes of noncoding variants, particularly in conserved regulatory regions, where sufficient alignment of the human and primate genomes is required to definitively determine if the variants are concordant.

[0436] Of the 504 known non-human primate species, approximately 60% are threatened with extinction due to illegal poaching and widespread habitat loss. Declining population sizes and potential extinctions represent irreparable losses of genetic diversity, urging urgent conservation efforts worldwide to benefit rare and irreplaceable species, as well as humanity itself.

[0437] Data generation and comparison

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

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

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

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

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

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

[0444]

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

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

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

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

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

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

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

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

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

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

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

[0456] Deep learning network architecture

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

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

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

[0460] Classifier performance on a holdout test set of 10,000 primate variants

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

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

[0463] Analysis of novel variants in DDD studies

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

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

[0466] Novel candidate gene discovery

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

[0468] ClinVar classification accuracy

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

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

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

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

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

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

[0475] In one embodiment, for the practical application of primateAI score, compared with the control example (such as Figure 21A ), based on the enrichment of new variant cases, in variants with dominant inheritance patterns, the preferred threshold is greater than 0.8, which may be pathogenic, less than 0.6, which may be benign, and between 0.6 and 0.8, which may be intermediate. In genes with recessive inheritance patterns, the preferred threshold is greater than 0.7, which may be pathogenic, and less than 0.5, which may be benign.

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

[0477] Figure 3 Schematic diagram of primateAI, a deep learning network architecture for pathogenicity classification. The model inputs include 51 amino acids (AA) flanking the reference sequence and the variant-substituted sequence, three 51-AA-length position weight matrices representing conservation from primate, mammalian, and vertebrate alignments, and pretrained secondary structure and solvent accessibility network inputs (also 51 AA in length).

[0478] Figure 4A 、 4B Tables 4 and 4C are Supplementary Table 16, showing the architecture details of an example deep learning model for pathogenicity prediction, primateAI. Shape specifies the shape of the output tensor for each layer of the model, and activation specifies the activations given to the neurons in that layer. Model inputs include a position-specific frequency matrix of the flanking amino acid sequence surrounding the variant (51 AA long, 20 deep), a one-hot-encoded human reference sequence and variable sequence (51 AA long, 20 deep), and a secondary structure and solvent accessibility model (51 AA long, 40 deep).

[0479] The example uses one-dimensional convolution. In other embodiments, the model can use different types of convolution, such as 2D convolution, 3D convolution, dilated or dilated convolution, transposed convolution, separable convolution, and depthwise convolution.

[0480] Certain layers also use batch normalization (Ioffe and Szegedy 2015). In batch normalization, during training, the distribution of each layer in a convolutional neural network (CNN) is altered, making the distribution different from one layer to another. This slows down the convergence of the optimization algorithm. Batch normalization is a technique that addresses this issue. The input to a batch normalization layer is annotated with x and its output with z. Batch normalization applies the following transformation to x:

[0481]

[0482] Batch Normalization uses μ and σ to mean-variance normalize the inputs and uses γ and β to linearly scale and shift the inputs. The normalization parameters μ and σ for the current layer of the training set are calculated using an exponential moving average. In other words, these parameters are not trainable. In contrast, gamma and beta are trainable. During training, the values ​​of μ and σ are used in the forward method during inference.

[0483] Figure 5 and 6 A deep learning network for predicting protein secondary structure and solvent accessibility was demonstrated. The model inputs were either a conserved position-weighted matrix generated by RaptorX software (for training on protein database sequences) or 99 vertebrate alignments (for training and perturbation on human protein sequences). The outputs of the second to final layers, each 51 aa, served as input to a deep learning network for pathogenicity classification.

[0484] Figure 7A and 7B Table 11 is a supplementary table showing the details of the example model architecture for a deep learning (DL) model for 3-state secondary structure prediction. The shape specifies the output tensor of each layer of the model.

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

[0486] Figure 20 Figure 2 shows the predicted pathogenicity score for each amino acid position in the SCN2A gene, with annotations for key functional domains. The bottleneck primateAI score for missense substitutions at each amino acid position is plotted along the gene.

[0487] In one embodiment, for the practical application of primateAI score, compared with the control example (such as Figure 21A ), based on the enrichment of new variant cases, in variants with dominant inheritance patterns, the preferred threshold is greater than 0.8, which may be pathogenic, less than 0.6, which may be benign, and between 0.6 and 0.8, which may be intermediate. In genes with recessive inheritance patterns, the preferred threshold is greater than 0.7, which may be pathogenic, and less than 0.5, which may be benign.

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

[0489] Figure 3 Schematic diagram of primateAI, a deep learning network architecture for pathogenicity classification. The model inputs include 51 amino acids (AA) flanking the reference sequence and the variant-substituted sequence, three 51-AA-length position weight matrices representing conservation from primate, mammalian, and vertebrate alignments, and pretrained secondary structure and solvent accessibility network inputs (also 51 AA in length).

[0490] Figure 4A 、 4B Tables 4 and 4C are Supplementary Table 16, showing the architecture details of an example deep learning model for pathogenicity prediction, primateAI. Shape specifies the shape of the output tensor for each layer of the model, and activation specifies the activations given to the neurons in that layer. Model inputs include a position-specific frequency matrix of the flanking amino acid sequence surrounding the variant (51 AA long, 20 deep), a one-hot-encoded human reference sequence and variable sequence (51 AA long, 20 deep), and a secondary structure and solvent accessibility model (51 AA long, 40 deep).

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

[0492] Certain layers also use batch normalization (Ioffe and Szegedy 2015). In batch normalization, during training, the distribution of each layer in a convolutional neural network (CNN) is altered, making the distribution different from one layer to another. This slows down the convergence of the optimization algorithm. Batch normalization is a technique that addresses this issue. The input to a batch normalization layer is annotated with x and its output with z. Batch normalization applies the following transformation to x:

[0493]

[0494] Batch Normalization uses μ and σ to mean-variance normalize the inputs and uses γ and β to linearly scale and shift the inputs. The normalization parameters μ and σ for the current layer of the training set are calculated using an exponential moving average. In other words, these parameters are not trainable. In contrast, gamma and beta are trainable. During training, the values ​​of μ and σ are used in the forward method during inference.

[0495] Figure 5 and 6 A deep learning network for predicting protein secondary structure and solvent accessibility was demonstrated. The model inputs were either a conserved position-weighted matrix generated by RaptorX software (for training on protein database sequences) or 99 vertebrate alignments (for training and perturbation on human protein sequences). The outputs of the second to final layers, each 51 aa, served as input to a deep learning network for pathogenicity classification.

[0496] Figure 7A and 7B Table 11 is supplemented with details of an example model architecture for a deep learning (DL) model for 3-state secondary structure prediction. Shape specifies the shape of the output tensor for each layer of the model, and Activation is the activation given to the neurons in that layer. The model input is a position-specific frequency matrix (55 AA long, 20 deep) for the flanking amino acid sequence surrounding the variant.

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

[0498] Figure 20 Figure 2 shows the predicted pathogenicity score for each amino acid position in the SCN2A gene, with annotations for key functional domains. The bottleneck primateAI score for missense substitutions at each amino acid position is plotted along the gene.

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

[0500] Figure 21B The distribution of primateAI prediction scores for de novo missense variants in individuals with developmental disorders decipherable (DDD) compared to unaffected individuals with corresponding Wilcoxon rank sum values ​​P is illustrated.

[0501] Figure 21C Figure 2. Comparison of classifiers for separation of de novo missense variants in DDD cases compared to controls. Wilcoxon rank sum test P values ​​for each classifier are shown.

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

[0503] Figure 22C Figure 5. Depiction of the classifiers comparing separation of novel missense variants in individuals with DDD versus controls across 605 genes. The y-axis represents the Wilcoxon rank sum test P value for each classifier.

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

[0505] Figure 22E The classification accuracy and AUC of each classifier are explained. Figure 22A The enrichment shown indicates that the classifier predicts an equal number of pathogenic and benign variants as expected. The classification accuracy shown is the average of the true positive rate and true negative rate. The dashed line indicates the maximum achievable AUC for a perfect classifier, considering that 33% of the de novo missense variants in DDD represent background.

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

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

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

[0509] Figure 24 Correction for the effect of sequencing depth on identifying common primate variants has been demonstrated. The probability of observing a given variant in a non-human primate species is negatively correlated with the sequencing depth at that position in the ExAC / gnomAD exome dataset. In contrast, lower gnomAD read depth does not affect the probability of observing a common human variant at that position (less than 0.1% allele frequency), because the large number of sequenced human exomes ensures that most common variants are identified. When selecting matching variants for each primate variant used in training the network, the probability of selecting the variant is adjusted for the effect of sequencing depth, and matching is performed on the trinucleotide context to control for mutation rate and gene conversion.

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

[0511] Figure 26 This is a line graph showing the effect of each position in and around the variant on the deep learning prediction score. We systematically zeroed the inputs for the adjacent amino acids surrounding the variant (positions -25 to +25) and measured the change in the neural network's predicted pathogenicity. The graph shows the average change in predicted pathogenicity score for 5,000 randomly selected variants at each adjacent amino acid position.

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

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

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

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

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

[0517] Figure 30A 、 30B Figures 30A and 30C compare the accuracy of the primateAI network and other classifiers in separating pathogenic and benign variants in 605 disease-associated genes. The scatter plot in Figure 30A shows the performance of each classifier on DDD cases versus controls (y-axis) and the accuracy of benign prediction on the retained primate dataset (x-axis). Figure 30BDifferent classifiers were compared for separating de novo missense variants in cases versus controls in 605 genes, as represented by receiver operating characteristic (ROC) curves, with the area under the curve (AUC) representing each classifier. Figure 30C The classification accuracy and AUC of the primateAI network are shown in Figures 28A, 28B, and 28C. 20 classifiers are listed. Figure 22A The enrichment shown indicates that the classifier predicts an equal number of pathogenic and benign variants. The classification accuracy shown is the average of the true positive rate and true negative rate. Assuming that 67% of de novo missense variants in DDD cases are pathogenic and 33% are diamond variants, and that 100% of de novo missense variants in controls are benign, the dashed line represents the maximum achievable AUC for a perfect classifier.

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

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

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

[0521] Figure 33 Supplementary Table 15 shows the performance comparison of deep learning networks using secondary structure markers annotated from the DSSP database and deep learning networks using predicted secondary structure markers.

[0522] Figure 34 Table 17 is supplemented with accuracy values ​​for the 10,000 retained primate variants evaluated and P values ​​for novel variants in DDD cases versus controls for each of the 20 classifiers. The primateAI model with only human data was a deep learning network trained using a labeled benign training dataset consisting of only common human variants (83.5K variants representing the 0.1% of the population), while the primateAI model with both human and primate data was a deep learning network trained on the full set of 385K labeled benign variants, including both common human and primate variants.

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

[0524] Figure 49A 、 49B , 49C, 49D, and 49E describe the missense / synonymous mutation rates across the human allele spectrum. Figure 49A This figure shows missense and synonymous variants observed in 123,136 human individuals in the ExAC / gnomAD database, grouped into four categories by allele frequency. Shaded gray bars indicate the number of synonymous variants in each category, while dark green bars indicate missense variants. The height of each bar is reduced to the number of synonymous variants in each category by allele frequency. After adjusting for mutation rate, the numbers of missense / synonymous variants and the missense / synonymous mutation rates are displayed. Figure 49B and 49C Describes variants common in chimpanzees (e.g. Figure 49B ) and chimpanzee single variants (e.g. Figure 49C Allelic spectrum of human missense and synonymous variants with consistent IBS (IBS). Red boxes indicate the depletion of chimpanzee missense variants at common allele frequencies (greater than 0.1%) relative to rare human allele frequencies (less than 0.1%) and also indicate the chi-squared (χ2) test value P.

[0525] Figure 49DRepresents human variants observed in at least one non-human primate species. Relative to ClinVar variants in a panel of 30 human samples sampled from ExAC / gnomAD allele frequencies (middle row) and variants observed in primates (bottom row), Figure 49E The number of benign and pathogenic missense variants in the entire ClinVar database (top row) is shown. Conflicting benign and pathogenic variants identified and variants annotated only as undetermined sense were excluded.

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

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

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

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

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

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

[0532] Figure 53 、 54 The missense:synonymous mutation rates of human variants consistent with the six primate states were elucidated in Figure 5 and 55. For ExAC / gnomAD variants consistent with the variant states, the missense:synonymous mutation rate pattern of the human allele spectrum was present in chimpanzees, bonobos, gorillas, orangutans, rhesus monkeys, and marmosets.

[0533] Figure 56 We simulated the saturation of new common missense variants discovered by increasing the size of the human cohort investigated. In the simulations, the genotypes of each sample were sampled according to the gnomAD allele frequencies. The fraction of discovered gnomAD common variants was averaged across 100 simulations for each sample size ranging from 10 to 10,000 variants.

[0534] Figure 57Figure 2. PrimateAI accuracy across different conservation profiles in the genome. The x-axis represents the percent alignment of the surrounding 51 amino acids with 99 vertebrate alignments. The y-axis represents the classification performance of PrimateAI accuracy for variants within each conserved region, based on a test dataset of 10,000 conserved primate variants.

[0535] Figure 58 Table 5 is supplemented with the impact of a benign training dataset on the labeling of common human variants and variants present in non-human primates.

[0536] Figure 59 The effect of allele frequency on the expected missense:synonymous mutation rate is shown in Supplementary Table 8. The expected number of synonymous and missense variants was calculated based on the allele frequency spectrum of variants within the intragenic region at least 20–30 nm away using the trinucleotide context for mutation rate and gene conversion bias.

[0537] Figure 60 ClinVar analysis is shown in Supplementary Table 9. According to one embodiment, missense variants with conflicting annotations and variants of uncertain significance were filtered out from variants downloaded from the October 2017 build of the ClinVar database, leaving 17,775 benign variants and 24,853 pathogenic variants.

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

[0539] Figure 62 Table 1 shows one embodiment of 14 additional candidate genes found in intellectual disability that did not previously meet the genome-wide significance threshold in the original DDD study.

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

[0541] Data Generation

[0542] All coordinates used in this article refer to the human genome construct UCSC hg19 / GRCh37, including coordinates for variants in other species, and were mapped to hg19 using multiple sequence alignments using the methods described in this section. Human protein-coding DNA sequences and multiple sequence alignments with 99 vertebrate genomes were downloaded from the UCSC Genome Browser for use in the hg19 construct (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.

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

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

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

[0546] To map variants to humans, we used a multispecies alignment of 99 vertebrates to ensure 1:1 mapping with orthologous human protein-coding regions. Mapping variants using orthologous multispecies alignments is crucial to eliminate artifacts caused by pseudogenes or retrotransposons, which occur when mapping SNPs directly between species using tools such as liftOver. If the genome structure of the species in dbSNP does not match that of the species in the 99 vertebrate multiple sequence alignment, we use liftOver to update the variant to the genome structure used in the multiple sequence alignment. If the variant appears in a reference / alternative manner, we accept the variant as state-identical. For example, if the human reference sequence has a G and the alternative allele has an A, the variant is considered to have the same state in other species with an A in the reference sequence and a G in the alternative allele. To ensure that variants have the same predicted protein-coding outcome in humans and other species, we require that the other two nucleotides in the codon be identical across species for both missense and synonymous variants. The polymorphisms for each species included in the analysis are listed in Supplementary Data File 1 and detailed metrics are shown in Supplementary Table 1 .

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

[0548] Corrected allele frequency for missense / synonymous ratios, mutation rates, genetic drift, and GC-based genes The impact of conversion

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

[0550] To calculate the expected missense / synonymous mutation ratio within 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 sufficiently distant to avoid the influence of extended splicing motifs. Because these regions are located near the edge of the exome capture sequence of the ExAC / gnomAD exomes, we removed all chrX regions to ensure unambiguous variant identification and excluded regions with an average read depth of <30. Each variant and its immediately upstream and downstream nucleotides belonged to one of 64 trinucleotide contexts. If we mutated the intervening nucleotide to three other bases, a total of 64 × 3 = 192 trinucleotide configurations were possible. Because trinucleotide configurations and their reverse complements are equivalent, 96 trinucleotide contexts were valid. We observed a very strong effect of trinucleotide context on mutation rate, with a smaller effect on GC-biased gene conversion, making trinucleotide context a useful tool for constructing these variants.

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

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

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

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

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

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

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

[0558]

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

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

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

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

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

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

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

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

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

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

[0569] Generation of benign variants for model training

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

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

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

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

[0574] Generate unlabeled variants to supplement benign training datasets

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

[0576] Filtering of unlabeled variants

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

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

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

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

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

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

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

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

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

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

[0587] Prediction of secondary structure and solvent accessibility

[0588] The deep learning network for pathogenicity prediction contains 36 total convolutional layers, 19 of which are 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 structure of most human proteins is unknown, we trained two models to enable the network to learn protein structure from primary sequence. Figure 6 Both models use the same network architecture and input. The input to the secondary structure and solvent accessibility networks is a 51-by-20 amino acid position frequency matrix encoding conservation information from a multiple sequence alignment of humans with 99 other vertebrates.

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

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

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

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

[0593] %. / buildFeature-i 1u7lA.fasta-c 10-o. / TGT / 1u7lA.tgt

[0594] %. / CNFsearch -a 30 -q 1u7lA

[0595] For each amino acid position in the dataset, we extracted a window from the position frequency matrix corresponding to the 51 flanking amino acids and used it to predict a secondary structure or solvent accessibility signature for the central amino acid within the 51-amino acid sequence. Using DSSP software, we derived secondary structure and relative solvent accessibility signatures directly from the known 3D crystal structure of the protein, without the need for prediction from the primary sequence. To incorporate the secondary structure and solvent accessibility networks into the pathogenicity prediction network, we calculated a position frequency matrix from a human-based multiple sequence alignment of 99 vertebrates. Although the conservation matrices generated from these two methods were generally similar, we backpropagated them through the secondary structure and solvent accessibility models during the training process for pathogenicity prediction, allowing for fine-tuning of parameter weights.

[0596] Model architecture and training

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

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

[0599] The best test accuracy of 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 30. The best test accuracy of 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 the predictions of the neural network using DSSP-annotated structural labels of approximately 4000 human proteins with crystal structures with the standard PrimateAI model using only predicted structural labels. We did not improve the accuracy of pathogenicity predictions when using DSSP-annotated labels (Supplementary Table 15).

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

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

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

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

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

[0605] Semi-supervised learning

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

[0607] Another common challenge in bioinformatics is data imbalance. This occurs when the class to be predicted is underrepresented in the data, either because instances belonging to that class are rare (noticeable cases) or difficult to obtain. Ironically, those rare classes are often the most important to learn from, as they may be associated with unique situations.

[0608] One algorithm for handling imbalanced data distributions is based on classifier ensembles. Limited amounts of labeled data naturally lead to weaker classifiers, but an ensemble of weaker classifiers often surpasses the performance of any individual component classifier. Furthermore, an ensemble of classifiers often improves the predictive accuracy obtained from a single classifier by a factor that justifies the effort and cost associated with learning multiple models. Intuitively, aggregating several classifiers should result in better overfitting control because averaging the high variability of the individual classifiers also averages out their overfitting.

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

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

[0611] We trained a deep learning network to distinguish between a carefully matched set of labeled benign and unlabeled variants to eliminate bias. In 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 primates.

[0612] We sampled a set of unlabeled variants that needed to be matched with benign variants on a three-nucleotide background (to control for mutation rate, genetic drift, and gene conversion) and adjusted for the effects of comparability and sequence coverage on variant confirmation. Because the number of unlabeled variants greatly exceeds the number of labeled benign variants, we obtained consistent prediction results by training eight models using the same labeled benign variants and eight randomly sampled unlabeled variants and taking the average of these models.

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

[0614] To reduce sampling bias, we trained an ensemble of eight models using the same set of benign training variants and eight different sets of pathogenic variants. Initially, we randomly sampled unknown variants to represent pathogenic variants. Next, the ensemble of models was used to iteratively score a set of unknown variants that had not participated in the previous training cycle. The highest-scoring pathogenic variants were then used to replace 5% of the random unknown variants from the previous cycle. Note that we retained 25% more top-scoring pathogenic variants than needed, allowing us to sample eight sets of pathogenic variants with different scores to replace the 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 initially randomly sampled unknown variants had been replaced by high-confidence pathogenic variants predicted by the ensemble model. Figure 42 The iterative equalization sampling process is explained.

[0615] Balancing benign and unknown training sets

[0616] Sampling unknown variants that match benign variants reduces bias in our model training. When randomly sampling unknown variants, deep learning models often extract biased information and produce simplistic solutions. For example, if the amino acid substitution K->M occurs more frequently in unknown variants than in benign variants, deep learning models will tend to classify the K->M substitution as pathogenic. Therefore, it is important to balance the distribution of amino acid substitutions between the two training sets.

[0617] Highly expressed mutation types, such as CpG transitions, have a significant representation bias among common benign variants. Orthologous variants from other primates also follow human mutation rates, implying an enrichment of highly expressed mutation types throughout the benign training set. If the sampling process for unknown variants is not well controlled and balanced, deep learning models may be more likely to classify CpG transitions as benign than less well-represented classes, such as translocations or non-CpG transitions.

[0618] To prevent the deep learning model from converging to a trivial, non-biological solution, we consider a trinucleotide context that balances benign and unknown variants. A trinucleotide consists of the base preceding the variant, the reference base for the variant, and the base following the variant. The reference base for the variant can be changed to any of the other three nucleotides. In total, there are 64 x 3 trinucleotide contexts.

[0619] Iterative Equalization Sampling

[0620] Cycle 1

[0621] We sampled unknown variants to match the exact number of benign variants in each three-nucleotide context. In other words, in the first cycle, we reflected both the benign and pathogenic training sets from the three-nucleotide context of the variants. This sampling strategy is based on the fact that variants with the same mutation rate between the benign and unknown groups are represented identically. This prevents the model from converging to a trivial solution based on mutation rate.

[0622] Cycles 2-20

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

[0624] Objectively, because cycle 1 used a fully matched training set, the predicted most pathogenic variants generated were free of any mutation rate bias. Therefore, no bias in this group was considered. The remaining 95% of the data was still controlled by the trinucleotide background mutation rate to prevent the model from converging to a trivial solution.

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

[0626] Cycle 21

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

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

[0629] Balanced sampling of unlabeled variants is crucial to eliminate biases unrelated to the deleteriousness of a variant. Without proper controls for confounding effects, deep learning can easily inadvertently introduce biases that affect discrimination between classes. Variants commonly found in humans tend to be enriched for variants from highly mutating classes, such as those located in CpG promoter regions. Similarly, primate polymorphisms follow human mutation rates, meaning that variants from highly mutating classes are enriched throughout the benign training set. If the sampling of unlabeled variants is not well controlled and balanced, deep learning networks tend to rely on mutation rate biases to classify variants, making them more likely to classify CpG transitions as benign than less-represented classes, such as translocations or non-CpG transitions. We sampled exactly the same number of unlabeled variants as labeled benign variants in each of the 96 trinucleotide contexts (discussed previously).

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

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

[0632] Generation of benign and unknown variants

[0633] Common variants in the human population

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

[0635] Common polymorphisms in great apes

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

[0637] Generation of unknown variants

[0638] By replacing the nucleotide at each base position with three other nucleotides, all possible variants can be generated from each base position in the canonical coding region. When a new codon is formed, it potentially changes the amino acid at that position. Synonymous mutations are filtered out.

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

[0640] Filtering variants

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

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

[0643] Validation and test data

[0644] To validate and test the pathogenicity model, according to one embodiment, we randomly sampled two sets of 10,000 benign variants from a large pool of benign variants for validation and testing. The remainder of the benign variants was used to train the deep learning model. These variants were specifically sampled from orthologous primate variants to ensure fair comparison between methods, as is the case when training on common human variants. According to one embodiment, we also randomly sampled two sets of 10,000 unknown variants for validation and testing, respectively. We ensured that the number of unknown variants in the 192 trinucleotide context matched the number of benign variants used for validation and testing, respectively.

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

[0646] Deep learning network architecture

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

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

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

[0650] Model Overview

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

[0652] The use of residual blocks

[0653] Figure 17The residual block is shown. Our deep learning models for pathogenicity prediction and for predicting secondary structure and solvent accessibility both employ the definition of the residual block, which was first described in this paper. The structure of the residual block is shown in the figure below. The input layer is first batch normalized and then normalized using the nonlinear activation "ReLU". After the activation, it passes through a 1D convolution layer. The intermediate layer output from the 1D convolution layer is again batch normalized and activated with ReLU, and then another 1D convolution layer is performed. At the end of the second 1D convolution, its output is merged with the original input. In this architecture, the input remains in its original state, and the residual connection is not affected by the nonlinear activation of the model.

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

[0655] Novelty of our model

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

[0657] Training secondary structure and solvent accessibility models

[0658] Data Preparation

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

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

[0661] Model architecture and training

[0662] We trained two end-to-end deep convolutional neural network models to predict the three-state secondary structure and three-state solvent accessibility of proteins, respectively. These two models have similar configurations, consisting of two input channels: one for the protein sequence and the other for the protein conservation profile. Each input channel has a dimension of Lx20, where L represents the length of the protein.

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

[0664] The activations from layer 2 pass through a series of nine residual blocks (layers 3 to 11) as defined above. The activations from layer 3 are fed back to layer 4, which in turn is fed back to layer 5, and so on. There are also skip connections that simply sum the outputs of every third residual block (layers 5, 8, and 11). The combined activations are then fed into two 1D convolutions with ReLU activations (layers 12 and 13). The activations from layer 13 are given to a softmax readout layer. This softmax computes the probabilities of the three output classes for a given input.

[0665] For the best secondary structure model, the 1D convolution was implemented with a dilation ratio of 1. For the solvent accessibility model, the final three residual blocks (layers 9, 10, and 11) were implemented with a dilation ratio of 2 to increase kernel coverage. Protein secondary structure strongly depends on interactions between immediately adjacent amino acids. Therefore, models with higher kernel coverage slightly improved performance. On the other hand, solvent accessibility is affected by long-range interactions between amino acids. Consequently, models with high kernel coverage using dilated convolutions achieved over 2% higher accuracy than models with shorter coverage.

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

[0667]

[0668]

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

[0670]

[0671]

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

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

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

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

[0676] Train a model to predict pathogenicity of variants

[0677] Input features of pathogenicity prediction models

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

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

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

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

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

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

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

[0685] Deep learning model training

[0686] Figure 19 A block diagram provides an overview of the deep learning model workflow. The pathogenicity training model includes five direct inputs and four indirect inputs. The five direct input features include the reference sequence (1a), alternative sequence (1b), primate conserved sequence (1c), mammalian conserved sequence (1d), and vertebrate conserved sequence (1e). The indirect inputs include secondary s...

Claims

1. A method of using a variant pathogenicity classifier, comprising: executing a recurrent neural network on at least one processor coupled to the memory, wherein the recurrent neural network is trained to determine pathogenicity of a variant; wherein the recurrent neural network is trained using benign training examples and pathogenic training examples of protein sequence pairs including reference protein sequences and alternative protein sequences as input; wherein the alternative protein sequences are generated from benign variants and pathogenic variants; wherein the benign variants include common human missense variants and non-human primate missense variants that occur at alternative non-human primate codon sequences that share a reference codon sequence that matches a human; and The recurrent neural network is trained to capture long-range dependencies in a reference protein sequence and in an alternative protein sequence.

2. The method according to claim 1, wherein: The recurrent neural network is a long short-term memory (LSTM) network.

3. The method according to claim 2, wherein: The recurrent neural network is a gated recurrent unit (GRU) network.

4. The method according to claim 3, wherein: The recurrent neural network is a multidimensional recurrent neural network.

5. A method of using a variant pathogenicity classifier, comprising: executing a recurrent neural network on at least one processor coupled to the memory, wherein the recurrent neural network is trained to determine pathogenicity of a variant; wherein the recurrent neural network is trained using benign training examples and pathogenic training examples of protein sequence pairs including reference protein sequences and alternative protein sequences as input; wherein the alternative protein sequence is generated from a benign variant and a pathogenic variant; and Wherein the benign variants include common human missense variants and non-human primate missense variants, wherein the non-human primate missense variants occur at alternative non-human primate codon sequences that share a reference codon sequence that matches humans.

6. A system comprising: Memory; one or more processors; a variant pathogenicity classifier, running on at least one processor coupled to the memory and trained to determine pathogenicity of a variant; wherein the variant pathogenicity classifier is trained using benign training examples and pathogenic training examples of protein sequence pairs including a reference protein sequence and an alternative protein sequence as input; wherein the alternative protein sequences are generated from benign variants and pathogenic variants; as well as Wherein the benign variants include common human missense variants and non-human primate missense variants, wherein the non-human primate missense variants occur at alternative non-human primate codon sequences that share a reference codon sequence that matches humans.

7. The system according to claim 6, wherein: The variant pathogenicity classifier is a recurrent neural network.

8. The system according to claim 7, wherein: The variant pathogenicity classifier captures long-range dependencies in the reference protein sequence and in alternative protein sequences.

9. The system according to claim 8, wherein: The variant pathogenicity classifier is a long short-term memory (LSTM) network.

10. The system according to claim 9, wherein: The variant pathogenicity classifier is a gated recurrent unit (GRU) network.

11. The system according to claim 10, wherein: The variant pathogenicity classifier is a multidimensional recurrent neural network.

12. The system according to claim 11, wherein: The variant pathogenicity classifier is a convolutional neural network.

13. The system according to claim 12, wherein: The variant pathogenicity classifier is a deep spatiotemporal neural network.

14. The system according to claim 13, wherein: The variant pathogenicity classifier uses a convolutional autoencoder.

Images