Aberrant splicing detection using convolutional neural networks (CNNs)

Abstract

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

Description

[0001] appendix

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

[0003] Priority application

[0004] This application claims priority or the benefit of U.S. Provisional Patent Application No. 62 / 573,125, filed on October 16, 2017, entitled “Deep Learning-Based Splice Site Classification,” by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM1001-1 / IP-1610-PRV); U.S. Provisional Patent Application No. 62 / 573,131, filed on October 16, 2017, entitled “Deep Learning-Based Aberrant Splicing Detection,” by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM 1001-2 / IP-1614-PRV); U.S. Provisional Patent Application No. 62 / 573,132, filed on October 16, 2017, entitled “Deep Learning-Based Aberrant Splicing Detection,” by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM 1001-3 / IP-1614-PRV); Farh, Sofia Kyriazopoulou Panagiotopoulou and Jeremy Francis McRae, entitled “Aberrant Splicing Detection Using Convolutional Neural Networks,” filed on October 16, 2017 (Attorney Docket No. ILLM1001-3 / IP-1615-PRV); and Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou and Jeremy Francis McRae, entitled “Prediction of Splicing from Raw Sequences Using Deep Learning,” filed on August 31, 2018 (Attorney Docket No. ILLM1001-10 / IP-1749-PRV). The provisional applications are hereby incorporated by reference herein for all purposes.

[0005] Incorporation

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

[0007] PCT Patent Application No. PCT / US18 / 55915, entitled “Deep Learning Based Splice Site Classification,” filed concurrently on October 15, 2018 by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM 1001-7 / IP-1610-PCT) and subsequently published as PCT Publication No. WO______.

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

[0009] U.S. nonprovisional patent application titled “Deep Learning Based Splice Site Classification” filed concurrently by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM1001-4 / IP-1610-US).

[0010] U.S. nonprovisional patent application titled “Deep Learning Based Aberrant Splicing Detection” filed concurrently by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM 1001-5 / IP-1614-US).

[0011] U.S. nonprovisional patent application titled “ABNORMAL SPLICING DETECTION USING CONVOLUTIONAL NEURAL NETWORKS (CNN)” filed concurrently by Kishore Jaganathan, Kai-How Farh, Sofia Kyriazopoulou Panagiotopoulou, and Jeremy Francis McRae (Attorney Docket No. ILLM 1001-6 / IP-1615-US).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[0039] Three table files in ASCII text format are submitted with this application and are incorporated by reference. The names, creation dates and sizes of the files are:

[0040] Table_S4_Mutation_Rate.txt August 31, 2018 2,452KB

[0041] Table_S5_gene_enrichment.txt 2018-08-31 362KB

[0042] Table_S6_Verification.txt August 31, 2018 362KB Technical Field

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

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

[0045] Machine Learning

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

[0047] Y=f(X)+∈

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

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

[0050]

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

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

[0053] Neural Networks

[0054] The single-layer perceptron (SLP) is the simplest model of a neural network. Figure 1 As shown in Figure 1, it consists of an input layer and an activation function. The input is passed through a weighted graph. The function f uses the sum of the inputs as an independent variable and compares it with a threshold θ.

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

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

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

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

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

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

[0061] Leung MK, Delong A, Alipanahi B et al. Machine Learning In GenomicMedicine: A Review of Computational Problems and Data Sets 2016; and LibbrechtMW, Noble WS. Machine Learning Applications In Genetics and Genomics. NatureReviews Genetics 2015;16(6):321-32. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0063] Figure 1 A single layer perceptron (SLP) is shown.

[0064] Figure 2 One implementation of a feed-forward neural network with multiple layers is shown.

[0065] Figure 3 Describes one implementation of how a convolutional neural network works.

[0066] Figure 4 Depicted is a block diagram of training a convolutional neural network in accordance with one implementation of the disclosed technology.

[0067] Figure 5 An implementation of a ReLU nonlinearity layer according to one implementation of the disclosed technology is shown.

[0068] Figure 6 A dilated convolution is shown.

[0069] Figure 7 is an implementation of a subsampling layer (average / max pooling) according to an implementation of the disclosed technology.

[0070] Figure 8 Depicts an implementation of two convolution layers in a convolutional layer.

[0071] Fig. 9 Depicted is a residual connection that reinjects original information downstream via feature map addition.

[0072] Fig.10 Depicts one implementation of residual blocks and skip connections.

[0073] Fig.11 One implementation of a stacked dilated convolution is shown.

[0074] Fig.12 The batch normalization forward pass is shown.

[0075] Fig.13 Shown is the batch normalization transformation at test time.

[0076] Fig.14 The batch normalization backward pass is shown.

[0077] Fig.15 Depicts the use of batch normalization layers with convolutional or densely connected layers.

[0078] Fig.16 An implementation of 1D convolution is shown.

[0079] Fig.17 Shows how global average pooling (GAP) works.

[0080] Fig.18 One implementation of a computing environment with training servers and production servers that can be used to implement the disclosed techniques is shown.

[0081] Fig.19 Depicts one implementation of the Atrous Convolutional Neural Network (ACNN) architecture, referred to herein as “SpliceNet.”

[0082] Fig. 20One implementation of a residual block that can be used by ACNN and convolutional neural networks (abbreviated as CNN) is shown.

[0083] Fig.21 Depicts another implementation of the ACNN architecture, referred to in this paper as “SpliceNet80”.

[0084] Fig. 22 Yet another implementation of the ACNN architecture is depicted, referred to herein as “SpliceNet400”.

[0085] Fig.23 Yet another implementation of the ACNN architecture is depicted, referred to herein as “SpliceNet2000”.

[0086] Fig.24 Yet another implementation of the ACNN architecture is depicted, referred to herein as “SpliceNet10000”.

[0087] Fig.25 , Fig.26 and Fig. 27 Various types of inputs processed by ACNN and CNN are shown.

[0088] Fig.28 It is shown that ACNN can be trained on at least 8 million non-splicing sites, and CNN can be trained on at least 1 million non-splicing sites.

[0089] Fig.29 A one-hot encoder is shown.

[0090] Fig.30 Depicts the training of ACNN.

[0091] Fig.31 A CNN is shown.

[0092] Fig.32 The training, validation, and testing of ACNN and CNN are shown.

[0093] Fig.33 Reference and alternative sequences are depicted.

[0094] Fig.34 Aberrant splicing detection is shown.

[0095] Fig.35 Shown is the processing pyramid of SpliceNet10000 for splice site classification.

[0096] Fig.36 Depicted is the processing pyramid of SpliceNet10000 for anomalous splicing detection.

[0097] Fig.37A , Fig.37B , Fig.37C , Fig.37D , Fig.37E , Fig.37F , Figure 37G and Fig.37H An implementation of using deep learning to predict splicing from raw sequences is shown.

[0098] Fig.38A , Fig.38B , Fig.38C , Fig.38D , Fig.38E , Fig.38F and Figure 38G Depicted is an implementation of validation of rare cryptic splicing mutations in RNA sequencing data.

[0099] Fig.39A , Fig.39B and Fig.39C One implementation of cryptic splice variants that frequently create tissue-specific alternative splicing is shown.

[0100] Fig.40A , Fig.40B , Fig.40C , Fig.40D and Fig.40E Depicts one realization in which a predicted cryptic splice variant is strongly deleterious in the human population.

[0101] Fig.41A , Fig.41B , Fig.41C , Fig.41D , Fig.41E and Fig.41F Shown is one realization of de novo cryptic splicing mutations in patients with rare genetic diseases.

[0102] Fig.42A and Fig.42B Depicted is the evaluation of various splicing prediction algorithms on LincRNA.

[0103] Fig.43A and Fig.43B Position-dependent effects of the TACTAAC branch point and the GAAGAA exon-splicing enhancer motif are shown.

[0104] Fig.44A and Fig.44B Delineated the effects of nucleosome positioning on splicing.

[0105] Fig.45 An example of calculating the effect size of a splicing-disrupting variant with complex effects is shown.

[0106] Fig.46A , Fig.46B and Fig.46C Evaluation of the SpliceNet-10k model on singletons and common variants is shown.

[0107] Fig.47A and Fig.47B Validation rates and effect sizes of splice site creation variants split by the position of the variant are depicted.

[0108] Fig.48A , Fig.48B , Fig.49C and Figure 49D depicts the evaluation of the SpliceNet-10k model on training and test chromosomes.

[0109] Fig.49A , Fig.49B and Fig.49C De novo cryptic splicing mutations in patients with rare genetic diseases were shown only from synonymous, intronic or non-translational regions.

[0110] Fig.50A and Fig.50B Shown are cryptic splicing de novo mutations in ASD and the proportion that are pathogenic DNMs.

[0111] Fig.51A , Fig.51B , Fig.51C , Fig.51D , Fig.51E , Fig.51F , Figure 51G , Fig.51H , Fig.51I and Fig.51J Depicted is RNA sequencing validation of predicted cryptic splicing de novo mutations in ASD patients.

[0112] Fig.52A and Fig.52B Validation rate and sensitivity for RNA sequencing of a model trained only on canonical transcripts are shown.

[0113] Fig.53A , Fig.53B and Fig.53C It is shown that overall modeling improves the SpliceNet-10k performance.

[0114] Figures 54A and 54B show the evaluation of SpliceNet-10k in regions of different exon densities.

[0115] Fig.55 is Table S1, which depicts one implementation of GTEx sampling for demonstrating effect size calculations and tissue-specific splicing.

[0116] Fig.56 is Table S2, which depicts one implementation of the cutoff values ​​used to evaluate the validation rate and sensitivity of different algorithms.

[0117] Fig.57 One implementation of per-gene enrichment analysis is shown.

[0118] Fig.58 One implementation of genome-wide enrichment analysis is shown.

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

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

[0121] introduce

[0122] Convolutional Neural Networks

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

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

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

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

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

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

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

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

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

[0132] Training a Convolutional Neural Network

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

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

[0135] Δw i =x i δ

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

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

[0138]

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

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

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

[0142]

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

[0144]

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

[0146]

[0147]

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

[0149]

[0150]

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

[0152]

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

[0154]

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

[0156]

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

[0158] v nm ←w nm +αδ hm a n

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

[0160]

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

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

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

[0164] v t+1 =μv-α▽wQ(z t ,w t )

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

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

[0167] Convolutional Layer

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

[0169] is defined as:

[0170]

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

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

[0173] Non-linear layers

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

[0175]

[0176]

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

[0178]

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

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

[0181]

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

[0183] Dilated Convolution

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

[0185] Subsampling Layer

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

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

[0188]

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

[0190] In average pooling, the output of the pooling neuron is the average of the input values ​​residing in the set of input neurons, as described below:

[0191]

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

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

[0194] Convolution Example

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

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

[0197] Forward pass

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

[0199]

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

[0201]

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

[0203]

[0204] Backward Propagation

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

[0206]

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

[0208]

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

[0210]

[0211]

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

[0213]

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

[0215]

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

[0217]

[0218]

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

[0220]

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

[0222]

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

[0224]

[0225]

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

[0227]

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

[0229]

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

[0231]

[0232]

[0233] Residual Connection

[0234] Fig. 9 Depicted is the addition of residual connections that reinject initial information downstream via feature maps. Residual connections involve reinjecting previous representations into the downstream data stream by adding past output tensors to later output tensors, which helps prevent information loss along the data processing stream. Residual connections address two common problems that plague any large-scale deep learning model: vanishing gradients and representational bottlenecks. In general, adding residual connections to any model with more than 10 layers is likely to be beneficial. As discussed above, residual connections involve making the output of an earlier layer available as an input to a later layer, effectively creating shortcuts in a sequential network. Instead of being concatenated to a later activation, the earlier output is added to the later activation, assuming that both activations are the same size. If they are of different sizes, a linear transformation that reshapes the earlier activation to the target shape can be used.

[0235] Residual Learning and Skip Connections

[0236] Fig.10 Depicts one implementation of residual blocks and skip connections. The main idea of ​​residual learning is that residual mappings are easier to learn than the original mappings. Residual networks stack multiple residual units to mitigate the drop in training accuracy. Residual blocks utilize special additive skip connections to combat vanishing gradients in deep neural networks. At the beginning of a residual block, the data flow is split into two streams: the first stream carries the unchanged input of the block, while the second stream applies weights and nonlinearities. At the end of the block, the two streams are merged using an element-wise sum. The main advantage of this type of construction is that it allows gradients to flow through the network more easily.

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

[0238] WaveNet

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

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

[0241]

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

[0243] Batch Normalization

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

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

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

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

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

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

[0250] Forward pass

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

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

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

[0254] Backward pass

[0255] Since normalization is a differentiable operation, we can Fig.14 The backward pass is computed as depicted in .

[0256] 1D Convolution

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

[0258] Global Average Pooling

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

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

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

[0262] the term

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[0294] Sequencing process

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[0316] d.VEP cache file

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

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

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

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

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

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

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

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

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

[0326] Variant Identification Application

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

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

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

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

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

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

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

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

[0335] Deep Learning in Genomics

[0336] Genetic variation can help explain many diseases. Each person has a unique genetic code, and there are many genetic variants within a group of people. Most harmful genetic variants have been purged from the genome by natural selection. It is important to determine which genetic variants may be pathogenic or harmful. This will help researchers focus on possible disease-causing gene variants and speed up the diagnosis and cure of many diseases.

[0337] Modeling the properties and functional effects (e.g., pathogenicity) of variants is an important but challenging task in the field of genomics. Despite the rapid development of functional genomic sequencing technologies, interpreting the functional consequences of variants remains a huge challenge due to the complexity of cell type-specific transcriptional regulatory systems.

[0338] Regarding pathogenicity classifiers, a deep neural network is an artificial neural network that utilizes multiple nonlinear and complex transformation layers to continuously model high-level features. Deep neural networks provide feedback via backpropagation that carries the difference between the observed and predicted outputs to adjust parameters. Deep neural networks have continued to advance with the availability of large training datasets, the power of parallel and distributed computing, and sophisticated training algorithms. Deep neural networks have facilitated significant advances in many fields such as computer vision, speech recognition, and natural language processing.

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

[0340] The goal of training a deep neural network is to optimize the weight parameters of each layer so that simpler features are gradually merged into complex features in order to learn the most appropriate hierarchical representation from the data. A single cycle of the optimization process is organized as follows. First, given a training dataset, the forward method sequentially computes the output in each layer and propagates the function signal forward through the network. In the final output layer, the target loss function measures the error between the inferred output and the given label. To minimize the training error, the backward method uses the chain rule to backpropagate the error signal and calculate the gradient with respect to all weights in the entire neural network. Finally, the weight parameters are updated using an optimization algorithm based on stochastic gradient descent. Batch gradient descent performs parameter updates on each complete dataset, while stochastic gradient descent provides a stochastic approximation by performing updates on each small set of data examples. Stochastic gradient descent has spawned several optimization algorithms. For example, the Adagrad and Adam training algorithms perform stochastic gradient descent while adaptively modifying the learning rate based on the update frequency and gradient moment of each parameter.

[0341] Another core element in deep neural network training is regularization, which refers to strategies designed to avoid overfitting and thus achieve good generalization performance. For example, weight decay adds a penalty term to the objective loss function, causing the weight parameter to converge to a smaller absolute value. During training, dropout randomly removes hidden units from a neural network, which can be viewed as a collection of possible subnetworks. To enhance the capabilities of dropout, a new activation function, maxout, and a dropout variant of a recurrent neural network called rnnDrop were proposed. In addition, batch normalization provides a new regularization method by normalizing the scalar features of each activation in a mini-batch and learning each mean and variance as parameters.

[0342] Considering that sequence data is multi-dimensional and high-dimensional, deep neural networks hold great promise in bioinformatics research due to their wide applicability and enhanced predictive power. Convolutional neural networks have been applied to solve sequence-based problems in genomics, such as motif discovery, pathogenic variant identification, and gene expression inference. Convolutional neural networks use a weight-sharing strategy that is particularly useful for studying DNA because it can capture sequence motifs, which are short and recurring local patterns in DNA that are believed to have important biological functions. A hallmark of convolutional neural networks is the use of convolutional filters. Unlike traditional classification methods that are based on carefully designed and handcrafted features, convolutional filters perform adaptive learning of features, similar to the process of mapping raw input data to information representations of knowledge. In this sense, convolutional filters act as a series of motif scanners, as a set of such filters can identify relevant patterns in the input and update themselves during the training process. Recurrent neural networks can capture long-range dependencies in sequence data of varying lengths, such as protein or DNA sequences.

[0343] Therefore, a powerful computational model to predict pathogenicity of variants could bring tremendous benefits to both basic science and translational research.

[0344] Currently, only 25–30% of patients with rare diseases receive a molecular diagnosis by testing protein-coding sequences, suggesting that the remaining diagnostic yield may lie in the 99% of the noncoding genome. Here, we describe a novel deep learning network that accurately predicts splice junctions from arbitrary pre-mRNA transcript sequences, enabling precise prediction of the splicing-altering effects of noncoding variants. Synonymous and intragenic variants with predicted splicing-altering outcomes were validated with high rates on RNA sequencing (RNA-seq) and were highly deleterious in the population. De novo variants with predicted splicing-altering outcomes were significantly enriched in patients with autism and intellectual disability compared with healthy controls and validated in 21 of 28 patients based on RNA-sequencing data. We estimate that 9–11% of pathogenic variants in patients with rare genetic diseases are caused by this class of previously underappreciated disease variants.

[0345] Exome sequencing has transformed clinical diagnosis for patients and families with rare genetic diseases, greatly reducing the time and cost of the diagnostic process when used as a first-line test (Monroe et al., 2016; Stark et al., 2016; Tan et al., 2017). However, in rare genetic disease cohorts, the diagnostic yield of exome sequencing is approximately 25-30% at best, and even with combined exome and microarray testing, the majority of patients remain undiagnosed (Lee et al., 2014; Trujillano et al., 2017; Yang et al., 2014). Noncoding regions play an important role in gene regulation and account for 90% of causal disease loci identified in unbiased genome-wide association studies of complex human diseases (Ernst et al., 2011; Farh et al., 2015; Maurano et al., 2012), suggesting that exome noncoding variants may also be an important causal mutation burden in rare genetic diseases. In fact, despite being located outside the essential GT and AG splicing dinucleotides, exogenous noncoding variants that disrupt the normal pattern of mRNA splicing (often referred to as cryptic splice variants) have been considered to play an important role in rare genetic diseases (Cooper et al., 2009; Padgett, 2012; Scotti and Swanson, 2016; Wang and Cooper, 2007). However, in clinical practice, cryptic splicing mutations are often overlooked due to our incomplete understanding of the splicing code and the resulting difficulty in accurately identifying splicing-altering variants outside the essential GT and AG dinucleotides (Wang and Burge, 2008).

[0346] Recently, RNA sequencing has emerged as a promising assay for detecting splicing abnormalities in Mendelian diseases (Cummings et al., 2017; Kremer et al., 2017), but its utility in clinical settings has so far remained limited to a small number of cases where the relevant cell type is known and biopsy is available. High-throughput screens for potential splicing-altering variants (Soemedi et al., 2017) have expanded the characterization of splicing variation but are less practical for evaluating random de novo mutations in genetic disease due to the very large genomic space in which splicing-altering mutations may arise. General prediction of splicing from arbitrary pre-mRNA sequences would make it possible to accurately predict the splicing-altering consequences of noncoding variants, significantly improving diagnostics for patients with genetic diseases. To date, despite progress in specific applications such as modeling sequence features of core splicing motifs ( Yeo and Burge, 2004 ), characterizing exonic splicing enhancers and silencers ( Fairbrother et al., 2002 ; Wang et al., 2004 ), and predicting inclusion of cassette exons ( Barash et al., 2010 ; Jha et al., 2017 ; Xiong et al., 2015 ), a general predictive model for splicing from native sequences that is close to spliceosome specific remains elusive.

[0347] Splicing of long pre-mRNAs into mature transcripts is remarkable for its precision and the clinical severity of splicing-altering mutations, but the cellular mechanisms that determine their specificity remain incompletely understood. Here, we trained a deep learning network with near in silico spliceosomal accuracy to identify exon-intron boundaries from pre-mRNA sequences with 95% accuracy, with over 80% validation of functional cryptic splicing mutations predicted by RNA sequencing. Noncoding variants predicted to alter splicing are highly deleterious in humans, with 80% of de novo cryptic splicing mutations experiencing negative selection, similar to the effects of other types of protein-truncating variants. De novo cryptic splicing mutations in patients with autism and intellectual disability are consistent with genes that are recurrently mutated by protein-truncating mutations, allowing for the discovery of additional candidate disease genes. We estimate that up to 24% of exogenous causal mutations in patients with rare genetic diseases are due to this class of previously underappreciated disease variants, underscoring the need for improved interpretation of the noncoding genome for the 99% of clinical sequencing applications.

[0348] Clinical exome sequencing has revolutionized the diagnostic approach for patients and families with rare genetic diseases, greatly reducing the time and cost of the diagnostic process when used as a first-line test. However, exome sequencing has been reported to have a diagnostic rate of 25%-30% in multiple large cohorts of patients with rare diseases and their parents, and even with combined exome and microarray testing, most patients remain undiagnosed. The noncoding genome is highly active in gene regulation, and noncoding variants account for up to 90% of the GWAS sample counts for common diseases, suggesting that rare variants in the noncoding genome may also account for a significant proportion of causal mutations in exogenous diseases such as rare genetic diseases and neoplastic diseases. However, the difficulty in interpreting variants in the noncoding genome means that, outside of large structural variants, the noncoding genome currently offers little additional diagnostic value relative to the rare exogenous variants that have the greatest impact on clinical management.

[0349] The role of splicing-altering mutations beyond the canonical GT and AG splicing dinucleotides has long been appreciated. In fact, these cryptic splice variants are the most common mutations in several rare genetic diseases, such as glycogen storage disease XI (Pompe disease) and erythrocytic protoporphyria. The extended splicing motifs at the 5' and 3' ends of introns are highly degenerate, and equally good motifs occur frequently in the genome, making it difficult for existing methods to accurately predict which noncoding variants may lead to cryptic splicing.

[0350] To better understand how the spliceosome achieves its specificity, we trained a deep learning neural network to predict for each nucleotide in a pre-mRNA transcript whether it is a splice acceptor, splice donor, or neither, using only the transcript sequence as input. Fig.37A Using canonical transcripts on even chromosomes as the training set and transcripts on odd chromosomes for testing (excluding paralogous genes), the deep learning network called exon-intron boundaries with 95% accuracy, even transcripts larger than 100KB, such as CFTR, which are typically reconstructed to nucleotide precision ( Fig.37B ).

[0351] We next sought to understand the specificity determinants used by the network to identify exon–intron boundaries with such remarkable accuracy. Unlike previous classifiers that operate on statistical or hand-crafted features, deep learning learns features directly from sequence in a hierarchical manner, enabling the conferring of additional specificity from long sequence context. Indeed, we found that the accuracy of the network was highly dependent on the length of the sequence context flanking the nucleotides to be predicted provided as input to the network (Table 1), and that performance was only slightly better than existing statistical methods when we trained the deep learning model using only 40-nt sequences. This suggests that deep learning adds little over existing statistical methods for identifying single 9–23 nt splice motifs, but rather that broader sequence context is key to distinguishing functional splice sites from nonfunctional sites with motifs of equal strength. Requiring the network to predict exons in which sequence perturbations occur showed that disrupting the donor motif often also results in abrogation of the acceptor signal ( Fig.37C ), as is frequently observed for exon skipping events in vivo , indicating that a considerable degree of specificity can be achieved simply by requiring pairing of strong acceptor and donor motifs at acceptable distances.

[0352] While a large body of evidence suggests that experimental perturbations of exon length have strong effects on exon gain and exon skipping, this does not explain why the accuracy of the deep learning network continues to increase beyond 1000-nt context. To better distinguish between local splicing motif-driven specificity and long-range specificity determinants, we trained a local network that received only 100-nt context input. Using the local network to score known contexts, we found that both exons and introns have an optimal length (no more than 115 nt for exons and no more than 1000 nt for introns) at which motif strength is minimal ( Fig.37D ). This relationship does not exist in the 10,000-nt deep learning network ( Fig.37E ), indicating that the variation in intron and exon lengths is fully accounted for in the wide-context deep learning network. Notably, the deep learning model that never gave the wide-context intron and exon boundaries was able to deduce these distances by inferring the locations of introns and exons from the sequence alone.

[0353] Systematic searches of hexamer space also revealed that the deep learning network exploited motifs in exon-intron definitions, particularly the branch point motif TACTAAC from positions −34 to −14, the well-characterized exonic splicing enhancer GAAGAA near the end of the exon, and the polyuridylic acid motif ( Fig.21 , 22 , 23 and 24).

[0354] We extended deep learning networks to assess genetic variants that alter splicing function by predicting exon-intron boundaries on reference transcript sequences as well as alternative transcript sequences containing the variants, and looking for any changes in exon-intron boundaries. Recent aggregated exome data from 60,706 individuals allowed us to assess the impact of negative selection on variants predicted to alter splicing function by examining their distribution in the allele spectrum. We found that predicted cryptic splice variants were strongly under negative selection ( Fig.38A ), as evidenced by their relative depletion at high allele frequencies compared to expected counts, and their magnitude of depletion is comparable to that of AG or GT splice-disrupting variants and stop-gain variants. The impact of negative selection is larger when considering cryptic splice variants that may cause frameshifts over those that cause in-frame changes ( Fig.38B Based on the depletion of frameshift cryptic splice variants compared with other types of protein-truncating variants, we estimate that 88% of confidently predicted cryptic splicing mutations are functional.

[0355] While there is not as much aggregated genome-wide data available as exome data, limiting the ability to detect the effects of natural selection in deep intronic regions, we were also able to calculate observed versus expected counts of cryptic splicing mutations in regions far from exons. Overall, we observed a 60% loss of cryptic splicing mutations in regions >50 nt from exon-intron boundaries ( Fig.38C ). The attenuated signal may be a combination of smaller sample sizes, greater difficulty in predicting deep intronic variants with whole genome data compared to exomes.

[0356] We can also use the observed versus expected number of cryptic splice variants to estimate the number of cryptic splice variants to be selected and how this compares to other types of protein-truncating variants. Because cryptic splice variants may only partially abolish splicing function, we also evaluated the number of observed versus expected cryptic splice variants at a more relaxed threshold and estimated that there are approximately 3 times more deleterious rare cryptic splice variants than rare AG or GT splice-disrupting variants in the ExAC dataset ( Fig.38D Each individual carries a maximum of 20 rare cryptic splicing mutations, which is approximately equal to the number of protein-truncating mutations ( Fig.38E ), although not all of these variants completely abolish splicing function.

[0357] The latest release of GTEx data includes 148 individuals with both whole-genome sequencing and RNA sequencing from multiple tissue sites, allowing us to look for the impact of rare cryptic splicing variants directly in RNA sequencing data. To estimate the situation encountered in rare disease sequencing, we considered only rare variants (singleton in the GTEx cohort, allele frequency <1% in 1000 genomes) and paired them with splicing events unique to the individual with that variant. Although differences in gene expression and tissue expression and the complexity of splicing abnormalities make it difficult to assess the sensitivity and specificity of deep learning predictions, we found that at stringent specificity thresholds, more than 90% of rare cryptic splicing mutations were validated on RNA sequencing ( Fig.39A A large number of aberrant splicing events that emerged in RNA sequencing appeared to be associated with variants predicted by the deep learning classifier to have limited effect, suggesting that they only partially affect splicing function. At these more sensitive thresholds, approximately 75% of novel splice sites were predicted to result in aberrations in splicing function ( Fig.38B ).

[0358] The deep learning network successfully predicted cryptic splice variants that were highly deleterious on population sequencing data and had a high validation rate on RNA sequencing, suggesting that this approach could be used to identify additional diagnoses in rare disease sequencing studies. To test this hypothesis, we examined de novo variants in exome sequencing studies of autism and neurodevelopmental disorders and demonstrated that cryptic splicing mutations were significantly enriched in affected individuals compared to their healthy siblings ( Fig.40A ). In addition, the enrichment of cryptic splicing mutations was slightly lower than that of protein-truncating mutations, indicating that approximately 90% of our predicted cryptic splicing mutations are functional. Based on these values, approximately no more than 20% of pathogenic protein-truncating variants can be attributed to cryptic splicing mutations in exons and nucleotides immediately adjacent to exons ( Fig.40B ). Extrapolating this number to genome-wide studies that are able to interrogate entire intronic sequences, we estimate that 24% of causal mutations in rare genetic diseases are due to cryptic splicing mutations.

[0359] We estimated the probability of the emergence of a de novo cryptic splicing mutation for each individual gene, allowing us to estimate the enrichment of cryptic splicing mutations in candidate disease genes compared to chance. De novo cryptic splicing mutations are highly enriched in genes that have been previously disrupted by protein-truncating variants ( Fig.40C ), suggesting that they cause disease primarily through haploinsufficiency rather than other modes of action. Adding predicted cryptic splicing mutations to the list of protein-truncating variants allowed us to identify 3 additional autism disease genes and 11 additional intellectual disability disease genes compared to using only protein-truncating variants ( Fig.40D ).

[0360] To assess the feasibility of validating cryptic splicing mutations in patients for whom possible disease tissue is unavailable (in this case, brain), we performed deep RNA sequencing of 37 individuals with predicted de novo cryptic splicing mutations from Simon's simplex collection and looked for aberrant splicing events present in this individual but absent in all other individuals in this experiment and in the 149 individuals in the GTEx cohort. We found that NNs in the 37 patients showed unique aberrant splicing on RNA sequencing ( Fig.40E ), which could be explained by predicted cryptic splice variants.

[0361] In summary, we validated a deep learning model that accurately predicts cryptic splice variants with sufficient precision to identify causal pathogenic variants in rare genetic diseases. We estimate that a large fraction of rare disease diagnoses caused by cryptic splicing are missed by considering only protein-coding regions and highlight the need to develop methods to account for the effects of penetrant rare variants in the noncoding genome.

[0362] result

[0363] Accurately predict splicing from raw sequences using deep learning

[0364] Using only the genomic sequence of the pre-mRNA transcript as input, we constructed a deep residual neural network (He et al., 2016a) to predict whether each position in the pre-mRNA transcript is a splice donor, a splice acceptor, or neither ( Fig.37A and Fig.21 , 22 , 23, and 24). Because splice donors and acceptors may 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 identify sequence determinants across very large genomic distances. In contrast to previous approaches that only considered short nucleotide windows near exon-intron boundaries (Yeo and Burge, 2004), or relied on artificially designed features (Xing et al., 2015) or experimental data such as expression or splicing factor binding (Jha et al., 2017), our neural network evaluates 10,000 nucleotides of flanking context sequence to predict the splicing function of each position in the pre-mRNA transcript, thereby learning splicing determinants directly from the raw sequence.

[0365] We used GENCODE-annotated pre-mRNA transcript sequences on a subset of human chromosomes (Harrow et al., 2012) to train the parameters of the neural network and used transcripts on the remaining chromosomes excluding paralogous genes to test the network's predictions. For the pre-mRNA transcripts in the test dataset, the network predicted splice sites with 95% top-k accuracy, which is the proportion of correctly predicted splice sites at a threshold where the number of predicted sites equals the number of actual splice sites in the test dataset (Boyd et al., 2012; Yeo and Burge, 2004). Even genes larger than 100 kb, such as CFTR, were often perfectly reconstructed to nucleotide accuracy ( Fig.37B ). To confirm that the network was not solely dependent on exonic sequence bias, we also tested the network with long noncoding RNAs. Although incomplete annotation of noncoding transcripts may reduce our accuracy, the network predicted known splice junctions in lincRNAs with a top-k accuracy of 84% ( Fig.42A and 42B ), indicating that it can approximate the spliceosomal behavior of arbitrary sequences without protein-coding selective pressure.

[0366] For each GENCODE-annotated exon in the test dataset (excluding the first and last exons of each gene), we also examined whether the network's prediction score correlated with the proportion of reads supporting exon addition and exon skipping based on RNA sequencing data from the Gene and Tissue Expression Atlas (GTEx Consortium et al., 2015) ( Fig.37C Exons constitutively spliced ​​in GTEx tissues had prediction scores close to 1 or 0, respectively, whereas exons that underwent a greater degree of alternative splicing (10%-90% exon gain across samples on average) tended toward intermediate values ​​(Pearson correlation coefficient = 0.78, P ≈ 0).

[0367] Next, we sought to understand the sequence determinants that the network exploits to achieve its extraordinary accuracy. We performed systematic in silico substitutions of every nucleotide near annotated exons and measured the effect on the network's prediction score at adjacent splice sites ( Fig.37E ). We found that perturbing the sequence of the splice donor motif often resulted in the loss of the network-predicted upstream splice acceptor site, as observed for exon skipping events in vivo, suggesting that exon definition between a paired upstream acceptor motif and a downstream donor motif set at an optimal distance confers a considerable degree of specificity (Berget, 1995). Other motifs that contribute to splicing signals include well-characterized binding motifs for the SR-protein family and branch points ( Fig.43A and 43B)(Fairbrother et al., 2002; Reed and Maniatis, 1988). The effects of these motifs are highly dependent on their location within exons, suggesting that their roles include specifying the precise positioning of intron-exon boundaries by discriminating between competing acceptor and donor sites.

[0368] Training the network with different input sequence contexts significantly affects the accuracy of splicing predictions ( Fig.37E ), indicating that long-range sequence determinants up to 10,000 nt from the splice site are critical for identifying functional splice sites from a large number of non-functional sites with near-optimal motifs. To examine long-range and short-range specific determinants, we compared the scores assigned to annotated splice sites by a model trained on an 80 nt sequence context (SpliceNet-80nt) with the full model trained on a 10,000 nt context (SpliceNet-10k). The network trained on an 80 nt sequence context assigned lower scores to splice sites connected to exons or introns of typical length (150 nt for exons and up to 1000 nt for introns) ( Fig.37F ), which is consistent with previous observations that splice motifs tend to be weaker in splice sites with unusually long or short exons and introns compared to such sites (Amit et al., 2012; Gelfman et al., 2012; Li et al., 2015). In contrast, the network trained on the 10,000 nt sequence context showed a preference for introns and exons of average length despite their weaker splice motifs, as it could account for the long-range specificity conferred by exon or intron length. Skipping of weaker motifs in long non-stop introns is consistent with faster RNA polymerase II elongation observed experimentally in the absence of exonic pauses, which may give the spliceosome less time to recognize suboptimal motifs (Close et al., 2012; Jonkers et al., 2014; Veloso et al., 2014). Our findings suggest that on average splice junctions have well-characterized long-range sequence determinants that have considerable specificity, explaining the high degree of tolerance to sequence degradation at most splice motifs.

[0369] As 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 potentially be exploited by the network in cases where chromatin state is predictable from the original 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 for splice site prediction, we ran a pair of optimal acceptor and donor motifs separated by 150 nt (approximately the size of an average exon) through the genome and asked the network to predict whether the motif pair would result in exon gain at that locus ( Figure 37G We found that positions predicted to favor exon gain were associated with high nucleosome occupancy even in intergenic regions (Spearman correlation coefficient = 0.36, P ≈ 0), and this effect persisted after controlling for GC content ( Fig.44A ). These results suggest that the network has implicitly learned to predict the position of nucleosomes from the original sequence and use it as a specific determinant for exon definition. Similar to exons and introns of average length, exons located on nucleosomes have weaker local splicing motifs ( Fig.44B ), consistent with a greater tolerance for degenerate motifs in the presence of compensatory factors ( Spies et al., 2009 ).

[0370] Although several studies have reported correlations between exon and nucleosome occupancy, a causal role for nucleosome positioning in exon definition has not been clearly established. Using data from 149 individuals with RNA sequencing and whole-genome sequencing from the Genotype-Tissue Expression (GTEx Consortium et al., 2015) cohort, we discovered novel exons that were exclusive to a single individual and corresponded to private splice site-creating mutations. These private exon creation events were significantly associated with existing nucleosome positioning in K562 and GM12878 cells (P = 0.006 by permutation test, Fig.37H ), even though these cell lines most likely lack the corresponding private gene mutation. Our results suggest that genetic variants are more likely to trigger the generation of new exons if the resulting new exons overlap regions occupied by existing nucleosomes, supporting a causal role for nucleosome positioning in promoting exon definition.

[0371] Validation of predicted cryptic splicing mutations in RNA sequencing data

[0372] We extended the deep learning network to assess genetic variants with altered splicing functions by predicting exon-intron boundaries for a reference pre-mRNA transcript sequence and an alternative transcript sequence containing the variant and taking the difference between the scores (Δ score). Fig.38A Importantly, the network was trained only on reference transcript sequences and splice junction annotations and never saw variant data during training, making prediction of variant effects a challenging test of the network’s ability to accurately model splicing sequence determinants.

[0373] We investigated the impact of cryptic splice variants in RNA sequencing data from the GTEx cohort (GTEx Consortium et al., 2015), which includes 149 individuals with both whole-genome sequencing and RNA sequencing from multiple tissues. To mimic the situation encountered in rare disease sequencing, we first focused on rare private mutations (occurring in only one individual in the GTEx cohort). We found that private mutations predicted by the neural network to have functional consequences were strongly enriched at novel private splice junctions and at skipped exon boundaries of private exon skipping events ( Fig.38B ), indicating that a large proportion of these predictions are functional.

[0374] To quantify the impact of splice site generation variants on the relative production of normal and aberrant splice isoforms, we measured the number of reads supporting a novel splicing event as a fraction of the total number of reads covering that site ( Fig.38C )(Cummings et al., 2017). For splice site disrupting variants, we observed low baseline rates of exon skipping for many exons, and the effect of the variant was to increase the fraction of exon skipped reads. Therefore, we took the larger of the two effects, calculating the decrease in the fraction of reads spliced ​​at the disrupted splice site and the increase in the fraction of reads skipped exons ( Fig.45 and STAR method).

[0375] Confidently predicted cryptic splice variants (Δ score ≥ 0.5) were validated by RNA sequencing at a rate three-quarters of the disruption rate of the basal GT and AG splicing variants ( Fig.38D ). Both the validation rate and effect size of cryptic splice variants closely track their Δ scores ( Fig.38D and 38E ), indicating that the prediction score of the model can well reflect the splicing-altering potential of the variant. Validated variants, especially those with lower scores (Δ score < 0.5), are often not fully penetrant and lead to alternative splicing that produces a mixture of abnormal and normal transcripts in RNA sequencing data ( Fig.38EOur estimates of validation rates and effect sizes are conservative and likely underestimate the true value because they do not account for changes in splicing isoforms and nonsense-mediated decay, which can preferentially degrade aberrantly spliced ​​transcripts because they often introduce premature stop codons ( Fig.38C and Fig.45 This is demonstrated by the average effect size of variants disrupting the essential GT and AG splice dinucleotides, which is less than 50% of the expected value for a fully penetrant heterozygous variant.

[0376] For cryptic splice variants that generate an aberrant splice isoform in at least three out of ten observed mRNA transcript copies, the network had a sensitivity of 71% when the variant was a nearby exon and 41% when the variant was in deep intronic sequences (Δ score ≥ 0.5, Fig.38F ). These findings suggest that prediction of deep intronic variants is more challenging, likely because deep intronic regions contain fewer specificity determinants that are selected to be present near exons.

[0377] To benchmark the performance of our network against existing methods, we selected three popular classifiers cited in the literature for diagnosis of rare genetic diseases: GeneSplicer (Pertea et al., 2001), MaxEntScan (Yeo and Burge, 2004), and NNSplice (Reese et al., 1997), and plotted the RNA sequencing validation rate and sensitivity at different thresholds ( Figure 38G As has been the experience of others in the field (Cummings et al., 2017), we found that existing classifiers lacked specificity given the very large number of noncoding variants that could affect splicing genome-wide, presumably because they focus on local motifs and largely fail to account for long-range specificity determinants.

[0378] Given the large gap in performance with existing methods, we performed additional controls to rule out the possibility that the results in the RNA sequencing data were confounded by overfitting. First, we repeated the validation and sensitivity analyses separately for private variants and variants present in more than one individual in the GTEx cohort ( Fig.46A , 46B and 46C). As neither the splicing machinery nor the deep learning model has access to allele frequency information, validating that the network performs similarly across allele frequency spectra is an important control. We found no significant difference in the validation rates of private and common cryptic splice variants in RNA sequencing at the same Δ score threshold (P > 0.05, Fisher’s exact test), indicating that the network’s predictions are robust to allele frequency.

[0379] Second, to validate the predictive power of the model across 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 those that occur in more distal regions. We found that cryptic splice variants were roughly equally distributed among the three groups, and at the same Δ score threshold, there were no significant differences in validation rates or effect sizes between the groups (P>0.3χ 2 The consistency test and the Mann-Whitney test with P>0.3, Fig.47A and 47B ).

[0380] Third, we performed RNA sequencing validation and sensitivity analysis on the variants on the chromosomes used for training and on the variants on the remaining chromosomes, respectively ( Fig.48A and 48B ). Although the network was trained only on reference genome sequences and splicing annotations and was not exposed to variant data during training, we wanted to rule out the possibility that variant predictions were biased due to the fact that the network had seen reference sequences in the training chromosomes. We found that the network performed equally well on variants from the training and test chromosomes, with no significant differences in validation rate or sensitivity (P>0.05, Fisher's exact test), suggesting that the network's variant predictions are unlikely to be explained by overfitting to the training sequences.

[0381] From the results of our model and other splicing prediction algorithms (compare Fig.37E and Figure 38G ) As can be seen, predicting cryptic splice variants is much more difficult than predicting annotated splice junctions. One important reason is that the basic distribution of exon addition rates differs between the two types of analysis. The vast majority of exons with GENCODE annotations have strong specificity determinants, resulting in constitutive splicing and prediction scores close to 1 ( Fig.37C In contrast, most cryptic splice variants are only partially expressed ( Fig.38D and 38E ), have low or moderate prediction scores, and often result in alternative splicing that produces a mixture of normal and aberrant transcripts. This makes the problem of predicting the impact of cryptic splice variants inherently more difficult than identifying annotated splice sites. Other factors, such as nonsense-mediated decay, unexplained isoform changes, and limitations of RNA sequencing analysis, further reduce the validation rate of RNA sequencing ( Fig.38C and Fig.45 ).

[0382] Tissue-specific alternative splicing is often caused by weak cryptic splice variantsAlternative splicing is a major mode of gene regulation that increases transcript diversity across tissues and developmental stages, and its dysregulation has been implicated in disease processes (Blencowe, 2006; Irimia et al., 2014; Keren et al., 2010; Licatalosi and Darnell, 2006; Wang et al., 2008). Unexpectedly, we found that the relative usage of novel splice junctions generated by cryptic splicing mutations can vary widely across tissues ( Fig.39A Furthermore, variants leading to tissue-specific differences in splicing can be reproducibly found between multiple individuals ( Fig.39B ), suggesting that tissue-specific biology may underlie these differences rather than random effects. We found that 35% of cryptic splice variants with weak and intermediate prediction scores (Δ scores 0.35–0.8) exhibited significant differences in the normal and aberrant transcript scores generated across tissues (χ 2 test, Bonferroni-corrected P < 0.01, Fig.39C This is in contrast to variants with higher prediction scores (Δ score > 0.8), which were significantly less likely to have a tissue-specific effect (P = 0.015). Our findings are consistent with previous observations that variably spliced ​​exons tend to have intermediate prediction scores, compared to constitutively spliced ​​exons that have scores close to 1 or 0, respectively ( Fig.37C ).

[0383] These results support a model in which tissue-specific factors, such as chromatin context and binding of RNA-binding proteins, may generate competition between two splice sites with close affinity (Gelfman et al., 2013; Luco et al., 2010; Shukla et al., 2011; Ule et al., 2003). Strong cryptic splice variants may completely shift splicing from the normal isoform to the abnormal isoform regardless of the epigenetic context, while weaker variants shift splice site selection closer to the decision boundary, resulting in the use of alternative splice sites in different tissue types and cellular contexts. This highlights the unexpected role that cryptic splicing mutations play in generating novel alternative splicing diversity, as natural selection will have the opportunity to retain mutations that generate useful tissue-specific alternative splicing.

[0384] Predicted cryptic splice variants are strongly deleterious in human populations

[0385] Although predicted cryptic splice variants have a high validation rate in RNA sequencing, in many cases the effects are not fully penetrant and a mixture of normal and aberrant splice isoforms is generated, raising the possibility that a small fraction of cryptic splice-altering variants are functionally insignificant. To explore the signature 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.

[0386] To measure the extent to which negative selection is acting on predicted splice-altering variants, we counted the number of predicted splice-altering variants found at common allele frequency (≥0.1% in the population) and compared this to the number of predicted splice-altering variants found at monoallelic frequency in ExAC (i.e., 1 in 60,706 individuals). Singleton variants represent recently created mutations that have been minimally filtered by purifying selection due to recent exponential population growth (Tennessen et al., 2012). In contrast, common variants represent a subset of neutral mutations that have been filtered by purifying selection. Therefore, the depletion of predicted splice-altering variants in the common allele spectrum relative to singleton variants provides an estimate of the fraction of predicted splice-altering variants that are deleterious and therefore functional. To avoid confounding effects on protein-coding sequences, we restricted our analysis to synonymous variants and intronic variants outside the essential GT or AG dinucleotides and excluded missense mutations that are also predicted to have splice-altering effects.

[0387] At common allele frequencies, confidently predicted cryptic splice variants (Δ score ≥ 0.8) are under strong negative selection, as evidenced by relative attrition relative to expectation ( Fig.40A At this threshold, most variants are expected to be close to fully penetrant in RNA sequencing data ( Fig.38D ), predicted synonymous and intronic cryptic splicing variants were reduced in common allele frequency by 78%, which was comparable to the 82% depletion of frameshift, stop gain, and essential GT or AG splicing disruption variants ( Fig.40B The impact of negative selection is greater when considering cryptic splice variants that would result in a frameshift rather than cryptic splice variants that cause an in-frame change ( Fig.40C The depletion of cryptic splice variants with frameshift consequences was almost identical to that of other classes of protein-truncating variants, indicating that the vast majority of confidently predicted cryptic splice mutations in near-intronic regions (≤50 nt from known exon-intron boundaries) are functional and have severe deleterious effects in humans.

[0388] To extend this analysis to deep intronic regions >50 nt from known exon-intron boundaries, we used pooled whole-genome sequencing data from 15,496 individuals from the Genome Aggregation Database (gnomAD) cohort (Lek et al., 2016) to calculate observed and expected cryptic splicing mutation counts at common allele frequencies. Overall, we observed a 56% depletion of common cryptic splicing mutations (Δ score ≥ 0.8) at >50 nt from exon-intron boundaries. Fig.40D ), which is consistent with the greater difficulty in predicting the effects of deep intronic variants, as we observed in RNA-sequencing data.

[0389] Next, we estimated the likelihood that cryptic splicing mutations cause exogenous diseases relative to other types of protein-coding variation by measuring the number of rare cryptic splicing mutations per individual in the gnomAD cohort. The fraction of cryptic splicing mutations predicted to be under negative selection ( Fig.40A ), the average person carries up to 5 rare functional cryptic splicing mutations (allele frequency < 0.1%) and up to 11 rare protein truncating mutations ( Fig.40E ). The number of cryptic splice variants was approximately twice that of the essential GT or AG splice disrupting variants. We caution that a large proportion of these cryptic splice variants may not completely lose gene function, either because they produce in-frame changes or because they do not completely switch splicing to an aberrant isoform.

[0390] De novo cryptic splicing mutations are a major cause of rare genetic diseases

[0391] Large-scale sequencing studies of patients with autism spectrum disorder and severe intellectual disability have revealed that de novo protein-coding mutations (missense, nonsense, frameshift, and essential splicing dinucleotides) can disrupt genes in 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 variants that act through altered splicing, we applied neural networks to predict the effects of de novo variants 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 populations.

[0392] Compared with the healthy control group (Δ score ≥ 0.1, Fig.41A , Fig.43A and Fig.43B ), de novo mutations that disrupt splicing are predicted to be 1.51-fold enriched in intellectual disability (P = 0.000416) and 1.31-fold enriched in autism spectrum disorder (P = 0.0203). When only synonymous and intronic mutations are considered ( Fig.49A , Fig.49B and Fig.49C ), splicing-disrupting mutations were also significantly enriched compared with controls, ruling out the possibility that the enrichment is explained solely by mutations with dual protein-coding and splicing effects. Based on the excess of de novo mutations in affected and unaffected individuals, and adjusting for the expected proportion of mutations in regions lacking sequencing coverage or variant ascertainment in each study, cryptic splicing mutations are estimated to account for approximately 11% of pathogenic mutations in autism spectrum disorder and 9% in intellectual disability ( Fig.41B In affected individuals, the majority of predicted de novo cryptic splicing mutations had a Δ score < 0.5 ( Fig.41C , Fig.50A and Fig.50B ), and a mixture of normal and aberrant transcripts is expected to be generated based on similarly scored variants in the GTEx RNA sequencing dataset.

[0393] To estimate the enrichment of cryptic splicing mutations in candidate disease genes compared to chance, we calculated the probability of a de novo cryptic splicing mutation for each individual gene using trinucleotide context-adjusted mutation rates (Samocha et al., 2014) (Table S4). Combining cryptic splicing mutations with protein-coding mutations in novel gene discovery yielded five additional candidate genes associated with intellectual disability and two additional genes associated with autism spectrum disorder ( Fig.41D and Fig.45 ), which fall below the discovery threshold (FDR < 0.01) when only protein-coding mutations are considered (Kosmicki et al., 2017; Sanders et al., 2015).

[0394] Experimental validation of de novo cryptic splicing mutations in autism

[0395] We obtained peripheral blood-derived lymphoblastoid cell lines (LCLs) from 36 individuals in the Simons Simplex Collection that harbor predicted de novo cryptic splice variants in genes with at least minimal LCL expression levels (De Rubeis et al., 2014; Sanders et al., 2012); each individual was the only case of autism in a direct family. As with most rare genetic diseases, relevant tissues and cell types (presumably the developing brain) are not available. We therefore performed high-depth mRNA sequencing (up to 350 million x 150 bp individual reads per sample, approximately 10 times more than in GTEx) to compensate for the weak expression of many of these transcripts in LCLs. To ensure that we were validating a representative set of predicted cryptic splice variants rather than just the top predictions, we applied relatively lenient thresholds (Δ scores > 0.1 for splicing-loss variants and Δ scores > 0.5 for splicing-gain variants; STAR Methods) and experimentally validated all de novo variants that met these criteria.

[0396] After excluding eight individuals with insufficient RNA sequencing coverage on the genes of interest, we found unique aberrant splicing events associated with predicted de novo cryptic splicing mutations in 21 of 28 patients ( Fig.41E and Fig.51A , 51B , 51C, 51D, 51E, 51F, 51G, 51H, 51I, and 51J). The other 35 individuals for whom deep LCL RNA sequencing was obtained and the 149 individuals in the GTEx cohort did not have these aberrant splicing events. Among the 21 confirmed de novo cryptic splicing mutations, we observed 9 new splice junctions, 8 exon skipping, and 4 intron retention, as well as more complex splicing aberrations ( Fig.41F , Fig.46A , Fig.46B and Fig.46C ). Seven cases did not show aberrant splicing in LCLs despite adequate transcript expression. Although a subset of them may represent false-positive predictions, some cryptic splicing mutations may result in tissue-specific alternative splicing that cannot be observed in LCLs under these experimental conditions.

[0397] Despite the limitations of the RNA sequencing assay, the high validation rate of predicted cryptic splicing mutations in ASD patients (75%) suggests that most predictions are valid. However, there was a significant difference in the number of predicted cryptic splicing mutations compared to controls (1.5-fold in DDD and 1.3-fold in ASD). Fig.41A ), the enrichment of de novo cryptic splice variants was only 38% of the effect size observed for de novo protein-truncating variants (2.5-fold for DDD and 1.7-fold for ASD) (Iossifov et al., 2014; McRae et al., 2017; De Rubeis et al., 2014). This allowed us to quantify that the clinical penetrance of functional cryptic splicing mutations is approximately 50% of that of classical forms of protein-truncating mutations (stop gain, frameshift, and essential splicing dinucleotides), as many of these mutations only partially interfere with the production of normal transcripts. In fact, some of the best characterized cryptic splicing 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. An estimate of clinical penetrance was calculated for all de novo variants that met a relative permissive threshold (Δ score ≥ 0.1), and variants with higher prediction scores were expected to have correspondingly higher penetrance.

[0398] Based on the excess of de novo mutations in cases relative to controls in the ASD and DDD cohorts, 250 cases could be explained by de novo cryptic splicing mutations compared to 909 cases explained by de novo protein-truncating variants ( Fig.41B ). Once the reduced penetrance of cryptic splicing mutations is taken into account, this is consistent with our previous estimates of the average number of rare cryptic splicing mutations (no more than 5) compared with rare protein-truncating mutations in the general population (no more than 11) ( Fig.38A The widespread distribution of cryptic splicing mutations in the genome suggests that the proportion of cases in neurodevelopmental disorders explained by cryptic splicing mutations (9-11%, Fig.41B ) may be generally applicable to other rare genetic diseases where the primary disease mechanism is loss of functional protein. To facilitate the interpretation of splicing-altering mutations, we pre-computed Δ score predictions for all possible single nucleotide substitutions genome-wide and make them available as a resource to the scientific community. We believe that this resource will advance the understanding of this previously underappreciated source of genetic variation.

[0399] Specific embodiments

[0400] Described herein are systems, methods, and articles for the manufacture of splice sites in genomic sequences (such as nucleotide sequences or amino acid sequences) using trained hole convolutional neural networks. One or more features of an embodiment may be combined with a basic embodiment. Non-mutually exclusive embodiments are taught to be combinable. One or more features of an embodiment may be combined with other embodiments. The present disclosure periodically reminds the user of these options. The omission of certain embodiments of the description repeating these options should not be considered a limitation on the combination of the preceding parts, and these descriptions are incorporated forward into the following embodiments by reference.

[0401] This section uses the terms module and phase interchangeably.

[0402] System embodiments of the present technical disclosure include one or more processors coupled to a memory. The memory is loaded with computer instructions to train a splice site detector that identifies splice sites in a genomic sequence (eg, a nucleotide sequence).

[0403] As shown in Figure 30, the system trains a dilated convolutional neural network (ACNN) with at least 50,000 donor splice site training samples, at least 50,000 acceptor splice site training samples, and at least 100,000 non-splicing site training samples. Each training sample is a target nucleotide sequence having at least one target nucleotide, each of which has at least 20 nucleotides on each flank.

[0404] ACCN uses dilated convolutions that allow for large receptive fields with few trainable parameters. The kernel of a dilated convolution is applied to an area larger than its length by skipping input values ​​in steps, also known as the dilated convolution rate or dilation factor. Dilated convolutions increase the spacing between convolution filter / kernel elements so that more closely spaced adjacent input entries (e.g., nucleotides, amino acids) are considered when performing the convolution operation. This allows long-range contextual dependencies to be incorporated into the input. Dilated convolutions retain portions of the convolution computation for reuse when processing adjacent nucleotides.

[0405] like Fig.30 As shown, the system provides a target nucleotide sequence flanked by at least 40 upstream context nucleotides and at least 40 downstream context nucleotides as an ACNN input to evaluate a training sample using the ACNN.

[0406] like Fig.30 As shown, based on this evaluation, the ACNN produces a triple score as output, making it possible for each nucleotide in the target nucleotide to be a donor splice site, an acceptor splice site, or a non-splicing site.

[0407] The disclosed system embodiments and other systems optionally include one or more of the following features. The system may also include features related to the disclosed method. For the sake of brevity, optional combinations of system features are not listed one by one. For each statutory set of basic features, features applicable to manufacturing systems, methods and articles are not repeated. The reader should understand how the features identified in this section can be easily combined with basic features in other statutory categories.

[0408] like Fig.25 , 26 As shown in Figures 27, the input may include a target nucleotide sequence having a target nucleotide, and the target nucleotide is flanked by 2500 nucleotides on both sides. In such an embodiment, the target nucleotide sequence is flanked by 5000 upstream context nucleotides and 5000 downstream context nucleotides.

[0409] The input may include a target nucleotide sequence having a target nucleotide flanked by 100 nucleotides on each side. In such an embodiment, the target nucleotide sequence is flanked by 200 upstream context nucleotides and 200 downstream context nucleotides.

[0410] The input may include a target nucleotide sequence having a target nucleotide flanked by 500 nucleotides on either side. In such an embodiment, the target nucleotide sequence is flanked by 1000 upstream context nucleotides and 1000 downstream context nucleotides.

[0411] like Fig.28 As shown, the system trains ACNN with 150,000 donor splice site training samples, at least 150,000 acceptor splice site training samples, and at least 800,000,000 non-splice site training samples.

[0412] like Fig.19 As shown, the ACNN may include multiple groups of residual blocks arranged in sequence from the lowest to the highest. Each group of residual blocks is parameterized by multiple convolution filters in the residual block, the residual block convolution window size and the residual block hole convolution rate.

[0413] like Fig.21 As shown in ,22,23 and 24, in ACNN, the dilated convolution rate evolves non-exponentially from low residual block groups to high residual block groups.

[0414] like Fig.21 As shown in ,22,23 and 24, in ACNN, the convolution window sizes are different between different groups of residual blocks.

[0415] The ACNN can be configured to evaluate an input comprising a target nucleotide sequence flanked by 40 upstream context nucleotides and 40 downstream context nucleotides. In such an embodiment, the ACNN comprises a set of four residual blocks and at least one skip connection. Each residual block has 32 convolutional filters, a convolution window size of 11, and a dilated convolution rate of 1. Fig.21 As shown, the implementation of ACNN in this article refers to SpliceNet80.

[0416] The ACNN can be configured to evaluate an input comprising a target nucleotide sequence flanked by 200 upstream context nucleotides and 200 downstream context nucleotides. In such an embodiment, the ACNN comprises at least two groups of four residual blocks and at least two skip connections. Each residual block in the first group has 32 convolutional filters, a convolution window size of 11, and a dilated convolution rate of 1. Each residual block in the second group has 32 convolutional filters, a convolution window size of 11, and a dilated convolution rate of 4. Fig. 22 As shown, the implementation of ACNN in this article refers to SpliceNet400.

[0417] The ACNN can be configured to evaluate an input comprising a target nucleotide sequence flanked by 1000 upstream context nucleotides and 1000 downstream context nucleotides. In such an embodiment, the ACNN comprises 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, a convolution window size of 11, and a dilated convolution rate of 1. Each residual block in the second group has 32 convolution filters, a convolution window size of 11, and a dilated convolution rate of 4. Each residual block in the third group has 32 convolution filters, a convolution window size of 21, and a dilated convolution rate of 19. Fig.23 As shown, the implementation of ACNN in this article refers to SpliceNet2000.

[0418] The ACNN can be configured to evaluate an input comprising a target nucleotide sequence flanked by 5,000 upstream context nucleotides and 5,000 downstream context nucleotides. In such an embodiment, the ACNN comprises at least four groups of residual blocks (four in each group) and at least four jump connections. Each residual block in the first group has 32 convolution filters, a convolution window size of 11, and a dilated convolution rate of 1. Each residual block in the second group has 32 convolution filters, a convolution window size of 11, and a dilated convolution rate of 4. Each residual block in the third group has 32 convolution filters, a convolution window size of 21, and a dilated convolution rate of 19. Each residual block in the fourth group has 32 convolution filters, a convolution window size of 41, and a dilated convolution rate of 25. Fig.24 As shown, the implementation of ACNN in this article refers to SpliceNet10000.

[0419] The triplet scores for each nucleotide in the target nucleotides may be exponentially normalized to sum to 1. In such an embodiment, the system classifies each nucleotide in the target nucleotides as a donor splice site, an acceptor splice site, or a non-splice site based on the highest score of the respective triplet scores.

[0420] like Fig.35 As shown, the dimension of ACNN input can be defined as (C u +L+C d )x4, where C is the number of upstream context nucleotides, C d is the number of downstream context nucleotides, and L is the number of nucleotides in the target nucleotide sequence. In one embodiment, the input dimension is (5000+5000+5000) x 4.

[0421] like Fig.35 As shown, the dimension of the ACNN input can be defined as Lx3. In one embodiment, the input dimension is 5000x3.

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

[0423] like Fig.32 As shown, ACNN evaluates the training samples in batches within a time period. The training samples are randomly sampled in batches. Each batch has a predetermined batch size. ACNN evaluates the training samples iteratively within multiple time periods (e.g., 1-10 time periods).

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

[0425] The system includes a one-hot encoder (such as Fig.29 As shown in Figure 2, the one-hot encoder sparsely encodes the training samples and provides the one-hot encoding as input.

[0426] ACNN is parameterized by multiple residual blocks, multiple skip connections, and multiple residual connections.

[0427] ACNN may include dimension-transforming convolutional layers that transform the spatial and feature dimensions of the aforementioned input.

[0428] like Fig. 20 As shown, each residual block may include at least one batch normalization layer, at least one rectified linear layer (ReLU for short), at least one dilated convolutional layer, and at least one residual connection. In such an embodiment, each residual block includes two batch normalization layers, two ReLU nonlinear layers, two dilated convolutional layers, and one residual connection.

[0429] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above system actions. Other embodiments may also include methods for performing the above system actions.

[0430] Other system embodiments of the disclosed technology include a trained splice site predictor running on a large number of parallel processors, the processors being coupled to a memory. The system trains a dilated convolutional neural network (ACNN) running on a large number of processors using at least 5,000 donor splice site training samples, at least 5,000 acceptor splice site training samples, and at least 100,000 non-splicing site training samples. Each training sample used for training is a nucleotide sequence, the nucleotide sequence including a target nucleotide with at least 400 nucleotides on both sides.

[0431] The system includes an input stage of an ACNN, the ACNN running on at least one of the plurality of processors and providing an input sequence of at least 801 nucleotides to evaluate the target nucleotide. Each target nucleotide is flanked by at least 400 nucleotides. In other embodiments, the system includes an input module of an ACNN, the ACNN running on at least one of the plurality of processors and providing an input sequence of at least 801 nucleotides to evaluate the target nucleotide.

[0432] The system includes an input stage of an ACNN that runs on at least one of the plurality of processors and converts the analysis of the ACNN into a classification score, making it possible that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site. In other embodiments, the system includes an input stage of an ACNN that runs on at least one of the plurality of processors and converts the analysis of the ACNN into a classification score, making it possible that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site.

[0433] Each feature discussed in this specific implementation section for the first system implementation is also applicable to this system implementation.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

[0434] ACNN can be trained with 150,000 donor splice site training samples, 150,000 acceptor splice site training samples, and 800,000,000 non-splicing site training samples. In another embodiment of the system, ACNN includes multiple groups of residual blocks arranged in a sequence from lowest to highest. In another embodiment of the system, each group of residual blocks is parameterized by multiple convolution filters in the residual block, residual block convolution window size, and residual block void convolution rate.

[0435] ACNN may include multiple groups of residual blocks arranged in a sequence from lowest to highest. Each group of residual blocks is parameterized by a number of convolutional filters in the residual block, a residual block convolution window size, and a residual block dilation convolution rate.

[0436] In ACNN, the dilated convolution rate evolves non-exponentially from low residual block groups to high residual block groups. Similarly, in ACNN, the convolution window size is different between different groups of residual blocks.

[0437] like Fig.18 As shown, ACNN can be trained through one or more training servers.

[0438] like Fig.18 As shown, the trained ACNN can be deployed to one or more production servers that receive input sequences from requesting clients. Fig.18 As shown, in such an embodiment, the production server processes the input sequence through the input stage and output stage of the ACNN to produce an output, which is sent to the client. Fig.18 As shown, in other embodiments, the production server processes the input sequence through the input stage and the output stage of the ACNN to generate an output, which is sent to the client.

[0439] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above system actions. Other embodiments may also include methods for performing the above system actions.

[0440] Other method embodiments of the disclosed technology include training a splice site predictor that identifies splice sites in a genomic sequence (eg, a nucleotide sequence).

[0441] The method includes providing a dilated convolutional neural network (ACNN) and an input sequence of at least 801 nucleotides to evaluate a target nucleotide, wherein both sides of the target nucleotide have at least 400 nucleotides.

[0442] ACNN can be trained with 5000 donor splice site training samples, 5000 acceptor splice site training samples and 100000 non-splicing site training samples. Each training sample used for training is a nucleotide sequence, and the nucleotide sequence includes a target nucleotide with at least 400 nucleotides on both sides.

[0443] The method also includes converting the ACNN analysis into a classification score, making it possible for each target nucleotide to be a donor splice site, an acceptor splice site, or a non-splicing site.

[0444] Each feature discussed in this specific implementation section for the first system implementation is also applicable to this system implementation.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

[0445] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above-mentioned system actions. Other embodiments may also include a system for performing the method, the system including a memory and one or more processors operable by the instructions stored in the memory.

[0446] Described herein are systems, methods, and articles for the manufacture of detecting abnormal splicing in genomic sequences (such as nucleotide sequences) using trained atrous convolutional neural networks. One or more features of an embodiment may be combined with a basic embodiment. Mutually non-exclusive embodiments are taught to be combinable. One or more features of an embodiment may be combined with other embodiments. The present disclosure periodically reminds the user of these options. The omission of certain embodiment descriptions repeating these options should not be considered a limitation on the combination of the preceding parts, and these descriptions are incorporated forward into the following embodiments by reference.

[0447] The system embodiment of the present technical disclosure includes one or more processors coupled to a memory. The memory is loaded with computer instructions to implement an abnormal splicing detector running on a plurality of processors, the processors running in parallel and coupled to the memory.

[0448] like Fig.34 As shown, the system includes a trained atrous convolutional neural network (ACNN) running on a large number of processors. ACCN uses atrous / dilated convolutions that allow for a large receptive field with a small number of trainable parameters. The kernel of the atrous / dilated convolution is applied to an area greater than its length by skipping input values ​​with a certain step size, also known as the atrous convolution rate or dilation factor. Atrous / dilated convolution increases the spacing of convolution filter / kernel elements so that adjacent input entries (such as nucleotides, amino acids) with larger spacing can be considered when performing convolution operations. This allows long-range context dependencies to be incorporated into the input. Atrous convolution retains part of the convolution calculation for reuse when processing adjacent nucleotides.

[0449] like Fig.34 As shown, ACNN classifies the target nucleotides in the input sequence and assigns the splice site score, making it possible for each nucleotide in the target nucleotide to be a donor splice site, an acceptor splice site or a non-splicing site. The input sequence includes at least 801 nucleotides, each of which has at least 400 nucleotides on both sides.

[0450] like Fig.34 As shown, the system also includes a classifier running on at least one large number of processors, and the classifier processes the reference sequence and the variant sequence through ACNN, so that each nucleotide in the reference sequence and the variant sequence is a donor splice site, an acceptor splice site or a non-splicing site. Each reference sequence and variant sequence has at least 101 target nucleotides, and each target nucleotide has at least 400 nucleotides on both sides. Fig.33 Reference sequences and alternative, variant sequences are described.

[0451] like Fig.34 As shown, based on the difference in the splice site scores of the target nucleotides in the reference sequence and the variant sequence, the system also determines whether the variant generating the variant sequence causes abnormal splicing and is pathogenic.

[0452] Embodiments of the disclosed systems and other systems optionally include one or more of the following features. The systems may also include features related to those described in the disclosed methods. For the sake of brevity, selective combinations of system features are not listed separately. For each statutory set of basic features, features applicable to manufacturing systems, methods, and articles are not repeated. The reader should understand how the features identified in this section can be easily combined with basic features in other statutory categories.

[0453] like Fig.34 As shown, differences in splice site scores can be determined by the position of the target nucleotide in the reference sequence and the variant sequence.

[0454] like Fig.34 As shown, for at least one target nucleotide position, when the global maximum difference in splice site scores is greater than a preset threshold, ACNN classifies variants that cause abnormal splicing and are pathogenic.

[0455] like Fig.17 As shown, for at least one target nucleotide position, when the global maximum difference in splice site scores is less than a preset threshold, ACNN classifies benign variants that do not cause abnormal splicing.

[0456] The threshold can be determined by a plurality of candidate thresholds, including processing a first set of reference sequences and variant sequences generated by benign common variants to generate a first set of abnormal splicing detections, processing a second set of reference sequences and variant sequences generated by pathogenic rare variants to generate a second set of abnormal splicing detections, and selecting at least one threshold for the classifier, wherein the threshold maximizes the number of abnormal splicing detections of the second set and minimizes the number of abnormal splicing detections of the first set.

[0457] In one embodiment, ACNN identifies variants that cause autism spectrum disorder (ASD). In another embodiment, ACNN identifies variants that cause developmental delay disorder (DDD).

[0458] like Fig.36 As shown, each reference sequence and variant sequence has at least 101 target nucleotides, each target nucleotide being flanked by at least 500 nucleotides.

[0459] like Fig.36 As shown, the splice site score of the target nucleotide in the reference sequence can be encoded in the first output of the ACNN, and the splice site score of the target nucleotide in the variant sequence can be encoded in the second output of the ACNN. In another embodiment, the first output is encoded as a first 101x3 matrix, and the second output is encoded as a second 101x3 matrix.

[0460] like Fig.36 As shown, in such an embodiment, each row of the first 101x3 matrix represents only the splice site score, making it possible for the target nucleotide in the reference sequence to be a donor splice site, an acceptor splice site, or a non-splice site.

[0461] Likewise, if Fig.36 As shown, in such an embodiment, each row of the second 101x3 matrix represents only the splice site score, making it possible that the target nucleotide in the variant sequence is a donor splice site, an acceptor splice site, or a non-splice site.

[0462] like Fig.36 As shown, in some embodiments, the splice site scores for each row of the first 101x3 matrix and the second 101x3 matrix can be exponentially normalized to sum to one.

[0463] like Fig.36 As shown, the classifier can compare the first 101x3 matrix and the second 101x3 matrix row by row and determine the change in the distribution of the splice site score based on the row. For at least one row-by-row comparison instance, when the distribution change is greater than a predetermined threshold, the ACNN classifies the variant that causes abnormal splicing and is pathogenic.

[0464] The system includes a one-hot encoder (such as Fig.29As shown in Figure 2, a one-hot encoder performs sparse encoding on the reference and variant sequences.

[0465] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

[0466] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above system actions. Other embodiments may also include methods for performing the above system actions.

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

[0468] The method includes classifying each nucleotide in the target subsequence into a donor splice site, an acceptor splice site and a non-splicing site, and processing the reference sequence through a trained atrous convolutional neural network (ACNN) to detect differential splicing patterns in the target subsequence of the input sequence.

[0469] Based on the processing, the method includes detecting differential splicing patterns in a reference target subsequence by classifying each nucleotide in the target subsequence into a donor splice site, an acceptor splice site, and a non-splice site.

[0470] The method includes processing a variant sequence by ACNN. The variant sequence differs from the reference sequence by at least one variant nucleotide, and the variant nucleotide is located in the variant target subsequence.

[0471] Based on the processing, the method includes detecting a second differential splicing pattern in the variant target subsequence by classifying each nucleotide in the variant target subsequence as a donor splice site, an acceptor splice site, and a non-splice site.

[0472] On a nucleotide by nucleotide basis, the method includes determining the difference between the first aberrant splicing pattern and the second aberrant splicing pattern by comparing the splice site classifications of the reference target subsequence and the variant target subsequence.

[0473] When the difference is greater than a predetermined threshold, the method includes classifying and storing in a memory the variant that causes abnormal splicing and is pathogenic.

[0474] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

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

[0476] The reference target subsequence and variant target subsequences can be aligned at nucleotide positions and differ by at least one variant nucleotide.

[0477] Each reference target subsequence and variant target subsequence has at least 40 target nucleotides, and each reference target subsequence and variant target subsequence is flanked by at least 40 nucleotides.

[0478] Each reference target subsequence and variant target subsequence has at least 101 target nucleotides, and each reference target subsequence and variant target subsequence is flanked by at least 5000 nucleotides.

[0479] A variant target subsequence may include two variants.

[0480] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above-mentioned system actions. Other embodiments may also include a system for performing the method, the system including a memory and one or more processors operable by the instructions stored in the memory.

[0481] Described herein are systems, methods, and articles for the manufacture of splice sites and abnormal splicing in genomic sequences (such as nucleotide sequences) using trained convolutional neural networks. One or more features of an embodiment may be combined with a basic embodiment. Mutually non-exclusive embodiments are taught to be combinable. One or more features of an embodiment may be combined with other embodiments. The present disclosure periodically reminds the user of these options. The omission of certain embodiment descriptions repeating these options should not be considered a limitation on the combination of the preceding parts, and these descriptions are incorporated forward into the following embodiments by reference.

[0482] The system embodiment of the present technical disclosure includes one or more processors coupled to a memory. The memory is loaded with computer instructions to train a splice site detector for identifying splice sites in a genomic sequence (eg, a nucleotide sequence).

[0483] The system trains a convolutional neural network (CNN) with at least 50,000 donor splice site training samples, at least 50,000 acceptor splice site training samples, and at least 100,000 non-splicing site training samples. Each training sample is a target nucleotide sequence having at least one target nucleotide, each of which has at least 20 nucleotides on each flank.

[0484] The system provides a target nucleotide sequence flanked by at least 40 upstream context nucleotides and at least 40 downstream context nucleotides as a CNN input for evaluating training samples using the CNN.

[0485] Based on this evaluation, the CNN then produces a triplet score as output, making it possible for each nucleotide in the target to be a donor splice site, an acceptor splice site, or a non-splicing site.

[0486] The disclosed system embodiments and other systems optionally include one or more of the following features. The system may also include features related to the disclosed method. For the sake of brevity, optional combinations of system features are not listed one by one. For each statutory set of basic features, features applicable to manufacturing systems, methods and articles are not repeated. The reader should understand how the features identified in this section can be easily combined with basic features in other statutory categories.

[0487] The input may include a target nucleotide sequence having a target nucleotide flanked by 100 nucleotides on each side. In such an embodiment, the target nucleotide sequence is flanked by 200 upstream context nucleotides and 200 downstream context nucleotides.

[0488] like Fig.28 As shown, the system trains CNN with at least 150,000 donor splice site training samples, at least 150,000 acceptor splice site training samples, and at least 100,000 non-splicing site training samples.

[0489] like Fig.31 As shown in FIG, CNN is parameterized through multiple convolutional layers, multiple convolutional filters, and multiple subsampling layers (such as maximum pooling and average pooling).

[0490] like Fig.31 As shown, a CNN may include one or more fully connected layers and a terminal classification layer.

[0491] CNNs may include dimension-transforming convolutional layers that transform the spatial and feature dimensions of the aforementioned inputs.

[0492] The triplet scores for each nucleotide in the target nucleotides may be exponentially normalized to sum to 1. In such an embodiment, the system classifies each nucleotide in the target nucleotides as a donor splice site, an acceptor splice site, or a non-splice site based on the highest score of the respective triplet scores.

[0493] like Fig.32 As shown, CNN evaluates the training samples in batches over a period of time. The training samples are randomly sampled in batches. Each batch has a predetermined batch size. CNN evaluates the training samples iteratively over multiple time periods (e.g., 1-10 time periods).

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

[0495] The system includes a one-hot encoder (such as Fig.29 As shown in Figure 2, the one-hot encoder sparsely encodes the training samples and provides the one-hot encoding as input.

[0496] The CNN is parameterized by multiple residual blocks, multiple skip connections, and multiple residual connections.

[0497] Each residual block may include at least one batch normalization layer, at least one rectified linear unit (ReLU) layer, at least one dimension change layer, and at least one residual connection. Each residual block may include two batch normalization layers, two ReLU nonlinear layers, two dimension change layers, and one residual connection.

[0498] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

[0499] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above system actions. Other embodiments may also include methods for performing the above system actions.

[0500] Other system embodiments of the disclosed technology include a splice site predictor trained on a large number of parallel processors coupled to a memory. The system trains a convolutional neural network (ACNN) running on a large number of processors using at least 5,000 donor splice site training samples, at least 5,000 acceptor splice site training samples, and at least 100,000 non-splicing site training samples. Each training sample used for training is a nucleotide sequence, and the nucleotide sequence includes a target nucleotide with at least 400 nucleotides on both sides.

[0501] The system includes an input stage of an ACNN, the CNN running on at least one of the plurality of processors and providing an input sequence of at least 801 nucleotides to evaluate the target nucleotide. Each target nucleotide has at least 400 nucleotides on either side. In other embodiments, the system includes an input module of a CNN, the ACNN running on at least one of the plurality of processors and providing an input sequence of at least 801 nucleotides to evaluate the target nucleotide.

[0502] The system comprises an input stage of a CNN, and the ACNN runs on at least one of the plurality of processors and converts the analysis of the CNN into a classification score, making it possible that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site. In other embodiments, the system comprises an input stage of a CNN, and the ACNN runs on at least one of the plurality of processors and converts the analysis of the CNN into a classification score, making it possible that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site.

[0503] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

[0504] The CNN was trained with 150,000 donor splice site training samples, 150,000 acceptor splice site training samples, and 800,000,000 non-splice site training samples.

[0505] CNN can be trained by one or more training servers.

[0506] The trained CNN can be deployed to one or more production servers that receive input sequences from requesting clients. In one such embodiment, the production server processes the input sequence through the input and output stages of the CNN to produce outputs, which are sent to the client. In other embodiments, the production server processes the input sequence through the input and output stages of the ACNN to produce outputs, which are sent to the client.

[0507] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above system actions. Other embodiments may also include methods for performing the above system actions.

[0508] Method embodiments of the disclosed technology include training a splice site predictor for identifying splice sites in a genomic sequence (e.g., a nucleotide sequence). The method includes providing a convolutional neural network (CNN) and an input sequence having at least 801 nucleotides to evaluate a target nucleotide, with at least 400 nucleotides flanking each side of the target nucleotide.

[0509] The CNN is trained by at least 50,000 donor splice site training samples, at least 50,000 acceptor splice site training samples and at least 100,000 non-splicing site training samples. Each training sample used for training is a nucleotide sequence, and the nucleotide sequence includes a target nucleotide with at least 400 nucleotides on both sides.

[0510] The method also includes converting the CNN analysis into a classification score, making it possible to determine whether each target nucleotide is a donor splice site, an acceptor splice site, or a non-splicing site.

[0511] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

[0512] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above-mentioned system actions. Other embodiments may also include a system for performing the method, the system including a memory and one or more processors operable by the instructions stored in the memory.

[0513] Another system embodiment of the present technical disclosure further includes one or more processors coupled to a memory. The memory is loaded with computer instructions to implement an abnormal splicing detector running on a plurality of processors, the processors running in parallel and coupled to the memory.

[0514] The system includes a trained convolutional neural network (CNN for short) running on a large number of processors.

[0515] like Fig.34 As shown, CNN classifies the target nucleotides in the input sequence and assigns the splice site score, making it possible for each nucleotide in the target nucleotide to be a donor splice site, an acceptor splice site or a non-splicing site. The input sequence includes at least 801 nucleotides, each of which has at least 400 nucleotides on both sides.

[0516] like Fig.34 As shown, the system also includes a classifier running on at least one large number of processors, and the classifier processes the reference sequence and the variant sequence through CNN, so that each nucleotide in the reference sequence and the variant sequence is a donor splice site, an acceptor splice site or a non-splicing site. Each reference sequence and variant sequence has at least 101 target nucleotides, and each target nucleotide has at least 400 nucleotides on both sides.

[0517] like Fig.34 As shown, based on the difference in the splice site scores of the target nucleotides in the reference sequence and the variant sequence, the system also determines whether the variant generating the variant sequence causes abnormal splicing and is pathogenic.

[0518] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not repeated herein and are considered repeated by reference.

[0519] The difference in splice site scores can be determined by the position of the target nucleotide in the reference sequence and the variant sequence.

[0520] For at least one target nucleotide position, when the global maximum difference in splice site scores is greater than a preset threshold, CNN classifies variants as causing aberrant splicing and pathogenic.

[0521] For at least one target nucleotide position, when the global maximum difference in splice site scores is less than a preset threshold, ACNN classifies benign variants that do not cause aberrant splicing.

[0522] The threshold can be determined by a plurality of candidate thresholds, including processing a first set of reference sequences and variant sequences generated by benign common variants to generate a first set of abnormal splicing detections, processing a second set of reference sequences and variant sequences generated by pathogenic rare variants to generate a second set of abnormal splicing detections, and selecting at least one threshold for the classifier, wherein the threshold maximizes the number of abnormal splicing detections of the second set and minimizes the number of abnormal splicing detections of the first set.

[0523] In one embodiment, ACNN identifies variants that cause autism spectrum disorder (ASD). In another embodiment, ACNN identifies variants that cause developmental delay disorder (DDD).

[0524] Each reference sequence and variant sequence has at least 101 target nucleotides, each target nucleotide being flanked by at least 1000 nucleotides.

[0525] The splice site scores of the target nucleotides in the reference sequence can be encoded in the first output of the CNN, and the splice site scores of the target nucleotides in the variant sequence can be encoded in the second output of the CNN. In one embodiment, the first output is encoded as a first 101x3 matrix, and the second output is encoded as a second 101x3 matrix.

[0526] In such an embodiment, each row of the first 101x3 matrix represents only the splice site score, making it possible that the target nucleotide in the reference sequence is a donor splice site, an acceptor splice site, or a non-splice site.

[0527] Likewise, in such an embodiment, each row of the second 101x3 matrix represents only the splice site score, making it possible that the target nucleotide in the variant sequence is a donor splice site, an acceptor splice site, or a non-splice site.

[0528] In some embodiments, the splice site scores for each row of the first 101x3 matrix and the second 101x3 matrix can be exponentially normalized to sum to unity.

[0529] The classifier can compare the first 101x3 matrix and the second 101x3 matrix row by row and determine the change in the distribution of the splice site score based on the row. For at least one row by row comparison instance, when the distribution change is greater than a predetermined threshold, the CNN classifies the variant that causes abnormal splicing and is pathogenic.

[0530] The system includes a one-hot encoder (such as Fig.29 As shown in Figure 2, a one-hot encoder performs sparse encoding on the reference and variant sequences.

[0531] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above system actions. Other embodiments may also include methods for performing the above system actions.

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

[0533] The method includes classifying each nucleotide in the target sequence into a donor splice site, an acceptor splice site and a non-splicing site, and processing the reference sequence through a trained atrous convolutional neural network (CNN) to detect differential splicing patterns in the target sequence of the input sequence.

[0534] Based on the processing, the method includes detecting differential splicing patterns in a reference target subsequence by classifying each nucleotide in the target subsequence into a donor splice site, an acceptor splice site, and a non-splice site.

[0535] The method includes processing a variant sequence by CNN. The variant sequence differs from the reference sequence by at least one variant nucleotide, and the variant nucleotide is located in the variant target subsequence.

[0536] Based on the processing, the method includes detecting a second differential splicing pattern in the variant target subsequence by classifying each nucleotide in the variant target subsequence as a donor splice site, an acceptor splice site, and a non-splice site.

[0537] On a nucleotide by nucleotide basis, the method includes determining the difference between the first aberrant splicing pattern and the second aberrant splicing pattern by comparing the splice site classifications of the reference target subsequence and the variant target subsequence.

[0538] When the difference is greater than a predetermined threshold, the method includes classifying and storing in a memory the variant that causes abnormal splicing and is pathogenic.

[0539] Each feature discussed in this specific embodiment section for other system and method embodiments is also applicable to this system embodiment.As mentioned above, all system features are not described herein again, and all system features are considered to be repeated by reference.

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

[0541] The reference target subsequence and variant target subsequences can be aligned at nucleotide positions and differ by at least one variant nucleotide.

[0542] Each reference target subsequence and variant target subsequence has at least 40 target nucleotides, and each reference target subsequence and variant target subsequence is flanked by at least 40 nucleotides.

[0543] Each reference target subsequence and variant target subsequence has at least 101 target nucleotides, and each reference target subsequence and variant target subsequence is flanked by at least 1000 nucleotides.

[0544] A variant target subsequence may include two variants.

[0545] Other embodiments may include a non-transitory computer-readable storage medium storing instructions, and a processor executes the instructions to perform the above-mentioned system actions. Other embodiments may also include a system for performing the method, the system including a memory and one or more processors operable by the instructions stored in the memory.

[0546] The above description enables the disclosed technology to be made and used. Various modifications of the disclosed embodiments are apparent, and the basic principles defined herein can be applied to other embodiments and applications without departing from the spirit and scope of the disclosed technology. Therefore, the disclosed technology is not limited to the described embodiments, but conforms to the widest scope consistent with the principles and features disclosed herein. The scope of the disclosed technology is defined by the attached claims.

[0547] Per gene enrichment analysis

[0548] Fig.57It is an embodiment of per gene enrichment analysis. In one embodiment, the abnormal splicing detector is further configured to perform per gene enrichment analysis, and per gene enrichment analysis is used to determine the pathogenicity of variants that have been determined to cause abnormal splicing. For a specific gene sampled from a group of individuals with a gene defect, per gene enrichment analysis includes applying the trained ACNN to identify candidate variants in a specific gene that causes abnormal splicing, determining the baseline number of mutations in a specific gene based on the sum of the observed trinucleotide mutation rates of the candidate variants and the product of the sum and the number of transmissions and the size of the group, applying the trained ACNN to identify new variants in a specific gene that causes abnormal splicing and comparing the baseline number of mutations with the number of new variants. Based on the comparison results, per gene enrichment analysis determines that a specific gene is related to a gene defect, and the new variant is pathogenic. In some embodiments, the gene defect is autism spectrum disorder (abbreviated as ASD). In other embodiments, the gene defect is developmental delay disorder (abbreviated as DDD).

[0549] like Fig.57 As shown in the embodiment, the abnormal splicing detector classifies 5 candidate variants in a specific gene that causes abnormal splicing. The three-nucleotide mutation rates of these 5 candidate variants are observed to be 10 -8 , 10 -2 , 10 -1 , 10 5 and 10 1 The baseline number of mutations for a particular gene was determined to be 10 based on the sum of the observed trinucleotide mutation rates of the five candidate variants and the product of the sum and the number of transmissions / number of chromosomes (2) and the size of the cohort (1000). -5 . Then compare it with the number of new variants (3).

[0550] In certain embodiments, the aberrant splicing detector is further configured to perform the comparison using a statistical test that generates a p-value as an output.

[0551] In other embodiments, the aberrant splicing detector is further configured to compare the baseline number of mutations to the number of de novo variants, and to determine, based on the comparison, that a particular gene is not associated with the gene profile and that the de novo variant is benign.

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

[0553] In another embodiment, at least some of the candidate variants are missense variants.

[0554] Genome-wide enrichment analysis

[0555] Fig.58It is an embodiment of whole genome enrichment analysis. In another embodiment, the abnormal splicing detector is further configured to perform whole genome enrichment analysis, and whole genome enrichment analysis is used to determine the pathogenicity of variants that have been determined to cause abnormal splicing. Whole genome enrichment analysis includes applying the trained ACNN to identify a first group of new variants of multiple genes sampled from a healthy individual group that cause abnormal splicing, applying the trained ACNN to identify a second group of new variants of multiple genes sampled from an individual group with a gene defect, and comparing the number of the first group and the second group of new variants and determining based on the comparison result that the second group of new variants is enriched in an individual group with a gene defect and the second group of new variants is pathogenic. In some embodiments, the gene defect is autism spectrum disorder (abbreviated as ASD). In other embodiments, the gene defect is developmental delay disorder (abbreviated as DDD).

[0556] In certain embodiments, 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 embodiment, the comparison can be further parameterized using the respective group sizes.

[0557] In certain embodiments, the aberrant splicing detector is further configured to compare the respective numbers of the first group and the second group, and determine based on the comparison result that the second group of novel variants is not enriched in the group of individuals with the gene defect and that the second group of novel variants is benign.

[0558] like Fig.58 As shown, in addition to the amount of mutation per individual (4), the mutation rates of the healthy group (0.001) and the affected group (0.004) are also depicted.

[0559] discuss

[0560] Despite the limited diagnostic yield of exome sequencing in patients with severe genetic defects, clinical sequencing focuses on rare coding mutations, and largely ignores noncoding genomic variation due to its difficulty in interpretation. Here we introduce a deep learning network that accurately predicts splicing of the original nucleotide sequence, thereby identifying noncoding mutations that disrupt the normal pattern of exons and introns with serious consequences on the resulting protein. We show that RNA sequencing technology can validate predicted cryptic splicing mutations at high rates, which are highly toxic to humans and are the main cause of rare genetic diseases.

[0561] By using a deep learning model as a computer model of the spliceosome, we are able to reconstruct specificity determinants with great precision for the in vivo accuracy of the spliceosome. We reiterate numerous findings from studies of the splicing machinery made over the past four decades, indicating that the spliceosome incorporates a large number of short-term and long-term specificity determinants in its decision making. In particular, we find that the presence of short-term and long-term determinants, such as exon and intron length and nucleosome positioning, explains much of the perceived degradation of splicing patterns, more than necessary to compensate and render additional specificity at the graph level. Our findings demonstrate the promise of deep learning models to provide biological insights, rather than simply acting as black-box classifiers.

[0562] Deep learning is a relatively new technique in biology and is not without potential trade-offs. By learning to automatically extract features from sequences, deep learning models can exploit novel sequence determinants that are not well described by human experts, but there is also a risk that the model may incorporate features that do not reflect the true behavior of the spliceosome. These irrelevant features can increase the apparent accuracy of predicting annotated exon-intron boundaries and reduce the accuracy of predicting the splicing change effects of any sequence change caused by genetic variation. Because accurate prediction of variants provides the strongest evidence that the model generalizes to real biology, we validate the predicted splicing change variants using three completely orthogonal approaches: RNA sequencing, natural selection in human populations and comparison with control cohorts, and enrichment for de novo variants. While the incorporation of irrelevant features into the model is not completely excluded, the resulting model appears to be sufficiently faithful to the real splicing biology to have important practical applications, such as identifying cryptic splicing mutations in patients with genetic diseases.

[0563] Compared to other classes of protein-truncating mutations, cryptic splicing mutations are particularly interesting due to incompletely penetrant splice-altering variants, as the phenomenon of ubiquitous alternative splicing tends to attenuate canonical splice sites relative to alternative splice sites, resulting in a mixture of variant and normal transcripts in RNA sequencing. The observation that these variants frequently drive tissue-variant alternative splicing highlights the unexpected role of cryptic splicing mutations in giving rise to novel alternative splicing diversity. Training deep learning models on splice junction annotations from RNA sequencing of relevant tissues would be a potential direction for development in order to derive tissue-specific models of alternative splicing. Training networks on annotations directly derived from RNA sequencing data could also help fill gaps in GENCODE annotations, improving model performance in variant prediction (e.g., Fig.52A and 52B shown).

[0564] Our understanding of how mutations in the noncoding genome contribute to human disease is far from complete. The discovery of exogenous de novo cryptic splicing mutations in childhood neurodevelopmental disorders suggests that improved interpretation of the noncoding genome could directly benefit patients with severe genetic defects. Cryptic splicing mutations also play an important role in cancer (Jung et al., 2015; Sanz et al., 2010; Supek et al., 2014), and recurrent somatic mutations in splicing factors have been shown to cause extensive changes in splicing specificity. Much work remains to be done to understand splicing regulation in different tissue and cellular contexts, especially in the context of mutations that directly affect proteins in the spliceosome. With advances in oligonucleotide therapeutics that can potentially target splicing defects with sequence specificity (Finkel et al., 2017), a better understanding of the regulatory mechanisms that govern this remarkable process could pave the way for new candidates for therapeutic intervention.

[0565] Fig.37A , 37B, 37C, 37D, 37E, 37F, 37G and 37H describe an embodiment of predicting splicing from raw sequences with deep learning.

[0566] like Fig.37A For each position of the pre-mRNA transcript, SpliceNet-10k uses 10,000 nucleotides of flanking sequences as input and predicts whether the position is a splice acceptor or a splice donor or whether the position is neither a splice acceptor nor a splice donor.

[0567] Fig.37B The full pre-mRNA transcripts for CFTR gene scoring using MaxEntScan (top) and SpliceNet-10k (bottom) are shown, along with the predicted acceptor (red arrows) and donor (green arrows) sites and actual locations of exons (black boxes). For each method, we applied a threshold that made the number of predicted sites equal the total number of actual sites.

[0568] like Fig.37C In this paper, we measured the inclusion rate of each exon in RNA sequencing technology and displayed the distribution of SpliceNet-10k scores of exons at different inclusion rates. The maximum values ​​of exon acceptor and donor scores are indicated.

[0569] like Fig.37D , Computer simulation of the impact of each nucleotide around exon 9 in the U2SURP gene. The vertical dimension of each nucleotide indicates a decrease in the strength of the predicted acceptor site (black arrow) when the nucleotide is mutated (Δ score).

[0570] like Fig.37F, the size of the input sequence context affects the accuracy of the network. Top-k accuracy is the score of correctly predicted splice sites at the threshold, and the number of predicted sites is equal to the actual number of existing sites. PR-AUC is the area under the precision-recall curve. We also show the Top-k accuracy and PR-AUC of three other algorithms for splice site detection.

[0571] like Fig.37F , the relationship between exon / intron length and strength of adjacent splice sites predicted by SpliceNet-80nt (local motif score) and SpliceNet-10k. The genome-wide distribution of exon length (yellow) and intron length (pink) is shown in context. The x-axis is logarithmic.

[0572] like Figure 37G , a pair of splice acceptor and splice donor motifs separated by 150 nt are shown along the HMGCR gene. The K562 nucleosome signal at each position and the probability of exon pair formation at this position are shown, as predicted by SpliceNet-10k.

[0573] like Fig.37H , shows the average K562 and GM12878 nucleosome signals near private mutations predicted by the SpliceNet-10k model for the creation of novel exons in the GTEx cohort. P values ​​obtained by permutation tests are shown.

[0574] Fig.38A , 38B, 38C, 38D, 38E, 38F and 38G describe an implementation method for validating rare cryptic splicing mutations in RNA sequencing technology data.

[0575] like Fig.38A For rs397515893, a pathogenic cryptic splice variant in an intron of MYBPC3 associated with cardiomyopathy, SpliceNet-10k predicted the acceptor and donor scores for each position in the gene with or without mutation in the pre-mRNA sequence to assess the splicing change effect of the mutation. The mutation delta score was the maximum change in the splicing prediction score within 50 nt of the variant.

[0576] like Fig.38B, we scored private genetic variants observed in one of 149 individuals in the GTEx cohort using the SpliceNet-10k model. Enrichment of private variants predicted to alter splicing (Δ score greater than 0.2, blue) or near private exon skipping sites (top) or private acceptor and donor sites (bottom) that do not affect splicing (Δ score less than 0.01, red) is shown. The y-axis indicates the number of times private splicing events and nearby private genetic variants occur together in the same individual, compared to the expected number from permutation.

[0577] like Fig.38C , shows an example of a heterozygous synonymous variant in PYGB that creates a novel donor site with incomplete penetrance. The RNA-seq range, number of junction reads, and junction locations are shown for variant and control individuals (blue and grey arrows). The effect size was calculated as the difference in new junction usage between individuals with the variant and those without the variant. In the lower overlay bar chart, we show the number of reads for the reference allele or the alternating allele that used the annotated or new junction (no splicing and new junction, respectively). The total number of reference allele reads is significantly different from the total number of alternating allele reads, indicating that 60% of transcripts are missing splicing at the new junctions in the RNA-seq data, presumably due to nonsense-mediated decay (NMD).

[0578] like Fig.38D , shows cryptic splicing mutations predicted by the SpliceNet-10k model that were partially validated using GTEx RNA sequencing data. The validation rate of disruptions of essential acceptor or donor dinucleotides (shown by dashed lines) is less than 100% based on coverage and nonsense-mediated decay.

[0579] like Fig.38E , shows the distribution of effect sizes for cryptic splicing predictions used for validation. The dashed line (50%) corresponds to the expected effect size for a fully penetrant heterozygous variant. Effect sizes for essential acceptor or donor dinucleotide disruptions measured were less than 50% based on nonsense-mediated decay or unexplained splice form changes.

[0580] like Fig.38F , shows the sensitivity of SpliceNet-10k for detecting splicing change private variants in the GTEx cohort at different Δ score cutoffs. Variants are split into deep intragenic variants (greater than 50 nt from exons) and near exonic variants (coinciding with exons or ≤ 50 nt from exon-intron boundaries).

[0581] like Figure 38G, showing the validation rate and sensitivity of SpliceNet10-k for different confidence cutoffs of splice sites and three other methods. The three points on the SpliceNet10-k curve represent the performance of SpliceNet10-k at Δ score cutoffs: 0.2, 0.5, and 0.8. For the other three algorithms, the three points on the curve indicate the performance at thresholds that predict the same number of cryptic splice variants as SpliceNet10-k at Δ score cutoffs: 0.2, 0.5, and 0.8.

[0582] Fig.39A , 39B and 39C represent an embodiment of frequently creating cryptic splice variants of tissue-specific alternative splicing.

[0583] like Fig.39A , shows an example of a heterozygous exon variant in CDC25B that creates a novel donor site. The variant belongs to only one individual in the GTEx cohort and shows tissue-specific alternative splicing that favors a large fraction of novel splice forms in muscle compared to fibroblasts (P = 0.006 by Fisher's exact test). In muscle and fibroblasts, the RNA sequencing technology range, junction read counts, and junction locations are shown for variant individuals and controls (blue and gray arrows).

[0584] like Fig.39B , shows an example of a heterozygous exogenous receptor-creating variant in FAM229B that exhibits consistent tissue-specific effects in all three individuals in the GTEx cohort with the variant. Arterial and lung RNA sequencing techniques for three individuals with the variant and a control individual are shown.

[0585] like Fig.39C , A subset of splice site creation variants in GTEx that were significantly associated with novel junctions across expressed tissues, as assessed by chi-square test. Validated cryptic splice variants with medium and low Δ scores were more likely to cause tissue-specific alternative splicing ((P = 0.015 by Fisher's exact test).

[0586] Fig.40A , 40B , 40C, 40D and 40E describe an embodiment of a predicted cryptic splice variant that is highly toxic in a human population.

[0587] like Fig.40A Synonymous mutations and intragenic variants with reliable predicted splicing change effects (Δ score ≥ 0.8) in the population are greatly depleted at general allele frequency relative to rare variants observed in only 60,706 individuals. -127) showed that 78% of the recently emerged predicted cryptic splice variants were sufficiently toxic to be removed by natural selection.

[0588] like Fig.40B , showing protein truncating variants and predicted synonymous mutations and intragenic cryptic splice variants in a subset of deleterious ExAC datasets, which are calculated in (A).

[0589] like Fig.40C , showing the subset of synonymous mutations and intragenic cryptic splicing-gained variants in the ExAC dataset that are partially deleterious (Δ score ≥ 0.8), split based on whether the variant causes a frameshift.

[0590] like Fig.40D , showing protein-truncating variants and predicted deep intragenic (greater than 50 nt from known exon-intron boundaries) cryptic splice variants in a subset of deleterious gnomAD datasets.

[0591] like Fig.40E , showing the average number of rare (allele frequency less than 0.1%) protein-truncating variants and rare functional cryptic splice variants per individual human genome. The number of cryptic splice variants predicted to be functional was assessed based on the deleterious prediction score. The total number of predictions is higher.

[0592] Fig.41A , 41B, 41C, 41D, 41E and 41F represent one embodiment of de novo cryptic splicing mutations in patients with rare inherited diseases.

[0593] like Fig.41A , shows predicted cryptic splicing de novo variants in patients from each coding developmental disorder cohort, individuals with autism spectrum disorder from the Autism Sequencing Consortium, and healthy controls. Variants identified between the DDD and ASD cohorts were adjusted for enrichment over healthy controls. Error bars represent 95% confidence intervals.

[0594] like Fig.41B , shows the estimated proportion of pathogenic de novo mutations in the functional categories of the DDD and ASD cohorts compared to healthy controls, based on the enrichment of each category.

[0595] like Fig.41C , showing the enrichment and excess of cryptic splicing de novo mutations in the DDD and ASD cohorts compared to healthy controls at different Δ score thresholds.

[0596] like Fig.41D , shows a list of novel candidate disease genes enriched for de novo mutations in the DDD and ASD cohorts (FDR < 0.01) when both predicted cryptic splicing mutations and protein-coding mutations are included in the enrichment analysis. Phenotypes present in multiple individuals are shown.

[0597] like Fig.41E , three examples of predicted de novo cryptic splicing mutations in autism patients validated by RNA sequencing technology, causing intron retention, exon skipping, and exon expansion, respectively. In each example, the RNA sequencing technology range and number of nodes of the affected monomer are shown at the top, and the unmutated control monomer is shown at the bottom. Sequences related to gene transcription are shown on the positive strand. Blue arrows and gray arrows distinguish the junction positions of individuals with variants and controls, respectively.

[0598] Fig.41F Validation status of 36 predicted cryptic splice sites for experimental validation by RNA sequencing technology.

[0599] Experimental model and subject details

[0600] Iosifov et al. previously published detailed information on 36 autistic patients in Nature 2014 (Table S1), which can be cross-referenced using the anonymous identifiers in column 1 of Table S4 in our paper.

[0601] Method details

[0602] I. Deep Learning for Splicing Prediction

[0603] SpliceNet architecture

[0604] We trained several ultra-deep convolutional neural network-based models to computationally predict the splicing of pre-mRNA nucleotide sequences. We designed four structures, namely, SpliceNet-80nt, SpliceNet-400nt, SpliceNet-2k, and SpliceNet-10k, which take as input 40, 200, 1,000, and 5,000 nucleotides on each side of the target position, respectively, and output the probability of the position being a splice acceptor and a splice donor. More precisely, the input of the model is a one-hot encoded nucleotide sequence, where A, C, G, and T (or equivalently U) are 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 (summing to 1), corresponding to the probability that the target position is a splice acceptor, a splice donor, and neither.

[0605] The basic unit of the SpliceNet structure 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 way ( Fig.21 , 22, 23 and 24). Residual blocks are often 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 difficult to train due to the exploding / vanishing gradient problem (Glorot and Bengio, 2010), and increasing the depth of such neural networks often led to large training errors (He et al., 2016a). Through a series of computational experiments, it was shown that a structure consisting of many residual blocks stacked one after another overcomes these problems (He et al., 2016a).

[0606] Fig.21 , 22 , 23 and 24 provide the complete SpliceNet architecture. The architecture consists of K stacked residual blocks connecting the input layer to the penultimate layer and convolutional units with softmax activation connecting the penultimate layer to the output layer. The residual blocks are stacked so 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 third residual block is added to the input of the penultimate layer. Such “skip connections” are often used in deep neural networks to increase the convergence rate during training (Oord et al., 2016).

[0607] 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 Koltum, 2016). Since a convolution kernel with window size W and dilation rate D can extract features across (W-1)D adjacent positions, a residual block with hyperparameters N, W, and D can extract features across 2(W-1)D adjacent positions. Therefore, The total adjacent span of the SpliceNet structure is obtained, where N i , W i and D i is the hyperparameter of the ith residual block. For SpliceNet-80nt, SpliceNet-400nt, SpliceNet-2k, and SpliceNet-10k structures, the number of residual blocks and the hyperparameters of each residual block are selected so that S is equal to 80, 400, 2,000, and 10,000, respectively.

[0608] In addition to the convolutional units, the SpliceNet structure has only normalization and nonlinear activation units. As a result, the model can be used in a sequence-to-sequence mode with variable sequence lengths (Oord et al., 2016). For example, the input of the SpliceNet-10k model is a one-hot encoded nucleotide sequence of length S / 2+l+S / 2, and the output is a lx 3 matrix corresponding to the three scores of the l-center position in the input, i.e., the position remaining after excluding the first and last S / 2 nucleotides. This feature can be exploited to save huge computations during training and testing. This is because most of the computations on positions close to each other are common, and when the shared computation is used in sequence-to-sequence mode, it only needs to be performed once by the model.

[0609] Our model uses a residual block structure, which has been widely adopted due to its success in image classification. The residual module consists of repeated units of convolutions interspersed with skip connections that allow information from earlier layers to skip the residual block. In each residual block, the input layer is first batch normalized, and then the activation layer is normalized using a rectified linear unit (ReLU). The activations then pass through a 1D convolutional layer. The intermediate output of the 1D convolutional layer is again batch normalized and activated with a ReLU, followed by another 1D convolutional layer. At the end of the second 1D convolution, its output is summed with the original input and added to the residual block, acting as a skip connection by allowing the original input information to bypass the residual block. In this structure, which the authors call a deep residual learning network, the input remains in its original state and the residual connections are not subject to the nonlinear activations of the model, allowing for efficient training of deeper networks.

[0610] After the residual block, the softmax layer calculates the probability of the three states for each amino acid, where the maximum softmax probability determines the state of the amino acid. The model is trained using the ADAM optimizer with the cumulative categorical cross entropy loss function of the entire protein sequence.

[0611] Dilated / atrous convolutions allow for large receptive fields with few trainable parameters. A dilated / atrous convolution is a convolution that applies a kernel to an area larger than its length by skipping input values ​​by a certain step size, also called the dilated convolution rate or dilation factor. Dilated / atrous convolutions add spacing between elements of the convolution filter / kernel so that more spaced neighboring input items (e.g., nucleotides, amino acids) can be considered when performing the convolution operation. This makes it possible to incorporate long-range contextual dependencies into the input. Dilated convolutions reserve part of the convolution computation for reuse when processing neighboring nucleotides.

[0612] The examples shown use 1D convolutions. In other embodiments, the model may use different types of convolutions, such as 2D convolutions, 2D convolutions, dilated or dilated convolutions, transposed convolutions, separable convolutions, and depthwise separable convolutions. Some layers also use ReLU activation functions, which greatly speed up the convergence of stochastic gradient descent compared to saturating nonlinearities such as sigmoid or hyperbolic tangent. Other examples of activation functions used by the disclosed technology include parametric ReLU, leaky ReLU, and exponential linear unit (ELU).

[0613] Some layers also use batch normalization (Ioffe and Szegedy 2015). Regarding batch normalization, the distribution of each layer in a convolutional neural network (CNN) changes during training and the distribution of each layer is different from each other. This reduces the convergence rate of the optimization algorithm. Batch normalization is a technique to overcome this problem. Denote the input of the batch normalization layer by x and its output by z. Batch normalization applies the following transformation to x:

[0614]

[0615] Batch Normalization applies mean variance normalization to the input x by using μ and σ and linear scaling, and shifts it by using γ and β. The normalization parameters μ and σ of the current layer are calculated 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 μ and σ values ​​calculated during training are used in the forward pass during inference.

[0616] Model training and testing

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

[0618] We used the following procedure to train and test the model in sequence-to-sequence mode with a block size of l = 5,000. For each gene, the mRNA transcript sequence between the canonical transcription start and stop sites was extracted from the hg19 / GRCh37 assembly. The input mRNA transcript sequence was one-hot encoded with the following sequence: A, C, G, T / U 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 to a multiple of 5,000 in length and then further zero-padded at the start and end with a flanking sequence length of S / 2, where S equals 80, 400, 2,000, and 10,000 for the SpliceNet-80nt, SpliceNet-400nt, SpliceNet-2k, and SpliceNet-10k models, respectively. The padded nucleotide sequence is then split into blocks of length S / 2+5,000+S / 2, such that the i-th block consists of nucleotide positions from 5,000(i-1)-S / 2+1 to 5,000i+S / 2. Similarly, the splicing output label sequence is one-hot encoded as follows: the non-splicing site, splice acceptor (corresponding to the first nucleotide of the exon), and splice donor (corresponding to the last nucleotide of the exon) are mapped to [1,0,0], [0,1,0], and [0,0,1], respectively. The one-hot encoded splicing output label sequence is zero-padded until the length is a multiple of 5,000, and then split into blocks of length 5,000 again, such that the i-th block consists of positions from 5,000(i-1)+1 to 5,000i. The one-hot encoded nucleotide sequence and the corresponding one-hot encoded label sequence are used as the input and target output of the model, respectively.

[0619] The models were trained 10 times on two NVIDIA GeForce GTX 1080Ti-GPUs with a batch size of 12. The Adam optimizer (Kingma and Ba, 2015) was used to minimize the categorical cross entropy loss between the target and predicted outputs during training. The learning rate of the optimizer was set to 0.001 for the first 6 times and then reduced by half for each subsequent time. For each architecture, we repeated the training process 5 times and obtained 5 trained models ( Fig.53A and 53B ). During testing, all five trained models are used to evaluate each input and the average of their outputs is used as the predicted output. We use the model used for the analysis in Figure 37A and other related figures.

[0620] For identification of variants involving splicing alterations Fig.38A-G, 39A-C, 40A-E, and 41A-F, we augmented the training set of GENCODE annotations to 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). This increased the number of splice junction annotations in the training set by approximately 50%. Training the network on the combined dataset improved the sensitivity for detecting splice alteration variants in RNA sequencing data, particularly for predicting deep intronic splice alteration variants, compared to networks trained on GENCODE annotations alone ( Fig.52A and 52B ), we use this network for analyses involving variant evaluation ( Fig.38A -G, 39A-C, 40A-E, and 41A-F and related figures). To ensure that there was no overlap between training and evaluation in the GTEx RNA-seq dataset, we only included junctions present in more than 5 individuals in the training dataset and only evaluated the performance of the network on variants present in 4 or fewer individuals. Details of novel splice junction identification are described in “Splice junction detection” in the GTEx analysis section of the Methods.

[0621] Top-k Accuracy

[0622] Since most positions are not splice sites, accuracy measures like the percentage of positions correctly classified are largely ineffective. Instead, we evaluate the model using two metrics that work well in this context, namely top-k accuracy and the area under the precision-recall curve. The top-k accuracy for a particular class is defined as: Suppose the test set has k positions that belong to that class. We choose the threshold so that we can accurately predict whether k test set positions belong to that class. The fraction of these k predicted positions that actually belong to that class is reported as the top-k accuracy. In practice, this is equivalent to the accuracy when the threshold is chosen so that precision and recall have the same value.

[0623] Model evaluation of lincRNA

[0624] 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 of the longest total exonic sequence for each lincRNA gene and referred to it as the canonical transcript of the gene. Since the reliability of lincRNA annotations is lower than that of protein-coding gene annotations and this misannotation affects our estimates of top-k accuracy, we used GTEx data to eliminate lincRNAs with potential annotation issues (see the “GTEx Dataset Analysis” section below for details on these data). For each lincRNA, we counted all split reads mapped across the length of the lincRNA across all GTEx samples (see “Splice Junction Detection” below for details). This is an estimate of the total junction-spanning reads for the lincRNA using either an annotated junction or a new junction. We also counted the number of reads that spanned the junction of the canonical transcript. We only considered lincRNAs for which at least 95% of the junction-spanning reads across all GTEx samples corresponded to the canonical transcript. We also required that all junctions of canonical transcripts be observed at least once in the GTEx cohort (excluding junctions spanning introns with length <10 nt). To calculate top-k accuracy, we only considered junctions of canonical transcripts of lincRNAs that passed the above filters (781 transcripts, 1047 junctions).

[0625] Identification of splice junctions from pre-mRNA sequences

[0626] exist Fig.37B In this paper, we compared the performance of MaxEntScan and SpliceNet-10k on identifying the canonical exon boundaries of genes from their sequences. We studied the CFTR gene, which is present in our test set and has 26 canonical splice acceptors and donors, as a case study and obtained the acceptor and donor scores for each of the 188,703 positions from the canonical transcription start site (chr7:117120017) to the canonical transcription end site (chr7:117308719) using MaxEntScan and SpliceNet-10k. A position was classified as a splice acceptor or donor if its corresponding score was greater than a chosen threshold while evaluating the top-k accuracy. MaxEntScan predicted 49 splice acceptors and 22 splice donors, of which 9 acceptors and 5 donors were true splice acceptors and donors, respectively. For better visualization, we show the pre-log scores of MaxScan (clipped to a maximum value of 2,500). SplitNet-10k predicted 26 splice acceptors and 26 splice donors, all of which were correct. Fig.42B , we repeated the analysis using the LINC00467 gene.

[0627] Assessment of exon inclusion at GENCODE-annotated splice junctions

[0628] We calculated the inclusion rate of all GENCODE-annotated exons based on GTEx RNA sequencing data ( Fig.37C ). For each exon excluding the first and last exon of each gene, we calculated the inclusion ratio as:

[0629]

[0630] Where L is the total read count of the junction from the previous canonical exon to the exon under consideration in 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 jump junction from the previous canonical exon to the next canonical exon.

[0631] The significance of each nucleotide in splice site recognition

[0632] exist Fig.37D In , we identified nucleotides that SplitNet-10k considered important for classifying a position as a splice acceptor. To do this, we considered the splice acceptor at chr3:142740192 in the U2SUPR gene present in our test set. The “importance score” of a nucleotide with respect to a splice acceptor is defined as follows: ref Denote the acceptor score of the considered splice acceptor. Recalculate the acceptor score by replacing the considered nucleotide with A, C, G, and T. Let these scores be represented by S A , S C , S G and S T The importance score of the nucleotide is estimated as:

[0633]

[0634] This process is often referred to as in-silico mutagenesis (Zhou and Troyanskaya, 2015). We mapped the 127 nucleotides from chr3:142740137 to chr3:142740263 so that the height of each nucleotide is its importance score with respect to the splice acceptor at chr3:142740192. The plotting function was adapted from DeepLIFT (Shrikumar et al., 2017) software.

[0635] Effects of TACTAAC and GAAGAA motifs on splicing

[0636] To investigate the effect of the position of the branch point sequence on receptor strength, we first used SplitNet-10k to obtain the receptor scores of 14,289 test set splice acceptors, y ref y represents the vector containing these scores. For each value of i from 0 to 100, we did the following: For each test set splice acceptor, we replaced the nucleotides at positions i to i-6 before the splice acceptor with TACTAAC and recalculated the acceptor scores using SplitNet-10k. The vector containing these scores is y alt,i express.

[0637] We are Fig.43A The following quantities are plotted as a function of i:

[0638] mean(y alt,i -y ref )

[0639] for Fig.43B , we repeated the same procedure using the SR-protein motif GAAGAA. In this case, we also investigated the impact of the motif when present after the splice acceptor and the effect on the donor strength. Based on a comprehensive study of k-mer space, GAAGAA and TACTAAC were the motifs with the greatest impact on both acceptor and donor strength.

[0640] The role of exon and intron length in splicing

[0641] To investigate the effect of exon length on splicing, we filtered out the 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 in order of increasing length. For each, we calculated a splicing score using SplicNet-80nt by averaging the acceptor score at the splice acceptor site and the donor score at the splice donor site. Fig.37F In , we plotted the splicing scores as a function of exon length. Before plotting, we applied the following smoothing procedure: x represents a vector containing the exon lengths and y represents a vector containing their corresponding splicing scores. We smoothed x and y using an average window of size 2,500.

[0642] We repeated this analysis by computing the splicing scores using SplicNet-10k. In the background, we show the length histogram of the 12,637 exons used for this analysis. We applied a similar analysis to investigate the effect of intron length on splicing, with the main difference that the first and last exons did not need to be excluded.

[0643] The role of nucleosomes in splicing

[0644] We downloaded nucleosome data of K562 cell line from UCSC genome browser. We used HMGR gene in the test set as an anecdotal example to demonstrate the influence of nucleosome positioning on SplicNet-10k score. For each position p in the gene, we calculated its "cultivated splicing score" as follows:

[0645] • Eight nucleotides from position p+74 to p+81 were replaced by the donor motif AGGTAAGG.

[0646] • Four nucleotides from p-78 to p-75 are replaced by the acceptor motif TAGG.

[0647] · 20 nucleotides from position p-98 to p-79 were replaced by the polypyrimidine tract CCTCCTTCGCCCTC.

[0648] · Seven nucleotides from p-105 to p-99 were replaced by the branch point sequence CACTAAC.

[0649] The average of the acceptor score at p-75 and the donor score at p-75 predicted by SplitNet-10k was taken as the implanted splicing score.

[0650] Figure 37G K562 nucleosome signals and embedded splicing scores at 5,000 positions from chr5:74652154 to chr5:74657153 are shown.

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

[0652] 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 in order of increasing nucleosome enrichment and calculated their acceptor scores using SplitNet-80nt. Fig.44BIn , the receptor score is a function of nucleosome enrichment. Before plotting, the Fig.37F We repeated this analysis using SplitNet-10k for 14,289 test splice donors.

[0653] Enrichment of nucleosome signals at novel exons

[0654] exist Fig.37H In this study, we wanted to look at the nucleosome signal around predicted novel exons. To ensure that we were looking at high-confidence novel exons, we only selected single variants (variants present in a single GTEx individual) where the predicted gained junction was completely private to the individual with the variant. Furthermore, to eliminate the confounding effects of nearby exons, we only looked at intronic variants that were at least 750 nt away from annotated exons. We downloaded nucleosome signals for GM12878 and K562 cell lines from the UCSC browser and extracted 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 carriers for variants that overlapped genes on the minus strand. We shifted the signal by 70 nt to the right from the acceptor site and 70 nt to the left from the donor site. After the shift, the nucleosome signals for both the acceptor and donor were centered in the middle of an ideal exon of 140 nt in length, which is the median length of exons in the GENCODE v19 annotation. Finally, we averaged the signals for all shifts and smoothed the resulting signal by calculating the mean within an 11 nt window centered at each position.

[0655] To test for association, we selected random single SNVs that were at least 750 nt away from annotated exons and predicted by the model to have no effect on splicing (△ score < 0.01). We created 1000 such random samples of SNVs, each with the same number of SNVs as used for Fig.37H The number of SNVs in the splice site gain set (128 sites) was the same. For each random sample, we calculated a smoothed average signal as described above. Since random SNVs are not predicted to generate new exons, we centered the nucleosome signal for each SNV on the SNV itself and randomly shifted it 70 nt to the left or 70 nt to the right. We then Fig.37H The nucleosome signal at the middle base was compared to the signal obtained from 1000 simulations at that base. An empirical p-value was calculated as the fraction of simulation sets with a value greater than or equal to the middle of the observed values ​​for the splice site gain variable.

[0656] Robustness of the network to differences in exon density

[0657] To investigate the generalizability of the network's predictions, we evaluated SpliNet-10k on regions of varying exon density. We first divided the test set positions into five categories based on the number of canonical exons that appear in a 10,000 nucleotide window (5,000 nucleotides on each side) (Figure 54). To ensure that the exon counts were integer values ​​for each position, we used the exon start in the window as a surrogate. For each category, we calculated the top-k accuracy and the area under the precision-recall curve. The number of positions and k values ​​varied across categories (see table below for details).

[0658]

[0659] The robustness of the network for each of the 5 models

[0660] Training multiple models and averaging their predictions as output is a common strategy in machine learning to achieve better predictive performance, known as ensemble learning. Fig.53A In , we show the top-k accuracy and area under the precision-recall curve of the 5 SplitNet-10k models we trained to construct the population. The results clearly demonstrate the stability of the training process.

[0661] We also calculated the Pearson correlation between their predictions. Since most positions in the genome are not splice sites, the correlation between the predictions of most models will be close to 1, making the analysis meaningless. To overcome this problem, we only considered those positions in the test set that were assigned an acceptor or donor score greater than or equal to 0.01 by at least one model. 53,272 positions (approximately the same number of splice sites and non-splice sites) met this criterion. Fig.53B The results are summarized in . The very high Pearson correlations between the model predictions further illustrate their robustness.

[0662] exist Fig.53C In , we show the effect of the number of models used to build the ensemble on performance. The results show that as the number of models increases, performance improves, but at a diminishing return.

[0663] II. Analysis of the GTEx RNA Sequencing Database

[0664] Delta scores for single nucleotide variants

[0665] We quantified the splicing changes caused by single nucleotide variants as follows: We first calculated the acceptor and donor scores for 101 positions (50 positions on each side) surrounding the variant using the reference nucleotide. Assume that these scores are represented by the vector a ref and d refWe then use the alternative nucleotide and recalculate the acceptor and donor scores. These scores are represented by the vector a alt and d alt express.

[0666] We estimate the following four quantities:

[0667] Δ score receptor gain) = max(a alt -a ref )

[0668] Δ score receptor loss) = max(a ref -a alt )

[0669] Δ score donor gain) = max(d alt -d ref )

[0670] Δ score donor gain) = max(d ref -d alt )

[0671] The maximum of these four scores is called the variant's delta score.

[0672] Standards for quality control and variant filtering

[0673] We downloaded GTEx VCF and RNA sequencing data from dbGaP (study accession number phs000424.v6.p1; https: / / www.ncbi.nlm.nih.gov / projects / gap / cgi-bin / study.cgi?study_id=phs000424.v6.p1).

[0674] We evaluated the performance of SplitNet on the autosomal SNVs of the 4 most frequently occurring individuals in the GTEx cohort. In particular, a variant was considered if it met the following criteria in at least one individual A:

[0675] 1. Unfiltered variants (the FILTER field of the VCF is PASS).

[0676] 2. The variant is not marked as MULTI_ALLELIC in the INFO field of individual A's VCF, and the VCF contains a single allele in the ALT field.

[0677] 3. Individual A is heterozygous for the variant.

[0678] 4. The ratio of alt_depth / (alt_depth+ref_depth) is between 0.25 and 0.75, where alt_depth and ref_depth are the interpretations of the supporting alternative alleles and reference alleles in a single A, respectively.

[0679] 5. In the VCF of individual A, the total depth alt_depth+ref_depth is between 20 and 300.

[0680] 6. The variant overlaps a gene body region. The gene body is defined as the region between the transcription start and stop of the canonical transcript of GENCODE (V24lift37).

[0681] 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 called variants that occurred in a single individual singletons and variants that occurred in 2-4 individuals common variants. We did not evaluate variants that occurred in more than 5 individuals to prevent overlap with the training dataset.

[0682] RNA sequencing interpretation and comparison

[0683] We mapped the reads of GTEx samples to the hg19 reference using OLego (Wu et al., 2013) such that the edit distance between the query read and the reference did not exceed 4 (parameter -M 4). Note that OLego can operate completely de novo and does not require any gene annotation. Since OLego looks for the presence of splice motifs at the ends of spliced ​​reads, its alignments can be biased towards or away from the reference around SNVs that destroy or create splice sites, respectively. To eliminate this bias, we further created an alternative reference sequence for each GTEx individual by inserting all SNVs of the individual with the PASS filter into the hg19 reference. We mapped all samples from each individual to the alternative reference sequence of that individual using OLego with the same parameters. For each sample, we combined the two alignment sets (against the hg19 reference and the individual's alternative reference) by picking the best alignment for each read pair. To select the best alignment for read pair P, we used the following steps:

[0684] 1. If both reads P were not mapped in either alignment, we randomly selected an alternative alignment of hg19 or P.

[0685] 2. If P had more unmapped ends in one set of alignments (e.g., both ends of P mapped to the alternative reference, but only one end mapped to hg19), we selected the alignment in which both ends of P were mapped.

[0686] 3. If both ends of P are mapped in both sets of alignments, we select the alignment with the least total mismatches, or randomly select an alignment with the same number of mismatches.

[0687] Detection of splice junctions in aligned RNA sequencing data

[0688] We used leafcutter_cluster, a utility in the leafcutter package (Li et al., 2018) to detect and count splice junctions in each sample. We required a single spliced ​​read to support a junction and assumed a maximum intron length of 500 kb (parameter -m 1-1500000). To obtain a high-confidence junction set for training deep learning models, we compiled the union of all leafcutter junctions in all samples and then removed from consideration junctions that met any of the following criteria:

[0689] 1. Either end of a junction that overlaps with an ENCODE blacklist region (table WgEncodedCapabilityConsensuesInclused in hg19 in the UCSC genome browser) or a simple repeat sequence (from the simple repeat track in hg19 in the UCSC genome browser).

[0690] 2. Both ends of the junction are located on non-canonical exons (based on the canonical transcript of GENCODE version V24lift37).

[0691] 3. The two ends of the junction are located on different genes, or either end is located in a non-gene region.

[0692] 4. Lack of essential GT / AG dinucleotides at both ends.

[0693] Junctions present in more than 5 individuals were used to augment the list of GENCODE annotated splice junctions for analysis of variant prediction ( Fig.38A -G, 39A-C, 40A-E, and 41A-F). Links to files containing lists of splice junctions used for training the models are provided in the Key Resources table.

[0694] Although we augmented the training dataset with junctions detected by leafcutter, we noticed that leafcutter filtered many junctions with good support in the RNA-sequencing data despite the relaxation parameters used. This artificially reduced our validation rate. Therefore, for the GTEx RNA-sequencing validation analysis ( Fig.38A-G and 39A-C), we recalculated a set of junctions and junction counts directly from the RNA sequencing read data. We counted all non-repeated split-mapped reads where the junction was flanked by at least 10 and 5 nt of MAPQ. A read was allowed to span more than two exons, in which case the read was counted for each junction with at least 5 nt of mapped sequence on either side.

[0695] Definition of Private Contact

[0696] A node is considered private in A if it satisfies at least one of the following conditions:

[0697] 1. In at least one sample from A, the junction has at least 3 reads and has never been observed in any other individual.

[0698] 2. At least two organizations meet both of the following conditions:

[0699] a. The average read count of junctions in samples from individual A in the tissue is at least 10.

[0700] b. On average, the normalized read counts for individual A are at least twice as high as any other individual in the tissue. Here, the normalized read count for a junction in a sample is defined as the number of reads for a junction normalized by the total number of reads at all junctions for the corresponding gene.

[0701] Tissues from other individuals (not individual A) with fewer than 5 samples were ignored by the test.

[0702] Private Contacts Around Singleton SNV enrichment

[0703] If one end of a private junction was accurately annotated based on the GENCODE annotation, we considered it a candidate for acceptor or donor gain and searched for singleton SNVs (SNVs that occur in a single GTEx individual) that are unique to the same individual within 150 nt of the unannotated end. If both ends of a private junction were annotated, it was considered a candidate for a private exon skipping event if it skipped at least one but no more than 3 identical exons of the same gene according to the GENCODE annotation. We then searched for singleton SNVs within 150 nt starting from the end of each skipped exon. Private junctions that were missing both ends in the GENCODE exon annotation were ignored because a large proportion of them had alignment errors.

[0704] To calculate the enrichment of singleton SNVs around novel private acceptors or donors ( Fig.38B, bottom), we summed the counts of singleton SNVs at each position relative to the private junction. If the overlapping gene was on the negative strand, the relative position was flipped. We divided the SNVs into two groups: SNVs private to an individual with a private junction and SNVs private to a different individual. To smooth the resulting signal, we averaged the counts in a 7nt window centered at each position. We then calculated the ratio of the smoothed counts of the first group (private to the same individual) to the smoothed counts of the second group (private to another individual). For the new private exon skipping term ( Fig.38B , top), we followed a similar procedure and aggregated counts of singleton SNVs around the ends of the skipped exons.

[0705] Validation of model predictions in GTEx RNA sequencing data

[0706] For 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 predictions of the deep learning model for the reference and other alleles and calculated a delta score. We also obtained the locations where the model predicted aberrant (new or disrupted) junctions. We then sought to determine whether there was evidence in the RNA sequencing data to support aberrant splicing in individuals with the variant at the predicted location. In many cases, models can be used to predict multiple effects of the same variant, e.g., a variant that disrupts an annotated splice donor may also increase the use of suboptimal donors, such as Fig.45 . In this case, the model predicts donor losses at annotated splice sites as well as donor gains at suboptimal sites. However, for validation purposes, we only considered the effect with the highest predicted Δ score for each variant. Therefore, for each variant, we considered the predicted effects of splice site creation and splice site destruction separately. Note that junctions that appeared in less than five individuals were excluded during model training to avoid evaluating new junctions when training the model.

[0707] Validation of cryptic splicing mutations based on private splice site prediction

[0708] For each private variant predicted to cause new junction formation, we used the network to predict the location of the newly created aberrant splice junction and looked at the RNA sequencing data to verify that this new junction only occurred in individuals with the SNV and no other GTEx. Similarly, for variants predicted to cause splice site loss affecting the splice site of exon X, we looked for new exon skipping events from the previous canonical exon (upstream of X according to GENCODE annotation) to the next canonical exon (downstream of X) that only occurred in individuals with the variant and no other GTEx. We excluded expected losses if the splice site predicted by the model to be lost was not annotated in GENCODE or had never been observed in GTEx individuals without the variant. We also excluded expected gains if the expected splice site to be gained was already annotated in GENCODE. To extend this analysis to common variants (present in two to four individuals), we also verified new junctions that were present in at least half of the individuals with the variant and absent in all individuals without the variant.

[0709] We could validate 40% of expected high-scoring (Δ score 0.5) acceptor and donor gains using the requirement that the expected aberrant splicing event be private to the individual with the variant, but only 3.4% of predicted high-scoring losses and 5.6% of essential GT or AG disruptions (with a substitution-based false positive rate of <0.2% – see the “Estimating false positive rates” section). There are two reasons for the difference in gain and loss confirmation rates. First, unlike gains, exon skipping events are rarely completely private to the individual with the variant because exons are often skipped at low baseline levels that can be observed with sufficiently deep RNA sequencing. Second, splice site losses can have other effects besides increased exon skipping, such as increased intron retention or increased usage of other suboptimal splice sites. For these reasons, we did not rely solely on the use of new private junctions to validate the model’s predictions, but instead validated variants based on quantitative evidence for increasing or decreasing junction usage in individuals with the variant that were predicted to be affected.

[0710] Validation of predicted cryptic splicing mutations by quantitative criteria

[0711] For the joint j in sample s, we obtain the normalized joint number cjs:

[0712]

[0713] Here, r jsis the raw junction count for junction j in the sample, and the sum in the denominator replaces all other junctions between annotated recipients and donors of the same gene as j (using the annotations from GENCODE v19). The asinh transformation is defined as It is similar to the logarithmic transformation commonly used to transform RNA sequencing data (Lonsdale et al., 2013), but defines it to be zero, thus eliminating the need for pseudocounts, which would be severely distorted since many junctions, especially new junctions, have low or zero counts. For large values, the asinh transformation behaves like a logarithmic transformation, but for small values, it is close to linear. Therefore, it is often used in datasets with a large number of near-zero values ​​(such as RNA-seq or ChIP-seq datasets) to prevent a few large values ​​from dominating the signal (Azad et al., 2016; Herring et al., 2018; Hoffman et al., 2012; Kasowski et al., 2013; SEQC / MAQC-III Consortium, 2014). As described below, in the "Consideration Criteria for Validation" section, samples with a denominator less than 200 in equation (1) were excluded from all validation analyses to avoid numerical problems.

[0714] For each junction j predicted to be gained or lost due to a SNV in a group of individuals I, we calculated the following Z-score in each tissue t separately:

[0715]

[0716] Among them A t is a sample set of individuals from I in organization t, and U t is the set of samples from all other individuals in tissue t. Note that for the same individual and tissue, there may be multiple samples in the GTEx dataset. As mentioned before, c js is the count of junction j in the sample. For predicted losses, we also calculated similar z scores for junction k that skipped the putative exon:

[0717]

[0718] Note that the loss that causes jumps results in a relative decrease in lost connections and a relative increase in jumps. This proves that jt and z kt The difference in molecular recovery is such that both scores tend to be negative for actual splice site loss.

[0719] Finally, we calculated the median Z score for all considered tissues. For losses, we calculated the median of each Z score according to equations (2) and (3), respectively. A loss prediction for a recipient or donor was considered valid if any of the following conditions were true:

[0720] 1. The median z-score in equation (2), which quantifies the relative loss of junctions, is less than the fifth percentile of the corresponding values ​​in the permuted data (-1.46), and the median z-score in equation (3), which quantifies the relative change in jumps, is nonpositive (zero, negative, or absent if no jump junctions were observed in any individual). In other words, there is strong evidence of reduced use of the affected junctions, but no evidence of reduced jumps in the affected individuals.

[0721] 2. The median of the z-scores from equation (3) is less than the fifth percentile value (-0.74) of the corresponding values ​​in the permuted data, and the median of the z-scores from equation (3) is non-positive.

[0722] 3. The median z score in equation (2) is less than the first percentile of the corresponding value in the permuted data (-2.54).

[0723] 4. The median z score in equation (3) is less than the first percentile of the corresponding value in the permuted data (-4.08).

[0724] 5. Junction skipping of the affected exon is observed in at least half of the individuals with the variant but not in the others (as described above in the “Validation of Predicted Cryptic Splicing Mutations Based on Private Splice Junctions” section).

[0725] The permutations used to obtain the above critical values ​​are described in the "Estimating False Validation Rate" section.

[0726] Empirically, we observed that we need to apply more stringent validation criteria for losses than for gains because, as described in the “Validation of cryptic splicing mutations based on private splice junction predictions” section, losses induce larger confounding effects than gains. The probability of observing a new junction near a private SNV by chance is small, so even slight evidence of a junction is sufficient for validation. In contrast, most predicted losses weaken existing junctions, and such weakening is harder to detect than changes caused by gains and is more likely to be attributed to noise values ​​in the RNA sequencing data.

[0727] Inclusion criteria for validation analysis

[0728] To avoid calculating z-scores in cases with low counts or poor coverage, we used the following criteria to filter variants for validation analysis:

[0729] 1. Only when the sample expresses the gene (in equation (1) (∑ g r gs >200)), it is considered for the above z-score calculation.

[0730] 2. Tissues were not considered for loss or gain z-score calculations if the mean counts for loss or 'reference' junctions in individuals without the variant were less than 10, respectively. The 'reference' junction is the canonical junction used before new junctions based on the GENCODE annotation (see the Effect Size Calculation section for details). The intuition is that we should not attempt to validate splice variants that affect a junction that is not expressed in control individuals. Similarly, we should not attempt to validate splice variants if the control individuals do not adequately express the transcript spanning the affected site.

[0731] 3. In the case of predicted splice site losses, samples from individuals without the variant were considered only if the sample had at least 10 lost junction counts. In the case of expected recipient or donor gains, samples from control individuals were considered only if the control individual had a “reference” junction count of at least 10. The intuition is that even in tissues where the affected junctions are expressed at high average levels (i.e., pass criterion 2), different samples may have drastically different sequencing depths, so only control samples with sufficient expression should be included.

[0732] 4. A tissue was considered only if at least one sample met the above criteria for an individual with the variant and at least 5 samples met the criteria for at least 2 different control individuals.

[0733] Variants without tissues that met the above criteria were considered unconfirmable and were excluded when calculating the confirmation rate. For splicing gain variants, we filtered those variants that were already present at GENCODE-annotated splice sites. Similarly, for splicing loss variants, we only considered those variants that would reduce the score of an existing GENCODE-annotated splice site. Overall, 55% and 44% of high-scoring (Δ score ≥ 0.5) predicted gains and losses, respectively, were considered to be confirmable and used for validation analysis.

[0734] Estimating the false validation rate

[0735] To ensure that the above procedure had a reasonable true validation rate, we first looked at SNVs that appeared in 1-4 GTEx individuals and disrupted essential GT / AG dinucleotides. We believed that such mutations almost certainly affect splicing, and therefore their validation rate should be close to 100%. Of the above disruptions, 39% were ascertainable according to the above criteria, and of those that were ascertainable, the validation rate was 81%. To estimate the false validation rate, we permuted the individual labels of the SNV data. For each SNV that appeared in k GTEx individuals, we selected a random subset of k GTEx individuals and assigned the SNV to them. We created 10 such random datasets and repeated the validation on them. The validation rates in the permuted datasets were: 1.7-2.1% for gains and 4.3-6.9% for losses, with medians of 1.8% and 5.7%, respectively. The higher false validation rate for losses and the lower validation rate for disruptions is due to the difficulty in validating splice site losses, as emphasized in the section “Validation of cryptic splicing mutations based on private splice site prediction”.

[0736] Calculating effect sizes of cryptic splice variants in RNA sequencing data

[0737] We define the “effect size” of a variant as the score of affected gene transcripts that have a changed splicing pattern due to the variant (e.g., a switch to a new acceptor or donor). As a reference example of a predicted splicing gain variant, consider Fig.38C For a predicted acquired donor A, we first determine the junction (AC) that is closest to the annotated acceptor C. We then determine the “reference” junction (BC), where B≠A is the annotated donor closest to A. We then calculate the relative usage of the new junction (AC) to the reference junction (BC) in each sample s:

[0738]

[0739] Here, r (AC)s is the raw explained count of junctions (AC) in sample s. For each tissue, we calculated the change in junction (AC) usage between individuals with the variant and all other individuals:

[0740]

[0741] Among them 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 is calculated as the median of the above differences across all tissues. The calculation is similar in the case of gain of receptors or in the case where the splice site generating the variant is intronic. A simplified form of the effect size calculation (assuming a single sample from individuals with and without the variant) is Fig.38C Shown in.

[0742] For the predicted losses, we first calculated the fraction of transcripts that skipped the affected exon. Fig.45 To predict the loss of donor C, we determined the junction (CE) of the next downstream annotated exon and the junction (AB) from the upstream exon to one of the putatively affected exons. We quantified the proportion of transcripts that skipped the affected exon as follows:

[0743]

[0744] As for the gain, we calculated the change in jump score between samples of individuals with the variant and samples of individuals without the variant:

[0745]

[0746] The fraction of skipped transcripts calculated as above cannot fully capture the effects of acceptor or donor loss, as such disruptions may also result in increased levels of intron retention or use of suboptimal splice sites. To account for some of these effects, we also calculated the usage of the missed junction (CE) relative to other junctions with the same acceptor E:

[0747]

[0748] Here, ∑r (·E)s is the sum of all junctions between any (annotated or novel) acceptor and donor E. This includes affected junctions (CE), skipped junctions (AE), and potential junctions from other suboptimal donors that compensated for the loss of C, such as Fig.45 We then calculated the change in relative usage of the affected joints:

[0749]

[0750] Note that unlike (5) and (7), which measure the increase in usage of junctions gained or skipped in individuals with the variant, in (9) we measure the decrease in usage of junctions lost, thus recovering the difference between the two components. For each tissue, the effect size is calculated as the maximum of (7) and (9). As for gains, the final effect size of the variant is the median effect size across the tissue.

[0751] Inclusion criteria for effect size analysis

[0752] Effect sizes for variants were considered for calculation only if the variant was deemed validated based on the criteria described in the previous section. To avoid calculating scores for very small aberrant transcripts, we only considered samples where both aberrant and reference junction counts were at least 10. Because most cryptic splice variants are in introns, effect sizes cannot be calculated directly by counting the number of reference and alternative reads that overlap with the variant. Therefore, effect sizes for losses were calculated indirectly based on the reduction in relative usage of normal splice junctions. For effect sizes for novel junction gains, aberrant transcripts may be subject to nonsense-mediated degradation, thereby attenuating the observed effect sizes. Despite the limitations of these measures, we observed a consistent trend toward smaller effect sizes for lower-scoring cryptic splice variants in both gain and loss events.

[0753] Expected effect size for fully penetrant heterozygous private SNVs

[0754] For a fully penetrant splice site creating variant that causes all transcripts of the variant haplotype in individuals carrying the variant to switch to the new junction, and assuming that the new junction does not occur in control individuals, then according to equation (5), the expected effect size would be 0.5.

[0755] Similarly, if a heterozygous SNV causes a new exon skipping event, and all transcripts of the affected haplotype switch to the skipping junction, the expected effect size in equation (7) is 0.5. If all transcripts from individuals with the variant switch to different junctions (skipping junctions or another compensation point), then in equation (8), the ratio of samples from individuals with the variant is 0.5, and the ratio from other samples is 1, so the difference in equation (9) is 0.5. It is assumed that there is no skipping or other junctions of receptor E in individuals without the variant. It is also assumed that splice site destruction does not trigger intron retention. In practice, at least low levels of intron retention are often associated with splice site destruction. In addition, even in the absence of alternative splicing variants, exon skipping is also very extensive. This explains why even for variants that destroy essential GT / AG dinucleotides, the measured effect sizes are all below 0.5.

[0756] For an expected effect size of 0.5 for a fully penetrant heterozygous variant, we also assumed that the variant did not trigger nonsense-mediated degradation (NMD). In the presence of NMD, both the numerator and denominator of equations (4), (6), and (8) would decrease, thereby reducing the observed effect size.

[0757] Fraction of transcripts subjected to nonsense-mediated degradation

[0758] for Fig.38C , since the variant is an exon, we can count the number of reads that span the variant and have the reference or alternative allele ("Ref(unspliced)" and "Alt(unspliced)", respectively). We also count the number of reads that are spliced ​​at a new splice site and the number of reads that may have an alternative allele ("Alt(new junction)"). Fig.38C In the example of, and other observed cases, we observed that the total number of reads from haplotypes with alternate alleles (the sum of “Alt(no splicing)” and “Alt(new junction)”) was less than the number of reads from the reference allele (“Ref(no splicing)”). Since we believe we have eliminated the reference bias in the mapping of reads, by mapping to reference and alternative haplotypes, and assuming that the number of reads is proportional to the number of transcripts for each allele, we expect the reference allele to account for half of the number of reads at the variant locus. We assume that the “missing” alternative allele reads correspond to transcripts from the alternative allele haplotype that are spliced ​​at the new junction and degraded by nonsense-mediated degradation (NMD). We call this group “Alt(NMD)”.

[0759] To determine whether the observed difference between the reference and alternative reads was significant, we calculated the probability of Alt(no splicing)+Alt(new junction) (or a small number) reads using the total number of trials with a success probability of 0.5 and Alt(no splicing)+Alt(new junction)+Ref(no splicing) under a binomial distribution. This is a conservative p-value because we underestimated the total number of "trials" by not counting transcripts that may have been degraded. Fig.38C The fraction of NMD transcripts was calculated as the number of “Alt(NMD)” reads out of the total number of reads spliced ​​at a new junction (Alt(NMD)+Alt(New Junction)).

[0760] Sensitivity of the network in detecting cryptic splicing junctions

[0761] To evaluate the sensitivity of the SpliceNet model ( Fig.38F ), we used SNVs within 20 nt of an affected splice site (i.e., a novel or disrupted acceptor or donor) and that did not overlap with an annotated exonic essential GT / AG dinucleotide and had an estimated effect size of at least 0.3 (see the section “Effect size calculation”). In all sensitivity plots, SNVs were defined as “nearby exonic” if they overlapped an annotated exon or were within 50 nt of an annotated exon boundary. All other SNVs were considered “deeply intronic”. Using a real dataset of strongly supported cryptic splice sites, we evaluated our model at varying Δ score thresholds and report the proportion of cryptic splice sites in the real dataset predicted by the model at that cutoff.

[0762] Comparison with existing splicing prediction models

[0763] We performed a head-to-head comparison of SpliceNet-10k, MaxEntScan (Yeo and Burge, 2004), GeneSplicer (Pertea et al., 2001), and NNSplice (Reese et al., 1997) on various metrics. We downloaded MaxEntScan and GeneSplicer software from http: / / genes.mit.edu / burgelab / maxent / download / and http: / / www.cs.jhu.edu / ~genomics / GeneSplicer / , respectively. NNSplice was not available as downloadable software, so we downloaded the training and test sets from http: / / www.fruitfly.org / data / seq_tools / datasets / Human / GENIE_96 / splicesets / and trained the model with the best performing architecture (Reese et al., 1997). As a sanity check, we replicated the test set metrics reported in (Reese et al., 1997).

[0764] To evaluate the Topk accuracy and area under the precision-recall curve of these algorithms, we scored each algorithm for all positions in the test set genes and lincRNAs ( Fig.37D ).

[0765] MaxEntScan and GeneSplicer outputs correspond to log odds ratios, while NNSplice and SpliceNet-10k outputs correspond to probabilities. To ensure that we give MaxEntScan and GeneSplicer the best chance, we use their default outputs and their transformed outputs (we first transform their outputs to make them correspond to probabilities) for the Δ score calculation. More precisely, the default output of MaxEntScan corresponds to

[0766]

[0767] In the process of transformation After that, it corresponds to the desired number. We compiled the GeneSplicer software twice, once with the RETURN_TRUE_PROB flag set to 0 and once with it set to 1. We selected the output strategy with the highest validation rate for RNA-Seq data (MaxEntScan: transformed output, GeneSplicer: default output).

[0768] In order to compare the confirmation rate and sensitivity of various algorithms ( Figure 38G ), we found a critical value at which all algorithms predicted the same number of gains and losses across the genome. That is, for each critical value of the SpliceNet-10k score, we found that each competing algorithm made the same number of gain predictions and the same number of loss predictions as SpiceNet-10k at that critical value. The selected critical values ​​are shown in Table S2.

[0769] Comparison of singleton variant and common variant prediction

[0770] We analyzed singleton SNVs and SNVs that appeared in 2-4 GTEx individuals. Fig.46A , 46B, and 46C) were validated and analyzed for sensitivity (as described in “Sensitivity analysis” and “Validation of model predictions”). To test whether there was a significant difference in confirmation rates between singleton and common variants, we performed Fisher’s exact test comparing confirmation rates for each Δ score group (0.2–0.35, 0.35–0.5, 0.5–0.8, 0.8–1) and each predicted effect (gain or loss of acceptor or donor). All P values ​​were greater than 0.05 after Bonferroni correction for 16 tests. We similarly tested sensitivity fo...

Claims

1. A computer-implemented method for predicting the likelihood of splicing sites in a pre-mRNA genomic sequence, the method comprising: A convolutional neural network (CNN) is trained on training examples of the pre-mRNA nucleotide sequence, the training examples comprising at least 50,000 training examples of donor splicing sites, at least 50,000 training examples of acceptor splicing sites, and at least 100,000 training examples of non-splicing sites, wherein the training comprises: Input training examples of one-hot encoded pre-mRNA nucleotide sequences, wherein each nucleotide sequence comprises at least one target nucleotide flanked by at least 100 upstream context nucleotides and at least 100 downstream context nucleotides; and Adjusting the parameters of the filter in the CNN by back-propagation to predict a score for the likelihood that each target nucleotide in the pre-mRNA nucleotide sequence is a donor splice site, an acceptor splice site, or a non-splicing site; The trained CNN is thus configured to accept as input a pre-mRNA nucleotide sequence that is one-hot encoded and includes at least one target nucleotide flanked by at least 100 nucleotides on each side, and to score the likelihood of each target nucleotide being a donor splice site, an acceptor splice site, or a non-splicing site.

2. The computer-implemented method of claim 1 , wherein: Each target nucleotide is flanked by 200 upstream context nucleotides and 200 downstream context nucleotides.

3. The computer-implemented method of claim 1 , further comprising training the CNN on at least 150,000 training examples of donor splice sites, 150,000 training examples of acceptor splice sites, and 1,000,000 training examples of non-splice sites.

4. The computer-implemented method of claim 1 , wherein: The CNN is parameterized by the number of convolutional layers, the number of convolutional filters, and the number of subsampling layers.

5. The computer-implemented method of claim 1 , wherein: The CNN includes one or more fully connected layers and a terminal classification layer.

6. The computer-implemented method of claim 1 , wherein: The CNN includes a dimensionality change layer that reshapes the spatial and feature dimensions of the previous input.

7. The computer-implemented method of claim 1 , wherein: The CNN is parameterized by the number of residual blocks, the number of skip connections, and the number of residual connections.

8. The computer-implemented method of claim 7, wherein: Each residual block includes at least one batch normalization layer, at least one rectified linear unit (ReLU) layer, at least one dimension change layer, and at least one residual connection.

9. The computer-implemented method of claim 7, wherein: Each residual block consists of two batch normalization layers, two rectified linear unit (ReLU) nonlinear layers, two dimension change layers, and a residual connection.

10. The computer-implemented method of claim 1, wherein: The CNN evaluates the training instances in batches during an epoch.

11. The computer-implemented method of claim 1 , wherein: The training instances are randomly sampled into batches, where each batch has a predetermined batch size.

12. The computer-implemented method of claim 10, wherein: The CNN iterates the evaluation of the training examples for at least 10 epochs.

13. The computer-implemented method of claim 1, wherein: The input includes a target nucleotide sequence having two adjacent target nucleotides.

14. The computer-implemented method of claim 13, wherein: The two adjacent target nucleotides are adenine, abbreviated as A, and guanine, abbreviated as G.

15. The computer-implemented method of claim 13, wherein: The two adjacent target nucleotides are guanine, abbreviated as G, and uracil, abbreviated as D.

16. A system for predicting the likelihood of splice sites in a pre-mRNA genomic sequence, comprising one or more processors coupled to a memory, the memory being loaded with computer instructions that, when executed on the one or more processors, perform the following actions: A convolutional neural network (CNN) is trained on training examples of the pre-mRNA nucleotide sequence, the training examples comprising at least 50,000 training examples of donor splicing sites, at least 50,000 training examples of acceptor splicing sites, and at least 100,000 training examples of non-splicing sites, wherein the training comprises: Input training examples of one-hot encoded pre-mRNA nucleotide sequences, wherein each nucleotide sequence comprises at least one target nucleotide flanked by at least 100 upstream context nucleotides and at least 100 downstream context nucleotides; and Adjusting the parameters of the filter in the CNN by back-propagation to predict a score for the likelihood that each target nucleotide in the pre-mRNA nucleotide sequence is a donor splice site, an acceptor splice site, or a non-splicing site; The trained CNN is thus configured to accept as input a pre-mRNA nucleotide sequence that is one-hot encoded and includes at least one target nucleotide flanked by at least 100 nucleotides on each side, and to score the likelihood of each target nucleotide being a donor splice site, an acceptor splice site, or a non-splicing site.

17. The system of claim 16, wherein: The CNN was trained on 150,000 training instances of donor splice sites, 150,000 training instances of acceptor splice sites, and 1,000,000 training instances of non-splice sites.

18. A system for abnormal splicing determination, comprising a plurality of processors operating in parallel coupled to a memory, the memory being loaded with computer instructions for determining abnormal splicing, the system comprising: A trained convolutional neural network (CNN) for processing pre-mRNA sequences running on the plurality of processors, comprising: Input module, which receiving a variant sequence of at least 101 target nucleotides flanked on each side by at least 400 nucleotides; and accessing a reference sequence comprising at least 101 target nucleotides flanked on each side by at least 400 nucleotides; a layer that convolves target nucleotides in the variant sequence and the reference sequence and assigns a splice site score to the target nucleotides to give the likelihood that each target nucleotide is a donor splice site, an acceptor splice site, or a non-splice site; and An output module determines whether a variant in the variant sequence causes abnormal splicing of any target nucleotide in the variant sequence and is therefore pathogenic based on the positional differences in the splice site scores of the target nucleotides in the variant sequence and the reference sequence.

19. The system of claim 18, wherein: The output module is configured to determine that a variant in the variant sequence causes abnormal splicing when the difference in the splice site score for at least one target nucleotide position is higher than a predetermined threshold.

20. The system of claim 18, wherein: The maximum difference in splice site scores at any position between the target nucleotides is below a predetermined threshold, further configured to classify the variant sequence as not causing aberrant splicing and therefore being benign.

21. The system according to any one of claims 19-20, wherein: The predetermined threshold is determined by: For multiple candidate thresholds, processing a first set of reference and variant sequence pairs generated from benign common variants to generate a first set of aberrant splicing detections; as well as processing a second set of reference and variant sequence pairs generated by the pathogenic rare variant to generate a second set of aberrant splicing assays; as well as At least one threshold is selected for use by the output module that maximizes the count of aberrant splicing detections in the second set of reference and variant sequence pairs and minimizes the count of aberrant splicing detections in the first set of reference and variant sequence pairs.

22. The system of claim 18, further configured to identify variants causing autism spectrum disorder (ASD).

23. The system of claim 18, further configured to identify variants causing developmental delay disorder DDD.

24. The system of claim 18, wherein: Each sequence of at least 101 consecutive target nucleotides is flanked on each side by an extended context of at least 1000 nucleotides.

25. The system of claim 18, wherein: The splice site scores of the target nucleotides in the reference sequence are encoded in a first output of the trained CNN, and the splice site scores of the target nucleotides in the variant sequences are encoded in a second output of the trained CNN.

26. The system of claim 25, wherein: The first output is encoded as a first 101x3 matrix, and the second output is encoded as a second 101x3 matrix.

27. The system of claim 26, wherein: Each row in the first 101×3 matrix uniquely represents a splice site score to give the likelihood that a target nucleotide in the reference sequence is a donor splice site, an acceptor splice site, or a non-splice site.

28. The system of claim 27, wherein: Each row in the second 101 x 3 matrix uniquely represents a splice site score to give the likelihood that a target nucleotide in the variant sequence is a donor splice site, an acceptor splice site, or a non-splice site.

29. The system according to any one of claims 26-28, wherein: The splice site scores in each row of the first 101×3 matrix and the second 101×3 matrix are exponentially normalized and summed to unity.

30. The system according to any one of claims 26-28, wherein: The output module performs a row-to-row comparison of the first 101×3 matrix and the second 101×3 matrix and determines changes in the splice site score distribution row by row.

31. The system of claim 30, wherein: For at least one instance of the row-to-row comparison, the corresponding change in the distribution is above a predetermined threshold, further configured to classify the variant as causing aberrant splicing and therefore being pathogenic.

32. The system of claim 18, wherein: The reference sequence and the variant sequences are one-hot encoded.

33. A method for detecting abnormal splicing, the method comprising: Accepting as input to a trained convolutional neural network (CNN) a pre-mRNA nucleotide sequence comprising at least 801 nucleotides, the pre-mRNA nucleotide sequence comprising at least 101 target nucleotides flanked on each side by at least 400 nucleotides; Processing a reference sequence and a variant sequence through the trained CNN to generate a splice site score to give a likelihood that each target nucleotide in the reference sequence and the variant sequence is a donor splice site, an acceptor splice site, or a non-splicing site, wherein the reference sequence and the variant sequence each have at least 101 target nucleotides flanked by at least 400 nucleotides on each side, upstream and downstream; and Based on the difference in the splice site scores of the target nucleotides in the reference sequence and the variant sequence, it is determined whether the variant in the variant sequence causes aberrant splicing and is therefore pathogenic.

34. A non-transitory computer readable storage medium having printed thereon computer program instructions for detecting abnormal splicing, wherein when the instructions are executed on a processor, a method comprising: Accepting as input to a trained convolutional neural network (CNN) pre-mRNA nucleotides in an input sequence comprising at least 801 nucleotides, the input sequence comprising at least 101 target nucleotides flanked on each side by at least 400 nucleotides; Processing a reference sequence and a variant sequence through the trained CNN to generate a splice site score to give a likelihood that each target nucleotide in the reference sequence and the variant sequence is a donor splice site, an acceptor splice site, or a non-splicing site, wherein the reference sequence and the variant sequence each have at least 101 target nucleotides and each target nucleotide is flanked on each side by at least 400 nucleotides; as well as Based on the difference in the splice site scores of the target nucleotides in the reference sequence and the variant sequence, it is determined whether the variant in the variant sequence causes aberrant splicing and is therefore pathogenic.

35. A non-transitory computer-readable storage medium having computer program instructions stored thereon, which when executed on a processor, implement the computer-implemented method according to any one of claims 1-15.

Images