Aberrant splicing detection using convolutional neural networks (CNNS)

Abstract

To construct a convolutional neural network-based classifier for variant classification, and in particular, train a convolutional neural network-based classifier with training data using a backpropagation-based gradient update technique that progressively match outputs of the classifier with corresponding ground truth labels.SOLUTION: The convolutional neural network-based classifier comprises groups of residual blocks, each group of residual blocks being parameterized by the number of convolution filters in the residual blocks, the convolution window size of the residual blocks, and the atrous convolution rate of the residual blocks, where the size of convolution window differs among groups of residual blocks, and the atrous convolution rate differs among groups of residual blocks. The training data includes benign training examples and pathogenic training examples of translated sequence pairs generated from benign variants and pathogenic variants.SELECTED DRAWING: Figure 31

Description

[Technical Field]

[0001] Additional notes This appendix contains a bibliography of possibly relevant references listed in the paper by the inventors. The subject matter of this paper is described in U.S. provisional applications to which this application claims priority / benefit. These references may be consulted by attorney upon request or via the Global Dossier.

[0002] Priority Application This application is a continuation of U.S. Provisional Patent Application No. 62 / 573,125, filed October 16, 2017, by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae, entitled "Deep Learning-Based Splice Site Classification" (Docket No. ILLM 1001-1 / IP-1610-PRV), U.S. Provisional Patent Application No. 62 / 573,131, filed October 16, 2017, by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae, entitled "Deep Learning-Based Aberrant Splicing Detection" (Docket No. ILLM 1001-2 / IP-1614 ... This application claims priority to or benefit of U.S. Provisional Patent Application No. 62 / 573,135, entitled "Aberrant Splicing Detection Using Convolutional Neural Networks (CNNs)," by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Docket No. ILLM 1001-3 / IP-1615-PRV), and U.S. Provisional Patent Application No. 62 / 726,158, entitled "Predicting Splicing from Primary Sequence with Deep Learning," by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Docket No. ILLM 1001-10 / IP-1749-PRV), filed August 31, 2018. The provisional applications are incorporated herein by reference for all purposes.

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

[0004] PCT Patent Application No. PCT / US18 / 55195, entitled "Deep Learning-Based Splice Site Classification," by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae, filed October 15, 2018 (Docket No. ILLM 1001-7 / IP-1610-PCT), subsequently published as PCT Publication No. WO2019 / 79198.

[0005] PCT Patent Application No. PCT / US18 / 55919, entitled "Deep Learning-Based Aberrant Splicing Detection," by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae, filed October 15, 2018 (Docket No. ILLM 1001-8 / IP-1614-PCT), subsequently published as PCT Publication No. WO2019 / 79200.

[0006] A concurrently filed U.S. non-provisional patent application entitled "Deep Learning-Based Splice Site Classification" by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Docket No. ILLM 1001-4 / IP-1610-US).

[0007] A concurrently filed U.S. non-provisional patent application entitled "Deep Learning-Based Aberrant Splicing Detection" by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Docket No. ILLM 1001-5 / IP-1614-US).

[0008] A concurrently filed U.S. non-provisional patent application entitled "Aberrant Splicing Detection Using Convolutional Neural Networks (CNNs)" by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Docket No. ILLM 1001-6 / IP-1615-US).

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

[0010] Reference 2-SO 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,

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

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

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

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

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

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

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

[0018] Reference 10-LC Piqueras, “AUTOREGRESSIVE MODEL BASED ON A DEEP CONVOLUTIONAL NEURAL NETWORK FOR AUDIO GENERATION”, Tampere University of Technology, 2016,

[0019] Reference 11-J. Wu, “Introduction to Convolutional Neural Networks”, Nanjing University, 2017,

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

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

[0022] Literature 1 describes a deep convolutional neural network architecture that accepts an input sequence and generates an output sequence that scores entries in the input sequence using a group of residual blocks including convolution filters with the same convolution window size, a batch normalization layer, a rectified linear unit (ReLU) layer, a dimensionality transformation layer, an atrous convolution layer with an exponentially increasing atrous convolution rate, skip connections, and a softmax classification layer. The disclosed technology uses the neural network components and parameters described in Literature 1. In one implementation, the disclosed technology modifies parameters of the neural network components described in Literature 1. For example, unlike Literature 1, the atrous convolution rate in the disclosed technology increases non-exponentially from a lower residual block group to a higher residual block group. In another example, unlike Literature 1, the convolution window size in the disclosed technology varies between groups of residual blocks.

[0023] Reference 2 provides a detailed description of the deep convolutional neural network architecture described in Reference 1.

[0024] Reference 3 describes atrous convolution used by the disclosed technology. In this specification, atrous convolution is also referred to as "dilated convolution." Atrous / dilated convolution enables large receptive fields with few trainable parameters. Atrous / dilated convolution is a convolution in which the kernel is applied over a region larger than its length by skipping input values ​​using a step, also called the atrous convolution rate or dilation factor. Atrous / dilated convolution adds spacing between elements of the convolution filter / kernel, allowing nearby input entries (e.g., nucleotides, amino acids) at a larger interval to be considered when the convolution operation is performed. This allows long-range compositional dependencies to be incorporated into the input. Atrous convolution saves partial convolution calculations so that they can be reused when neighboring nucleotides are processed.

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

[0026] Reference 5 describes the skip connection used by the disclosed technology. In this specification, the skip connection is also referred to as a "highway network."

[0027] Reference 6 describes the tightly coupled convolutional network architecture used by the disclosed technique.

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

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

[0030] Reference 9 also describes the atrous / dilated convolution used by the disclosed technique.

[0031] Reference 10 describes various architectures of deep neural networks that can be used with the disclosed techniques, including convolutional neural networks, deep convolutional neural networks, and deep convolutional neural networks with atrous / dilated convolutions.

[0032] Reference 11 provides a detailed description of convolutional neural networks that can be used with the disclosed techniques, including algorithms for training convolutional neural networks that include subsampling layers (e.g., pooling) and fully connected layers.

[0033] Reference 12 provides a detailed description of various convolution operations that can be used with the disclosed technique.

[0034] Reference 13 describes various architectures of convolutional neural networks that can be used with the disclosed technique.

[0035] Incorporation by reference table submitted electronically with the application Three table files in ASCII text format have been submitted with this application and are incorporated by reference. The file names, creation dates, and sizes are as follows:

[0036] table_S4_mutation_rates.txt August 31, 2018 2,452KB

[0037] table_S5_gene_enrichment.txt August 31, 2018 362KB

[0038] table_S6_validation.txt August 31, 2018 362KB

[0039] The disclosed technology relates to artificial intelligence-type computers and digital data processing systems and corresponding data processing methods and products (i.e., knowledge-based systems, inference systems, and knowledge acquisition systems) for emulating intelligence, including systems for reasoning with uncertainty (e.g., fuzzy logic systems), adaptive systems, machine learning systems, and artificial neural networks. In particular, the disclosed technology relates to training deep convolutional neural networks using deep learning-based techniques. [Background technology]

[0040] The subject matter described in this section should not be assumed to be prior art merely as a result of the reference to prior art in this section. Similarly, problems associated with subject matter mentioned in this section or presented as background art should not be assumed to have been recognized in the prior art. The subject matter in this section merely represents different approaches, which may themselves correspond to implementations of the claimed technology.

[0041] Machine Learning In machine learning, input variables are used to predict output variables. The input variables are often called features and are denoted as X=(X1, X2, ..., X k ), where each X i , i∈1,…,k are features. The output variables are often called response or dependent variables, and are the variables Y i The relationship between Y and the corresponding X can be expressed by the following general formula: Y=f(X)+∈

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

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

[0044]

number

[0045] In the above formula,

number

number

[0046] function

number

number

[0047] Neural Networks The Single Layer Perceptron (SLP) is the simplest model of a neural network. It has one input layer and one activation function, as shown in Figure 1. The inputs are passed through a weighted graph. The function f takes the sum of the inputs as its argument and compares it with a threshold θ.

[0048] Figure 2 shows one implementation of a fully connected neural network containing multiple layers. A neural network is a system of interconnected artificial neurons (e.g., a1, a2, a3) that exchange messages with each other. The illustrated neural network 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(·). Each connecting line has a numerical weight (e.g., w ) that is adjusted 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 will respond correctly when fed an image to be recognized. The input layer processes the raw input, and the hidden layer processes the output from the input layer based on the weights of the connections between the input and hidden layers. The output layer takes the output from the hidden layer and processes it based on the weights of the connections between the hidden and output layers. The network contains multiple layers of feature detection neurons, each with a number of neurons that respond to different combinations of inputs from the previous layer. These layers are structured as follows: the first layer detects a set of primitive patterns in the input image data, the second layer detects patterns of patterns, and the third layer detects patterns of those patterns.

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

[0050] [Patent Document 1] International Publication No. 07 / 010252 [Patent Document 2] International Application No. 2007 / 003798 [Patent Document 3] US Patent Application Publication No. 2009 / 0088327 [Patent Document 4] US Patent Application Publication No. 2016 / 0085910 [Patent Document 5] US Patent Application Publication No. 2013 / 0296175 [Patent Document 6] International Publication No. 04 / 018497 [Patent Document 7] U.S. Patent No. 7,057,026 [Patent Document 8] International Publication No. 91 / 06678 [Patent Document 9] International Publication No. 07 / 123744 [Patent Document 10] U.S. Patent No. 7,329,492 [Patent Document 11] U.S. Patent No. 7,211,414 [Patent Document 12] U.S. Patent No. 7,315,019 [Patent Document 13] U.S. Patent No. 7,405,281 [Patent Document 14] US Patent Application Publication No. 2008 / 0108082 [Patent Document 15] U.S. Patent No. 5,641,658 [Patent Document 16] US Patent Application Publication No. 2002 / 0055100 [Patent Document 17] U.S. Patent No. 7,115,400 [Patent Document 18] US Patent Application Publication No. 2004 / 0096853 [Patent Document 19] US Patent Application Publication No. 2004 / 0002090 [Patent Document 20] US Patent Application Publication No. 2007 / 0128624 [Patent Document 21] US Patent Application Publication No. 2008 / 0009420 [Patent Document 22] US Patent Application Publication No. 2007 / 0099208 [Patent Document 23] International Publication No. 04 / 018497 [Patent Document 24] US Patent Application Publication No. 2007 / 0166705 [Patent Document 25] U.S. Patent No. 7,057,026 [Patent Document 26] US Patent Application Publication No. 2008 / 0280773 [Patent Document 27] U.S. Patent Application No. 13 / 018255 [Patent Document 28] International Publication No. 00 / 4018497 [Patent Document 29] US Patent Application Publication No. 2007 / 0166705 [Patent Document 30] U.S. Patent No. 7,057,026 [Patent Document 31] International Patent Application No. 2013 / 030867 [Patent Document 32] International Publication No. 2014 / 142831 [Non-patent literature]

[0051] [Non-Patent Document 1] T. Ching et al., Opportunities And Obstacles For Deep Learning In Biology And Medicine, www.biorxiv.org:142760, 2017 [Non-patent document 2] Angermueller C, Parnamaa T, Parts L, Stegle O. Deep Learning For Computational Biology. Mol Syst Biol. 2016;12:878 [Non-patent document 3] Park Y, Kellis M. 2015 Deep Learning For Regulatory Genomics. Nat. Biotechnol. 33, pp. 825-826. (doi:10.1038 / nbt.3313) [Non-patent document 4] Min, S., Lee, B., and Yoon, S. Deep Learning in Bioinformatics. Brief. Bioinform. bbw068 (2016) [Non-patent document 5] Leung MK, Delong A, Alipanahi B, et al. Machine Learning In Genomic Medicine: A Review of Computational Problems and Data Sets 2016 [Non-patent document 6] Libbrecht MW, Noble WS. Machine Learning Applications In Genetics and Genomics. Nature Reviews Genetics 2015;16(6):321-32 [Non-Patent Document 7] Bentley et al., Nature 456:53-59 (2008) [Non-patent document 8] Lizardi et al., Nat. Genet. 19:225-232 (1998) [Non-Patent Document 9] Dunn, Tamsen & Berry, Gwenn & Emig-Agius, Dorothea & Jiang, Yu & Iyer, Anita & Udar, Nitin & Stromberg, Michael. (2017). Pisces: An Accurate and Versatile Single Sample Somatic and Germline Variant Caller. 595-595. 10.1145 / 3107411.3108203 [Non-Patent Document 10] T Saunders, Christopher & Wong, Wendy & Swamy, Sajani & Becq, Jennifer & J Murray, Lisa & Cheetham, Keira. (2012). Strelka: Accurate somatic small-variant calling from sequenced tumor-normal sample pairs. Bioinformatics (Oxford, England). 28. 1811-7. 10.1093 / bioinformatics / bts271 [Non-Patent Document 11] Kim, S., Scheffler, K., Halpern, AL, Bekritsky, MA, Noh, E., Kallberg, M., Chen, X., Beyter, D., Krusche, P., and Saunders, CT (2017). Strelka2: Fast and accurate variant calling for clinical sequencing applications [Non-Patent Document 12] Stromberg, Michael & Roy, Rajat & Lajugie, Julien & Jiang, Yu & Li, Haochen & Margulies, Elliott. (2017). Nirvana: Clinical Grade Variant Annotator. 596-596. 10.1145 / 3107411.3108204 [Non-Patent Document 13] Iossifov et al., Nature 2014 Summary of the Invention [Means for solving the problem]

[0052] In the drawings, like reference characters generally refer to like parts throughout the different views. Further, the drawings are not necessarily drawn to scale, emphasis instead generally being placed on illustrating the principles of the disclosed technology. In the following description, various implementations of the disclosed technology are described with reference to the following drawings: [Brief explanation of the drawings]

[0053] [Figure 1] FIG. 1 illustrates a Single Layer Perceptron (SLP). [Figure 2] FIG. 1 illustrates one implementation of a feedforward neural network that includes multiple layers. [Figure 3] FIG. 1 illustrates one implementation of the functionality of a convolutional neural network. [Figure 4] FIG. 1 is a block diagram of training a convolutional neural network according to one implementation of the disclosed technology. [Figure 5] FIG. 1 illustrates one implementation of a ReLU nonlinear layer in accordance with one implementation of the disclosed technology. [Figure 6] FIG. 1 illustrates dilated convolution. [Figure 7] FIG. 1 illustrates an implementation of a subsampling layer (average / max pooling) according to one implementation of the disclosed technology. [Figure 8] FIG. 1 illustrates an implementation of two-layer convolution of a convolution layer. [Figure 9] FIG. 10 illustrates residual connections that reinject prior information downstream via feature map addition. [Figure 10] FIG. 1 illustrates one implementation of residual blocks and skip connections. [Figure 11] FIG. 1 illustrates an implementation of stacked dilated convolution. [Figure 12] FIG. 1 illustrates a batch normalization forward pass. [Figure 13] FIG. 1 illustrates a batch normalization transformation during testing. [Figure 14]FIG. 1 illustrates a batch normalization backward pass. [Figure 15] FIG. 1 illustrates the use of batch normalization layers with convolutional or densely connected layers. [Figure 16] FIG. 1 illustrates one implementation of 1D convolution. [Figure 17] FIG. 1 illustrates how Global Average Pooling (GAP) works. [Figure 18] FIG. 1 illustrates one implementation of a computing environment including a training server and a production server that can be used to implement the disclosed techniques. [Figure 19] FIG. 1 illustrates one implementation of an Atrous Convolutional Neural Network (ACNN) architecture, referred to herein as "SpliceNet." [Figure 20] FIG. 1 illustrates one implementation of a residual block that can be used by ACNNs and convolutional neural networks (abbreviated CNNs). [Figure 21] FIG. 1 illustrates another implementation of the architecture of an ACNN, referred to herein as "SpliceNet80." [Figure 22] FIG. 1 illustrates yet another implementation of an ACNN architecture, referred to herein as “SpliceNet400.” [Figure 23] FIG. 1 illustrates yet another implementation of the architecture of an ACNN, referred to herein as “SpliceNet2000.” [Figure 24] FIG. 1 illustrates yet another implementation of the architecture of an ACNN, referred to herein as “SpliceNet10000.” [Figure 25] FIG. 1 illustrates an ACNN and the different types of inputs processed by the CNN. [Figure 26] FIG. 1 illustrates an ACNN and the different types of inputs processed by the CNN. [Figure 27] FIG. 1 illustrates an ACNN and the different types of inputs processed by the CNN. [Figure 28] FIG. 1 shows an ACNN that can be trained on at least 8 million non-splice sites and a CNN that can be trained on at least 1 million non-splice sites. [Figure 29] FIG. 1 illustrates a one-hot encoder. [Figure 30] FIG. 1 illustrates the training of ACNN. [Figure 31] FIG. 1 is a diagram illustrating a CNN. [Figure 32] FIG. 1 illustrates training, validation, and testing of ACNN and CNN. [Figure 33] FIG. 1 shows reference and alternative sequences. [Figure 34] FIG. 1 shows aberrant splicing detection. [Figure 35] FIG. 1 illustrates the SpliceNet10000 processing pyramid for splice site classification. [Figure 36] FIG. 1 illustrates the SpliceNet10000 processing pyramid for aberrant splice site detection. [Figure 37A] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37B] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37C] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37D] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37E] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37F] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37G] FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 37H]FIG. 1 illustrates one implementation of deep learning to predict splicing from primary sequence. [Figure 38A] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 38B] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 38C] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 38D] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 38E] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 38F] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 38G] FIG. 1 shows one implementation of validation of rare cryptic splice mutations in RNA-seq data. [Figure 39A] FIG. 1 illustrates an implementation in which potential splice variants frequently form tissue-specific alternative splicing. [Figure 39B] FIG. 1 illustrates an implementation in which potential splice variants frequently form tissue-specific alternative splicing. [Figure 39C] FIG. 1 illustrates an implementation in which potential splice variants frequently form tissue-specific alternative splicing. [Figure 40A] FIG. 1 shows one implementation in which predicted potential splice variants have a strong adverse effect in the human population. [Figure 40B] FIG. 1 shows one implementation in which predicted potential splice variants have a strong adverse effect in the human population. [Figure 40C]FIG. 1 shows one implementation in which predicted potential splice variants have a strong adverse effect in the human population. [Figure 40D] FIG. 1 shows one implementation in which predicted potential splice variants have a strong adverse effect in the human population. [Figure 40E] FIG. 1 shows one implementation in which predicted potential splice variants have a strong adverse effect in the human population. [Figure 41A] FIG. 1 shows one implementation of de novo cryptic splice mutations in patients with rare genetic diseases. [Figure 41B] FIG. 1 shows one implementation of de novo cryptic splice mutations in patients with rare genetic diseases. [Figure 41C] FIG. 1 shows one implementation of de novo cryptic splice mutations in patients with rare genetic diseases. [Figure 41D] FIG. 1 shows one implementation of de novo cryptic splice mutations in patients with rare genetic diseases. [Figure 41E] FIG. 1 shows one implementation of de novo cryptic splice mutations in patients with rare genetic diseases. [Figure 41F] FIG. 1 shows one implementation of de novo cryptic splice mutations in patients with rare genetic diseases. [Figure 42A] FIG. 1 shows evaluation of various splicing prediction algorithms for lincRNAs. [Figure 42B] FIG. 1 shows evaluation of various splicing prediction algorithms for lincRNAs. [Figure 43A] Figure 1 shows the position-dependent effects of the TACTAAC branchpoint and the GAAGAA intra-exonic splice enhancer motif. [Figure 43B] Figure 1 shows the position-dependent effects of the TACTAAC branchpoint and the GAAGAA intra-exonic splice enhancer motif. [Figure 44A] FIG. 1 shows the influence of nucleosome positioning on splicing. [Figure 44B]FIG. 1 shows the influence of nucleosome positioning on splicing. [Figure 45] FIG. 1 shows an example of calculating effect sizes for splice break variants with compound effects. [Figure 46A] Figure 1 shows the evaluation of the SpliceNet-10k model on singletons and common variants. [Figure 46B] Figure 1 shows the evaluation of the SpliceNet-10k model on singletons and common variants. [Figure 46C] Figure 1 shows the evaluation of the SpliceNet-10k model on singletons and common variants. [Figure 47A] FIG. 1 shows validation rates and effect sizes of splice site-forming variants split by variant position. [Figure 47B] FIG. 1 shows validation rates and effect sizes of splice site-forming variants split by variant position. [Figure 48A] FIG. 10 shows the evaluation of the SpliceNet-10k model on training and test chromosomes. [Figure 48B] FIG. 10 shows the evaluation of the SpliceNet-10k model on training and test chromosomes. [Figure 48C] FIG. 10 shows the evaluation of the SpliceNet-10k model on training and test chromosomes. [Figure 48D] FIG. 10 shows the evaluation of the SpliceNet-10k model on training and test chromosomes. [Figure 49A] FIG. 1 shows de novo cryptic splice variants in patients with rare genetic diseases from synonymous region sites, intron region sites, or untranslated region sites only. [Figure 49B] FIG. 1 shows de novo cryptic splice variants in patients with rare genetic diseases from synonymous region sites, intron region sites, or untranslated region sites only. [Figure 49C]FIG. 1 shows de novo cryptic splice variants in patients with rare genetic diseases from synonymous region sites, intron region sites, or untranslated region sites only. [Figure 50A] FIG. 1 shows cryptic splice de novo mutations in ASD as a proportion of pathogenic DNMs. [Figure 50B] FIG. 1 shows cryptic splice de novo mutations in ASD as a proportion of pathogenic DNMs. [Figure 51A] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51B] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51C] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51D] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51E] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51F] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51G] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51H] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51I] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 51J] FIG. 1 shows RNA-seq validation of predicted potential splice de novo mutations in ASD patients. [Figure 52A]FIG. 1 shows the validation rate and sensitivity to RNA-seq of models trained only on canonical transcripts. [Figure 52B] FIG. 1 shows the validation rate and sensitivity to RNA-seq of models trained only on canonical transcripts. [Figure 53A] FIG. 10 shows that ensemble modeling improves SpliceNet-10k performance. [Figure 53B] FIG. 10 shows that ensemble modeling improves SpliceNet-10k performance. [Figure 53C] FIG. 10 shows that ensemble modeling improves SpliceNet-10k performance. [Figure 54A] FIG. 1 shows the evaluation of SpliceNet-10k in regions of varying exon density. [Figure 54B] FIG. 1 shows the evaluation of SpliceNet-10k in regions of varying exon density. [Figure 55] Table S1 shows one implementation of the GTEx samples used to demonstrate effect size calculations and tissue-specific splicing. [Figure 56] Table S2 shows one implementation of the cutoffs used to assess the validation rate and sensitivity of each different algorithm. [Figure 57] FIG. 1 shows one implementation of per-gene enrichment analysis. [Figure 58] FIG. 1 illustrates one implementation of genome-wide enrichment analysis. [Figure 59] FIG. 1 is a simplified block diagram of a computer system that can be used to implement the disclosed techniques. DETAILED DESCRIPTION OF THE INVENTION

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

[0055] Introduction Convolutional Neural Networks Convolutional neural networks are a special kind 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, i.e., in the case of images, patterns seen 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 a spatial hierarchy of patterns.

[0056] Regarding the first point, after learning a pattern in the upper right corner of a photo, a convolutional layer can recognize this pattern anywhere, for example, in the upper left corner. A densely connected network needs to re-learn this pattern if it appears in a new location. This makes convolutional neural networks data-efficient, as they require fewer training samples to learn a representation and have the ability to generalize.

[0057] Regarding the second point, the first convolutional layer can learn small local patterns such as edges, the second convolutional layer learns larger patterns composed of features from the first layer, and so on. This allows convolutional neural networks to efficiently learn increasingly complex and abstract visual concepts.

[0058] A convolutional neural network learns highly nonlinear mappings by interconnecting multiple layers of artificial neurons arranged in different layers with dependent activation functions. It contains one or more convolutional layers interspersed with one or more submapping and nonlinear layers, typically followed by one or more fully connected layers. Each element in a convolutional neural network receives input from a set of features in the previous layer. Each convolutional neural network learns simultaneously because neurons in the same feature map have identical weights. These locally shared weights reduce the network's complexity, allowing convolutional neural networks to avoid the complexity of data reconstruction during feature extraction and regression or classification processes when multidimensional input data enters the network.

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

[0060] For example, the first convolutional layer takes a feature map of size (28, 28, 1) and outputs a feature map of size (26, 26, 32), computing 32 filters on that input. Each of these 32 output channels contains a 26x26 grid of values, and the 26x26 grid is a response map of the filter on the input, showing the response of that filter pattern at different locations in the input. That's what is meant by the term feature map: every 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.

[0061] Convolution is defined by two main parameters: (1) the size of the patches extracted from the input—these are typically 1x1, 3x3, or 5x5—and (2) the depth of the output feature map—the number of filters computed by the convolution. Often, these start at depth 32, continue to depth 64, and end at depth 128 or 256.

[0062] Convolution works by sliding these windows of size 3x3 or 5x5 over the 3D input feature map, stopping at every position and extracting a 3D patch of surrounding features (shape (window_height, window_width, input_depth)). Each such 3D patch is then converted (via a tensor product with the same learned weight matrix, called the convolution kernel) into a 1D vector of shape (output_depth, ). All of these vectors are then spatially reassembled into a 3D output map of shape (height, width, output_depth). Every spatial location in the output feature map corresponds to the same location in the input feature map (e.g., the bottom right corner of the output captures information about the bottom right corner of the input). For example, in a 3x3 window, the vector output [i, j, :] comes from the 3D patch input [i-1:i+1, j-1:J+1, :]. The complete process is shown in detail in Figure 3.

[0063] A convolutional neural network (CN) contains a convolutional layer that performs convolution operations between input values ​​and a convolutional filter (a matrix of weights) that is learned through repeated gradient updates during training. If (m, n) is the filter size and W is the weight matrix, the convolutional layer convolves the input X with W 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 takes as it slides across the input is called the stride, and the filter area (m×n) is called the receptive field. The same convolutional filter is applied across various positions in the input, reducing the number of weights to be learned. This also enables position-invariant learning: if a significant pattern is present in the input, the convolutional filter will learn that pattern regardless of where it is located in the array.

[0064] Training a convolutional neural network 4 shows a block diagram of training a convolutional neural network according to one implementation of the disclosed technology. The convolutional neural network is tuned or trained so that input data results in a particular output estimate. The convolutional neural network is tuned using backpropagation based on comparing the output estimates with the ground truth until the output estimates gradually match or approach the ground truth.

[0065] Convolutional neural networks are trained by adjusting the weights between neurons based on the difference between the ground truth and the actual output, which can be mathematically written as follows, where δ = (ground truth) - (actual output):

[0066]

number

[0067] In one implementation, the training rules are W nm ←W nm +α(t m -φ m )a n It is defined as follows.

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

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

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

[0071]

number

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

[0073]

number

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

[0075]

number

[0076] This is done for all hidden layers, resulting in activations written as follows:

[0077]

number

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

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

[0080]

number

[0081] The weights of the output layer are v mk ←v mk +αδ ok φ m It will be updated as follows.

[0082] The weights of the hidden layer are learned using a learning rate α. v nm ←w nm +αδ hm a n It will be updated as follows.

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

number

number

number

number

number

number

[0084]

number

[0085] where α is the learning rate. Furthermore, the loss is calculated as the average of a set of n data pairs. This calculation terminates when the learning rate α becomes low enough at linear convergence. In one implementation, the gradient is calculated by using only selected data pairs sent to the Nesterov accelerated gradient and adaptive gradient to increase computational efficiency.

[0086] In one implementation, the convolutional neural network calculates the cost function using stochastic gradient descent (SGD), which calculates the gradients of the weights in the loss function by a single randomized set of data vs. z t This is approximated by calculating only from the v t +1=μv-α▽ w q(z t , w t ) w t +1=w t +v t +1

[0087] In the above equation, α is the learning rate, μ is the momentum, and t is the current weight state before the update. The convergence rate of SGD is approximately O(1 / t) when the learning rate α is reduced sufficiently both fast and slow. In other implementations, convolutional neural networks use various loss functions, such as Euclidean loss and softmax loss. In a further implementation, the Adam stochastic optimizer is used by convolutional neural networks.

[0088] Convolutional Layer The convolutional layer of a convolutional neural network acts as a feature extractor. It acts as an adaptive feature extractor that can learn and decompose input data into hierarchical features. In one implementation, a convolutional layer takes two images as input and produces a third image as output. In such an implementation, convolution operates on two images in two dimensions (2D), where one image is the input image and the other image, called the "kernel," is applied as a filter on the input image to produce the output image. Thus, for an input vector f of length n and a kernel g of length m, the convolution of f and g, f*g, is defined as follows:

[0089]

number

[0090] The convolution operation involves sliding a kernel over the input image. For each position of the kernel, the convolution value of the kernel and the input image is multiplied and the results are added. The sum of the products is the value of the output image at the point in the input image where the kernel is centered. The different outputs obtained from multiple kernels are called feature maps.

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

[0092] Nonlinear Layer FIG. 5 illustrates an implementation of a nonlinear layer according to one implementation of the disclosed technology. The nonlinear layer uses different nonlinear trigger functions to clearly identify likely features on each hidden layer. The nonlinear layer uses various specific functions to implement nonlinear triggering, including rectified linear unit (ReLU), hyperbolic tangent, absolute hyperbolic tangent, sigmoid, and continuous trigger (nonlinear) functions. In one implementation, the ReLU activation implements the function y=max(x, 0), keeping the input and output sizes of the layer the same. The advantage of using ReLU is that convolutional neural networks train many times faster. ReLU is a nonlinear, nonsaturating activation function that is linear with respect to the input when the input value is greater than zero and zero otherwise. Mathematically, the ReLU activation function is written as follows:

[0093]

number

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

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

[0096] In yet another implementation, the convolutional neural network uses a sigmoid unit activation function, which is a continuous saturating function described by the logistic function:

[0097]

number

[0098] where β = 1. The sigmoid unit activation function does not produce negative activations and is antisymmetric only with respect to the y-axis.

[0099] Dilated convolution Figure 6 shows a dilated convolution. Dilated convolution is sometimes called atrous convolution, which literally means "having holes." The French name comes from the algorithm a trous, which computes the fast dyadic wavelet transform. In these types of convolution layers, the inputs corresponding to each field of the filter are not neighboring points. This is illustrated in Figure 6. The distance between the inputs depends on the dilation factor.

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

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

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

[0103] In average pooling, the output of the pooling neuron is the average value of the input values ​​together with the set of input neurons, as described by the following equation:

[0104]

number

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

[0106] In Figure 7, the input has size 4x4. In 2x2 subsampling, the 4x4 image is divided into non-overlapping matrices of size 2x2. In average pooling, the average of the four values ​​is an integer output. In max pooling, the maximum of the four values ​​in the 2x2 matrix is ​​an integer output.

[0107] Convolution Example FIG. 8 shows an implementation of a two-layer convolution. In FIG. 8, an input of 2048 dimensions is convolved. In Convolution 1, the input is convolved by a convolution layer with two channels of 16 kernels of size 3×3. The resulting 16 feature maps are then normalized by a ReLU activation function in ReLU1 and then pooled in Pool1 by average pooling using a 16-channel pooling layer with a 3×3 kernel. In Convolution 2, the output of Pool1 is then convolved by another convolution layer consisting of 16 channels of 30 kernels of size 3×3. This is followed by another ReLU2 and average pooling in Pool2 with a kernel size of 2×2. The convolution layers use various strides and padding, for example, zero, one, two, and three. The resulting feature vector has 512 dimensions according to one implementation.

[0108] In other implementations, convolutional neural networks use 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, for example, one, two, or three fully connected layers and 100 to 200 neurons per layer. In another implementation, the convolutional neural network is a deep network with more layers and fewer neurons per layer, for example, five, six, or eight fully connected layers and 30 to 50 neurons per layer.

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

[0110]

number

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

[0112]

number

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

[0114]

number

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

[0116]

number

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

[0118]

number

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

[0120]

number

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

[0122]

number

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

[0124]

number

[0125] The weight and bias variance for row x, column y in the mth feature map of the previous layer, which receives input from k neurons in the hidden layer, is determined by the following equation:

[0126]

number

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

[0128]

number

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

[0130]

number

[0131] The weight and bias variations at row x and column y in the mth feature map of subsample layer S and convolutional layer C are determined by the following equations:

[0132]

number

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

[0134]

number

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

[0136]

number

[0137] The weight and bias variance in row r, column c of the mth convolution core of the kth feature map of the lth convolution layer C is as follows:

[0138]

number

[0139] Residual Connections Figure 9 shows residual connections, which reinject prior information downstream via feature map addition. Residual connections involve reinjecting previous representations downstream by appending past output tensors to later output tensors, helping to prevent information loss along the data processing flow. Residual connections address two common problems that arise in any large-scale deep learning model: vanishing gradients and representation bottlenecks. In general, adding residual connections to any model with more than 10 layers is likely to be beneficial. As mentioned above, residual connections involve making the output of a previous layer available as the input of a later layer, effectively creating a shortcut in a sequential network. Rather than being concatenated to form a later activation, the previous output is added to the later activation, assuming both activations are the same size. If the activations have different sizes, a linear transformation can be used to reshape the previous activation to the target shape.

[0140] Residual Learning and Skip Connections Figure 10 shows one implementation of residual blocks and skip connections. The main idea of ​​residual learning is that residual mappings are much easier to learn than the original mappings. Residual networks stack many residual units to mitigate the degradation of training accuracy. Residual blocks utilize special additive skip connections to address gradient vanishing in deep neural networks. At the beginning of the residual block, the data flow is split into two streams: the first stream holds the block's unmodified input, while the second stream applies weights and nonlinearities. At the end of the block, the two streams are merged using element-wise summation. The main advantage of such a configuration is that gradients flow easily through the network.

[0141] Deep convolutional neural networks (CNNs) benefit from residual networks, are easily trainable, and have achieved improved accuracy in image classification and object detection. A convolutional feedforward network connects the output of the lth layer as the input to the (l+1)th layer, thereby generating the following layer transitions x l =H l (x l-1 ) The residual block generates the discriminant function x l =H l (x l-1 )+x l-1 The advantage of the residual block is that the gradient can flow directly from the later layer to the earlier layer via the discriminant function. However, the discriminant function and H l The outputs of are combined by addition, which may hinder the information flow in the network.

[0142] WaveNet WaveNet is a deep neural network for generating raw audio waveforms. WaveNet is distinguished from other convolutional neural networks by its ability to achieve a relatively large "field of view" at low cost. Furthermore, it can add signal conditioning locally and globally, allowing WaveNet to be used as a text to speech (TTS) engine with multiple voices, where the TTS provides local conditioning and a specific voice provides global conditioning.

[0143] The key building block of WaveNet is the causal dilated convolution. As an extension to causal dilated convolutions, WaveNet also allows stacking of these convolutions, as shown in Figure 11. To obtain the same receptive field with the dilated convolutions in this figure, another dilation layer is required. This stack is an iteration of dilated convolutions, connecting the outputs of the dilated convolution layers to a single output. This allows WaveNet to obtain a large "field of view" of one output node at a relatively low computational cost. By comparison, to obtain a field of view of 512 inputs, a fully convolutional network (FCN) requires 511 layers. For a dilated convolutional network, eight layers are required. Stacked dilated convolutions require only seven layers with two stacks or six layers with four stacks. To grasp the difference in computational power required to cover the same field of view, the following table shows the number of weights required in a network with one filter per layer and a filter width of 2. Furthermore, it is assumed that the network uses 8-bit binary encoding.

[0144] [Table 1]

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

[0146] Batch normalization Batch normalization is a method for accelerating deep network training by making data normalization an integral part of the network architecture. Batch normalization can adaptively normalize data during training, even as the mean and variance change over time. Batch normalization works by internally maintaining an exponential moving average of the batchwise mean and variance of the data seen during training. The primary effect of batch normalization is that it aids gradient propagation—as well as residual connections—and thus enables deep networks. Some very deep networks can only be trained if they include multiple batch normalization layers.

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

[0148] Batch normalization provides a definition for feeding forward inputs and computing gradients with respect to parameters and themselves with respect to the inputs through a backward pass. In practice, a batch normalization layer is inserted after a convolutional or fully connected layer, but before the output is fed into the activation function. In a convolutional layer, different elements—i.e., activations—of the same feature map at different positions are normalized in the same way to follow the convolutional property. Therefore, all activations in a mini-batch are normalized at every position, not just activation-by-activation.

[0149] Internal covariate shift is the reason deep architectures are notoriously slow to train, as deep networks not only have to learn new representations at each layer, but also account for changes in the network distributions.

[0150] Covariate shift is a well-known problem in the deep learning field in general and occurs frequently in real-world problems. A common covariate shift problem is a difference in distribution between the training set and the test set, which leads to suboptimal generalization performance. This problem is usually addressed by a standardization or whitening preprocessing step. However, the whitening operation in particular is computationally expensive and therefore impractical in an online setting, especially when the covariate shift occurs across different layers.

[0151] Internal covariate shift is a phenomenon in which 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 functional relationships remain the same. To avoid the costly computation of the covariate matrix and to decorrelate and whiten the data at every layer and step, we normalize the distribution of each input feature at each layer across each mini-batch to have zero mean and standard deviation 1.

[0152] Forward Pass During the forward pass, the mini-batch mean and variance are calculated. For 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 by the learned scale and shift parameters. Batch normalization forward pass f BN is shown in FIG.

[0153] In Fig. 12, β is the batch mean,

number

[0154] Because normalization is a differentiable transformation, errors are propagated within these learned parameters, and thus, the expressive power of the network can be restored by learning a discriminative transformation. Conversely, by learning scale and shift parameters that are identical to the corresponding batch statistics, the batch normalization transformation does not affect the network if that were the optimal behavior to perform. At test time, since the input does not depend on other samples from the mini-batch, the batch mean and variance are replaced by their respective population statistics. Another approach is to maintain running averages of the batch statistics during training and use these to calculate the network output at test time. At test time, the batch normalization transformation can be represented as illustrated in Figure 13. In Figure 13, μ D and

number

[0155] Backward Pass Since normalization is a differentiable operation, the backward pass can be calculated as shown in FIG.

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

[0157] Global Average Pooling Figure 17 shows how global average pooling (GAP) works. By taking and recording the spatial average of features in the previous layer, global average pooling can be used to replace fully connected (FC) layers for classification. This reduces the training load and bypasses the overfitting problem. Global average pooling applies structure ahead of the model and is equivalent to a linear transformation with predefined weights. Global average pooling reduces the number of parameters and eliminates fully connected layers. Fully connected layers are typically the layers that rely most heavily on parameters and connections, and global average pooling constitutes a much cheaper way to achieve similar results. The key idea of ​​global average pooling is to generate an average value from each feature map in the previous layer as a recorded confidence coefficient and feed it directly into a softmax layer.

[0158] Global average pooling has three advantages: (1) there are no extra parameters in the global average pooling layer, so overfitting is avoided in the global average pooling layer; (2) the output of global average pooling is the average of the entire feature map, so global average pooling is more robust to spatial transformations; (3) the number of parameters in the fully connected layer is very large, typically accounting for 50% of all parameters in the entire network; replacing the fully connected layer with global average pooling can significantly reduce the model size, so global average pooling is very useful in model compression.

[0159] Global average pooling is significant because stronger features in the previous layer are expected to have higher average values. In some implementations, global average pooling can be used as a proxy for classification scores. The feature map under global average pooling can be interpreted as a confidence map and a force correspondence between the feature map and the category. Global average pooling is particularly effective when the features in the previous layer have sufficient abstraction for direct classification, but global average pooling alone is not sufficient when multi-level features are to be combined as a group, such as a part model, and this combination is best performed by adding a simple fully connected layer or other classifier after global average pooling.

[0160] term All literature and similar material 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 material, is expressly incorporated by reference in its entirety. In the event that one or more of the incorporated literature and similar material differs from or conflicts with this application, including but not limited to defined terms, term usage, techniques described, etc., this application controls.

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

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

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

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

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

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

[0167] The term "chromosome" refers to the gene carrier of a living cell, which is derived from a chromatin strand containing DNA and protein components (especially histones). The traditional internationally recognized system of numbering individual human genome chromosomes is used herein.

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

[0169] As used herein, the term "sample" generally refers to a biological fluid, cell, tissue, organ, or sample from an organism containing a nucleic acid or mixture of nucleic acids containing at least one nucleic acid sequence whose sequence and / or phase is to be determined. Such samples include, but are not limited to, saliva / oral fluid, amniotic fluid, blood, blood fractions, fine needle biopsy samples (e.g., surgical biopsies, fine needle biopsies, etc.), urine, ascites, pleural fluid, tissue explants, organ culture fluid, and any other tissue or cell preparation, or fractions or derivatives thereof or fractions or derivatives isolated therefrom. While samples are often obtained from human subjects (e.g., patients), samples can be obtained from any organism that possesses chromosomes, including, but not limited to, dogs, cats, horses, goats, sheep, cows, pigs, etc. Samples may be used directly as obtained from a biological source or after pretreatment to modify the sample's properties. For example, such pretreatment may include preparing plasma from blood, diluting viscous fluids, etc. Pretreatment methods may include, but are not limited to, filtration, precipitation, dilution, distillation, mixing, centrifugation, freezing, lyophilization, concentration, amplification, nucleic acid fragmentation, inactivation of interfering components, addition of reagents, lysis, and the like.

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

[0171] The term " paired-end sequencing " refers to the sequencing method that determines the sequence of both ends of target fragment.Paired-end sequencing can facilitate the detection of genome rearrangement and repeat fragment, as well as gene fusion and new transcription product.Methods for paired-end sequencing are described in PCT Publication No. WO07010252, PCT Application No. PCTGB2007 / 003798 and US Patent Application Publication No. US2009 / 0088327, each of which is incorporated herein by reference. In one example, the sequence of operations may be performed as follows: (a) generating clusters of nucleic acids, (b) linearizing the nucleic acids, (c) hybridizing a first sequencing primer and performing repeated cycles of extension, scanning, and deblocking as described above, (d) "flipping" the target nucleic acid on the flow cell surface by synthesizing a complementary copy, (e) linearizing the resynthesized strand, and (f) hybridizing a second sequencing primer and performing repeated cycles of extension, scanning, and deblocking as described above. The flipping operation may be performed by delivering reagents as described above for a single cycle of bridge amplification.

[0172] The term "reference genome" or "reference sequence" refers to any particular known genomic sequence, whether partial or complete, of an organism that can be used to reference identified sequences from a subject. For example, reference genomes used for human subjects, as well as many other organisms, can be found at the National Center for Biotechnology Information at ncbi.nlm.nih.gov. "Genome" refers to the complete genetic information of an organism or virus, represented by nucleic acid sequence. A genome includes both genes and non-coding DNA sequences. A reference sequence may be larger than the reads to which it is aligned. For example, it may 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 10 times larger. 5 times larger, or at least about 10 6 times larger, or at least about 10 7The reference genome sequence may be 100 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 implementations, the reference chromosome is a chromosome sequence from the human genome version hg19. Such sequences may be referred to as chromosome reference sequences, and the term reference genome is intended to cover such sequences. Other examples of reference sequences include genomes of other species, and even chromosomes of any species, partial chromosomal regions (e.g., strands), etc. In various implementations, the reference sequence is a consensus sequence derived from multiple individuals, or other combinations. However, in certain applications, the reference sequence may be obtained from a specific individual.

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

[0174] Next-generation sequencing methods include, for example, sequencing-by-synthesis technology (Illumina), pyrosequencing (454), ion semiconductor technology (Ion Torrent sequencing), single-molecule real-time sequencing (Pacific Biosciences), and sequencing by ligation (SOLiD sequencing). Depending on the sequencing method, the length of each read can vary from about 30 bp to over 10,000 bp. For example, Illumina's sequencing method using a SOLiD sequencer generates nucleic acid reads of about 50 bp. In another example, Ion Torrent sequencing generates nucleic acid reads of up to 400 bp, and 454 pyrosequencing generates nucleic acid reads of about 700 bp. In yet another example, single-molecule real-time sequencing can generate reads of 10,000 bp to 15,000 bp. Thus, in some implementations, the nucleic acid sequence reads have a length of 30-100 bp, 50-200 bp, or 50-400 bp.

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

[0176] The term "raw fragment" refers to sequence data for a portion of a genomic sequence of interest that at least partially overlaps a specified location or secondary location of interest within a sample read or sample fragment. Non-limiting examples of raw fragments include duplex stitched fragments, simplex stitched fragments, duplex un-stitched fragments, and simplex un-stitched fragments. The term "raw" is used to indicate that raw fragments contain sequence data that bear some relationship to the sequence data in the sample read, regardless of whether they represent supporting variants that correspond to, authenticate, or confirm potential variants in the sample read. The term "raw fragment" does not indicate that a fragment necessarily contains supporting variants that validate against a variant call in the sample read. For example, when a sample read is determined by a variant calling application to represent a first variant, the variant calling application may determine that one or more raw fragments lack a corresponding type of "support" that could otherwise be expected to occur given the variant in the sample read.

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

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

[0179] The term "variant" refers to a nucleic acid sequence that differs from a reference nucleic acid. Typical nucleic acid sequence variants include, but are not limited to, single nucleotide polymorphisms (SNPs), short insertion / deletion polymorphisms (indels), copy number variations (CNVs), microsatellite markers or tandem repeats, and structural variations. Somatic variant calling is the activity of identifying variants present at low frequency in a DNA sample. Somatic variant calling is noteworthy in the context of cancer treatment. Cancer is caused by the accumulation of mutations in DNA. DNA samples from tumors typically have heterogeneity, containing some normal cells, some cells in the early stages of cancer progression (with relatively few mutations), and some late-stage cells (with relatively many mutations). Because of this heterogeneity, somatic mutations often appear at low frequency when sequencing tumors (e.g., from FFPE samples). For example, an SNV may be found in only 10% of reads covering a given base. A variant to be classified as somatic or germline by a variant classifier is also referred to herein as a "variant under test."

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

[0181] The term "variant frequency" refers to the relative frequency of an allele (gene variant) at a particular locus within a population, expressed as a proportion or percentage. For example, the proportion or percentage may be the proportion of all chromosomes in the population that have that allele. For example, a sample variant frequency refers to the relative frequency of an allele / variant at a particular locus / position along a gene sequence of interest across a "population" corresponding to the number of reads and / or samples obtained for the gene sequence of interest from individuals. As another example, a baseline variant frequency refers to the relative frequency of an allele / variant at a particular locus / position along one or more baseline gene sequences in a "population" corresponding to the number of reads and / or samples obtained for one or more baseline gene sequences from a population of normal individuals.

[0182] "Variant allele frequency (VAF)" refers to the percentage of observed sequenced reads that match a variant divided by the total coverage at the target position. VAF is a measure of the proportion of sequenced reads that carry the variant.

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

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

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

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

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

[0188] "Read depth" (traditionally a number followed by an "x") refers to the number of sequenced reads with overlapping alignments at the target position. It is often expressed as an average or percentage above a cutoff across a set of intervals (such as exons, genes, or panels). For example, a clinical report might state that the panel average coverage was 1,105x, with 98% of the target bases covered >100x.

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

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

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

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

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

[0194] Nucleic acids can be prepared to contain a known primer sequence adjacent to an unknown target sequence. To initiate the first SBS sequencing cycle, one or more differently labeled nucleotides and a DNA polymerase can be flowed into or through a flow cell by a fluid flow subsystem. A single type of nucleotide can be added at a time, or the nucleotides used in the sequencing procedure can be specifically designed with reversible termination properties, allowing each cycle of the sequencing reaction to occur simultaneously in the presence of several types of labeled nucleotides (e.g., A, C, T, G). The nucleotides can contain detectable labeling moieties, such as fluorophores. When the four nucleotides are mixed, the polymerase can select the correct base to incorporate, and each sequence is extended by one base. Unincorporated nucleotides can be washed away by flowing a wash solution through the flow cell. One or more lasers can excite the nucleic acids, inducing fluorescence. The fluorescence emitted from the nucleic acids is based on the fluorophores of the incorporated bases, and different fluorophores can emit radiation at different wavelengths. A deblocking reagent is added to the flow cell, which removes the reversible terminator group from the extended and detected DNA strand. The deblocking reagent can then be washed away by flushing a wash solution through the flow cell. The flow cell is then ready for another cycle of sequencing, beginning with the introduction of labeled nucleotides as described above. The fluidic and detection operations can be repeated several times to complete the sequencing run. Exemplary sequencing methods are described, for example, in Bentley et al., Nature 456:53-59 (2008), International Publication No. WO04 / 018497, U.S. Patent No. 7,057,026, International Publication No. WO91 / 06678, International Publication No. WO07 / 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, which are incorporated herein by reference.

[0195] In some implementations, nucleic acids can be attached to a surface and amplified before or during sequencing.For example, amplification can be performed using bridge amplification to form nucleic acid clusters on 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, which are incorporated herein by reference. Another useful method for amplifying nucleic acids on a surface is rolling circle amplification (RCA), described, for example, in Lizardi et al., Nat. Genet. 19:225-232 (1998) and U.S. Patent Application Publication No. 2007 / 0099208A1, which are incorporated herein by reference.

[0196] One exemplary SBS protocol utilizes modified nucleotides with removable 3' blocks, as described, for example, in International Publication No. WO04 / 018497, U.S. Patent Application Publication No. 2007 / 0166705A1, and U.S. Patent No. 7,057,026, which are incorporated herein by reference. For example, repeated cycles of SBS reagents can be delivered to a flow cell to which target nucleic acids are attached, for example, as a result of a bridge amplification protocol. Nucleic acid clusters can be converted into a single stranded format using a linearization solution. The linearization solution can include, for example, a restriction endonuclease that can cleave one strand of each cluster. Other methods of cleavage can be used as an alternative to restriction or cleavage enzymes, including, inter alia, chemical degradation (e.g., cleavage of diol linkages with periodate), cleavage of abasic sites by exposure to heat or alkali, cleavage with endonucleases (e.g., "USER" as supplied by NEB, Ipswich, Mass., USA, part number M5505S), cleavage of ribonucleotides incorporated into amplification products otherwise consisting of deoxyribonucleotides, photochemical cleavage, or cleavage of peptide linkers. After the linearization operation, a sequencing primer can be delivered to a flow cell under conditions for hybridization of the sequencing primer to the target nucleic acid to be sequenced.

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

[0198] Although the sequencing operations are described above with respect to a particular SBS protocol, it will be understood that other protocols for sequencing any of a wide variety of other molecular analyses can be implemented as desired.

[0199] One or more processors of the system then receive the sequencing data for subsequent analysis. The sequencing data can be formatted in various ways, such as a .BAM file. The sequencing data can include, for example, a large number of sample reads. The sequencing data can include multiple sample reads with corresponding sample sequences of nucleotides. Although only one sample read is described, it should be understood that the sequencing data can include, for example, hundreds, thousands, hundreds of thousands, or millions of sample reads. Different sample reads can have different numbers of nucleotides. For example, a sample read can have from 10 to about 500 or more nucleotides. The sample reads can span the entire genome of the source. As an example, the sample reads are directed toward predetermined genetic loci, such as those genetic loci with suspected STRs or suspected SNPs.

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

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

[0202] Such a variant calling application can include four sequentially executed modules:

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

[0204] (2) Pisces Variant Caller: Calls small SNVs, insertions, and deletions. Pisces includes a variant folding algorithm for fusing variants separated by read boundaries, a basic filtering algorithm, and a simple Poisson-based variant confidence scoring algorithm. The output is a VCF.

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

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

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

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

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

[0210] b. CNVs / SVs (Yu). Annotations for copy number variations and structural variants may be provided. Similar to annotations for small variants, transcripts overlapping with SVs and previously reported structural variants may be annotated in online databases. Unlike small variants, not all overlapping transcripts need to be annotated, as there are too many transcripts overlapping with large SVs. Instead, all overlapping transcripts belonging to partially overlapping genes may be annotated. In particular, for these transcripts, the affected introns, exons, and consequences caused by the structural variant may be reported. An option is available that allows output of all overlapping transcripts, but basic information for these transcripts, such as gene symbols and flags indicating whether they are canonical overlaps or partially overlapping transcripts, may be reported. For each SV / CNV, it is also important to know whether these variants have been studied and their frequency in different populations. Therefore, we reported overlapping SVs in external databases such as 1000 Genomes, DGV, and ClinGen. To avoid using an arbitrary cutoff to determine which SVs are overlapping, instead, all overlapping transcripts can be used and the mutual overlap can be calculated, i.e., the overlap length can be divided by the minimum of the lengths of these two SVs.

[0211] c. Reporting Supplementary Annotations. Supplementary annotations are of two types: minor variants and structural variants (SVs). SVs are modeled as intervals, and overlapping SVs can be identified using interval arrays as described above. Minor variants are modeled as points and matched by location and (optionally) allele. As such, they are searched using a binary search-like algorithm. Because the supplementary annotation database can be quite large, a fairly small index is created that maps chromosome locations to file locations where the supplementary annotations are located. The index is a sorted array of objects (consisting of chromosome location and file location) that can be binary searched using the location. To keep the index size small, multiple locations (up to a certain maximum count) are compressed into a single object that stores only the value for the first location and the delta for subsequent locations. We use a binary search, so the execution time is O(lg n), where n is the number of items in the database.

[0212] d. VEP cache file

[0213] e. Transcript Database. The transcript cache and supplementary database (SAdb) files are serialized dumps of data objects such as transcripts and supplementary annotations. We use the Ensemble VEP cache as our data source for the cache. To create the cache, all transcripts are inserted into an interval array, and the final state of the array is stored in the cache file. Therefore, during annotation runs, we only need to load the precomputed interval array and perform lookups on it. Because the cache is loaded in memory and lookups are very fast (as explained above), finding overlapping transcripts is extremely fast in Nirvana (with a profile of less than 1% of the total run time).

[0214] f. Supplementary Database. Data sources for the SAdb are listed under Supplementary Material. The SAdb for a minor variant is created by a k-way merge of all data sources, such that each object in the database (identified by reference name and location) retains all associated supplementary annotations. Issues encountered while parsing the data source files are detailed on the Nirvana homepage. To limit memory usage, only the SA index is loaded into memory. This index allows for fast lookup of file locations for supplementary annotations. However, adding supplementary annotations has been identified as Nirvana's largest bottleneck (profiled at ~30% of total execution time), since the data must be fetched from disk.

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

[0216] Such variant annotation tools can include preprocessing. For example, Nirvana included numerous annotations from external data sources, such as ExAC, EVS, the 1000 Genomes Project, dbSNP, ClinVar, Cosmic, DGV, and ClinGen. To fully utilize these databases, we must remove irrelevant information from them. We implemented different strategies to handle conflicts present from different data sources. For example, for multiple dbSNP entries for the same position and alternative alleles, we combine all IDs into a comma-separated list of IDs, and if there are multiple entries with different CAF values ​​for the same allele, we use the first CAF value. For conflicting ExAC and EVS entries, we consider the number of sample counts, and the entry with the higher sample count is used. For the 1000 Genomes Project, we removed the allele frequencies of conflicting alleles. Another issue is information inaccuracy. Although we extracted allele frequency information exclusively from the 1000 Genomes Project, we noticed that for GRCh38, the allele frequencies reported in the info field did not exclude samples for which genotypes were unavailable, resulting in reduced frequencies for unavailable variants for all samples. To ensure the accuracy of our annotation, we calculate true allele frequencies using all of the individual-level genotypes. As we know, the same variant can have different representations based on different alignments. To ensure that we can accurately report information for already identified variants, we must preprocess variants from different resources to ensure consistent representation. For all external data sources, we trimmed alleles to remove overlapping nucleotides in both the reference and alternative alleles.In ClinVar, we directly parsed the xml file and performed a 5' alignment for all variants, which is often used in vcf files. Different databases can contain the same set of information. To avoid unnecessary duplication, we removed some redundant information. For example, we removed variants in the DGV that have the 1000 Genomes Project as their data source because we have already reported these variants in 1000 Genomes with more detailed information.

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

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

[0219] The one or more processors then perform validation operations in association with the potential variant calls. The validation operations may be based on quality scores and / or a hierarchy of stepped tests, as described below. When the validation operations authenticate or verify the potential variant call, the validation operations pass the variant call information (from the variant caller application) to the sample report generator. Alternatively, when the validation operations invalidate or deem the potential variant call unsuitable, the validation operations pass a corresponding indication (e.g., a negative indicator, no call indicator, an invalid call indicator) to the sample report generator. The validation operations may also pass a confidence score related to the confidence that the variant call or invalid call designation is correct.

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

[0221] The sample report may indicate that a variant call has been determined but has not been validated. When a variant call is determined to be invalid, the sample report may indicate additional information regarding the criteria for the decision not to validate the variant call. For example, the additional information in the report may include a description of the raw fragments and the degree to which the raw fragments support or refute the variant call (e.g., counts). Additionally or alternatively, the additional information in the report may include a quality score obtained by an implementation described herein.

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

[0223] The assignment operation may include analyzing a series of n nucleotides of the identification sequence to determine whether the series of n nucleotides of the identification sequence effectively matches one or more of the selected sequences. In certain implementations, the assignment operation may include analyzing the first n nucleotides of the sample sequence to determine whether the first n nucleotides of the sample sequence effectively matches one or more of the selected sequences. The number n may have a wide variety of values ​​that may be programmed into the protocol or entered by a user. For example, the number n may be defined as the number of nucleotides of the shortest selected sequence in the database. This number may be a predetermined number. The predetermined number of nucleotides may be, for example, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, or 30. However, other implementations may use fewer or more nucleotides. The number n may be selected by an individual, such as a user of the system. The number n can be based on one or more conditions. For example, the number n can be defined as the number of nucleotides in the shortest primer sequence in the database or a specified number, whichever is smaller. In some implementations, a minimum value for n, such as 15, can be used so that primer sequences with fewer than 15 nucleotides can be designated as exceptions.

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

[0225] If the sample read effectively matches one of the selected sequences in the database (i.e., an exact match or a near match as described above), the sample read is assigned or assigned to a genetic locus correlated with the selected sequence. This may be referred to as locus calling or tentative locus calling, in which the sample read is called for a genetic locus correlated with the selected sequence. However, as described above, the sample read may be called for multiple genetic loci. In such implementations, further analysis may be performed to call or assign the sample read to only one of the potential genetic loci. In some implementations, the sample read compared to the database of reference sequences is the first read from paired-end sequencing. When performing paired-end sequencing, a second read (representing an unprocessed fragment) correlated to the sample read is obtained. After assignment, subsequent analysis performed on the assigned read may be based on the type of genetic locus being called for the assigned read.

[0226] The sample reads are then analyzed to identify potential variant calls. Among other things, the results of the analysis identify potential variant calls, sample variant frequencies, a reference sequence, and the location within the genome sequence of interest where the variants occurred. For example, if a gene locus is known to contain SNPs, the assigned reads called for the gene locus may be analyzed to identify the SNPs in the assigned reads. If a gene locus is known to contain polymorphic repetitive DNA elements, the assigned reads may be analyzed to identify or characterize the polymorphic repetitive DNA elements within the sample reads. In some implementations, if the assigned reads effectively match the STR locus and the SNP locus, a warning or flag may be assigned to the sample read. The sample read may be designated as both an STR locus and an SNP locus. Analyzing may include aligning the assigned reads according to an alignment protocol to determine the sequence and / or length of the assigned reads. Alignment protocols may include methods described in International Patent Application No. PCT / US2013 / 030867, filed March 15, 2013 (International Publication No. WO2014 / 142831), which is incorporated herein by reference in its entirety.

[0227] The one or more processors then analyze the raw fragments to determine whether a supporting variant is present at the corresponding position within the raw fragment. Various types of raw fragments may be identified. For example, the variant caller may identify a type of raw fragment that represents a variant that validates the original variant call. For example, the type of raw fragment may represent a double-stitched fragment, a single-stitched fragment, a double-unstitched fragment, or a single-unstitched fragment. Optionally, other raw fragments may be identified instead of or in addition to the foregoing examples. In connection with identifying each type of raw fragment, the variant caller also identifies the position within the raw fragment at which the supporting variant occurred, as well as a count of the number of raw fragments that displayed the supporting variant. For example, the variant caller may output an indication indicating that 10 reads of the raw fragment have been identified as representing a double-stitched fragment with a supporting variant at a specific position X. The variant caller may also output an indication indicating that 5 reads of the raw fragment have been identified as representing a single-unstitched fragment with a supporting variant at a specific position Y. The variant caller may also output a number of raw fragments that corresponded to the reference sequence and therefore did not contain supporting variants that would otherwise provide evidence to validate potential variant calls in the genome sequence of interest.

[0228] Next, the count of raw fragments containing the supporting variant, as well as the position where the supporting variant occurred, are retained. Additionally or alternatively, the count of raw fragments that did not contain the supporting variant at the position of interest (relative to the position of the potential variant call within the sample read or sample fragment) can be retained. Additionally or alternatively, the count of raw fragments that correspond to the reference sequence and do not authenticate or confirm the potential variant call can be retained. The determined information is output to the variant call validation application, including the count and type of raw fragments that support the potential variant call, the position of the supporting variant within the raw fragment, the count of raw fragments that do not support the potential variant call, and the like.

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

[0230] Deep Learning in Genomics Genetic variation can help explain many diseases. Every human has a unique genetic code, and within groups of individuals, there are many genetic variants. Most genetic variants with adverse effects have been depleted from the genome by natural selection. It is important to identify which genetic variations are likely to be pathogenic or have adverse effects. This will help researchers focus on genetic variants with high pathogenic potential, accelerating the pace of diagnosis and treatment for many diseases.

[0231] Modeling the properties and functional effects (e.g., pathogenicity) of variants is an important yet challenging task in the field of genomics. Despite rapid advances in functional genome sequencing technologies, interpreting the functional consequences of variants remains a major problem due to the complexity of cell-type-specific transcriptional regulatory systems.

[0232] With respect to pathogenicity classifiers, deep neural networks are a type of artificial neural network that uses multiple nonlinear complex transformation layers to sequentially model high-level features. Deep neural networks provide feedback via backpropagation conveying the difference between observed and predicted outputs to adjust parameters. The development of deep neural networks is due to the availability of large training datasets, the power of parallel and distributed computing, and advanced training algorithms. Deep neural networks have driven significant advances in numerous fields, including computer vision, speech recognition, and natural language processing.

[0233] Convolutional neural networks (CNNs) and recurrent neural networks (RNNs) are building blocks of deep neural networks. Convolutional neural networks have been particularly successful in image recognition, using architectures with convolutional layers, nonlinear layers, and pooling layers. Recurrent neural networks are designed to exploit sequential information in input data with cyclic connections among their constituent units, such as perceptrons, long- and short-term memory units, and gated recurrent units. In addition, many other emerging deep neural networks have been proposed for limited configurations, such as deep spatiotemporal neural networks, multidimensional recurrent neural networks, and convolutional autoencoders.

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

[0235] Another core element in training deep neural networks is regularization, which refers to a strategy intended to avoid overfitting and thus achieve good generalization performance. For example, weight decay adds a penalty term to the objective loss function to encourage weight parameters to converge to smaller absolute values. Dropout removes hidden units from a neural network during training, which may be viewed as an ensemble of possible subnetworks. To enhance the capabilities of dropout, new activation functions have been proposed: maxout and a variant of dropout for recurrent neural networks called rnnDrop. Furthermore, batch normalization provides a new regularization method through learning scalar feature normalization and the mean and variance of each activation as parameters within a mini-batch.

[0236] Given the multidimensional and high-dimensional nature of sequencing data, deep neural networks hold great promise for bioinformatics research due to their broad applicability and high predictive capabilities. Convolutional neural networks have been adapted to solve sequence-based problems in genomics, such as motif discovery, pathogenic variant identification, and gene expression inference. Convolutional neural networks use a weight-sharing strategy that is particularly useful in studying DNA because they can capture sequence motifs and recursively identify local patterns in DNA that are short and likely to have significant biological functions. The hallmark of convolutional neural networks is the use of convolutional filters. Unlike traditional classification approaches based on meticulously designed, handcrafted features, convolutional filters perform adaptive feature learning, a process similar to mapping raw input data to information representations of knowledge. In this sense, convolutional filters act as a series of motif scanners, because a set of such filters can recognize relevant patterns in the input and update itself during the training procedure. Recurrent neural networks can capture long-range dependencies in sequential data of variable length, such as protein or DNA sequences.

[0237] Therefore, powerful computational models for predicting the pathogenicity of variants could be of enormous benefit to both basic science and translational research.

[0238] Currently, only 25–30% of patients with rare diseases receive a molecular diagnosis from examination of protein-coding sequences, suggesting that the remaining diagnostic yield may lie within the non-coding 99% of the genome. Here, we describe a novel deep learning network that accurately predicts splice junctions from any pre-mRNA transcript sequence, enabling accurate prediction of the splice alteration effects of non-coding variants. Synonymous and intronic mutations with predicted splice alteration outcomes validate at high rates on RNA-seq and have significant deleterious effects within the human population. De novo mutations with predicted splice alteration outcomes are significantly enriched in patients with autism and intellectual disability compared to healthy controls, and are validated against RNA-seq data in 21 of 28 of these patients. We estimate that 9–11% of pathogenic mutations in patients with rare genetic disorders are caused by this previously underappreciated class of disease variants.

[0239] Exome sequencing is transforming clinical diagnosis for patients and families with rare genetic diseases, significantly reducing the time and cost of the endless journey to seek a diagnosis when adopted as a first-line test (Monroe et al., 2016; Stark et al., 2016; Tan et al., 2017). However, the diagnostic yield of exome sequencing is approximately 25-30% in rare genetic disease cohorts, and the majority of patients remain undiagnosed even after combined exome and microarray testing (Lee et al., 2014; Trujillano et al., 2017; Yang et al., 2014). Noncoding regions play a significant role in gene regulation, accounting for 90% of causal disease loci discovered in unbiased genome-wide association studies of complex human diseases (Ernst et al., 2011; Farh et al., 2015; Maurano et al., 2012), suggesting that penetrant noncoding variants may also contribute to a significant burden of causal mutations in rare genetic diseases. Indeed, penetrant noncoding variants, often referred to as cryptic splice variants, that disrupt the normal pattern of mRNA splicing despite being outside the essential GT and AG splice dinucleotides, have long been recognized to play a significant role in rare genetic diseases (Cooper et al., 2009; Padgett, 2012; Scotti and Swanson, 2016; Wang and Cooper, 2007). However, cryptic splice mutations are often overlooked in clinical practice due to our incomplete understanding of the splicing code and the resulting difficulty in accurately identifying splice-altering variants outside of the essential GT and AG dinucleotides ( Wang and Burge, 2008 ).

[0240] Recently, RNA-seq has emerged as a promising assay for detecting splicing abnormalities in Mendelian diseases (Cummings et al., 2017; Kremer et al., 2017). However, its clinical utility remains limited to a small number of cases where the relevant cell type is known and biopsies are accessible. High-throughput screening assays for potential splice-altering variants (Soemedi et al., 2017) have expanded the characterization of splicing variations, but the genomic space in which splice-altering mutations can occur is extremely large, making it less practical for assessing random de novo mutations in genetic diseases. General prediction of splicing from any pre-mRNA sequence would potentially enable accurate prediction of the splice-altering consequences of non-coding variants, substantially improving the diagnosis of patients with genetic diseases. To date, a general predictive model of splicing from unprocessed sequences that approaches spliceosome specificity remains elusive, despite progress in specific applications such as modeling sequence properties of core splicing motifs ( Yeo and Burge, 2004 ), characterizing exon splice enhancers and silencers ( Fairbrother et al., 2002 ; Wang et al., 2004 ), and predicting cassette exon inclusion ( Barash et al., 2010 ; Jha et al., 2017 ; Xiong et al., 2015 ).

[0241] Splicing long pre-mRNAs into mature transcripts is remarkable for its accuracy and the clinical malignancy of splice-altering mutations; moreover, our understanding of the cellular mechanisms that determine its specificity remains incomplete. Here, we train a deep learning network that approaches the precision of the spliceosome in silico, identifying exon-intron boundaries from pre-mRNA sequences with 95% accuracy and predicting functional cryptic splice mutations with a validation rate of over 80% on RNA-seq. Non-coding variants predicted to alter splicing have a strong negative impact in human populations, with 80% of newly generated cryptic splice mutations undergoing negative selection, similar to the effects of other classes of protein-truncating variants. De novo cryptic splice mutations in patients with autism and intellectual disability target the same genes that are recurrently mutated by protein-truncating mutations, enabling the discovery of additional candidate disease genes. We estimate that up to 24% of penetrant causal mutations in patients with rare genetic diseases may be due to this previously underappreciated class of disease variants, highlighting the need to improve interpretation of the 99% of the genome that is non-coding for clinical sequencing applications.

[0242] Clinical exome sequencing has revolutionized diagnostics for patients and families with rare genetic diseases, significantly reducing the time and cost of the endless journey to diagnosis when adopted as a first-line test. However, diagnostic yields for exome sequencing have been reported at 25–30% in several large cohorts of rare disease patients and their parents, with the majority of patients remaining undiagnosed even after combined exome and microarray testing. The non-coding genome is highly active in genetic regulation, with non-coding variants accounting for ~90% of GWAS hits for common diseases, suggesting that rare variants in the non-coding genome may also account for a significant proportion of causal mutations in rare genetic disorders and penetrant diseases such as oncology. However, the difficulty of interpreting variants within the non-coding genome means that, outside of large structural variants, the non-coding genome currently offers little additional diagnostic benefit for rare penetrant variants, which have the greatest impact on clinical management.

[0243] The role of splice-altering mutations outside the canonical GT and AG splice dinucleotides has long been evaluated in rare diseases. Indeed, these cryptic splice variants are the most common mutations for several rare genetic disorders, such as glycogen storage disease XI (Pompe disease) and erythropoietic protoporphyria. Extended splice motifs at the 5' and 3' ends of introns are highly variable, and equally favorable motifs occur frequently in genomes, making accurate prediction of which non-coding variants may cause cryptic splicing impractical using existing methods.

[0244] To better understand how the spliceosome achieves its specificity, we trained a deep learning neural network to predict, for each nucleotide in the pre-mRNA, whether it is a splice acceptor, a splice donor, or neither, using only the transcript sequence as its input (Figure 37A). Using canonical transcripts on even chromosomes as the training set and transcripts on odd chromosomes (excluding paralogs) as the test set, the deep learning network called exon-intron boundaries with 95% accuracy, often perfectly reconstructing them to nucleotide accuracy, even for transcripts over 100 kb long, such as CFTR (Figure 37B).

[0245] We next sought to understand the specificity determinants used by the network to recognize exon-intron boundaries with such remarkable accuracy. In contrast to previous classifiers that operate based on statistically or ergonomically designed features, deep learning learns features directly from sequences in a hierarchical manner, allowing additional specificity to be imparted from long-range sequence context. Indeed, we found that the accuracy of the network was highly dependent on the length of the sequence context adjacent to the nucleotide under prediction provided as input into the network (Table 1), and when we trained a deep learning model using only 40-nt sequences, performance only modestly exceeded that of existing statistical methods. This indicates that deep learning adds little beyond existing statistical methods for recognizing individual 9- to 23-nt splicing motifs, but that broader sequence context is key to distinguishing functional splice sites from nonfunctional sites with equally strong motifs. Asking the network to predict where on the exon the sequence would be perturbed showed that truncating the donor motif also caused the acceptor signal to disappear (Figure 37C), a finding frequently observed in exon skipping events in vivo, indicating that remarkable specificity is conferred simply by requiring pairing between strong acceptor and donor motifs at an acceptable distance.

[0246] Although a large body of evidence indicates that experimental perturbations of exon length have a strong effect on exon inclusion versus exon skipping, this does not explain why the accuracy of deep learning networks continues to increase beyond 1,000 nt of configuration. To better distinguish between local splice motif-driven specificity and long-range specificity determinants, we trained local networks that take only 100 nt of configuration as input. Using the local networks to score known junctions, we find that both exons and introns have optimal lengths (~115 nt for exons and ~1,000 nt for introns) at which motif strength is minimized (Figure 37D). This relationship is absent in the 10,000-nt deep learning network (Figure 37E), indicating that intron and exon length variation is already fully incorporated within the wide-configuration deep learning network. Notably, intron and exon boundaries were never fed into the deep learning model in a broad configuration, which indicated that these distances could be derived by inferring the positions of exons and introns from sequence alone.

[0247] Systematic exploration of hexamer space also reveals that deep learning networks exploit motifs in exon-intron definition, particularly the branchpoint motif TACTAAC at positions −34 to −14, the well-characterized exon splice enhancer GAAGAA near the end of the exon, and poly-U motifs that are typically part of polypyrimidine tracts but also appear to act as exon splice silencers ( Figures 21 , 22 , 23 , and 24 ).

[0248] We extend our deep learning network to assess genetic variants for splice-altering function by predicting exon-intron boundaries in both reference and variant-containing transcript sequences and then searching for altered exon-intron boundaries. The recent availability of aggregated exome data from 60,706 individuals allows us to assess the impact of negative selection on predicted variants and altered splice function by examining their distribution within the allele frequency spectrum. We find that predicted cryptic splice variants are under strong negative selection (Figure 38A), evidenced by their relative depletion at high allele frequencies compared to expected counts, with the magnitude of depletion comparable to that of AG or GT splice-breaking variants and stop-gain variants. The impact of negative selection is even greater when considering cryptic splice variants that cause frameshifts over variants that cause in-frame changes (Figure 38B). Based on the depletion of frameshift cryptic splice variants compared with other classes of protein-truncating variations, we confidently estimate that 88% of predicted cryptic splice mutations are functional.

[0249] Although whole-genome data as aggregated as exome data are not available, restricting our ability to detect the effects of natural selection in deep intronic regions also allowed us to calculate the observed versus expected counts of cryptic splice mutations distant from exon regions. Overall, we observe a 60% depletion of cryptic splice mutations at distances >50 nt from exon-intron boundaries (Figure 38C). The signal attenuation is likely due to a combination of the smaller sample size of whole-genome data compared with exomes and the greater difficulty in predicting the effects of deep intronic variants.

[0250] We can also use the observed versus expected number of cryptic splice variants to estimate the number of cryptic splice variants under selection and how this compares to other classes of protein-truncating variants. Because cryptic splice variants only partially abolish splice function, we also evaluated the number of observed versus expected cryptic splice variants at a more relaxed threshold and estimated that there are approximately three times more deleterious rare cryptic splice variants compared with rare AG or GT splice-cleavage variants in the ExAC dataset (Figure 38D). Each individual has approximately 20 rare cryptic splice mutations, roughly equal to the number of protein-truncating variants (Figure 38E), although not all of these variants completely abolish splice function.

[0251] The recent release of the GTEx data, including 148 individuals with both whole-genome sequencing and RNA-seq from multiple tissue sites, allows us to explore the effects of rare cryptic splice variants directly within the RNA-seq data. To approximate the scenario encountered in rare disease sequencing, we considered only rare variants (singleton variants in the GTEx cohort and <1% allele frequency in 1000 genomes) and paired these with splicing events that were unique to the individual carrying the variant. Although differences in gene and tissue expression and the complexity of splice aberrations make it difficult to assess the sensitivity and specificity of deep learning predictions, we found that, at a stringent specificity threshold, >90% of rare cryptic splice mutations validated on RNA-seq (Figure 39A). Many aberrant splicing events present in RNA-seq appear to be associated with variants predicted by the deep learning classifier to have modest effects, suggesting that they only partially affect splice function. At these more sensitive thresholds, approximately 75% of novel junctions are predicted to cause abnormalities in splicing function (FIG. 38B).

[0252] The success of our deep learning network in predicting cryptic splice variants that have a strong negative impact on population sequencing data and a high rate of validation on RNA-seq suggests that this method could be used to identify additional diagnoses in rare disease sequencing studies. To test this hypothesis, we investigated de novo variants in exome sequencing studies for autism and neurodevelopmental disorders and demonstrated that cryptic splice variants were significantly enriched in affected individuals versus their healthy siblings (Figure 40A). Furthermore, enrichment for cryptic splice variants was slightly lower than that for protein-truncating variants, indicating that approximately 90% of our predicted cryptic splice variants are functional. Based on these values, we can attribute approximately 20% of disease-causing protein-truncating variants to cryptic splice variants in exons and their immediate neighboring nucleotides (Figure 40B). Extrapolating this figure to whole-genome studies where entire intron sequences can be interrogated, we estimate that 24% of causal mutations in rare genetic diseases are due to cryptic splice mutations.

[0253] We estimated the probability of calling de novo cryptic splice mutations for each individual gene and compared it to chance, allowing us to estimate the enrichment of cryptic splice mutations in candidate disease genes. De novo cryptic splice mutations were strongly enriched in genes previously hit by protein-truncating variants but not by missense variants (Figure 40C), indicating that most may cause disease through haploinsufficiency rather than other mechanisms of action. By adding predicted cryptic splice mutations to the list of protein-truncating variants, we were able to identify three additional disease genes in autism and 11 additional disease genes in intellectual disability compared to using protein-truncating variants alone (Figure 40D).

[0254] To assess the feasibility of validating cryptic splice mutations in patients for whom likely disease tissue was unavailable (in this case, brain), we performed deep RNA-seq on 37 individuals with predicted de novo cryptic splice mutations from Simon's Simplex Collection, looking for aberrant splicing events present in those individuals but absent in all other individuals in the experiment and 149 individuals from the GTEx cohort. We found that NN of the 37 patients exhibited unique aberrant splicing on RNA-seq (Figure 40E) that was explained by the predicted cryptic splice variants.

[0255] In summary, we demonstrate a deep learning model that accurately predicts cryptic splice variants with sufficient accuracy to help identify causal disease mutations in rare genetic diseases. We estimate that a substantial portion of rare disease diagnoses caused by cryptic splicing are currently missed by considering only protein-coding regions, and we emphasize the need to develop methods for interpreting the effects of penetrant rare variants in non-coding genomes.

[0256] result Accurate prediction of splicing from primary sequence using deep learning We constructed a deep residual neural network (He et al., 2016a) that predicts whether each position in a pre-mRNA transcript is a splice donor, a splice acceptor, or neither (Figure 37A, Figure 21, Figure 22, Figure 23, and Figure 24) using only the genomic sequence of the pre-mRNA transcript as input. Because splice donors and splice acceptors can be separated by tens of thousands of nucleotides, we employed a novel network architecture consisting of 32 dilated convolutional layers (Yu and Koltun, 2016) that can recognize sequence determinants spanning very large genomic distances. In contrast to previous methods that consider only short nucleotide windows flanking exon-intron boundaries ( Yeo and Burge, 2004 ) or rely on ergonomically designed features ( Xing et al., 2015 ) or experimental data such as expression or splice factor binding ( Jha et al., 2017 ), our neural network learns splicing determinants directly from the primary sequence by evaluating 10,000 nucleotides of the flanking constituent sequences to predict the splice function of each position in the pre-mRNA transcript.

[0257] We used GENCODE-annotated pre-mRNA transcript sequences (Harrow et al., 2012) on a subset of human chromosomes to train the neural network parameters and tested the network's predictions using transcripts on the remaining chromosomes, excluding paralogs. For pre-mRNA transcripts in the test dataset, the network predicted splice junctions with a Top-k accuracy of 95%, the percentage of correctly predicted splice sites at a threshold, where the number of predicted sites is equal to the actual number of splice sites present in the test dataset (Boyd et al., 2012; Yeo and Burge, 2004). Even-length genes over 100 kb, such as CFTR, are often perfectly reconstructed with nucleotide accuracy (Figure 37B). To confirm that the network is not simply dependent on exon sequence bias, we tested the network on long noncoding RNAs. Despite incomplete non-coding transcript annotation, which would be expected to reduce our accuracy, the network predicted known splice junctions in lincRNAs with 84% top-k accuracy (Figures 42A and 42B), demonstrating its ability to approximate the behavior of the spliceosome on arbitrary sequences without protein-coding selective pressure.

[0258] For each GENCODE-annotated exon in the test dataset (excluding the first and last exons of each gene), we also investigated whether the network's prediction score correlated with the proportion of reads supporting exon inclusion versus exon skipping based on RNA-seq data from the Gene and Tissue Expression atlas (GTEx) (The GTEx Consortium et al., 2015) (Figure 37C). Exons that were constitutively spliced ​​in or out across GTEx tissues had prediction scores close to 1 or 0, respectively, whereas exons that underwent a substantial degree of alternative splicing (between 10% and 90% exon inclusion averaged over the samples) tended toward intermediate scores (Pearson correlation = 0.78, P = 0).

[0259] We next sought to understand the sequence determinants utilized by the network to achieve its remarkable accuracy. We performed systematic in silico substitutions of each nucleotide near annotated exons and measured the effect on the network's prediction score at adjacent splice sites (Figure 37E). We found that truncating the sequence of a splice donor motif frequently caused the network to predict that the upstream splice acceptor site was also lost, as observed in exon skipping events in vivo, indicating that significant specificity is conferred by exon definition between a pair of upstream acceptor motifs and a set of downstream donor motifs at optimal distances (Berget, 1995). Additional motifs contributing to splicing signals include well-characterized binding motifs and branch points of the SR protein family (Figures 43A and 43B) (Fairbrother et al., 2002; Reed and Maniatis, 1988). The effect of these motifs depends heavily on their position within the exon, suggesting that their role involves specifying the precise positioning of intron-exon boundaries by distinguishing between competing acceptor and donor sites.

[0260] Training the network with variable input sequence configurations had a striking effect on the accuracy of splice prediction (Figure 37E), indicating that long-range sequence determinants up to 10,000 nt away from the splice site are essential for distinguishing functional splice junctions from the numerous nonfunctional sites with near-optimal motifs. To examine long- and short-range specificity determinants, we compared the scores assigned to junctions annotated by a model trained with an 80-nt sequence configuration (SpliceNet-80nt) with the full model trained with a 10,000-nt configuration (SpliceNet-10k). Networks trained on 80-nt sequence contexts assign lower scores to junctions adjacent to exons or introns of typical lengths (150 nt for exons and ~1000 nt for introns) (Figure 37F), consistent with previous observations that such sites tend to have weaker splice motifs compared to unusually long or short exon and intron splice sites (Amit et al., 2012; Gelfman et al., 2012; Li et al., 2015). In contrast, networks trained on 10,000-nt sequence contexts show a preference for average-length introns and exons despite their relatively weaker splice motifs, accounting for the long-range specificity conferred by the exon or intron length. Skipping of weaker motifs in long, uninterrupted introns is consistent with the faster RNA polymerase II elongation observed experimentally in the absence of exon posing, which may allow the spliceosome to recognize suboptimal motifs more quickly (Close et al., 2012; Jonkers et al., 2014; Veloso et al., 2014). Our findings suggest that the average splice junction possesses favorable long-range sequence determinants that confer substantial specificity, accounting for the large degree of sequence degeneracy tolerated at most splice motifs.

[0261] Because splicing occurs co-transcriptionally (Cramer et al., 1997; Tilgner et al., 2012), the interplay between chromatin state and co-transcriptional splicing may also guide exon definition (Luco et al., 2011) and may potentially be exploited by the network to the extent that chromatin state is predictable from primary sequence. In particular, genome-wide studies of nucleosome positioning have shown that nucleosome occupancy is higher in exons (Andersson et al., 2009; Schwartz et al., 2009; Spies et al., 2009; Tilgner et al., 2009). To test whether the network uses sequence determinants of nucleosome positioning in splice site prediction, we examined optimal acceptor motif and donor chief pairs separated by 150 nt (approximately the size of an average exon) across the genome and asked the network to predict whether the motif pair would result in exon inclusion at that locus (Figure 37G). We found that positions predicted to be favorable for exon inclusion correlated with positions with high nucleosome occupancy (Spearman correlation = 0.36, P = 0), even in intergenic regions, and this effect persisted after controlling for GC content (Figure 44A). These results suggest that the network implicitly learns to predict nucleosome positioning from primary sequence and uses it as a specificity determinant in exon definition. Similar to average-length exons and introns, exons positioned on nucleosomes have weak local splice motifs (Figure 44B), consistent with a higher tolerance for degenerate motifs in the presence of compensatory factors (Spies et al., 2009).

[0262] Although multiple studies have reported a correlation between exon and nucleosome occupancy, the causal role for nucleosome positioning during exon definition remains uncertain. Using data from 149 individuals with both RNA-seq and whole-genome sequencing from the Genotype-Tissue Expression (GTEx) cohort (The GTEx Consortium et al., 2015), we identified novel exons that were private to a single individual and corresponded to private splice-site-generating gene mutations. These private exon generation events were significantly associated with pre-existing nucleosome positioning in K562 and GM12878 cells (P = 0.006 by permutation test, Figure 37H), even though these cell lines most likely lacked the corresponding private gene mutations. Our results show that genetic variants are more likely to trigger the generation of novel exons when the resulting novel exons overlap existing nucleosome-occupied regions, supporting a causal role for nucleosome positioning in promoting exon definition.

[0263] Validation of predicted cryptic splice mutations in RNA-seq data We extended the deep learning network to assess genetic variants for splice alteration function by predicting exon-intron boundaries on both the reference pre-mRNA transcript sequence and the alternative transcript sequence containing the variant and taking the difference between the scores (Δ-score). Importantly, the network was only trained on the reference transcript sequence and splice junction annotations, never seeing variant data during training, making prediction of variant effects a challenging test of the network's ability to accurately model sequence determinants of splicing.

[0264] We explored the effects of cryptic splice variants in RNA-seq data in the GTEx cohort (The GTEx Consortium et al., 2015), which includes 149 individuals with both whole-genome sequencing and RNA-seq from multiple tissues. To approximate scenarios encountered in rare disease sequencing, we first focused on rare private mutations (present in only one individual in the GTEx cohort). We found that private mutations predicted by neural networks to have functional consequences were strongly enriched at private de novo splice junctions and at the boundaries of skipped exons in private exon-skipping events, suggesting that the majority of these predictions are functional.

[0265] To quantify the effect of splice-site generation variants on the relative production of normal and aberrant splice isoforms, we measured the number of reads supporting a novel splice event as a percentage of the total number of reads covering the site (Figure 38C) (Cummings et al., 2017). For splice-site truncation variants, we observed that many exons have a low baseline rate of exon skipping, and the effect of the variant is to increase the proportion of exon-skipping reads. Therefore, we calculated both the decrease in the proportion of reads spliced ​​at the break junction and the increase in the proportion of reads that skipped the exon, and took the larger of the two effects (Figure 45 and STAR methodology).

[0266] Potential splice variants predicted with high confidence (Δ score ≥ 0.5) validate on RNA-seq at three-quarters the rate of intrinsic GT or AG splice cleavage (Figure 38D). Both the validation rate and effect size of potential splice variants closely track their Δ score (Figure 38D and Figure 38E), demonstrating that the model's predicted score is a good proxy for a variant's splice alteration potential. Validated variants, especially those with low scores (Δ score < 0.5), often have incomplete penetrance, resulting in alternative splicing generating a mixture of both aberrant and normal transcripts in the RNA-seq data. Our estimates of validation rate and effect size are conservative and may underestimate their true value due to both unexplained splice isoform changes and nonsense-mediated decay, which frequently introduces premature stop codons and therefore preferentially degrades aberrantly spliced ​​transcripts (Figure 38C and Figure 45). This is evidenced by the fact that the average effect size of variants truncating essential GT and AG splice dinucleotides is less than the 50% expected for fully penetrant heterozygous variants.

[0267] For cryptic splice variants that generate aberrant splice isoforms in at least 3 / 10 observed copies of an mRNA transcript, the network has a sensitivity of 71% when the variant is close to an exon and 41% when the variant is within a deep intron sequence (Δ score ≥ 0.5, Figure 38F). These findings indicate that deep intronic variants are more difficult to predict because deep intronic regions contain fewer specificity determinants that are sometimes selected to reside near exons.

[0268] To benchmark the performance of our network against existing methods, we selected three popular classifiers referenced in the rare genetic disease diagnostic literature: GeneSplicer (Pertea et al., 2001), MaxEntScan (Yeo and Burge, 2004), and NNSplice (Reese et al., 1997), and plotted the RNA-seq validation rate and sensitivity at varying thresholds (Figure 38G). Similar to the experience of others in the field (Cummings et al., 2017), we find that existing classifiers have insufficient specificity when presented with a large number of genome-wide noncoding variants that may potentially affect splicing, likely because they focus on local motifs and are largely unrelated to long-range specificity determinants.

[0269] Given the significant performance gap compared to existing methods, we performed additional controls to eliminate the possibility that our results on RNA-seq data were confounded by overfitting. First, we repeated the validation and sensitivity analysis separately for private variants and variants present in multiple individuals within the GTEx cohort (Figures 46A, 46B, and 46C). Because neither the splicing machinery nor the deep learning model have access to allele frequency information, verifying that the network has similar performance across the allele frequency spectrum is an important control. We found that at the same delta-score threshold, private and common potential splice variants did not show significant differences in validation rates in RNA-seq (P > 0.05, Fisher's exact test), indicating that the network's predictions are robust to allele frequency.

[0270] Second, to validate the model's predictions on different types of cryptic splice variants that can create novel splice junctions, we separately evaluated variants that generate novel GT or AG dinucleotides, those that affect extended acceptor or donor motifs, and variants that occur within more distant regions. We find that cryptic splice variants are roughly equally distributed among the three groups, and that at the same delta score threshold, there are no significant differences in validation rates or effect sizes between groups (χ for uniformity of P > 0.3, respectively). 2 P > 0.3 by Mann-Whitney U test (Figures 47A and 47B).

[0271] Third, we performed RNA-seq validation and sensitivity analyses separately for variants on the chromosomes used for training and for variants on the rest of the chromosomes (Figure 48A and Figure 48B). Although the network was trained only on the reference genome sequence and splice annotations and was not exposed to variant data during training, we wanted to rule out the possibility that bias in variant predictions could arise from the fact that the network was looking at the reference sequence in the chromosome on which it was training. We found that the network performed equally well on variants from the training and test chromosomes, with no significant difference in validation rate or sensitivity (P > 0.05, Fisher's exact test), indicating that the network's variant predictions could not be explained by overfitting the training sequence.

[0272] Predicting potential splice variants is a more challenging problem than predicting annotated splice junctions, as reflected by the results of our model and other splice prediction algorithms (compare Figure 37E and Figure 38G). An important reason is the difference in the underlying distribution of exon coverage between the two analyses. The vast majority of GENCODE-annotated exons have strong specificity determinants, resulting in constitutive splicing and prediction scores close to 1 (Figure 37C). In contrast, most potential splice variants are only partially penetrant (Figure 38D and Figure 38E), have low-to-moderate prediction scores, and frequently cause alternative splicing that generates a mixture of both normal and aberrant transcripts. This makes the latter problem of predicting the effects of potential splice variants inherently more challenging than identifying annotated splice sites. Additional factors such as nonsense-mediated decay, unexplained isoform variation, and limitations of RNA-seq assays further contribute to lowering the RNA-seq validation rate (Figure 38C and Figure 45).

[0273] Tissue-specific alternative splicing frequently arises from weak cryptic splice variants Alternative splicing is a major mode of gene regulation used to increase transcript diversity in different tissues and developmental stages, and its dysregulation is associated with disease processes (Blencowe et al., 2006; Irimia et al., 2014; Keren et al., 2010; Licatalosi and Darnell, 2006; Wang et al., 2008). Unexpectedly, we find that the relative utilization of novel splice junctions formed by cryptic splice mutations can vary substantially between tissues (Figure 39A). Furthermore, variants causing tissue-specific differences in splicing are reproducible in multiple individuals (Figure 39B), indicating that tissue-specific biology, rather than stochastic effects, underlies these differences. We find that 35% of cryptic splice variants with weak and intermediate predicted scores (Δ scores 0.35–0.8) exhibit significant differences in the proportion of normal and abnormal transcripts produced between tissues (χ 2 (Bonferroni correlation P<0.01 for the Δ test, Figure 39C). This contrasts with variants with high prediction scores (Δ score >0.8), which were significantly less likely to produce tissue-specific effects (P=0.015). Our findings are consistent with previous observations that alternatively spliced ​​exons tend to have intermediate prediction scores (Figure 37C) compared with constitutively spliced-in or -spliced-out exons, which have scores close to 1 or 0, respectively.

[0274] These results support a model in which tissue-specific factors, such as chromatin organization and binding of RNA-binding proteins, can sway the competition between two closely favored splice junctions (Gelfman et al., 2013; Luco et al., 2010; Shukla et al., 2011; Ule et al., 2003). Strong cryptic splice variants are likely to completely shift splicing from the normal to the aberrant isoform regardless of epigenetic organization, whereas weaker variants drive splice junction selection closer to the decision boundary, resulting in the use of alternative junctions in different tissue types and cellular contexts. This highlights the unexpected role played by cryptic splice mutations in generating novel alternative splicing diversity; natural selection then has the opportunity to preserve mutations that form useful tissue-specific alternative splicing variants.

[0275] Predicted cryptic splice variants have strong adverse effects in the human population Although predicted cryptic splice variants have a high validation rate in RNA-seq, in many cases their effects are not fully penetrant, generating a mixture of both normal and aberrant splice isoforms, increasing the probability that a proportion of these cryptic splice-altering variants will not be functionally significant. To examine the signatures of natural selection on predicted cryptic splice variants, we scored each variant present in 60,706 human exomes from the Exome Aggregation Consortium (ExAC) database (Lek et al., 2016) and identified variants predicted to alter exon-intron boundaries.

[0276] To measure the degree of negative selection acting on predicted splice-altering variants, we counted the number of predicted splice-altering variants found at common allele frequencies (≥0.1% in the human population) and compared it to the number of predicted splice-altering variants at singleton allele frequencies in ExAC (i.e., 1 in 60,706 individuals). Due to the recent exponential expansion of the human population size, singleton variants represent recently formed mutations that have been minimally filtered out by purifying selection (Tennessen et al., 2012). In contrast, common variants represent a subset of neutral mutations that have been passed through the sieve of purifying selection. Thus, depletion of predicted splice-altering variants in the common allele frequency spectrum relative to singleton variants provides an estimate of the proportion of predicted splice-altering variants that have deleterious effects and are therefore functional. To avoid confounding effects on protein-coding sequences, we restricted our analysis to synonymous and intronic variants located outside essential GT and AG dinucleotides and excluded missense mutations that are also predicted to have splice-altering effects.

[0277] At common allele frequencies, confidently predicted cryptic splice variants (Δ score ≥ 0.8) are under strong negative selection, as evidenced by their relative depletion compared to expectations (Figure 40A). At this threshold, where most variants are predicted to be close to full penetrance in the RNA-seq data (Figure 38D), predicted synonymous and intronic cryptic splice mutations are only 78% depleted at common allele frequencies, comparable to the 82% depletion of frameshift, stop-gain, and essential GT or AG splice break variants (Figure 40B). The impact of negative selection is even greater when considering frameshift cryptic splice variants over variants causing in-frame changes (Figure 40C). The depletion of cryptic splice variants with frameshift consequences was nearly identical to that of other classes of protein-truncating variations, indicating that the majority of confidently predicted cryptic splice mutations within near-intronic regions (≤50 nt from known exon-intron boundaries) are functional and have strong adverse effects in the human population.

[0278] To extend this analysis to deep intronic regions >50 nt from known exon-intron boundaries, we aggregated whole-genome sequencing data from 15,496 individuals from the Genome Aggregation Database (gnomAD) cohort (Lek et al., 2016) and calculated the observed and expected counts of cryptic splice mutations at common allele frequencies. Overall, we observed a 56% depletion (Δ score ≥ 0.8) of common cryptic splice mutations at distances >50 nt from exon-intron boundaries (Figure 40D), consistent with the significant difficulty in predicting the impact of deep intronic variants, as observed in the RNA-seq data.

[0279] We next sought to estimate the potential for cryptic splice mutations to contribute to penetrant genetic disease relative to other types of protein-coding variation by measuring the number of rare cryptic splice mutations per individual in the gnomAD cohort. Based on the predicted proportion of cryptic splice mutations under negative selection (Figure 40A), the average human carries ~5 rare functional cryptic splice mutations (allele frequency <0.1%) compared with ~11 rare protein-truncating variants (Figure 40E). Cryptic splice variants outnumber essential GT or AG splice truncation variants by approximately 2:1. We note that a significant proportion of these cryptic splice variants may not completely abolish gene function by forming in-frame alterations or by not completely shifting splicing to aberrant isoforms.

[0280] De novo cryptic splice mutations are a major cause of rare genetic disorders Large-scale sequencing studies of patients with autism spectrum disorder and severe intellectual disability have demonstrated the central role of de novo protein-coding mutations (missense, nonsense, frameshift, and essential splice dinucleotides) that truncate genes within neurodevelopmental pathways ( Fitzgerald et al., 2015 , Iossifov et al., 2014 , McRae et al., 2017 , Neale et al., 2012 , De Rubeis et al., 2014 , Sanders et al., 2012 ). To assess the clinical impact of noncoding mutations acting through altered splicing, we applied neural networks to predict the effects of de novo mutations in 4,293 individuals with intellectual disability from the Deciphering Developmental Disorders cohort (DDD) (McRae et al., 2017), 3,953 individuals with autism spectrum disorder (ASD) from the Simons Simplex Collection (De Rubeis et al., 2014; Sanders et al., 2012; Turner et al., 2016) and the Autism Sequencing Consortium, and 2,073 unaffected sibling controls from the Simons Simplex Collection. To control for differences in de novo variant ascertainment across studies, we normalized the expected number of de novo variants so that the number of synonymous mutations per individual was the same across cohorts.

[0281] De novo mutations predicted to break splices are enriched 1.51-fold (P = 0.000416) and 1.30-fold (P = 0.0203) in patients with intellectual disability and autism spectrum disorder compared to healthy controls (Δ score ≥ 0.1, Figures 41A, 43A, and 43B). Splice-breaking mutations are also significantly enriched in patients versus controls when considering only synonymous and intronic mutations, except that enrichment can only be explained by mutations with dual protein-coding and splicing effects (Figures 49A, 49B, and 49C). Based on the excess of de novo mutations in affected versus unaffected individuals, cryptic splice mutations are estimated to comprise approximately 11% of pathogenic mutations in patients with autism spectrum disorder and 9% in patients with intellectual disability (Figure 41B), after adjusting for the expected proportion of mutations in regions lacking sequencing coverage or variant confirmation in each study. Most de novo predicted cryptic splice mutations in affected individuals had a delta score <0.5 (Figure 41C, Figure 50A, and Figure 50B) and would be expected to generate a mixture of normal and abnormal transcripts based on variants with similar scores in the GTEx RNA-seq dataset.

[0282] To estimate the enrichment of cryptic splice mutations within candidate disease genes relative to chance, we calculated the probability of calling a de novo cryptic splice mutation for each individual gene by adjusting for mutation rate using trinucleotide composition (Samocha et al., 2014) (Table S4). Combining both cryptic splice mutations and protein-coding mutations in de novo gene discovery yields five additional candidate genes associated with intellectual disability and two additional genes associated with autism spectrum disorder (Figures 41D and 45) that are below the discovery threshold (FDR < 0.01) when considering only protein-coding mutations (Kosmicki et al., 2017; Sanders et al., 2015).

[0283] Experimental validation of de novo cryptic splice mutations in patients with autism We obtained peripheral blood-derived lymphoblastoid cell lines (LCLs) from 36 individuals from the Simons Simplex Collection who harbored predicted de novo cryptic splice mutations in genes with at least minimal levels of LCL expression (De Rubeis et al., 2014; Sanders et al., 2012), with each individual representing only autism cases within their immediate family. As is the case with most rare genetic diseases, relevant tissues and cell types (likely the developing brain) were not accessible. Therefore, we performed high-depth mRNA sequencing (~350 million 150-bp single reads per sample, roughly 10x the coverage of GTEx) to compensate for the weak representation of many of these transcripts in LCLs. To ensure that we were validating a representative set of predicted potential splice variants, rather than simply the top predictions, we applied relatively permissive thresholds (Δ score > 0.1 for splice loss variants and Δ score > 0.5 for splice gain variants, STAR method) and performed experimental validation on all de novo variants that met these criteria.

[0284] After excluding eight individuals with insufficient RNS-seq coverage in the gene of interest, we identified unique aberrant splicing events associated with predicted de novo cryptic splice mutations in 21 of 28 patients (Figure 41E and Figure 51A, Figure 51B, Figure 51C, Figure 51D, Figure 51E, Figure 51F, Figure 51G, Figure 51H, Figure 51I, and Figure 51J). These aberrant splicing events were absent in the other 35 individuals for whom deep LCL RNA-seq was obtained, as well as in 149 individuals from the GTEx cohort. We observed nine cases of de novo splice generation, eight cases of exon skipping, and four cases of intron retention, as well as more complex splicing abnormalities, among the 21 confirmed de novo cryptic splice mutations (Figure 41F, Figure 46A, Figure 46B, and Figure 46C). Seven cases showed no aberrant splicing in LCLs despite sufficient expression of the transcript. Although a subset of these may represent false-positive predictions, some cryptic splice mutations may result in tissue-specific alternative splicing that is not observable in LCLs under these experimental conditions.

[0285] The high validation rate (75%) of predicted cryptic splice variants in patients with autism spectrum disorder indicates that most predictions are functional, despite the limitations of the RNA-seq assay. However, the enrichment of de novo cryptic splice variants in cases compared with controls (1.5-fold in DDD and 1.3-fold in ASD; Figure 41A) is only 38% of the effect size observed for de novo protein-truncating variants (2.5-fold in DDD and 1.7-fold in ASD) (Iossifov et al., 2014; McRae et al., 2017; De Rubeis et al., 2014). This allows us to quantify that functional cryptic splice variants have roughly 50% of the clinical penetrance of classic forms of protein-truncating mutations (stop-gain, frameshift, and essential splice dinucleotide), as many only partially disrupt the production of normal transcripts. Indeed, some of the most well-characterized cryptic splice mutations in Mendelian diseases, such as c.315-48T>C in FECH (Gouya et al., 2002) and c.-32-13T>G in GAA (Boerkoel et al., 1995), are hypomorphic alleles associated with milder phenotypes or later age of onset. Clinical penetrance estimates are calculated for all de novo variants that meet a relatively permissive threshold (Δ score ≥ 0.1), and variants with stronger predictive scores would be expected to have correspondingly higher penetrance.

[0286] Based on the excess of de novo mutations in cases versus controls across the ASD and DDD cohorts, 250 cases could be explained by de novo cryptic splice mutations compared to 909 cases that could be explained by de novo protein-truncating variants (Figure 41B). This is consistent with our previous estimate of the average number of rare cryptic splice mutations (~5) compared to rare protein-truncating variants (~11) per person in the general population, after incorporating the reduced penetrance of cryptic splice mutations. The wide distribution of cryptic splice mutations across the genome suggests that the proportion of cases explained by cryptic splice mutations in neurodevelopmental disorders (9–11%, Figure 41B) is likely to generalize to other rare genetic disorders in which the primary disease mechanism is loss of functional protein. To facilitate the interpretation of splice-altering mutations, we precalculate Δ-score predictions for all possible single-nucleotide substitutions genome-wide and provide them as a resource to the scientific community. We believe this resource will advance our understanding of this previously underappreciated source of genetic variation.

[0287] Specific Implementations We describe manufacturing systems, methods, and articles of manufacture for using trained atrous convolutional neural networks to detect splice sites in genomic sequences (e.g., nucleotide sequences or amino acid sequences). One or more features of one implementation may be combined with a base implementation. Implementations that are not mutually exclusive are taught as combinable. One or more features of one implementation may be combined with other implementations. The present disclosure periodically informs users of these options. The omission from some implementations of repeating references to these options should not be considered as limiting the combinations taught in the preceding sections, and these references are incorporated by reference into each of the following implementations in turn.

[0288] In this section, the terms module and stage are used interchangeably.

[0289] A system implementation of the disclosed technology includes one or more processors coupled to a memory that is loaded with computer instructions for training a splice site detector to identify splice sites within a genomic sequence (e.g., a nucleotide sequence).

[0290] As shown in Figure 30, the system trains an Atrous Convolutional Neural Network (ACNN) on at least 50,000 training examples of donor splice sites, at least 50,000 training examples of acceptor splice sites, and at least 100,000 training examples of non-splicing sites. Each training example is a target nucleotide sequence having at least one target nucleotide flanked by at least 20 nucleotides on each side.

[0291] ACNN is a convolutional neural network that uses atrous / dilated convolution, which allows for large receptive fields with few trainable parameters. Atrous / dilated convolution is a convolution in which the kernel is applied over a region larger than its length by skipping input values ​​using a step, also called the atrous convolution rate or dilation factor. Atrous / dilated convolution adds spacing between elements of the convolution filter / kernel, so that neighboring input entries (e.g., nucleotides, amino acids) at a larger interval are considered when the convolution operation is performed. This allows long-range compositional dependencies to be incorporated into the input. Atrous convolution saves partial convolution calculations so that they can be reused when neighboring nucleotides are processed.

[0292] As shown in Figure 30, to evaluate training examples using an ACNN, the system provides a target nucleotide sequence as input to the ACNN, which is further flanked by at least 40 upstream constituent nucleotides and at least 40 downstream constituent nucleotides.

[0293] Based on this evaluation, the ACNN then generates as output a triplet score for the likelihood that each nucleotide in the target nucleotide sequence is a donor splice site, an acceptor splice site, or a non-splicing site, as shown in Figure 30.

[0294] This system implementation and other disclosed systems optionally include one or more of the following features. The system may also include features described in connection with the disclosed methods. For brevity, alternative combinations of system features are not individually listed. Features applicable to manufacturing systems, manufacturing methods, and articles of manufacture are not repeated for each set of statutory classes of base features. The reader will understand how the features specified in this section can be readily combined with base features in other statutory classes.

[0295] As shown in Figures 25, 26, and 27, the input can include a target nucleotide sequence having a target nucleotide flanked on each side by 2500 nucleotides, in such implementations, the target nucleotide sequence is further flanked by 5000 upstream constituent nucleotides and 5000 downstream constituent nucleotides.

[0296] The input can include a target nucleotide sequence having a target nucleotide flanked on each side by 100 nucleotides, hi such implementations, the target nucleotide sequence is further flanked by 200 upstream and 200 downstream constituent nucleotides.

[0297] The input can include a target nucleotide sequence having a target nucleotide flanked on each side by 500 nucleotides, hi such implementations, the target nucleotide sequence is further flanked by 1000 upstream and 1000 downstream constituent nucleotides.

[0298] As shown in Figure 28, the system can train an ACNN on 150,000 training examples of donor splice sites, 150,000 training examples of acceptor splice sites, and 800,000,000 training examples of non-splicing sites.

[0299] As shown in Figure 19, an ACNN can include groups of residual blocks arranged from lowest to highest, where each group of residual blocks is parameterized by the number of convolution filters in the residual block, the convolution window size of the residual block, and the atrous convolution rate of the residual block.

[0300] As shown in Figures 21, 22, 23, and 24, in ACNN, the atrous convolution rate increases non-exponentially from the lower residual block group to the higher residual block group.

[0301] As shown in Figures 21, 22, 23, and 24, in ACNN, the convolution window size varies between groups of residual blocks.

[0302] The ACNN can be configured to evaluate an input containing a target nucleotide sequence flanked by 40 upstream and 40 downstream nucleotides. In one such implementation, the ACNN includes a group of four residual blocks and at least one skip connection. Each residual block has 32 convolution filters, 11 convolution window sizes, and 1 atrous convolution rate. This implementation of the ACNN is referred to herein as "SpliceNet80" and is shown in FIG. 21.

[0303] The ACNN can be configured to evaluate an input containing a target nucleotide sequence flanked by 200 upstream and 200 downstream constituent nucleotides. In one such implementation, the ACNN includes at least two groups of four residual blocks and at least two skip connections. Each residual block in the first group has 32 convolution filters, 11 convolution window sizes, and 1 atrous convolution rate. Each residual block in the second group has 32 convolution filters, 11 convolution window sizes, and 4 atrous convolution rate. This implementation of the ACNN is referred to herein as "SpliceNet400" and is shown in FIG. 22.

[0304] The ACNN can be configured to evaluate an input containing a target nucleotide sequence flanked by 1,000 upstream and 1,000 downstream constituent nucleotides. In one such implementation, the ACNN includes at least three groups of four residual blocks and at least three skip connections. Each residual block in the first group has 32 convolution filters, an 11 convolution window size, and an atrous convolution rate of 1. Each residual block in the second group has 32 convolution filters, an 11 convolution window size, and an atrous convolution rate of 4. Each residual block in the third group has 32 convolution filters, an 21 convolution window size, and an atrous convolution rate of 19. This implementation of the ACNN is referred to herein as "SpliceNet2000" and is shown in FIG. 23.

[0305] The ACNN can be configured to evaluate an input containing a target nucleotide sequence flanked by 5,000 upstream constituent nucleotides and 5,000 downstream constituent nucleotides. In one such implementation, the ACNN includes at least four groups of four residual blocks and at least four skip connections. Each residual block in the first group has 32 convolution filters, an 11 convolution window size, and an atrous convolution rate of 1. Each residual block in the second group has 32 convolution filters, an 11 convolution window size, and an atrous convolution rate of 4. Each residual block in the third group has 32 convolution filters, a 21 convolution window size, and an atrous convolution rate of 19. Each residual block in the fourth group has 32 convolution filters, a 41 convolution window size, and an atrous convolution rate of 25. This implementation of the ACNN is referred to herein as "SpliceNet10000" and is shown in Figure 24.

[0306] Each triplet score for each nucleotide in the target nucleotide sequence may be exponentially normalized to sum to 1. In one such implementation, the system classifies each nucleotide in the target nucleotide sequence as a donor splice site, an acceptor splice site, or a non-splicing site based on the highest score in the respective triplet scores.

[0307] As shown in Figure 35, the input dimension of ACNN is (C u +L+C d ) × 4, and C u is the number of upstream nucleotides, and C d is the number of downstream constituent nucleotides, and L is the number of nucleotides in the target nucleotide sequence. In one implementation, the dimension of the input is (5000+5000+5000)×4.

[0308] As shown in Figure 35, the dimension of the output of the ACNN may be defined as L x 3. In one implementation, the dimension of the output is 5000 x 3.

[0309] As shown in Figure 35, each group of residual blocks can generate an intermediate output by processing the preceding input. The dimension of the intermediate output may be defined as (I-[{(W-1)*D}*A]) x N, where I is the dimension of the preceding input, W is the convolution window size of the residual block, D is the atrous convolution rate of the residual block, A is the number of atrous convolution layers in the group, and N is the number of convolution filters in the residual block.

[0310] As shown in Figure 32, ACNN evaluates training examples in epochs in a batch-wise manner. Training examples are randomly sampled into batches. Each batch has a predetermined batch size. ACNN repeats the evaluation of training examples over multiple epochs (e.g., 1 to 10).

[0311] The input may include a target nucleotide sequence having two adjacent target nucleotides. The two adjacent target nucleotides may be adenine (abbreviated A) and guanine (abbreviated G). The two adjacent target nucleotides may be guanine (abbreviated G) and uracil (abbreviated U).

[0312] The system comprises a one-hot encoder (shown in Figure 29) that sparsely encodes the training examples and provides the one-hot encoding as input.

[0313] The ACNN may be parameterized by the number of residual blocks, the number of skip connections, and the number of residual connections.

[0314] ACNNs can include dimension-transforming convolutional layers that reshape the spatial and feature dimensions of the preceding input.

[0315] 20, each residual block may include at least one batch normalization layer, at least one normalized linear layer (abbreviated as ReLU), at least one atrous convolutional layer, and at least one residual connection. In one such implementation, each residual block includes two batch normalization layers, two ReLU nonlinear layers, two atrous convolutional layers, and one residual connection.

[0316] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the operations of the systems described above. Yet other implementations may include methods for performing the operations of the systems described above.

[0317] Another system implementation of the disclosed technology includes a trained splice site predictor running on multiple processors operating in parallel and coupled to a memory. The system trains an atrous convolutional neural network (ACNN) running on the multiple processors on at least 50,000 training examples of donor splice sites, at least 50,000 training examples of acceptor splice sites, and at least 100,000 training examples of non-splicing sites. Each training example used in training is a nucleotide sequence containing a target nucleotide flanked on each side by at least 400 nucleotides.

[0318] The system includes an input stage of the ACNN executing on at least one of the multiple processors, providing an input sequence of at least 801 nucleotides for evaluation of a target nucleotide, each target nucleotide being flanked on each side by at least 400 nucleotides. In another implementation, the system includes an input module of the ACNN executing on at least one of the multiple processors, providing an input sequence of at least 801 nucleotides for evaluation of a target nucleotide.

[0319] The system includes an ACNN output stage that runs on at least one of the multiple processors and translates the results of the ACNN analysis into a classification score for the likelihood that each of the target nucleotides is a donor splice site, an acceptor splice site, or a non-splicing site. In another implementation, the system includes an ACNN output module that runs on at least one of the multiple processors and translates the results of the ACNN analysis into a classification score for the likelihood that each of the target nucleotides is a donor splice site, an acceptor splice site, or a non-splicing site.

[0320] Each of the features described in this specific implementation section for the first system implementation applies equally to this system implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0321] The ACNN may be trained on 150,000 training examples of donor splice sites, 150,000 training examples of acceptor splice sites, and 800,000,000 training examples of non-splicing sites. In another implementation of the system, the ACNN includes groups of residual blocks arranged in order from lowest to highest. In yet another implementation of the system, each group of residual blocks is parameterized by the number of convolution filters in the residual block, the convolution window size of the residual block, and the atrous convolution rate of the residual block.

[0322] An ACNN can include groups of residual blocks arranged from lowest to highest, where each group of residual blocks is parameterized by the number of convolutional filters in the residual block, the convolution window size of the residual block, and the atrous convolution rate of the residual block.

[0323] In ACNN, the atrous convolution rate increases non-exponentially from the lower residual block group to the higher residual block group. Also in ACNN, the convolution window size varies between groups of residual blocks.

[0324] An ACNN can be trained on one or more training servers, as shown in FIG. 18.

[0325] The trained ACNN may be deployed on one or more production servers that receive input sequences from requesting clients, as shown in Figure 18. In one such implementation, the production server processes the input sequences through the input and output stages of the ACNN, as shown in Figure 18, to generate output that is transmitted to the clients. In another implementation, the production server processes the input sequences through the input and output modules of the ACNN, as shown in Figure 18, to generate output that is transmitted to the clients.

[0326] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the operations of the systems described above. Yet other implementations may include methods for performing the operations of the systems described above.

[0327] A method implementation of the disclosed technology involves training a splice site detector to identify splice sites within a genomic sequence (e.g., a nucleotide sequence).

[0328] The method involves feeding an input sequence of at least 801 nucleotides to an Atrous Convolutional Neural Network (abbreviated ACNN) for evaluation of target nucleotides each flanked by at least 400 nucleotides on each side.

[0329] The ACNN is trained on at least 50,000 training examples of donor splice sites, at least 50,000 training examples of acceptor splice sites, and at least 100,000 training examples of non-splicing sites. Each of the training examples used in training is a nucleotide sequence containing a target nucleotide flanked on each side by at least 400 nucleotides.

[0330] The method further includes translating the results of the ACNN analysis into a classification score for the likelihood that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site.

[0331] Each of the features described in this specific implementation section for the first system implementation applies equally to this method implementation. As indicated above, all system features are not repeated here and should be considered repeated by reference.

[0332] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the methods described above. Yet another implementation may include a system comprising a memory and one or more processors operable to execute instructions stored in the memory to perform the methods described above.

[0333] We describe manufacturing systems, methods, and articles of manufacture for using trained atrous convolutional neural networks to detect aberrant splicing in genomic sequences (e.g., nucleotide sequences). One or more features of one implementation may be combined with a base implementation. Implementations that are not mutually exclusive are taught as combinable. One or more features of one implementation may be combined with other implementations. The present disclosure periodically informs users of these options. The omission from some implementations of repeated references to these options should not be considered as limiting the combinations taught in the preceding sections, and these references are incorporated by reference into each of the following implementations in turn.

[0334] A system implementation of the disclosed technology includes one or more processors coupled to a memory that is loaded with computer instructions that implement an abnormal splicing detector operating in parallel and running on multiple processors coupled to the memory.

[0335] As shown in Figure 34, the system includes a trained atrous convolutional neural network (ACNN) running on multiple processors. ACNN is a convolutional neural network that uses atrous / dilated convolution, which enables large receptive fields with few trainable parameters. Atrous / dilated convolution is a convolution in which the kernel is applied over a region larger than its length by skipping input values ​​using a step, also called the atrous convolution rate or dilation factor. Atrous / dilated convolution adds spacing between elements of the convolution filter / kernel, allowing nearby input entries (e.g., nucleotides, amino acids) at a larger interval to be considered when the convolution operation is performed. This allows long-range compositional dependencies to be incorporated into the input. Atrous convolution saves partial convolution calculations for reuse when neighboring nucleotides are processed.

[0336] As shown in Figure 34, ACNN classifies target nucleotides in an input sequence and assigns a splice site score for the likelihood that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site. The input sequence contains at least 801 nucleotides, and each target nucleotide is flanked by at least 400 nucleotides on each side.

[0337] As shown in Figure 34, the system also includes a classifier running on at least one of the multiple processors that processes the reference and variant sequences through the ACNN and generates a splice site score for the likelihood that each target nucleotide in the reference and variant sequences is a donor splice site, an acceptor splice site, or a non-splicing site. The reference and variant sequences each have at least 101 target nucleotides, with each target nucleotide flanked by at least 400 nucleotides on each side. Figure 33 shows the reference and alternative / variant sequences.

[0338] The difference in splice site scores of the target nucleotides in the reference and variant sequences then determines whether the variant that generated the variant sequence causes aberrant splicing and is therefore pathogenic, as shown in Figure 34.

[0339] This implementation and other disclosed systems optionally include one or more of the following features. The systems may also include features described in connection with the disclosed methods. For brevity, alternative combinations of system features are not individually listed. Features applicable to manufacturing systems, manufacturing methods, and articles of manufacture are not repeated for each set of statutory classes of base features. The reader will understand how the features specified in this section can be readily combined with base features in other statutory classes.

[0340] As shown in Figure 34, the difference in splice site scores can be determined position-by-position between target nucleotides in the reference and variant sequences.

[0341] As shown in Figure 34, when the global maximum difference in splice site scores for at least one target nucleotide position is higher than a predetermined threshold, ACNN classifies the variant as causing aberrant splicing and therefore pathogenic.

[0342] As shown in Figure 17, when the global maximum difference in splice site scores for at least one target nucleotide position is below a predetermined threshold, ACNN classifies the variant as not causing aberrant splicing and therefore benign.

[0343] The threshold may be determined from among a plurality of candidate thresholds, including processing a first set of reference and variant sequence pairs generated by benign common variants to generate a first set of aberrant splicing detections, processing a second set of reference and variant sequence pairs generated by pathogenic rare variants to generate a second set of aberrant splicing detections, and selecting at least one threshold for use in the classifier that maximizes the count of aberrant splicing detections in the second set and minimizes the count of aberrant splicing detections in the first set.

[0344] In one implementation, the ACNN identifies variants that cause autism spectrum disorders (ASD), hi another implementation, the ACNN identifies variants that cause developmental delay disorders (DDD).

[0345] As shown in Figure 36, the reference sequence and the variant sequence can each have at least 101 target nucleotides, and each target nucleotide can be flanked by at least 5000 nucleotides on each side.

[0346] As shown in Figure 36, the splice site scores of the target nucleotides in the reference sequence may be encoded in a first output of the ACNN, and the splice site scores of the target nucleotides in the variant sequence may be encoded in a second output of the ACNN. In one implementation, the first output is encoded as a first 101 x 3 matrix, and the second output is encoded as a second 101 x 3 matrix.

[0347] As shown in Figure 36, in one such implementation, each row in the first 101 x 3 matrix uniquely represents a splice site score for the likelihood that a target nucleotide in the reference sequence is a donor splice site, an acceptor splice site, or a non-splicing site.

[0348] Also, as shown in Figure 36, in one such implementation, each row in the second 101 x 3 matrix uniquely represents a splice site score for the likelihood that a target nucleotide in the variant sequence is a donor splice site, an acceptor splice site, or a non-splicing site.

[0349] As shown in FIG. 36, in some implementations, the splice site scores in each row of the first 101×3 matrix and the second 101×3 matrix may be exponentially normalized to sum to one.

[0350] As shown in Figure 36, the classifier can perform row-to-row comparisons of the first 101 x 3 matrix and the second 101 x 3 matrix to determine, for each row, a change in the distribution of splice site scores. For at least one instance of the row-to-row comparison, when the change in distribution is higher than a predetermined threshold, the ACNN classifies the variant as causing aberrant splicing and therefore pathogenic.

[0351] The system comprises a one-hot encoder (shown in Figure 29) that sparsely encodes the reference and variant sequences.

[0352] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this system implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0353] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the operations of the systems described above. Yet other implementations may include methods for performing the operations of the systems described above.

[0354] Method implementations of the disclosed technology include detecting genomic variants that cause aberrant splicing.

[0355] The method involves processing a reference sequence through an Atrous Convolutional Neural Network (abbreviated ACNN) that is trained to detect differential splicing patterns within a target subsequence of an input sequence by classifying each nucleotide within the target subsequence as a donor splice site, an acceptor splice site, or a non-splicing site.

[0356] The method includes detecting a first differential splicing pattern within the reference target subsequence by classifying each nucleotide within the reference target subsequence as a donor splice site, an acceptor splice site, or a non-splicing site based on the processing.

[0357] The method includes processing a variant sequence through an ACNN, where the variant sequence and the reference sequence differ by at least one variant nucleotide located within the variant target subsequence.

[0358] The method includes detecting a second differential splicing pattern within the variant target subsequence by classifying each nucleotide within the variant target subsequence as a donor splice site, an acceptor splice site, or a non-splicing site based on the processing.

[0359] The method includes determining, on a nucleotide-by-nucleotide basis, differences between the first and second differential splicing patterns by comparing the splice site classifications of the reference and variant target subsequences.

[0360] When the difference is higher than a predetermined threshold, the method includes classifying the variant as causing aberrant splicing and therefore pathogenic, and storing the classification result in memory.

[0361] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this method implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0362] The differential splicing pattern can identify the positional distribution of splicing events within the target subsequence. Examples of splicing events include at least one of cryptic splicing, exon skipping, mutually exclusive exons, alternative donor sites, alternative acceptor sites, and intron retention.

[0363] The reference target subsequence and the variant target subsequence are aligned with respect to nucleotide position and can differ by at least one variant nucleotide.

[0364] The reference target subsequence and the variant target subsequence may each have at least 40 nucleotides and each be flanked by at least 40 nucleotides on each side.

[0365] The reference target subsequence and the variant target subsequence may each have at least 101 nucleotides and each be flanked by at least 5000 nucleotides on each side.

[0366] The variant target subsequence can include two variants.

[0367] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the methods described above. Yet another implementation may include a system comprising a memory and one or more processors operable to execute instructions stored in the memory to perform the methods described above.

[0368] We describe manufacturing systems, methods, and articles of manufacture for using trained convolutional neural networks to detect splice sites and aberrant splicing within genomic sequences (e.g., nucleotide sequences). One or more features of one implementation may be combined with a base implementation. Implementations that are not mutually exclusive are taught as combinable. One or more features of one implementation may be combined with other implementations. The present disclosure periodically informs users of these options. The omission from some implementations of repeated references to these options should not be considered as limiting the combinations taught in the preceding sections, and these references are incorporated herein by reference in each of the following implementations in turn.

[0369] A system implementation of the disclosed technology includes one or more processors coupled to a memory that is loaded with computer instructions for training a splice site detector to identify splice sites within a genomic sequence (e.g., a nucleotide sequence).

[0370] The system trains a convolutional neural network (CNN) on at least 50,000 training examples of donor splice sites, at least 50,000 training examples of acceptor splice sites, and at least 100,000 training examples of non-splicing sites, where each training example is a target nucleotide sequence having at least one target nucleotide flanked by at least 20 nucleotides on each side.

[0371] To evaluate training examples using a CNN, the system provides as input to the CNN a target nucleotide sequence that is further flanked by at least 40 upstream constituent nucleotides and at least 40 downstream constituent nucleotides.

[0372] Based on this evaluation, the CNN then generates as output a triplet score for the likelihood that each nucleotide in the target nucleotide sequence is a donor splice site, an acceptor splice site, or a no-splice site.

[0373] This system implementation and other disclosed systems optionally include one or more of the following features. The system may also include features described in connection with the disclosed methods. For brevity, alternative combinations of system features are not individually listed. Features applicable to manufacturing systems, manufacturing methods, and articles of manufacture are not repeated for each set of statutory classes of base features. The reader will understand how the features specified in this section can be readily combined with base features in other statutory classes.

[0374] The input can include a target nucleotide sequence having a target nucleotide flanked on each side by 100 nucleotides, hi such implementations, the target nucleotide sequence is further flanked by 200 upstream and 200 downstream constituent nucleotides.

[0375] As shown in Figure 28, the system can train a CNN on 150,000 training examples of donor splice sites, 150,000 training examples of acceptor splice sites, and 1,000,000 training examples of non-splicing sites.

[0376] As shown in FIG. 31, a CNN can be parameterized by the number of convolutional layers, the number of convolutional filters, and the number of subsampling layers (e.g., max pooling and average pooling).

[0377] As shown in FIG. 31, a CNN can include one or more fully connected layers and a terminal classification layer.

[0378] CNNs can include dimension-transforming convolutional layers that reshape the spatial and feature dimensions of the preceding input.

[0379] Each triplet score for each nucleotide in the target nucleotide sequence may be exponentially normalized to sum to 1. In one such implementation, the system classifies each nucleotide in the target nucleotide sequence as a donor splice site, an acceptor splice site, or a non-splicing site based on the highest score in the respective triplet scores.

[0380] As shown in Figure 32, the CNN evaluates training examples in epochs in a batch-wise manner. The training examples are randomly sampled into batches. Each batch has a predetermined batch size. The CNN repeats the evaluation of the training examples over multiple epochs (e.g., 1 to 10).

[0381] The input may include a target nucleotide sequence having two adjacent target nucleotides. The two adjacent target nucleotides may be adenine (abbreviated A) and guanine (abbreviated G). The two adjacent target nucleotides may be guanine (abbreviated G) and uracil (abbreviated U).

[0382] The system comprises a one-hot encoder (shown in Figure 32) that sparsely encodes the training examples and provides the one-hot encoding as input.

[0383] A CNN can be parameterized by the number of residual blocks, the number of skip connections, and the number of residual connections.

[0384] Each residual block can include at least one batch normalization layer, at least one normalized linear layer (abbreviated as ReLU), at least one dimension transformation layer, and at least one residual connection. Each residual block can include two batch normalization layers, two ReLU nonlinear layers, two dimension transformation layers, and one residual connection.

[0385] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this system implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0386] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the operations of the systems described above. Yet other implementations may include methods for performing the operations of the systems described above.

[0387] Another system implementation of the disclosed technology includes a trained splice site predictor running on multiple processors operating in parallel and coupled to a memory. The system trains a convolutional neural network (CNN) running on the multiple processors on at least 50,000 training examples of donor splice sites, at least 50,000 training examples of acceptor splice sites, and at least 100,000 training examples of non-splicing sites. Each training example used in training is a nucleotide sequence containing a target nucleotide flanked on each side by at least 400 nucleotides.

[0388] The system includes an input stage of a CNN executing on at least one of the multiple processors, providing an input sequence of at least 801 nucleotides for evaluation of a target nucleotide, each target nucleotide being flanked on each side by at least 400 nucleotides. In another implementation, the system includes an input module of a CNN executing on at least one of the multiple processors, providing an input sequence of at least 801 nucleotides for evaluation of a target nucleotide.

[0389] The system includes a CNN output stage that runs on at least one of the multiple processors and translates the results of the CNN analysis into a classification score for the likelihood that each of the target nucleotides is a donor splice site, an acceptor splice site, or a non-splicing site. In another implementation, the system includes a CNN output module that runs on at least one of the multiple processors and translates the results of the CNN analysis into a classification score for the likelihood that each of the target nucleotides is a donor splice site, an acceptor splice site, or a non-splicing site.

[0390] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this system implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0391] The CNN can be trained on 150,000 training examples of donor splice sites, 150,000 training examples of acceptor splice sites, and 800,000,000 training examples of non-splicing sites.

[0392] A CNN can be trained on one or more training servers.

[0393] The trained CNN can be deployed on one or more production servers that receive input sequences from requesting clients. In one such implementation, the production server processes the input sequences through the input and output stages of the CNN to generate output that is transmitted to the client. In another implementation, the production server processes the input sequences through the input and output stages of the CNN to generate output that is transmitted to the client.

[0394] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the operations of the systems described above. Yet other implementations may include methods for performing the operations of the systems described above.

[0395] A method implementation of the disclosed technology includes training a splice site detector to identify splice sites within a genomic sequence (e.g., a nucleotide sequence). The method includes providing an input sequence of at least 801 nucleotides to a convolutional neural network (abbreviated CNN) for evaluation of target nucleotides, each flanked on each side by at least 400 nucleotides.

[0396] The CNN is trained on at least 50,000 training examples of donor splice sites, at least 50,000 training examples of acceptor splice sites, and at least 100,000 training examples of non-splicing sites. Each of the training examples used in training is a nucleotide sequence containing a target nucleotide flanked on each side by at least 400 nucleotides.

[0397] The method further includes translating the results of the CNN analysis into a classification score for the likelihood that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site.

[0398] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this method implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0399] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the methods described above. Yet another implementation may include a system comprising a memory and one or more processors operable to execute instructions stored in the memory to perform the methods described above.

[0400] Yet another system implementation of the disclosed technology includes one or more processors coupled to a memory that is loaded with computer instructions implementing an abnormal splicing detector operating in parallel and running on multiple processors coupled to the memory.

[0401] The system includes a trained convolutional neural network (CNN) running on multiple processors.

[0402] As shown in Figure 34, the CNN classifies target nucleotides in an input sequence and assigns a splice site score for the likelihood that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site. The input sequence contains at least 801 nucleotides, and each target nucleotide is flanked by at least 400 nucleotides on each side.

[0403] As shown in Figure 34, the system also includes a classifier running on at least one of the multiple processors that processes the reference and variant sequences through a CNN and generates a splice site score for the likelihood that each target nucleotide in the reference and variant sequences is a donor splice site, an acceptor splice site, or a non-splicing site. The reference and variant sequences each have at least 101 target nucleotides, and each target nucleotide is flanked by at least 400 nucleotides on each side.

[0404] The difference in splice site scores of the target nucleotides in the reference and variant sequences then determines whether the variant that generated the variant sequence causes aberrant splicing and is therefore pathogenic, as shown in Figure 34.

[0405] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this system implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0406] The difference in splice site scores can be determined position by position between the target nucleotides in the reference and variant sequences.

[0407] The CNN classifies a variant as causing aberrant splicing and therefore pathogenic when the maximum global difference in splice site scores for at least one target nucleotide position is higher than a predetermined threshold.

[0408] When the global maximum difference in splice site scores for at least one target nucleotide position is below a predetermined threshold, the CNN classifies the variant as not causing aberrant splicing and therefore benign.

[0409] The threshold may be determined from among a plurality of candidate thresholds, including processing a first set of reference and variant sequence pairs generated by benign common variants to generate a first set of aberrant splicing detections, processing a second set of reference and variant sequence pairs generated by pathogenic rare variants to generate a second set of aberrant splicing detections, and selecting at least one threshold for use in the classifier that maximizes the count of aberrant splicing detections in the second set and minimizes the count of aberrant splicing detections in the first set.

[0410] In one implementation, the CNN identifies variants that cause autism spectrum disorders (ASD), hi another implementation, the CNN identifies variants that cause developmental delay disorders (DDD).

[0411] The reference sequence and the variant sequence each have at least 101 target nucleotides, and each target nucleotide can be flanked on each side by at least 1000 nucleotides.

[0412] The splice site scores of the target nucleotides in the reference sequence may be encoded in a first output of the CNN, and the splice site scores of the target nucleotides in the variant sequence may be encoded in a second output of the CNN. In one implementation, the first output is encoded as a first 101 x 3 matrix, and the second output is encoded as a second 101 x 3 matrix.

[0413] In one such implementation, each row in the first 101 x 3 matrix uniquely represents a splice site score for the likelihood that a target nucleotide in the reference sequence is a donor splice site, an acceptor splice site, or a non-splicing site.

[0414] Also, in one such implementation, each row in the second 101 x 3 matrix uniquely represents a splice site score for the likelihood that a target nucleotide in the variant sequence is a donor splice site, an acceptor splice site, or a non-splicing site.

[0415] In some implementations, the splice site scores in each row of the first 101×3 matrix and the second 101×3 matrix may be exponentially normalized to sum to one.

[0416] The classifier can perform row-to-row comparisons of the first 101 x 3 matrix and the second 101 x 3 matrix and determine, for each row, a change in the distribution of splice site scores. For at least one instance of the row-to-row comparison, when the change in distribution is higher than a predetermined threshold, the CNN classifies the variant as causing aberrant splicing and therefore pathogenic.

[0417] The system comprises a one-hot encoder (shown in Figure 29) that sparsely encodes the reference and variant sequences.

[0418] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the operations of the systems described above. Yet other implementations may include methods for performing the operations of the systems described above.

[0419] Method implementations of the disclosed technology include detecting genomic variants that cause aberrant splicing.

[0420] The method involves processing a reference sequence through an Atrous Convolutional Neural Network (abbreviated CNN) that has been trained to detect differential splicing patterns within a target subsequence of an input sequence by classifying each nucleotide within the target subsequence as a donor splice site, an acceptor splice site, or a non-splicing site.

[0421] The method includes detecting a first differential splicing pattern within the reference target subsequence by classifying each nucleotide within the reference target subsequence as a donor splice site, an acceptor splice site, or a non-splicing site based on the processing.

[0422] The method includes processing a variant sequence through a CNN, wherein the variant sequence and the reference sequence differ by at least one variant nucleotide located within the variant target subsequence.

[0423] The method includes detecting a second differential splicing pattern within the variant target subsequence by classifying each nucleotide within the variant target subsequence as a donor splice site, an acceptor splice site, or a non-splicing site based on the processing.

[0424] The method includes determining, on a nucleotide-by-nucleotide basis, differences between the first and second differential splicing patterns by comparing the splice site classifications of the reference and variant target subsequences.

[0425] When the difference is higher than a predetermined threshold, the method includes classifying the variant as causing aberrant splicing and therefore pathogenic, and storing the classification result in memory.

[0426] Each of the features described in this particular implementation section relative to other system and method implementations applies equally to this method implementation, and as indicated above, all system features are not repeated here and should be considered repeated by reference.

[0427] The differential splicing pattern can identify the positional distribution of splicing events within the target subsequence. Examples of splicing events include at least one of cryptic splicing, exon skipping, mutually exclusive exons, alternative donor sites, alternative acceptor sites, and intron retention.

[0428] The reference target subsequence and the variant target subsequence are aligned with respect to nucleotide position and can differ by at least one variant nucleotide.

[0429] The reference target subsequence and the variant target subsequence may each have at least 40 nucleotides and each be flanked by at least 40 nucleotides on each side.

[0430] The reference target subsequence and the variant target subsequence may each have at least 101 nucleotides and each be flanked by at least 1000 nucleotides on each side.

[0431] The variant target subsequence can include two variants.

[0432] Other implementations may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform the methods described above. Yet another implementation may include a system comprising a memory and one or more processors operable to execute instructions stored in the memory to perform the methods described above.

[0433] The preceding description is presented to enable making and using the disclosed technology. It will be apparent that various modifications can be made to the disclosed implementations, and the general principles defined herein may be applied to other implementations and applications without departing from the spirit or scope of the disclosed technology. Thus, the disclosed technology is not intended to be limited to the implementations shown, but is to be accorded the widest scope consistent with the principles and features disclosed herein. The scope of the disclosed technology is defined by the appended claims.

[0434] Per-gene enrichment analysis Figure 57 shows one implementation of per-gene enrichment analysis. In one implementation, the aberrant splicing detector is further configured to implement per-gene enrichment analysis to determine the pathogenicity of variants determined to cause aberrant splicing. For a particular gene sampled from a cohort of individuals suffering from a genetic disease, per-gene enrichment analysis includes applying a trained ACNN to identify candidate variants in the particular gene that cause aberrant splicing, determining a baseline number of mutations for the particular gene based on summing the observed trinucleotide mutation rates of the candidate variants and multiplying the sum by the transmission count and the size of the cohort, applying the trained ACNN to identify de novo variants in the particular gene that cause aberrant splicing, and comparing the baseline number of mutations with the count of de novo variants. Based on the output of the comparison, per-gene enrichment analysis determines that the particular gene is associated with the genetic disease and that the de novo variant is pathogenic. In some implementations, the genetic disease is autism spectrum disorder (ASD). In other implementations, the genetic disease is developmental delay disorder (abbreviated DDD).

[0435] In the example shown in Figure 57, five candidate variants in a particular gene are classified by the aberrant splicing detector as causing aberrant splicing. These five candidate variants are classified as causing aberrant splicing by the aberrant splicing detector. -8 , 10 -2 , 10 -1 , 10 5 , and 10 1 The baseline number of mutations for a particular gene is 10 based on summing the observed trinucleotide mutation rates for each of the five candidate variants and multiplying that sum by the transmission / chromosome count (2) and the cohort size (100). -5 This is then compared to the de novo variant count (3).

[0436] In some implementations, the aberrant splicing detector is further configured to perform the comparison using a statistical test that produces a p-value as an output.

[0437] In other implementations, the aberrant splicing detector is further configured to compare the baseline number of mutations with the count of de novo variants and, based on the output of the comparison, determine that the particular gene is not associated with a genetic disease and that the de novo variants are benign.

[0438] In one implementation, at least some of the candidate variants are protein truncation variants.

[0439] In another implementation, at least some of the candidate variants are missense variants.

[0440] Genome-wide enrichment analysis Figure 58 shows an implementation of genome-wide enrichment analysis. In another implementation, the aberrant splicing detector is further configured to implement genome-wide enrichment analysis to determine the pathogenicity of variants determined to cause aberrant splicing. The genome-wide enrichment analysis includes applying a trained ACNN to identify a first set of de novo variants that cause aberrant splicing in multiple genes sampled from a cohort of healthy individuals, applying the trained ACNN to identify a second set of de novo variants that cause aberrant splicing in multiple genes sampled from a cohort of individuals suffering from a genetic disease, comparing the counts of the first and second sets, and determining, based on the output of the comparison, that the second set of de novo variants is enriched in the cohort of individuals suffering from the genetic disease and therefore pathogenic. In some implementations, the genetic disease is autism spectrum disorder (abbreviation ASD). In other implementations, the genetic disease is developmental delay disorder (abbreviation DDD).

[0441] In some implementations, the aberrant splicing detector is further configured to perform the comparison using a statistical test that generates a p-value as an output. In one implementation, this comparison can be further parameterized by the respective cohort sizes.

[0442] In some implementations, the aberrant splicing detector is further configured to compare the respective counts of the first and second sets and, based on the output of the comparison, determine that the second set of de novo variants is not enriched in the cohort of individuals suffering from the genetic disease and is therefore benign.

[0443] In the example shown in Figure 58, the mutation rate in the healthy cohort (0.001) and the mutation rate in the affected cohort (0.004) are illustrated along with the mutation rate per individual (4).

[0444] Essay Despite the limited diagnostic yield of exon sequencing in patients with severe genetic disorders, clinical sequencing has focused on rare coding mutations, largely ignoring noncoding genomic variation due to interpretation challenges. Here, we introduce a deep learning network that accurately predicts splicing from primary nucleotide sequences, thereby identifying noncoding mutations that disrupt the normal patterning of exons and introns, which have profound consequences on the resulting protein. We demonstrate that the predicted cryptic splice mutations are highly validated by RNA-seq, have a strong deleterious effect in the human population, and are a major cause of rare genetic diseases.

[0445] Using deep learning networks as a computer-simulation-based model of the spliceosome, we were able to reconstruct the specificity determinants that enable the spliceosome to achieve remarkable precision in vivo. We reaffirm many of the discoveries made over the past 40 years of research into the splicing machinery and show that the spliceosome integrates multiple short- and long-range specificity determinants into its decisions. In particular, we find that the perceived degeneracy of most splice motifs is explained by the presence of long-range determinants, such as exon / intron length and nucleosome positioning, that more than compensate for and obviate the need for additional motif-level specificity. Our findings demonstrate the promise of deep learning models for providing biological insights, rather than simply acting as black-box classifiers.

[0446] Deep learning is a relatively new technique in biology and is not without potential trade-offs. By automatically learning to extract features from sequences, deep learning models can take advantage of novel sequence determinants that human experts cannot adequately account for, but they also run the risk of incorporating features that do not reflect the true behavior of the spliceosome. These irrelevant features may increase the apparent accuracy of predicting annotated exon-intron boundaries but decrease the accuracy of predicting the splice-altering effects of any sequence changes induced by genetic variation. Because accurate variant prediction provides the strongest evidence that a model can generalize to true biology, we provide validation of predicted splice-altering variants using three fully orthogonal methods: RNA-seq, natural selection in human populations, and de novo variant enrichment in case-control cohorts. While this does not completely preclude the incorporation of irrelevant features into the model, the resulting model appears to be sufficiently reliable for the true biology of splicing to be of great value for practical applications such as identifying potential splice mutations in patients suffering from genetic diseases.

[0447] Compared with other classes of protein-truncating mutations, a particularly intriguing aspect of cryptic splice mutations is the widespread phenomenon of alternative splicing due to splice-altering variants with incomplete penetrance, which tend to weaken canonical splice sites relative to alternative splice sites, resulting in the production of a mixture of both aberrant and normal transcripts in RNA-seq data. The observation that these variants frequently drive tissue-specific alternative splicing highlights the unexpected role that cryptic splice mutations play in generating novel alternative splicing diversity. A promising future direction would be to train deep learning models on splice junction annotations from RNA-seq of relevant tissues, thereby obtaining tissue-specific models of alternative splicing. Training networks on annotations derived directly from RNA-seq data would also help fill gaps in GENCODE annotations and improve model performance in variant prediction (Figures 52A and 52B).

[0448] Our understanding of how mutations in the noncoding genome cause disease in humans is still far from complete. The discovery of penetrant de novo cryptic splice mutations in childhood neurodevelopmental disorders demonstrates that improved interpretation of the noncoding genome can directly benefit patients suffering from serious genetic disorders. Cryptic splice mutations also play a major role in cancer (Jung et al., 2015; Sanz et al., 2010; Supek et al., 2014), and recurrent somatic mutations in splice factors have been shown to generate widespread alterations in splicing specificity (Graubert et al., 2012; Shirai et al., 2015; Yoshida et al., 2011). Much work remains to be done to understand the regulation of splicing in different tissues and cellular contexts, especially in the case of mutations that directly affect proteins in the spliceosome. In light of recent advances in oligonucleotide therapeutics that could potentially target splicing defects in a sequence-specific manner ( Finkel et al., 2017 ), a better understanding of the regulatory mechanisms governing this remarkable process may pave the way for novel candidates for therapeutic intervention.

[0449] Figures 37A, 37B, 37C, 37D, 37E, 37F, 37G, and 37H show one implementation of predicting splicing from primary sequence by deep learning.

[0450] With reference to Figure 37A, for each position within the pre-mRNA transcript, SpliceNet-10k uses 10,000 nucleotides of the flanking sequence as input and predicts whether the position is a splice acceptor, a splice donor, or neither.

[0451] With reference to Figure 37B, the complete pre-mRNA transcript for the CFTR gene scored using MaxEntScan (top) and SpliceNet-10k (bottom) is shown, along with the predicted acceptor (red arrow) and donor (green arrow) sites and the actual locations of the exons (black boxes). For each method, we applied a threshold that made the number of predicted sites equal the total number of actual sites.

[0452] For each exon, we measured the exon inclusion rate on RNA-seq and show the SpliceNet-10k score distribution for exons at different inclusion rates. Shown are the maximum acceptor and donor scores for the exon.

[0453] With reference to Figure 37D, the effect of computer-simulated mutation of each nucleotide around exon 9 in the U2SURP gene is shown. The vertical size of each nucleotide indicates the predicted decrease in strength (Δ score) of the acceptor site (black arrow) when that nucleotide is mutated.

[0454] With reference to Figure 37E, the effect of the size of the input sequence construct on the accuracy of the network is shown. Top-k precision is the proportion of correctly predicted splice sites at a threshold where the number of predicted sites equals the actual number of sites present. PR-AUC is the area under the precision-recall curve. We also show the top-k precision and PR-AUC of three other algorithms for splice site detection.

[0455] With reference to Figure 37F, the relationship between exon / intron length and the strength of adjacent splice sites as predicted by SpliceNet-80nt (local motif score) and SpliceNet-10k is shown. Genome-wide distributions of exon length (yellow) and intron length (pink) are shown in the background. The x-axis is on a logarithmic scale.

[0456] With reference to Figure 37G, pairs of splice acceptor and donor motifs are spaced 150 nt apart and run along the HMGCR gene. Illustrated at each position is the K562 nucleosome signal and the likelihood of the pair forming an exon at that position, as predicted by SpliceNet-10k.

[0457] For Figure 37H, the average K562 and GM12878 nucleosome signals near private mutations predicted by the SpliceNet-10k model to form novel exons in the GTEx cohort are shown. The p-values ​​from the permutation test are shown.

[0458] Figures 38A, 38B, 38C, 38D, 38E, 38F, and 38G show one implementation of validation of rare cryptic splice mutations in RNA sequencing data.

[0459] With reference to Figure 38A, to assess the splice-altering impact of a mutation, SpliceNet-10k predicts acceptor and donor scores at each position within the pre-mRNA sequence of a gene with and without the mutation, as shown here for rs397515893, a potentially pathogenic splice variant in the MYBPC3 intron associated with cardiomyopathy. The delta score value for a mutation is the maximum change in splice prediction score within 50 nt of the variant.

[0460] For Figure 38B, we scored private genetic variants (observed in 1 of 149 individuals in the GTEx cohort) using the SpliceNet-10k model. Shown is the enrichment for private variants predicted to alter splicing (Δ score > 0.2, blue) or not affect splicing (Δ score < 0.01, red) near private exon skipping junctions (top) or private acceptor and donor sites (bottom). The y-axis shows the number of times a private splice event and nearby private genetic variants co-occur in the same individual compared to the expected number obtained through substitution.

[0461] With reference to Figure 38C, an example of a heterozygous synonymous variant in PYGB that forms a novel donor site with incomplete penetrance is shown. RNA-seq coverage, junction read counts, and junction location (blue and gray arrows) are illustrated for individuals with the variant and control individuals. Effect sizes are calculated as the difference in novel junction (AC) usage between individuals with and without the variant. In the stacked bar graph below, we show the number of reads with reference or alternative alleles that used the annotated junction or novel junction ("No Splice" and "Novel Junction," respectively). The total number of reference reads was significantly different from the total number of alternative reads (P = 0.018, binomial test), suggesting that 60% of transcripts splicing at novel junctions are missing in the RNA-seq data, likely due to nonsense-mediated decay (NMD).

[0462] With reference to Figure 38D, the percentage of potential splice mutations predicted by the SpliceNet-10k model validated against GTEx RNA-seq data is shown. Validation rates for essential acceptor or donor dinucleotide (dashed line) cleavages are less than 100% due to coverage and nonsense-mediated decay.

[0463] With reference to Figure 38E, the distribution of effect sizes for validated cryptic splice predictions is shown. The dashed line (50%) corresponds to the expected effect size of a fully penetrant heterozygous variant. The measured effect size of essential acceptor or donor dinucleotide cleavage is less than 50% due to nonsense-mediated decay or unexplained isoform changes.

[0464] With reference to Figure 38F, the sensitivity of SpliceNet-10k in detecting splice-altering private variants within the GTEx cohort at different delta score cutoffs is shown. Variants are divided into deep intronic variants (>50 nt from an exon) and near-exon variants (≦50 nt from an overlapping exon or exon-intron boundary).

[0465] Referring to Figure 38G, the validation rate and sensitivity of SpliceNet-10k and three other methods for splice site prediction at different confidence cutoffs are shown. The three points on the SpliceNet-10k curve show the performance of SpliceNet-10k at ΔScore cutoffs of 0.2, 0.5, and 0.8. For the other three algorithms, the three points on the curve show performance at ΔScore cutoffs of 0.2, 0.5, and 0.8, thresholds that predict the same number of potential splice variants as SpliceNet-10k.

[0466] Figures 39A, 39B, and 39C show one implementation in which cryptic splice variants frequently form tissue-specific alternative splicing.

[0467] Referring to Figure 39A, an example of a heterozygous exonic variant in CDC25B forming a novel donor site is shown. This variant is private to a single individual in the GTEx cohort and indicates tissue-specific alternative splicing that favors a greater proportion of the novel splice isoform in muscle compared to fibroblasts (P = 0.006 by Fisher's exact test). RNA-seq coverage, splice read counts, and splice location (blue and gray arrows) are illustrated for individuals with the variant in both muscle and fibroblasts and for control individuals.

[0468] Referring to Figure 39B, an example of a heterozygous exon acceptor-forming variant in FAM229B is shown, showing consistent tissue-specific effects across all three individuals in the GTEx cohort who carry the variant. RNA-seq for arteries and lungs is shown for three individuals with the variant and a control individual.

[0469] With reference to Figure 39C, the proportion of splice-site-forming variants within the GTEx cohort associated with significantly heterogeneous use of novel junctions across expressing tissues, as assessed by a chi-squared test for homogeneity, is shown. Validated cryptic splice variants with low to moderate Δ score values ​​were more likely to result in tissue-specific alternative splicing (P=0.015, Fisher's exact test).

[0470] Figures 40A, 40B, 40C, 40D, and 40E show one implementation in which predicted potential splice variants have a strong adverse effect in the human population.

[0471] With reference to Figure 40A, synonymous and intronic variants (≤50 nt from a known exon-intron boundary and excluding essential GT and AG dinucleotides) with confidently predicted splice alteration effects (Δ score ≥ 0.8) are strongly depleted in common allele frequencies (≥ 0.1%) in the human population for rare variants observed only once in 60,706 individuals, with an odds ratio of 4.58 (P < 10 by chi-square test). -127 ) show that 78% of recently emerged predicted cryptic splice variants have sufficient deleterious effects to be eliminated by natural selection.

[0472] With respect to Figure 40B, the proportion of predicted synonymous and intronic cryptic splice variants in the ExAC dataset with protein truncating variants and deleterious effects, calculated as in (A), is shown.

[0473] With respect to Figure 40C, the proportion of synonymous and intronic potential splice gain variants in the ExAC dataset with adverse effects is shown, split based on whether the variant is predicted to cause a frameshift (Δ score ≥ 0.8).

[0474] With respect to Figure 40D, the percentage of predicted deep intronic (>50 nt from a known exon-intron boundary) potential splice variants in the gnomAD dataset that have protein truncation variants and adverse effects is shown.

[0475] With respect to Figure 40E, the average number of rare (gene frequency <0.1%) protein-truncating variants and rare functional cryptic splice variants per individual human genome is shown. The number of cryptic splice mutations predicted to be functional is estimated based on the proportion of predictions with adverse effects. The total number of predictions is higher.

[0476] Figures 41A, 41B, 41C, 41D, 41E, and 41F show one implementation of de novo cryptic splice mutations in patients with rare genetic diseases.

[0477] With reference to Figure 41A, predicted potential splice de novo mutations per individual are shown for patients from the Deciphering Developmental Disorders cohort (DDD), individuals with autism spectrum disorder (ASD) from the Simons Simplex Collection and the Autism Sequencing Consortium, and healthy controls. Enrichment in the DDD and ASD cohorts over healthy controls is illustrated, adjusting for variant ascertainment between cohorts. Error bars indicate 95% confidence intervals.

[0478] With reference to Figure 41B, the estimated proportion of pathogenic de novo mutations by functional category for the DDD and ASD cohorts based on enrichment for each category compared to healthy controls is shown.

[0479] With reference to FIG. 41C, the enrichment and excess of potential splice de novo mutations within the DDD and ASD cohorts compared to healthy controls at different delta score thresholds is shown.

[0480] With reference to Figure 41D, a list of novel candidate disease genes enriched for de novo mutations in the DDD and ASD cohorts (FDR<0.01) is shown when predicted cryptic splice mutations were included together with protein-coding mutations in the enrichment analysis. Phenotypes that were present in multiple individuals are depicted.

[0481] With reference to Figure 41E, three examples of predicted de novo cryptic splice mutations in autism patients, validated by RNA-seq, are shown, resulting in intron retention, exon skipping, and exon extension, respectively. For each example, the RNA-seq coverage and junction counts for affected individuals are shown at the top, and for control individuals without the mutation are shown at the bottom. Sequences are shown on the sense strand relative to the gene transcript. Blue and gray arrows indicate the location of the junctions in individuals with the variant and control individuals, respectively.

[0482] With reference to Figure 41F, the validation status for 36 predicted potential splice sites selected for experimental validation by RNA-seq is shown.

[0483] Experimental model and subject details Subject details for 36 autistic individuals have been previously published by Iossifov et al., Nature 2014 (Table S1) and can be cross-referenced using the anonymous identifiers in column 1 of Table S4 in our paper.

[0484] Method details I. Deep Learning for Splice Prediction SpliceNet Architecture We trained several ultra-deep convolutional neural network-based models to computationally predict splicing from pre-mRNA nucleotide sequences. We designed four architectures, SpliceNet-80nt, SpliceNet-400nt, SpliceNet-2k, and SpliceNet-10k, that use 40, 200, 1,000, and 5,000 nucleotides as input on each side of a position of interest, respectively, and output the probabilities that the position is a splice acceptor and donor. More precisely, the input to the model is a one-hot encoded sequence of nucleotides, with A, C, G, and T (or equivalently U) encoded as [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], and [0, 0, 0, 1], respectively, and the output of the model consists of three scores that add up to 1, corresponding to the probability that the position in question is a splice acceptor, a splice donor, or neither.

[0485] The basic unit of the SpliceNet architecture is the residual block (He et al., 2016b), which consists of a batch normalization layer (Ioffe and Szegedy, 2015), a rectified linear unit (ReLU), and convolutional units organized in a specific manner (Figures 21, 22, 23, and 24). Residual blocks are commonly used when designing deep neural networks. Before the development of residual blocks, deep neural networks consisting of many convolutional units stacked one after another were very difficult to train due to problems with exploding / vanishing gradients (Glorot and Bengio, 2010), and increasing the depth of such neural networks often resulted in higher training errors (He et al., 2016a). Through a comprehensive set of computational experiments, an architecture consisting of many residual blocks stacked one after another was shown to overcome these problems (He et al., 2016a).

[0486] The complete SpliceNet architecture is presented in Figures 21, 22, 23, and 24. The architecture consists of K stacked residual blocks connecting the input layer to the penultimate layer, and a convolutional unit with softmax activation connecting the penultimate layer to the output layer. The residual blocks are stacked such that the output of the i-th residual block is connected to the input of the i+1-th residual block. Furthermore, the output of every fourth residual block is appended to the input of the penultimate layer. Such "skip connections" are commonly used in deep neural networks to increase convergence speed during training (Oord et al., 2016).

[0487] Each residual block has three hyperparameters N, W, and D, where N represents the number of convolution kernels, W represents the window size, and D represents the dilation rate of each convolution kernel (Yu and Koltun, 2016). Since a convolution kernel with window size W and dilation rate D extracts features spanning (W-1)D neighboring positions, a residual block with hyperparameters N, W, and D extracts features spanning 2(W-1)D neighboring positions. Therefore, the total neighborhood coverage of the SpliceNet architecture is

number

[0488] The SpliceNet architecture has only convolutional units plus normalization and nonlinear activation units. As a result, the model can be used in sequence-to-sequence mode with variable sequence lengths (Oord et al., 2016). For example, the SpliceNet-10k model (S = 10,000) encodes one-hot encoded nucleotide sequences of length S / 2 + l + S / 2, and the output is an l × 3 matrix corresponding to the three scores of l central positions in the input, i.e., the positions remaining after excluding the first and last S / 2 nucleotides. This feature can be exploited to achieve enormous computational savings in training and even testing. This is due to the fact that most of the computations for positions close to each other are common, and the shared computations need only be performed once by the model when used in sequence-to-sequence mode.

[0489] Our model employs a residual block architecture, which has become widely used due to its success in image classification. The residual block contains repeated units of convolution interspersed with skip connections, which allow information from earlier layers to skip the residual block. In each residual block, the input layer is first batch normalized, followed by an activation layer using rectified linear units (ReLU). The activations are then passed through a 1D convolutional layer. This intermediate output from the 1D convolutional layer is again batch normalized and ReLU activated, followed by another 1D convolutional layer. At the end of the second 1D convolution, we sum the output with the original input within 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 the authors call a deep residual learning network, the input is kept in its original state, and the residual connections have no nonlinear activation from the model, allowing for the effective training of deeper networks.

[0490] Following the residual block, a softmax layer calculates the probability of three states for each amino acid, among which the maximum softmax probability determines the state of the amino acid. The model is trained using the ADAM optimizer with a cumulative multi-class cross-entropy loss function for all protein sequences.

[0491] Atrous / dilated convolution allows for large receptive fields with few trainable parameters. Atrous / dilated convolution is a convolution in which the kernel is applied over a region larger than its length by skipping input values ​​using a step, also called the atrous convolution rate or dilation factor. Atrous / dilated convolution adds spacing between elements of the convolution filter / kernel, so that neighboring input entries (e.g., nucleotides, amino acids) at a larger interval are considered when the convolution operation is performed. This allows long-range compositional dependencies to be incorporated into the input. Atrous convolution saves partial convolution calculations for reuse when neighboring nucleotides are processed.

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

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

[0494]

number

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

[0496] Training and testing the model We downloaded the GENCODE (Harrow et al., 2012) V24lift37 gene annotation table from the UCSC table browser, extracted annotations for 20,287 protein-coding genes, and selected the primary transcript when multiple isoforms were available. We removed genes that did not have splice junctions and divided the remaining genes into training and test set genes as follows: Genes belonging to chromosomes 2, 4, 6, 8, 10-22, X, and Y were used to train the model (13,384 genes, 130,796 donor-acceptor pairs). We randomly selected 10% of the training genes and used them to determine early stopping points during training, while the remainder were used to train the model. To test the model, we used genes from chromosomes 1, 3, 5, 7, and 9 that did not have paralogs (1,652 genes, 14,289 donor-acceptor pairs). For this purpose, we referred to the human gene paralog list from http: / / grch37.ensembl.org / biomart / martview.

[0497] We trained and tested a model in sequence-sequence mode with chunks of size l = 5,000 using the following procedure: For each gene, the mRNA transcript sequence between the canonical transcription start and end sites was extracted from the hg19 / GRCh37 assembly. The input mRNA transcript sequence was one-hot encoded as follows: A, C, G, and T / U were mapped to [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], and [0, 0, 0, 1], respectively. The one-hot encoded nucleotide sequence was zero-padded until its length was a multiple of 5,000, and then further zero-padded at the beginning and end with flanking sequences of length S / 2, where S is equal to 80, 400, 2,000, and 10,000 for the SpliceNet-80nt, SpliceNet-400nt, SpliceNet-2k, and SpliceNet-10k models, respectively. The padded nucleotide sequence was then divided into blocks of length S / 2 + 5,000 + S / 2, such that the i-th block consisted of nucleotide positions 5,000(i-1)-S / 2 + 1 to 5,000i + S / 2. Similarly, the splice output tag sequence was one-hot encoded as follows: The splice site, splice acceptor (the first nucleotide of the corresponding exon), and splice donor (the last nucleotide of the corresponding exon) were mapped to [1, 0, 0], [0, 1, 0], and [0, 0, 1], respectively. The one-hot-encoded splice output tagged sequence was zero-padded to a length that was a multiple of 5,000 and then divided into blocks of length 5,000, such that the i-th block consisted of positions 5,000(i-1) + 1 to 5,000i. The one-hot-encoded nucleotide sequence and the corresponding one-hot-encoded tagged sequence were used as the input to the model and the target output of the model, respectively.

[0498] The models were trained on two NVIDIA GeForce GTX 1080 Ti GPUs with a batch size of 12 for 10 epochs. The multicategory cross-entropy loss between the target output and the predicted output was minimized using the Adam optimizer (Kingma and Ba, 2015) during training. The optimizer's learning rate was set to 0.001 for the first six epochs and then reduced by a factor of 2 for each subsequent epoch. For each architecture, we repeated the training procedure five times and obtained five trained models (Figures 53A and 53B). During testing, each input was evaluated using all five trained models, and the average of their outputs was used as the predicted output. We used these models for the analysis in Figure 37A and other related figures.

[0499] For the analyses shown in Figures 38A-G, 39A-C, 40A-E, and 41A-F, which involve the identification of splice-altering variants, we augmented the training set of GENCODE annotations to also include novel splice junctions commonly observed in the GTEx cohort on chromosomes 2, 4, 6, 8, 10-22, X, and Y (67,012 splice donors and 62,911 splice acceptors), which increased the number of splice junction annotations in the training set by ~50%. Training a network on the combined dataset improved the sensitivity of detecting splice change variants in RNA-seq data compared to a network trained on GENCODE annotations alone (Figures 52A and 52B), particularly for predicting deep intron splice change variants, and we used this network for analyses involving variant evaluation (Figures 38A-G, 39A-C, 40A-E, and 41A-F and related figures). To ensure that the GTEx RNA-seq dataset did not contain overlap between training and evaluation, we included only junctions present in five or more individuals in the training dataset and only evaluated the network's performance on variants present in four or fewer individuals. Details of novel splice junction identification are described in the "Splice Junction Detection" section of the GTEx Analysis section of the Methods.

[0500] Top-k accuracy Accuracy metrics such as the percentage of correctly classified positions are largely ineffective due to the fact that the majority of positions are not splice sites. We instead evaluated our models using two metrics that are valid in such settings: top-k precision and the area under the precision-recall curve. Top-k precision for a particular class is defined as follows: Suppose the test set has k positions that belong to a class. We choose a threshold such that exactly k test set positions are predicted as belonging to that class. The proportion of these k predicted positions that truly belong to this class is reported as Top-k precision. In effect, this is equivalent to precision when the threshold is chosen so that precision and recall have the same value.

[0501] Model evaluation on lincRNA We obtained a list of all lincRNA transcripts based on the GENCODE V24lift37 annotation. Unlike protein-coding genes, lincRNAs are not assigned a primary transcript in the GENCODE annotation. To minimize redundancy in the validation set, we identified the transcript with the longest total exon sequence for each lincRNA gene and called this the canonical transcript for the gene. Because lincRNA annotations are expected to be less reliable than those for protein-coding genes, and because such misannotations affect our estimates of Top-k accuracy, we eliminated lincRNAs with potential annotation issues using GTEx data (see the "Analysis on the GTEx Dataset" section below for details on these data). For each lincRNA, we counted all split reads that mapped across the length of the lincRNA across all GTEx samples (see "Splice Junction Detection" below for details). This was an estimate of all junction-spanning reads of the lincRNA, using either annotated or novel junctions. We also counted the number of reads spanning canonical transcript junctions. We considered only lincRNAs for which at least 95% of junction-spanning reads across all GTEx samples corresponded to the canonical transcript. We also required that all junctions of the canonical transcript were observed at least once within the GTEx cohort (excluding intron-spanning junctions <10 nt in length). To calculate Top-k accuracy, we considered only canonical transcript junctions of lincRNAs that passed the above filters (781 transcripts, 1047 junctions).

[0502] Identifying splice junctions from pre-mRNA sequences Figure 37B compares the performance of MaxEntScan and SpliceNet-10k in identifying canonical exon boundaries of genes from sequences. We used the CFTR gene, which is in our test set and has 26 canonical splice acceptors and donors, as a case study. We used MaxEntScan and SpliceNet-10k to obtain acceptor and donor scores for each of 188,703 positions from the canonical transcription start site (chr7:117,120,017) to the canonical transcription end site (chr7:117,308,719). Positions were classified as splice acceptors or donors if their corresponding scores were greater than a threshold selected while assessing Top-k accuracy. MaxEntScan predicted 49 splice acceptors and 22 splice donors, of which 9 and 5 were true splice acceptors and donors, respectively. For better visualization, we show the MaxEntScan pre-log scores (clipped to a maximum of 2,500). SpliceNet-10k predicted 26 splice acceptor and 26 splice donor sites, all of which were correct. In Figure 42B, we repeated the analysis using the LINC00467 gene.

[0503] Estimation of exon inclusion at GENCODE-annotated splice junctions We calculated the inclusion rate of all GENCODE-annotated exons from the GTEx RNA-seq data (Figure 37C). For each exon, excluding the first and last exons of each gene, we calculated the inclusion rate as follows:

[0504]

number

[0505] where L is the total read count of the junction from the previous canonical exon to the exon under consideration across all GTEx samples, R is the total read count of the junction from the exon under consideration to the next canonical exon, and S is the total read count of the skipping junction from the previous canonical exon to the next canonical exon.

[0506] The significance of various nucleotides for splice site recognition In Figure 37D, we identify nucleotides considered important by SpliceNet-10k toward classifying a position as a splice acceptor. To do this, we considered the splice acceptor at chr3:142,740,192 in the U2SURP gene, which is in our test set. The "importance score" of a nucleotide with respect to a splice acceptor is defined as: s ref Let s denote the acceptor score of the splice acceptor under consideration. The acceptor score is recalculated by replacing the nucleotide under consideration with A, C, G, and T. These scores are denoted as s A , s C , s G , and s T The importance score of a nucleotide is estimated as follows:

[0507]

number

[0508] This procedure is often referred to as in silico mutagenesis (Zhou and Troyanskaya, 2015). We plotted 127 nucleotides from chr3:142,740,137 to chr3:142,740,263, with the height of each nucleotide representing the importance score for the splice acceptor at chr3:142,740,192. The plotting function was adapted from DeepLIFT (Shrikumar et al., 2017) software.

[0509] Effect of TACTAAC and GAAGAA motifs on splicing To study the effect of branchpoint sequence position on acceptor strength, we first obtained acceptor scores for 14,289 test set splice acceptors using SpliceNet-10k. ref Let denote a vector containing these scores. For each value of i ranging from 0 to 100, we did the following: For each test set splice acceptor, we replaced the nucleotides at positions i through i-6 before the splice acceptor with TACTAAC and recalculated the acceptor score using SpliceNet-10k. The vector containing these scores is denoted by y alt,i In Figure 43A we plot the following quantities as a function of i: mean(y alt,i -y ref )

[0510] In Figure 43B, we repeated the same procedure using the SR-protein motif GAAGAA. In this case, we also studied the effect of the motif when present after a splice acceptor, as well as its effect on donor strength. GAAGAA and TACTAAC were the motifs with the greatest effect on acceptor and donor strength, based on a comprehensive search within k-mer space.

[0511] The role of exon and intron length in splicing To examine the effect of exon length on splicing, we filtered out test set exons that were either the first or last exon. This filtering step removed 1,652 of the 14,289 exons. We sorted the remaining 12,637 exons by increasing length. For each of them, we calculated a splicing score by averaging the acceptor score at the splice acceptor site and the donor score at the splice donor site using SpliceNet-80nt. We plot the splicing score as a function of exon length in Figure 37F. Before plotting, we applied the following smoothing procedure: Let x represent a vector containing the exon length, and y represent a vector containing its corresponding splicing score. We smoothed both x and y using an averaging window of size 2,500.

[0512] We repeated this analysis by calculating splicing scores using SpliceNet-10k. In the background, we illustrate a histogram of the lengths of the 12,637 exons considered in this analysis. We applied a similar analysis to examine the effect of intron length on splicing, with the key difference being that we did not need to exclude the first and last exons.

[0513] The role of nucleosomes in splicing We downloaded nucleosome data for the K562 cell line from the UCSC genome browser. We used the HMGR gene in our test set as an example to demonstrate the impact of nucleosome positioning on the SpliceNet-10k score. For each position p within the gene, we calculated a "plant splicing score" as follows: · Eight nucleotides from positions p+74 to p+81 were replaced by the donor motif AGGTAAGG. · Four nucleotides at positions p-78 to p-75 were replaced by the acceptor motif TAGG. · 20 nucleotides from positions p-98 to p-79 were replaced by the polypyrimidine tract CCTCCTTTTTCCTCGCCCTC. Seven nucleotides at positions p-105 to p-99 were replaced by the branch point sequence CACTAAC. The average of the acceptor score at p-75 and the donor score at p+75 predicted by SpliceNet-10k was used as the plant splicing score.

[0514] The plant splicing scores for the K562 nucleosome signal as well as 5,000 positions from chr5:74,652,154 to chr5:74,657,153 are shown in Figure 37G.

[0515] To calculate the genome-wide Spearman correlation between these two tracks, we randomly selected 1,000,000 intergenic positions that were at least 100,000 nt away from all canonical genes. For each of these positions, we calculated the plant splicing score as well as the average K562 nucleosome signal (a window size of 50 was used for averaging). The correlation between these two values ​​across 1,000,000 positions is shown in Figure 37G. We further subclassified these positions based on GC content (estimated using the nucleotides between the plant acceptor and donor motifs) with a bin size of 0.02. We show the genome-wide Spearman correlation for each bin in Figure 44A.

[0516] For each of the 14,289 test set splice acceptors, we extracted nucleosome data within 50 nucleotides on each side and calculated its nucleosome enrichment as the average signal on the exon side divided by the average signal on the intron side. We sorted the splice acceptors by increasing nucleosome enrichment and calculated their acceptor scores using SpliceNet-80nt. The acceptor scores are plotted as a function of nucleosome enrichment in Figure 44B. Before plotting, the smoothing procedure used in Figure 37F was applied. We repeated this analysis for the 14,289 test set splice donors using SpliceNet-10k.

[0517] Enrichment of nucleosome signals in novel exons In Figure 37H, we wanted to look at the nucleosome signals around predicted novel exons. To ensure we were looking at reliable novel exons, we selected only singleton variants (variants present in a single GTEx individual) whose predicted gain-of-junction was completely private to the individual carrying the variant. In addition, to remove confounding effects from nearby exons, we only looked at intronic variants that were at least 750 nt away from the annotated exon. We downloaded the nucleosome signals for the GM12878 and K562 cell lines from the UCSC browser and extracted the nucleosome signals within 750 nt of each predicted novel acceptor or donor site. We averaged the nucleosome signals between the two cell lines and flipped the signal vectors for variants that overlapped the gene on the minus strand. We shifted the signal from the acceptor site by 70 nt to the right and the signal from the donor site by 70 nt to the left. After shifting, the nucleosome signals for both the acceptor and donor sites were centered in the middle of an idealized exon 140 nt in length, which is the median exon length in the GENCODE v19 annotation. Finally, we smoothed the resulting signals by averaging all shifted signals and calculating the mean within an 11 nt window centered at each position.

[0518] To test for association, we selected random singleton SNVs located at least 750 nt away from annotated exons and predicted by the model to have no effect on splicing (Δ score < 0.01). We created 1,000 random samples of such SNVs, each with the same number of SNVs (128 sites) as the set of splice site gain sites used in Figure 37H. For each random sample, we calculated the smoothed mean signal as described above. Because random SNVs were not predicted to form novel exons, we centered the nucleosome signal from each SNV on itself and randomly shifted it either 70 nt to the left or 70 nt to the right. We then compared the nucleosome signal at the middle base in Figure 37H with the signal obtained from 1,000 simulations at that base. Empirical p-values ​​were calculated as the proportion of simulated sets that had a median value equal to or greater than the value observed for splice site gain variants.

[0519] Robustness of the network to differences in exon density To examine the generalizability of the network's predictions, we evaluated SpliceNet-10k in regions with varying exon density. We first divided the test set positions into five categories according to the number of canonical exons present within a 10,000-nucleotide window (5,000 nucleotides on each side) (Figure 54). To ensure that the exon count was an integer value for each position, we used the number of exon starts present within the window as a proxy. For each category, we calculated the Top-k precision and the area under the precision-recall curve. The number of positions and the value of k were different for different categories (detailed in the table below).

[0520] [Table 2]

[0521] Robustness of the network for each of the five models in the ensemble Training multiple models and using the average of their predictions as the output is a common strategy in machine learning to obtain better predictive performance, called ensemble learning. In Figure 53A, we show the Top-k accuracy and area under the precision-recall curves of the five SpliceNet-10k models we trained to construct the ensemble. The results clearly demonstrate the stability of the training process.

[0522] We also calculated the Pearson correlation between predictions. Because most positions in the genome are not splice sites, the correlation between most model predictions would be close to 1, rendering the analysis meaningless. To overcome this problem, we considered only positions in the test set that were assigned an acceptor or donor score of 0.01 or greater by at least one model. This criterion was met for 53,272 positions (roughly equal numbers of splice and non-splice sites). These results are summarized in Figure 53B. The very high Pearson correlation between the model predictions further illustrates their robustness.

[0523] We show the effect of the number of models used to construct the ensemble on performance in Figure 53 C. These results show that performance improves as the number of models increases, but with diminishing returns.

[0524] II. Analysis on the GTEx RNA-seq dataset Delta score for single nucleotide variants We quantified splicing changes due to single-nucleotide variants as follows: We first used a reference nucleotide to calculate acceptor and donor scores for 101 positions around the variant (50 positions on each side). These scores are then respectively expressed as vectors a ref and d refWe then recalculated the acceptor and donor scores using the alternative nucleotides. These scores are represented by the vector a alt and d alt It is assumed that the sine wave is represented by the following formula: We evaluated four quantities: Δ score (acceptor gain) = max(a alt -a ref ) Δ score (acceptor loss) = max(a ref -a alt ) Δ score (donor gain) = max(d alt -d ref ) Δ score (donor loss) = max(d ref -d alt )

[0525] The maximum of these four scores is called the delta score of the variant.

[0526] Variant quality control and filtering criteria We downloaded GTEx VCF and RNA-seq data from dbGaP (study accession phs000424.v6.p1; https: / / www.ncbi.nlm.nih.gov / projects / gap / cgi-bin / study.cgi?study_id=phs000424.v6.p1).

[0527] We evaluated the performance of SpliceNet on autosomal SNVs that occurred in at most four individuals in the GTEx cohort. In particular, a variant was considered if it met the following criteria in at least one individual A: 1. The variant was not filtered (the FILTER field in the VCF was PASS). 2. The variant was not marked as MULTI_ALLELIC in the INFO field of individual A's VCF, and the VCF contained a single allele in the ALT field. 3. Individual A was heterozygous for the variant. 4. The ratio alt_depth / (alt_depth + ref_depth) is between 0.25 and 0.75, where alt_depth and ref_depth are the number of reads supporting the alternative and reference alleles in individual A, respectively. 5. The total depth, alt_depth + ref_depth, was between 20 and 300 in the VCF of individual A. 6. The variant overlapped the gene body region. The gene body was defined as the region between the start and end of transcription of the canonical transcript from GENCODE (V24lift37).

[0528] For variants that met these criteria in at least one individual, we considered all individuals in which the variant occurred (even if they did not meet the above criteria) to have the variant. We refer to variants that occur in a single individual as singletons and variants that occur in two to four individuals as common. We did not evaluate variants that occur in five or more individuals to avoid overlap with the training dataset.

[0529] RNA-seq read alignment We used OLego (Wu et al., 2013) to map the reads of GTEx samples to the hg19 reference, allowing an edit distance of at most 4 between the query read and the reference (parameter -M 4). Note that OLego can operate completely de novo and does not require gene annotation. Because OLego examines the presence of splicing motifs at the ends of split reads, its alignments may be biased toward or against the reference around SNVs that cut or create splice sites, respectively. To eliminate such bias, we further created alternative reference sequences for each GTEx individual by inserting all SNVs of the individual into the hg19 reference with a PASS filter. We used OLego with the same parameters to map all samples from each individual to that individual's alternative reference sequence. For each sample, we then combined the two sets of alignments (against the hg19 reference and against the individual's alternative reference) by picking the best alignment for each read pair. To select the best alignment for a read pair P, we used the following procedure. 1. If both reads of P are unmapped in both sets of alignments, we randomly choose an alternative alignment of hg19 or P. 2. If P had more unmapped ends in one set of alignments than the others (e.g., both ends of P were mapped to an alternative reference, but only one end was mapped to hg19), we chose the alignment in which both ends of P were mapped. 3. If both ends of P map in both sets of alignments, we choose the alignment with the smallest total mismatches, or a random one if the number of mismatches is the same.

[0530] Detecting splice junctions in aligned RNA-seq data We detected and counted splice junctions in each sample using leafcutter_cluster (Li et al., 2018), a utility in the leafcutter package. We required that a single split read support a junction and assumed a maximum intron length of 500 Kb (parameters -m 1 -l 500000). To obtain a high-confidence set of junctions for training the deep learning model, we compiled a union of all leafcutter junctions across all samples and removed from consideration junctions that met any of the following criteria: 1. Either end of the junction overlapped an ENCODE blacklist region (table wgEncodeDacMapabilityConsensusExcludable in hg19 from the UCSC genome browser) or a simple repeat (Simple Repeats track in hg19 from the UCSC genome browser). 2. Both ends of the junction were on non-canonical exons (based on the canonical transcript from GENCODE version V24lift37). 3. The two ends of the junction were on different genes, or one end was in a non-genic region. 4. Either end lacked the essential GT / AG dinucleotide.

[0531] Junctions that were present in five or more individuals were used to augment the list of GENCODE-annotated splice junctions for analysis of variant prediction (Figures 38A-G, 39A-C, 40A-E, and 41A-F). A link to a file containing the list of splice junctions used to train the model is provided in the Key Resources table.

[0532] Although we augmented the training dataset using junctions detected by leafcutter, we found that leafcutter filtered many junctions with good support in the RNA-seq data, even with relaxed parameters. This artificially reduced our validation rate. Therefore, for the GTEx RNA-seq validation analysis (Figures 38A-G and 39A-C), we recalculated the set of junctions and junction counts directly from the RNA-seq read data. We counted all non-overlapping split-mapped reads with at least 5 nt aligned on each side of the junction by MAPQ at least 10 times. Reads were allowed to span more than two exons; in that case, reads were counted toward each junction with at least 5 nt of sequence mapped on both sides.

[0533] Defining a Private Junction A mating was considered private in individual A if it met at least one of the following criteria: 1. The junction had at least three reads in at least one sample from A and was never observed in any other individual. 2. There were at least two organizations that met both of the following criteria: a. The average read count of the junction within samples from individual A within a tissue was at least 10. b. Individual A had, on average, at least 2x more normalized reads than any other individual in that tissue, where the normalized read count of a junction within a sample was defined as the number of reads at the junction normalized by the total number of reads across all junctions for the corresponding gene.

[0534] Tissues with fewer than five samples from other individuals (non-A) were ignored in this test.

[0535] Enrichment of singleton SNVs around private junctions If a private junction had exactly one end annotated, based on the GENCODE annotation, we considered it a candidate for acceptor or donor gain and searched for singleton SNVs (SNVs occurring in a single GTEx individual) that were private in the same individual within 150 nt of the unannotated end. If a private junction had both ends annotated, we considered it a candidate for a private exon skipping event if it skipped at least one but no more than three exons of the same gene based on the GENCODE annotation. We then searched for singleton SNVs within 150 nt of each end of the skipped exon. Private junctions that did not have both ends in the GENCODE exon annotation were ignored, as a substantial proportion of these were alignment errors.

[0536] To calculate the enrichment of singleton SNVs around a novel private acceptor or donor (Figure 38B, bottom), we tallied the counts of singleton SNVs at each position relative to the private junction. If the overlapping gene was on the minus strand, the relative positions were flipped. We divided the SNVs into two groups: SNVs that were private in the individual with the private junction and SNVs that were private in different individuals. To smooth the resulting signal, we averaged the counts in a 7-nt window centered at each position. We then calculated the ratio of the smoothed counts from the first group (private in the same individual) to the smoothed counts from the second group (private in different individuals). For novel private exon skipping (Figure 38B, top), we followed a similar procedure and tallied the counts of singleton SNVs around the end of the skipped exon.

[0537] Validation of model predictions within GTEx RNA-seq data For either private variants (occurring in one individual in the GTEx cohort) or common variants (occurring in two to four individuals in the GTEx cohort), we obtained the deep learning model's predictions for the reference and alternative alleles and calculated a delta score. We also obtained the configurations where the model predicted the aberrant (de novo or truncated) junction. We then sought to determine whether there was evidence in the RNA-seq data supporting splicing abnormalities in individuals carrying variants at the predicted configurations. In many cases, the model can predict multiple effects for the same variant; for example, a variant that truncates an annotated splice donor can also increase the utilization of a suboptimal donor, as in Figure 45. In that case, the model might predict both a donor loss at the annotated splice site and a donor gain at the suboptimal site. However, for validation purposes, we considered only the effect with the highest predicted delta score for each variant. Therefore, for each variant, we considered the predicted splice site formation and splice site cleavage effects separately. Junctions occurring in fewer than five individuals were excluded during model training to avoid evaluating the model on the novel junctions on which it was trained.

[0538] Validation of predicted cryptic splice mutations based on private splice junctions For each private variant predicted to cause de novo junction formation, we used the network to predict the location of the newly created aberrant splice junction and looked at the RNA-seq data to validate if such a de novo junction appeared only in the individual with the SNV and not in any other GTEx individuals. Similarly, for variants predicted to cause splice site loss affecting a splice site in exon X, we looked for de novo exon skipping events from the previous canonical exon (those upstream of X based on the GENCODE annotation) to the next canonical exon (those downstream of X) that appeared only in the individual with the variant and not in any other individuals in GTEx. We excluded predicted losses if the splice site predicted to be lost by the model was not annotated in GENCODE or was never observed in GTEx individuals without the variant. We also excluded predicted gains if the splice site predicted to be gained was already annotated in GENCODE. To extend this analysis to common variants (present in two to four individuals), we also validated novel junctions that were present in at least half of the individuals with the variant and absent in all individuals without the variant.

[0539] Using the requirement that predicted aberrant splice events be private to individuals carrying the variant, we were able to validate 40% of predicted high-score (Δ score ≥ 0.5) acceptor and donor gains, but only 3.4% of predicted high-score losses and 5.6% of essential GT or AG truncations (with a misvalidation rate of <0.2% based on substitutions—see the "Estimating Misvalidation Rates" section). The discrepancy in validation rates for gains and losses is twofold. First, unlike gains, exon skipping events are rarely globally private to individuals carrying the variant because exons are often skipped at a low baseline level, which can be observed with sufficiently deep RNA-seq. Second, splice site losses can have other effects in addition to enhancing exon skipping, such as increasing intron support or increasing the utilization of alternative suboptimal splice sites. For these reasons, we did not rely entirely on private de novo junctions to validate model predictions; we also validated variants based on quantitative evidence for increased or decreased usage of predicted affected junctions in individuals carrying the variant.

[0540] Validation of predicted cryptic splice mutations through quantitative criteria For a junction j from sample s, we calculate the normalized junction count c js was obtained.

[0541]

number

[0542] where r js is the raw junction count for junction j in sample s, and the sum in the denominator is taken over all other junctions between annotated acceptors and donors of the same gene as j (using annotations from GENCODE v19). The asinh transformation is

number

[0543] For each gain or loss junction j predicted to be caused by an SNV occurring in a set I of individuals, we calculated a z-score in each tissue t.

[0544]

number

[0545] where A t is the sample set from individuals in I in organization t, and U t is the set of samples from all other individuals in tissue t. Note that there may be multiple samples in the GTEx dataset for the same individual and tissue. As before, c jsis the count for junction j in sample s. For predicted losses, we also calculated a similar z-score for junction k that skips the exon predicted to be affected as follows:

[0546]

number

[0547] Note that losses that result in skipping will cause a relative decrease in loss junctions and a relative increase in skipping. This is due to the numerator z jt and z kt This justifies the reversal of the difference in , and therefore both of these scores tend to be negative for actual splice site losses.

[0548] Finally, we calculated the median z-score across all considered tissues. For losses, we calculated the median of each of the z-scores separately from equations (2) and (3). An acceptor or donor loss prediction was considered validated if any of the following was true: 1. The median z-score from equation (2), quantifying the relative loss of splicing, was below the 5th percentile (-1.46) of the corresponding value in the permutation data, and the median z-score from equation (3), quantifying the relative change in skipping, was non-positive (zero, negative, or missing, which is when the skipping splice was not observed in any individual). In other words, there was strong evidence for a reduction in the utilization of the affected splice and no evidence suggesting a reduction in skipping in affected individuals. 2. The median z-score from equation (3) was less than the 5th percentile (-0.74) of the corresponding values ​​in the permuted data, and the median z-score from equation (3) was non-positive. 3. The median z-score from equation (2) was less than the first percentile (-2.54) of the corresponding values ​​in the permuted data. 4. The median z-score from equation (3) was less than the first percentile (-4.08) of the corresponding values ​​in the permuted data. 5. Junctions skipping the affected exon are observed in at least half of the individuals with the variant and in none of the others (as described above in "Validation of predicted cryptic splice mutations based on private splice junctions").

[0549] A description of the permutations used to obtain the above cutoffs is provided in the section "Estimating Misvalidation Rates."

[0550] Empirically, we observed that, as explained in the section "Validating Predicted Cryptic Splice Mutations Based on Private Splice Junctions," more stringent validation criteria should be applied to losses compared to gains because losses tend to result in more mixed effects than gains. Observing a novel junction near a private SNV is unlikely to occur by chance, and therefore even faint evidence of a junction will be sufficient for validation. In contrast, most predicted losses result in weakening of existing junctions, and such weakening is more difficult to detect than on-off changes caused by gains and is more likely to be attributed to noise in the RNA-seq data.

[0551] Inclusion criteria for validation analyses To avoid calculating z-scores when counts were low or coverage was poor, we used the following criteria to filter variants for validation analysis: 1. A sample was considered for the calculation of the z-score above only if it expressed the gene (Σ g r gs >200). 2. Tissues were not considered for calculation of loss or gain z-scores if the mean count of loss or "reference" junctions, respectively, in individuals without the variant was less than 10. The "reference" junction is the canonical junction used before the gain of the novel junction, based on the GENCODE annotation (see the Effect Size Calculation section for details). The intuition is that we should not attempt to validate splice loss variants affecting junctions that are not expressed in control individuals. Similarly, we should not attempt to validate splice gain variants when control individuals did not sufficiently represent transcripts spanning the affected site. 3. In the case of predicted splice site loss, samples from individuals without the variant were considered only if they had at least 10 counts of the lost junction. In the case of predicted acceptor or donor site loss, samples from control individuals were considered only if they had at least 10 counts of the "reference" junction. The intuition is that even in tissues with large average representation of the affected junction (i.e., meeting criterion 2), different samples can have widely different sequencing depths, and therefore only control samples with sufficient representation should be included. 4. Tissues were considered only if there was at least one sample meeting the above criteria from an individual with the variant, as well as at least five samples meeting the above criteria from at least two distinct control individuals.

[0552] Variants for which no tissues met the criteria for consideration were considered non-confirmable and were excluded when calculating validation rates. For splice gain variants, we filtered those that appeared at existing GENCODE-annotated splice sites. Similarly, for splice loss variants, we considered only those that reduced the score of existing GENCODE-annotated splice sites. Overall, 55% and 44% of predicted gains and losses with high scores (Δ score ≥ 0.5), respectively, were considered confirmable and used in validation analyses.

[0553] Estimating the false validation rate To confirm that the above procedure has a reasonable true validation rate, we first looked at SNVs occurring in 1–4 GTEx individuals, truncating essential GT / AG dinucleotides. We argued that such mutations almost certainly affect splicing, resulting in a validation rate approaching 100%. Of such truncations, 39% were confirmable based on the criteria described above, and among those confirmable, the validation rate was 81%. To estimate the false validation rate, we permuted individual labels in the SNV data. For each SNV occurring in k GTEx individuals, we selected a random subset of k GTEx individuals and assigned the SNV to them. We created 10 such randomized datasets and repeated the validation process on them. Validation rates in the permuted datasets ranged from 1.7–2.1% for gains and 4.3–6.9% for losses, with medians of 1.8% and 5.7%, respectively. The high false validation rate for losses and the relatively low validation rate for essential truncations are due to the difficulty of validating splice site losses, as highlighted in the section "Validation of predicted cryptic splice mutations based on private splice junctions."

[0554] Calculating effect sizes of potential splice variants in RNA-seq data We defined the "effect size" of a variant as the proportion of transcripts of an affected gene whose splicing pattern was altered by the variant (e.g., the proportion that switched to a novel acceptor or donor). As a reference example for predicted splice gain variants, consider the variant in Figure 38C. For predicted gain donor A, we first identified the closest annotated junction (AC) to acceptor C. We then identified a "reference" junction (BC), where B≠A is the closest annotated donor to A. Then, in each sample s, we calculated the relative usage of the novel junction (AC) compared to the reference junction (BC) as follows:

[0555]

number

[0556] where r (AC)s is the raw read count of the junction (AC) in sample s. For each tissue, we calculated the change in junction (AC) usage between individuals carrying the variant and all other individuals as follows:

[0557]

number

[0558] where A t is the sample set from individuals with the variant in tissue t, and U tis the set of samples from other individuals in tissue t. The final effect size was calculat...

Claims

1. A system comprising: at least one processor; 1. A non-transitory computer-readable storage medium containing instructions that, when executed by the at least one processor, cause the system to: identifying a target nucleotide sequence comprising the target nucleotide and constituent nucleotides; providing the target nucleotide sequence as an input to an atrous convolutional neural network trained to generate splice site predictions for an input nucleotide sequence based on feature extraction, the atrous convolutional neural network including groups of residual blocks parameterized by the number of convolutional filters in the residual blocks of each group of residual blocks; generating a splice site prediction for the target nucleotide based on the constituent nucleotides in the target nucleotide sequence, using the atrous convolutional neural network to perform atrous convolution to extract learned features representative of a splicing pattern of the target nucleotide sequence; a non-transitory computer-readable storage medium for causing the Including, the system.

2. The system of claim 1, wherein the target nucleotide comprises a variant nucleotide.

3. The system described in claim 1, wherein the constituent nucleotides include at least 200 upstream nucleotides adjacent to the target nucleotide and at least 200 downstream nucleotides adjacent to the target nucleotide.

4. The system described in claim 1, wherein the constituent nucleotides include at least 5,000 upstream nucleotides adjacent to the target nucleotide and at least 5,000 downstream nucleotides adjacent to the target nucleotide.

5. The system of claim 1, further comprising instructions that, when executed by the at least one processor, cause the system to generate the splice site prediction for the target nucleotide by the atrous convolutional neural network by generating a splice site score for the target nucleotide and a non-splicing site score for the target nucleotide based on the constituent nucleotides in the target nucleotide sequence.

6. The system described in claim 5, further comprising instructions that, when executed by the at least one processor, cause the system to identify the highest score among the splice site score of the target nucleotide and the non-splicing site score of the target nucleotide.

7. The system of claim 1, wherein the group of residual blocks is further parameterized by a convolution window size of the residual block and an atrous convolution rate of the residual block.

8. The system described in claim 7, wherein the atrous convolution rate includes atrous convolution that stores partial convolution calculations so that they can be reused when adjacent nucleotides of the constituent nucleotide are processed by the atrous convolution neural network.

9. The system of claim 1, wherein the group of residual blocks is further parameterized by a number of residual blocks, a number of skip connections, and a number of residual connections.

10. Each group of residual blocks generates an intermediate output by processing a preceding input, and the dimension of the intermediate output is (I-[{(W-1)*D}*A])×N, where: I is the dimension of the preceding input, W is the convolution window size of the residual block for each group of residual blocks; D is the atrous convolution rate of the residual block of each group of residual blocks; A is the number of atrous convolutional layers in each group of residual blocks, The system of claim 1 , wherein N is the number of convolution filters in the residual block of each group of residual blocks.

11. A non-transitory computer-readable storage medium storing instructions that, when executed by at least one processor, cause a system to: identifying a target nucleotide sequence comprising the target nucleotide and constituent nucleotides; providing the target nucleotide sequence as an input to an atrous convolutional neural network trained to generate splice site predictions for an input nucleotide sequence based on feature extraction, the atrous convolutional neural network including groups of residual blocks parameterized by the number of convolutional filters in the residual blocks of each group of residual blocks; generating a splice site prediction for the target nucleotide based on the constituent nucleotides in the target nucleotide sequence, using the atrous convolutional neural network to perform atrous convolution to extract learned features representative of a splicing pattern of the target nucleotide sequence; A non-transitory computer-readable storage medium that causes 12. The non-transitory computer-readable storage medium of claim 11, wherein the target nucleotide comprises a variant nucleotide.

13. A non-transitory computer-readable storage medium as described in claim 11, wherein the constituent nucleotides include at least 200 upstream nucleotides adjacent to the target nucleotide and at least 200 downstream nucleotides adjacent to the target nucleotide.

14. A non-transitory computer-readable storage medium as described in claim 11, further storing instructions that, when executed by the at least one processor, cause the system to generate the splice site prediction for the target nucleotide by the atrous convolutional neural network by generating a splice site score for the target nucleotide and a non-splicing site score for the target nucleotide based on the constituent nucleotides in the target nucleotide sequence.

15. A non-transitory computer-readable storage medium as described in claim 14, further storing instructions that, when executed by the at least one processor, cause the system to identify the highest score among the splice site score of the target nucleotide and the non-splicing site score of the target nucleotide.

16. A computer-implemented method comprising: identifying a target nucleotide sequence comprising the target nucleotide and constituent nucleotides; providing the target nucleotide sequence as an input to an atrous convolutional neural network trained to generate splice site predictions for an input nucleotide sequence based on feature extraction, the atrous convolutional neural network including groups of residual blocks parameterized by the number of convolutional filters in the residual blocks of each group of residual blocks; generating splice site predictions for the target nucleotide based on the constituent nucleotides in the target nucleotide sequence using the atrous convolutional neural network that performs atrous convolution to extract learned features representative of the splicing pattern of the target nucleotide sequence; 20. A computer-implemented method comprising:

17. The computer-implemented method of claim 16, wherein the group of residual blocks is further parameterized by a convolution window size of the residual block and an atrous convolution rate of the residual block.

18. The computer-implemented method of claim 17, wherein the atrous convolution rate includes atrous convolution that stores partial convolution calculations so that they can be reused when adjacent nucleotides of the constituent nucleotide are processed by the atrous convolution neural network.

19. The computer-implemented method of claim 16, wherein the group of residual blocks is further parameterized by a number of residual blocks, a number of skip connections, and a number of residual connections.

20. Each group of residual blocks generates an intermediate output by processing a preceding input, and the dimension of the intermediate output is (I-[{(W-1)*D}*A])×N, where: I is the dimension of the preceding input, W is the convolution window size of the residual block for each group of residual blocks; D is the atrous convolution rate of the residual block of each group of residual blocks; A is the number of atrous convolutional layers in each group of residual blocks, 17. The computer-implemented method of claim 16, wherein N is the number of convolution filters in the residual block of each group of residual blocks.

Images