Covariate adjustments for temporal data from phenotypic measures for different drug use patterns

Abstract

A computer-implemented method for predicting phenotypic shifts in response to the use of multiple drugs for multiple phenotypes in a cohort of individuals with multiple confounding factors, the cohort of individuals having associated phenotypic measures, covariate measures, and drug use patterns for two separate time points, the phenotypic measures for the first and second time points being covariate-adjusted and drug use-adjusted through the use of biostatistics.

Description

[Technical field]

[0001] (Priority Application) This application claims the benefit and priority of the following: U.S. Provisional Patent Application No. 63 / 294,813, entitled “PERIODIC MASK PATTERN FOR REVELATION LANGUAGE MODELS,” filed on December 29, 2021 (Attorney Docket No. ILLM1063-1 / IP-2296-PRV); U.S. Provisional Patent Application No. 63 / 294,816, entitled “CLASSIFYING MILLIONS OF VARIANTS OF UNCERTAIN SIGNIFICANCE USING PRIMATE SEQUENCING AND DEEP LEARNING,” filed on December 29, 2021 (Attorney Docket No. ILLM1064-1 / IP-2297-PRV); U.S. Provisional Patent Application No. 63 / 294,820, entitled “IDENTIFYING GENES WITH DIFFERENTIAL SELECTIVE CONSTRAINT BETWEEN HUMANS AND NON-HUMAN PRIMATES,” filed on December 29, 2021 (Attorney Docket No. ILLM1065-1 / IP-2298-PRV); U.S. Nonprovisional Patent Application No. 63 / 294,827, entitled “DEEP LEARNING NETWORK FOR EVOLUTIONARY CONSERVATION,” filed on December 29, 2021 (Attorney Docket No. ILLM1066-1 / IP-2299-PRV); U.S. Provisional Patent Application No. 63 / 294,828, entitled “INTER-MODEL PREDICTION SCORE RECALIBRATION,” filed on December 29, 2021 (Attorney Docket No. ILLM1067-1 / IP-2301-PRV); U.S. Provisional Patent Application No. 63 / 294,830, entitled “SPECIES-DIFFERENTIABLE EVOLUTIONARY PROFILES,” filed on December 29, 2021 (Attorney Docket No. ILLM1068-1 / IP-2302-PRV); U.S. Provisional Patent Application No. 63 / 351,283, entitled “OPTIMIZED BURDEN TEST BASED ON NESTED T-TESTS THAT MAXIMIZE SEPARATION BETWEEN CARRIERS AND NON-CARRIERS,” filed on June 10, 2022 (Attorney Docket No. ILLM1070-1 / IP-2368-PRV); U.S. Provisional Patent Application No. 63 / 351,299, entitled “RARE VARIANT POLYGENIC RISK SCORES,” filed on June 10, 2022 (Attorney Docket No. ILLM1071-1 / IP-2378-PRV), and U.S. Provisional Patent Application No. 63 / 351,317, entitled “COVARIATE CORRECTION INCLUDING DRUG USE FROM TEMPORAL DATA,” filed on June 10, 2022 (Attorney Docket No. ILLM1073-1 / IP-2387-PRV).

[0002] The priority applications are incorporated by reference as if fully set forth herein.

[0003] FIELD OF THEINVENTION The disclosed technology relates to artificial intelligence based computers and digital data processing systems and corresponding data processing methods and products for mimicking intelligence (i.e., knowledge-based systems, inference systems, and knowledge acquisition systems), including systems for reasoning with uncertainty (e.g., fuzzy logic systems), adaptive systems, machine learning systems, and artificial neural networks. In particular, the disclosed technology relates to the use of deep convolutional neural networks to analyze ordered data.

[0004] (Related Applications) This application is related to U.S. Nonprovisional Patent Application No. 17 / 968,285 (Attorney Docket No. ILLM1070-2 / IP-2368-US), entitled “OPTIMIZED BURDEN TEST BASED ON NESTED T-TESTS THAT MAXIMIZE SEPARATION BETWEEN CARRIERS AND NON-CARRIERS,” filed on October 18, 2022. The related application is hereby incorporated by reference for all purposes.

[0005] This application is related to U.S. Nonprovisional Patent Application No. 17 / 968,723 (Attorney Docket No. ILLM1071-2 / IP-2378-US), entitled "RARE VARIANT POLYGENIC RISK SCORES," filed on October 18, 2022. The related application is hereby incorporated by reference for all purposes.

[0006] (Incorporated by reference) The following are incorporated by reference for all purposes as if fully set forth herein: Sundaram,L.et al.Predicting the clinical impact of human mutation with deep neural networks.Nat.Genet.50,1161-1170(2018), Jaganathan,K.et al.Predicting splicing from primary sequence with deep learning.Cell 176,535-548(2019), U.S. Patent Application No. 62 / 573,144, entitled “TRAINING A DEEP PATHOGENICITY CLASSIFIER USING LARGE-SCALE BENIGN TRAINING DATA,” filed on October 16, 2017 (Attorney Docket No. ILLM1000-1 / IP-1611-PRV); U.S. Patent Application No. 62 / 573,149, entitled “PATHOGENICITY CLASSIFIER BASED ON DEEP CONVOLUTIONAL NEURAL NETWORKS (CNNs),” filed on October 16, 2017 (Attorney Docket No. ILLM1000-2 / IP-1612-PRV); U.S. Patent Application No. 62 / 573,153, entitled “DEEP SEMI-SUPERVISED LEARNING THAT GENERATES LARGE-SCALE PATHOGENIC TRAINING DATA,” filed on October 16, 2017 (Attorney Docket No. ILLM1000-3 / IP-1613-PRV); U.S. Patent Application No. 62 / 582,898, entitled “PATHOGENICITY CLASSIFICATION OF GENOMIC DATA USING DEEP CONVOLUTIONAL NEURAL NETWORKS (CNNs),” filed on November 7, 2017 (Attorney Docket No. ILLM 1000-4 / IP-1618-PRV); U.S. Patent Application No. 16 / 160,903, entitled “DEEP LEARNING-BASED TECHNIQUES FOR TRAINING DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed on October 15, 2018 (Attorney Docket No. ILLM1000-5 / IP-1611-US); U.S. Patent Application Serial No. 16 / 160,986, entitled "DEEP CONVOLUTIONAL NEURAL NETWORKS FOR VARIANT CLASSIFICATION," filed on October 15, 2018 (Attorney Docket No. ILLM1000-6 / IP-1612-US); U.S. Patent Application Serial No. 16 / 160,968, entitled “SEMI-SUPERVISED LEARNING FOR TRAINING AN ENSEMBLE OF DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed on October 15, 2018 (Attorney Docket No. ILLM1000-7 / IP-1613-US); U.S. patent application Ser. No. 16 / 160,978, entitled “DEEP LEARNING-BASED SPLICE SITE CLASSIFICATION,” filed on October 15, 2018 (Attorney Docket No. ILLM1001-4 / IP-1680-US); U.S. Patent Application No. 16 / 407,149, entitled “DEEP LEARNING-BASED TECHNIQUES FOR PRE-TRAINING DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed on May 8, 2019 (Attorney Docket No. ILLM1010-1 / IP-1734-US); U.S. Patent Application No. 17 / 232,056, entitled “DEEP CONVOLUTIONAL NEURAL NETWORKS TO PREDICT VARIANT PATHOGENICITY USING THREE-DIMENSIONAL (3D) PROTEIN STRUCTURES,” filed on April 15, 2021 (Attorney Docket No. ILLM1037-2 / IP-2051-US); U.S. Patent Application No. 63 / 175,495, entitled “MULTI-CHANNEL PROTEIN VOXELIZATION TO PREDICT VARIANT PATHOGENICITY USING DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed on April 15, 2021 (Attorney Docket No. ILLM1047-1 / IP-2142-PRV); U.S. Patent Application No. 63 / 175,767, entitled “EFFICIENT VOXELIZATION FOR DEEP LEARNING,” filed on April 16, 2021 (Attorney Docket No. ILLM1048-1 / IP-2143-PRV); U.S. Patent Application No. 17 / 468,411, entitled “ARTIFICIAL INTELLIGENCE-BASED ANALYSIS OF PROTEIN THREE-DIMENSIONAL (3D) STRUCTURES,” filed on September 7, 2021 (Attorney Docket No. ILLM1037-3 / IP-2051A-US); U.S. Provisional Patent Application No. 63 / 253,122, entitled “PROTEIN STRUCTURE-BASED PROTEIN LANGUAGE MODELS,” filed on October 6, 2021 (Attorney Docket No. ILLM1050-1 / IP-2164-PRV); U.S. Provisional Patent Application No. 63 / 281,579, entitled “PREDICTING VARIANT PATHOGENICITY FROM EVOLUTIONARY CONSERVATION USING THREE-DIMENSIONAL (3D) PROTEIN STRUCTURE VOXELS,” filed on November 19, 2021 (Attorney Docket No. ILLM1060-1 / IP-2270-PRV); U.S. Provisional Patent Application No. 63 / 281,592, entitled “COMBINED AND TRANSFER LEARNING OF A VARIANT PATHOGENICITY PREDICTOR USING GAPED AND NON-GAPED PROTEIN SAMPLES,” filed on November 19, 2021 (Attorney Docket No. ILLM1061-1 / IP-2271-PRV); U.S. Provisional Patent Application No. 63 / 294,813, entitled “PERIODIC MASK PATTERN FOR REVELATION LANGUAGE MODELS,” filed on December 29, 2021 (Attorney Docket No. ILLM1063-1 / IP-2296-PRV); U.S. Provisional Patent Application No. 63 / 294,816, entitled “CLASSIFYING MILLIONS OF VARIANTS OF UNCERTAIN SIGNIFICANCE USING PRIMATE SEQUENCING AND DEEP LEARNING,” filed on December 29, 2021 (Attorney Docket No. ILLM1064-1 / IP-2297-PRV); U.S. Provisional Patent Application No. 63 / 294,820, entitled “IDENTIFYING GENES WITH DIFFERENTIAL SELECTIVE CONSTRAINT BETWEEN HUMANS AND NON-HUMAN PRIMATES,” filed on December 29, 2021 (Attorney Docket No. ILLM1065-1 / IP-2298-PRV); U.S. Nonprovisional Patent Application No. 63 / 294,827, entitled “DEEP LEARNING NETWORK FOR EVOLUTIONARY CONSERVATION,” filed on December 29, 2021 (Attorney Docket No. ILLM1066-1 / IP-2299-PRV); U.S. Provisional Patent Application No. 63 / 294,828, entitled “INTER-MODEL PREDICTION SCORE RECALIBRATION,” filed on December 29, 2021 (Attorney Docket No. ILLM1067-1 / IP-2301-PRV), and U.S. Provisional Patent Application No. 63 / 294,830, entitled “SPECIES-DIFFERENTIABLE EVOLUTIONARY PROFILES,” filed on December 29, 2021 (Attorney Docket No. ILLM1068-1 / IP-2302-PRV). [Background technology]

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

[0008] Genomics in the broad sense, also called functional genomics, aims to characterize the function of all genomic elements of an organism by using genome-scale assays such as genome sequencing, transcriptome profiling, and proteomics. Genomics has emerged as a data-driven science and operates not by testing preconceived models and hypotheses, but by discovering novel properties from the exploration of genome-scale data. Applications of genomics include finding associations between genotypes and phenotypes, discovering biomarkers for patient stratification, predicting gene functions, and mapping biochemically active genomic regions and residues such as transcriptional enhancers and single nucleotide polymorphisms (SNPs) using biostatistical analysis.

[0009] Genomics data is too large and complex to be mined solely by visual inspection of pairwise correlations. Instead, analytical tools are needed to support the discovery of unexpected relationships, derive novel hypotheses and models, and make predictions. Unlike some algorithms in which assumptions and domain expertise are hard-coded, machine learning algorithms are designed to automatically detect patterns in data. Thus, machine learning algorithms are well suited for data-driven science, especially genomics. However, the performance of machine learning algorithms can be highly dependent on how the data is represented, i.e., how each variable (also called a feature) is calculated. For example, to classify tumors as malignant or benign from a fluorescent microscopy image, a pre-processing algorithm can detect cells, identify cell types, and generate a list of cell counts for each cell type.

[0010] The machine learning model can take estimated cell counts, which are an example of hand-designed features, as input features to classify tumors. The central problem is that classification performance is highly dependent on the quality and relevance of these features. For example, relevant visual features such as cell morphology, distance between cells or localization within an organ are not captured in the cell counts, and this incomplete representation of the data can reduce classification accuracy.

[0011] Deep learning, a sub-discipline of machine learning, addresses this problem by embedding feature computation into the machine learning model itself, generating an end-to-end model. This result has been achieved through the development of deep neural networks, which are machine learning models that involve successive primitive operations that compute increasingly complex features by taking the results of previous operations as input. Deep neural networks can improve prediction accuracy by discovering relevant features of high complexity, such as cell morphology and spatial organization of cells in the above example. The construction and training of deep neural networks has been made possible by the explosion of data, advances in algorithms, and substantial increases in computing power, especially through the use of graphical processing units (GPUs).

[0012] The goal of supervised learning is to obtain a model that takes features as input and returns a prediction of a so-called target variable. An example of a supervised learning problem is the problem of predicting whether an intron will be spliced ​​or not (target), given features on the RNA such as the presence or absence of canonical splice site sequences, the location of splicing branch points, or the intron length. Training a machine learning model refers to learning its parameters, which generally involves minimizing a loss function on the training data with the goal of making accurate predictions on unknown data.

[0013] For many supervised learning problems in computational biology, the input data can be represented as a table with multiple columns or features, each of which contains numerical or categorical data that is potentially useful for making predictions. Some input data are naturally represented as tabular features (e.g., temperature or time), while other input data must first be transformed using a process called feature extraction to fit into a tabular representation (e.g., converting deoxyribonucleic acid (DNA) sequences to k-mer counts). For intron-splicing prediction problems, the presence or absence of canonical splice site sequences, the location of splicing branch points, and intron lengths can be preprocessed features collected in tabular form. Tabular data is the norm for a wide range of supervised machine learning models, ranging from simple linear models such as logistic regression to more flexible nonlinear models such as neural networks and many others.

[0014] Logistic regression is a binary classifier, i.e., a supervised learning model that predicts a binary target variable. Specifically, logistic regression predicts the probability of a positive class by calculating a weighted sum of input features that are mapped to the [0,1] interval using a sigmoid function, a type of activation function. The parameters of logistic regression, or other linear classifiers that use different activation functions, are the weights in the weighted sum. Linear classifiers fail when the weighted sum of input features, for example, cannot sufficiently distinguish the classes of whether an intron is spliced ​​or not. To improve prediction performance, new input features can be added manually by transforming or combining existing features in new ways, for example, by taking exponentiations or pairwise products.

[0015] Neural networks automatically learn these nonlinear feature transformations using hidden layers, each of which can be thought of as multiple linear models with outputs transformed by a nonlinear activation function such as a sigmoid function or the more common rectified-linear unit (ReLU). Together, these layers organize the input features into related complex patterns, facilitating the task of distinguishing between two classes.

[0016] Deep neural networks use many hidden layers, and when each neuron receives input from all neurons in the previous layer, the layer is said to be fully connected. Neural networks are generally trained using stochastic gradient descent, an algorithm suitable for training models on very large data sets. Implementation of neural networks using modern deep learning frameworks allows rapid prototyping with different architectures and data sets. Fully connected neural networks can be used for several genomics applications, including predicting the proportion of exons spliced ​​for a given sequence from sequence features such as the presence of splice factor binding motifs or sequence conservation, prioritizing potentially disease-causing genetic variants, and predicting cis-regulatory elements in a given genomic region using features such as chromatin marks, gene expression, and evolutionary conservation.

[0017] For effective prediction, local dependencies in spatial and longitudinal data must be considered. For example, shuffling of DNA sequences or image pixels severely disrupts information patterns. These local dependencies set spatial or longitudinal data apart from tabular data, where feature ordering is arbitrary. Consider the problem of classifying genomic regions as bound vs. unbound by a particular transcription factor, where binding regions are defined as high-confidence binding events in chromatin immunoprecipitation followed by sequencing (ChIP-seq) data. Transcription factors bind to DNA by recognizing sequence motifs. Fully connected layers based on sequence-derived features such as the number of k-mer instances in a sequence or position weight matrix (PWM) matches can be used for this task. Such models can generalize well to sequences with the same motif located at different positions, because k-mer or PWM instance frequencies are robust to shifting motifs in the sequence. However, they cannot recognize patterns where transcription factor binding depends on the combination of multiple motifs with distinct intervals. Furthermore, the number of possible k-mers grows exponentially with the k-mer length, which poses both conservation and overfitting challenges.

[0018] A convolutional layer is a special form of a fully connected layer in which the same fully connected layer is applied locally, for example within a 6 bp window, to all sequence positions. This approach can also be viewed as scanning the sequence using multiple PWMs, for example for the transcription factors GATA1 and TAL1. By using the same model parameters across positions, the total number of parameters is dramatically reduced and the network can detect motifs at positions not seen during training. Each convolutional layer scans the sequence with several filters by generating a scalar value at every position that quantizes the match between the filter and the sequence. As in a fully connected neural network, a nonlinear activation function (typically ReLU) is applied at each layer. A pooling operation is then applied, which aggregates the activations in successive bins across the position axis, typically taking the maximum or average activation for each channel. Pooling reduces the effective sequence length and coarsens the signal. Subsequent convolutional layers can construct the output of the previous layer and detect whether the GATA1 and TAL1 motifs were present within a certain distance range. Finally, the output of the convolutional layers can be used as input to a fully connected neural network to perform the final prediction task. Thus, different types of neural network layers (e.g., fully connected and convolutional layers) can be combined within a single neural network.

[0019] Convolutional neural networks (CNNs) can predict various molecular phenotypes based on DNA sequence alone. Applications include classification of transcription factor binding sites, as well as prediction of molecular phenotypes such as chromatin features, DNA contact maps, DNA methylation, gene expression, translation efficiency, RBP binding, and microRNA (miRNA) targets. In addition to predicting molecular phenotypes from sequences, convolutional neural networks can be applied to more technical tasks traditionally addressed by hand-designed bioinformatics pipelines. For example, convolutional neural networks can predict guide RNA specificity, denoise ChIP-seq, improve Hi-C data resolution, predict laboratory origin from DNA sequences, and call genetic variants. Convolutional neural networks have also been used to model long-range dependencies in genomes. Although interacting regulatory elements may be located far apart on unfolded linear DNA sequences, these elements are often proximal in real 3D chromatin conformations. Thus, modeling of molecular phenotypes from linear DNA sequences can be improved by allowing long-range dependencies, albeit with a crude approximation of chromatin, and allowing the model to implicitly learn aspects of 3D organization such as promoter-enhancer loops. This is achieved by using dilated convolutions with receptive fields of up to 32 kb. Dilated convolutions also allow splice sites to be predicted from sequences using receptive fields of 10 kb, thereby allowing integration of gene sequences over distances as long as a typical human intron (see Jaganathan, K. et al. Predicting splicing from primary sequence with deep learning. Cell 176, 535-548 (2019)).

[0020] Different types of neural networks can be characterized by their parameter sharing schemes. For example, fully connected layers have no parameter sharing, while convolutional layers impose translational invariance by applying the same filter at every position of their input. Recurrent neural networks (RNNs) are an alternative to convolutional neural networks for processing sequential data such as DNA sequences or time series that implement a different parameter sharing scheme. Recurrent neural networks apply the same operation to each sequence element. This operation takes as input the memory of the previous sequence element and the new input. It updates the memory and optionally emits an output that is either passed to subsequent layers or used directly as a model prediction. By applying the same model to each sequence element, recurrent neural networks are invariant to position index in the processed sequence. For example, recurrent neural networks can detect open reading frames in DNA sequences regardless of their position in the sequence. This task requires the recognition of a specific sequence of inputs, such as a start codon followed by an in-frame stop codon.

[0021] The main advantage of recurrent neural networks over convolutional neural networks is that, in theory, they can carry over information through infinitely long sequences via memory. Furthermore, recurrent neural networks can naturally handle sequences of widely varying length, such as mRNA sequences. However, convolutional neural networks combined with various tricks (such as dilated convolutions) can reach performance comparable to or even better than recurrent neural networks for sequence modeling tasks such as audio synthesis and machine translation. Recurrent neural networks can aggregate the output of convolutional neural networks to predict single-cell DNA methylation status, RBP binding, transcription factor binding, and DNA accessibility. Furthermore, recurrent neural networks apply sequential operations, so they cannot be easily parallelized and are therefore much slower computationally than convolutional neural networks.

[0022] Although most of the human genetic code is common to all humans, each human has a unique genetic code. In some cases, the human genetic code may contain outliers, called genetic variants, that may be common among a relatively small group of individuals in the human population. For example, a particular human protein may contain a particular sequence of amino acids, but variants of that protein may differ by one amino acid in an otherwise identical particular sequence.

[0023] Gene variants can be pathogenic and can result in disease. Most such gene variants have been depleted from the genome by natural selection, but the ability to identify which gene variants are likely to be pathogenic can help researchers focus on these gene variants to gain understanding of the corresponding diseases and their diagnosis, treatment, or cure. The clinical interpretation of millions of human gene variants remains unclear. Some of the most frequent pathogenic variants are single nucleotide missense mutations that change the amino acids of proteins. However, not all missense mutations are pathogenic.

[0024] Models that can predict molecular phenotypes directly from biological sequences can be used as in silico perturbation tools to investigate the association between genetic and phenotypic variations and have emerged as new methods for quantitative trait locus identification and variant prioritization. These approaches are of great importance considering that the majority of variants identified by genome-wide association studies of complex phenotypes are non-coding, which makes it difficult to estimate their effect and contribution to the phenotype. Furthermore, linkage disequilibrium results in blocks of co-inherited variants, which makes it difficult to pinpoint individual causal variants. Thus, sequence-based deep learning models that can be used as matching tools to assess the impact of such variants provide a promising approach to find potential drivers of complex phenotypes. One example is predicting the effect of non-coding single nucleotide variants and short insertions or deletions (indels) indirectly from the differences between two variants on transcription factor binding, chromatin accessibility or gene expression prediction. Another example is predicting novel splice site generation from the quantitative effect of gene variants on sequence or splicing.

[0025] To predict the pathogenicity of missense variants from protein sequence and sequence conservation data, an end-to-end deep learning approach for variant effect prediction is applied (see Sundaram, L. et al. Predicting the clinical impact of human mutation with deep neural networks. Nat. Genet. 50, 1161-1170 (2018), referred to herein as "PrimateAI"). PrimateAI uses a deep neural network trained on variants of known pathogenicity with data augmentation using cross-species information. In particular, PrimateAI uses wild-type and mutant protein sequences to compare differences and determine the pathogenicity of the variant using a trained deep neural network. Such an approach utilizing protein sequences for pathogenicity prediction is promising because it can avoid the circularity problem and overfitting to prior knowledge. However, the number of clinical data available in ClinVar is relatively small compared to the number of data sufficient to effectively train a deep neural network. To overcome this data scarcity, PrimateAI uses common human variants and primate-derived variants as benign data, and simulated variants based on trinucleotide context as unlabeled data.

[0026] PrimateAI outperforms conventional methods when trained directly on sequence alignments. PrimateAI learns important protein domains, conserved amino acid positions, and sequence dependencies directly from training data consisting of approximately 120,000 human samples. PrimateAI substantially outperforms other variant pathogenicity prediction tools in distinguishing benign and pathogenic de novo mutations in candidate developmental disorder genes and in reproducing prior knowledge in ClinVar. These results suggest that PrimateAI is an important step forward for variant classification tools that can reduce reliance on prior knowledge of clinical reports.

[0027] Central to protein biology is the understanding of how structural elements give rise to observed functions. The plethora of protein structural data enables the development of computational methods to systematically derive the rules governing structure-function relationships. However, the performance of these methods critically depends on the choice of protein structural representation.

[0028] Protein sites are microenvironments within a protein structure that are differentiated by their structural or functional role. Sites can be defined by a three-dimensional (3D) location and the local neighborhood around this location where the structure or function resides. Central to rational protein engineering is the understanding of how the structural arrangement of amino acids creates functional features within a protein site. Determination of the structural and functional roles of individual amino acids in a protein provides information to aid in the manipulation and modification of protein function. Identifying functionally or structurally important amino acids allows for focused engineering efforts such as site-directed mutagenesis to modify the functional properties of a target protein. Alternatively, this knowledge can help avoid engineering designs that disable desired functions.

[0029] Since it is established that structure is much more conserved than sequence, the increase in protein structural data provides an opportunity to systematically study the underlying patterns governing structure-function relationships using data-driven approaches. A fundamental aspect of any computational protein analysis is how protein structural information is represented. The performance of machine learning methods often depends more on the choice of data representation than on the machine learning algorithm used. A good representation efficiently captures the most important information, whereas a poor representation produces a noisy distribution lacking the underlying pattern.

[0030] The plethora of protein structures and the recent success of deep learning algorithms provide an opportunity to develop tools for automatically extracting task-specific representations of protein structures.

[0031] Of the more than 70,000,000 possible missense variants in the human genome, the majority have unknown clinical significance, and only about 0.1% have been annotated in clinical variant databases. Opportunities arise to understand the role of rare penetrant variants in common diseases. Accurately distinguishing harmful variants from those with benign outcomes could benefit both precision medicine and targeted drug development. [Brief description of the drawings]

[0032] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee. Color drawings may also be available in PAIR via the Supplemental Content tab.

[0033] In the drawings, like reference characters generally refer to like parts throughout the different views. Also, the drawings are not necessarily to scale, emphasis instead being placed upon illustrating the principles of the disclosed technology. In the following description, various embodiments of the disclosed technology are described with reference to the following drawings, in which: [Figure 1A] FIG. 1 is a schematic illustrating a method for determining weighted rare variant PRS. [Figure 1B] 1 illustrates an exemplary calculation of a rare variant polygenic risk score for a particular plurality of genes and a particular phenotype. [Diagram 2] 1 illustrates an exemplary computer system that can be used to implement the disclosed techniques. [Diagram 3] Illustrate examples of genetic variants. [Figure 4] This is an example illustrating phenotypic effects in response to genotypic changes. [Diagram 5] Illustrates multiple phenotypes corresponding to patient X with cardiovascular disease. [Figure 6] Graphical contrast of genome-wide association studies for common versus rare variants. [Figure 7] Each illustrates a genetic association test for specific individual common variants or specific aggregated rare variants at gene resolution. [Figure 8] FIG. 1 is a schematic diagram of a method for optimizing rare variant summation tests. [Figure 9] FIG. 1 is a flow diagram illustrating the process of a system for determining pathogenicity of a variant. [Figure 10] 1 illustrates an exemplary processing architecture of a pathogenicity classifier, according to one implementation of the disclosed technology. [Figure 11] 1 illustrates an exemplary computer system that can be used to implement the disclosed techniques. [Figure 12] 1 illustrates one embodiment of determining the final pathogenicity score. [Figure 13] FIG. 1 is a flow diagram of a process for correcting for multiple testing within each gene. [Figure 14] FIG. 1 is a schematic diagram of a method for correcting phenotypic values ​​for covariates. [Figure 15] FIG. 1 is a schematic depicting predicting phenotypic shift in response to the use of multiple drugs on multiple phenotypes. [Figure 16] FIG. 1 is a schematic diagram depicting the experimental set-up for obtaining drug use patterns and phenotype data for specific cohorts. [Figure 17] Illustrates multiple phenotypes corresponding to patient X with cardiovascular disease. [Figure 18] Graph quantifying the total number of significant gene-phenotype pairs identified for different types of summation tests. [Figure 19-1] A collection of graphs illustrating rare deleterious variants influencing disease severity and age of onset identified by the pathogenicity classifier PrimateAI-3D. [Figure 19-2] A collection of graphs illustrating rare deleterious variants influencing disease severity and age of onset identified by the pathogenicity classifier PrimateAI-3D. [Figure 19-3]A collection of graphs illustrating rare deleterious variants influencing disease severity and age of onset identified by the pathogenicity classifier PrimateAI-3D. [Figure 20] Graph of the mean absolute Spearman correlation between different pathogenicity scores and phenotypic values. [Figure 21] Heatmap comparison of rare deleterious variants with common genome-wide association study variants. [Figure 22] Further comparison of rare deleterious variants with common genome-wide association study variants. [Figure 23-1] FIG. 1 illustrates cholesterol pathways and total cholesterol distribution across all individuals in the UK Biobank cohort. [Figure 23-2] FIG. 1 illustrates cholesterol pathways and total cholesterol distribution across all individuals in the UK Biobank cohort. [Figure 24] Includes measurement of rare variant PRS performance. [Diagram 25] 13 is a graph of PRS outlier enrichment. [Figure 26] 13 is a graph of PRS outliers for quantitative phenotypes. [Figure 27] 1 includes graphs of normalized total cholesterol distribution from two separate cohorts. [Figure 28] Graphs comparing rare variant PRS outliers and phenotypes between two separate cohorts are included. [Figure 29] Includes graphs illustrating performance results by ethnicity. [Diagram 30] 1 is a table comparing effect sizes and frequencies for common PRS variants and rare PRS genes used for normal cholesterol levels. [Diagram 31] Graph of the average proportion of phenotypic variance explained by different pathogenicity scoring methods for a set of 34 gene-phenotype pairs, which were selected based on their enrichment for rare missense and LoF variants. [Diagram 32]Heatmap of enrichment of rare variants in GWAS genes for all pairwise comparisons between quantitative and clinical phenotypes. [Diagram 33] Includes a graph of the distribution of the absolute value of the ratio between the mean effect size of singleton LoF variants and the effect size of the most significant GWAS variant for the same gene. [Diagram 34] The number of variants per individual in any of the genes from the gene-phenotype pair that were significant at a 5% false discovery rate is shown. [Diagram 35] A comparison of effect sizes for training vs. testing data splits is shown. [Diagram 36] Included are graphs comparing common variant PRS subsets with rare variant PRS subsets. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0034] The following discussion is presented to enable those skilled in the art to make and use the disclosed technology and is provided in the context of a particular application and its requirements. Various modifications to the disclosed embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be applied to other embodiments 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 embodiments shown, but is to be accorded the widest scope consistent with the principles and features disclosed herein.

[0035] The detailed description of the various embodiments can be better understood when read in conjunction with the accompanying drawings. To the extent that the figures show diagrams of functional blocks of the various embodiments, the functional blocks do not necessarily show a division between hardware circuitry. Thus, for example, one or more of the functional blocks (e.g., a module, a processor, or a memory) may be implemented in a single piece of hardware (e.g., a general-purpose signal processor or a block of random access memory, a hard disk, etc.) or in multiple pieces of hardware. Similarly, a program may be a stand-alone program, may be incorporated as a subroutine in an operating system, may be a function in an installed software package, etc. It should be understood that the various embodiments are not limited to the arrangements and instrumentalities shown in the figures.

[0036] The processing engines and databases in the figures designated as modules can be implemented in hardware or software and need not be divided in exactly the same blocks as shown in the figures. Some modules may be implemented on different processors, computers or servers, or even spread among many different processors, computers or servers. In addition, it will be understood that some of the modules can be operated in parallel or in a different order than shown in the figures without affecting the functionality achieved. The modules in the figures can also be considered as flow chart steps in a method. Also, a module does not necessarily have to have all code located contiguously in memory. Some portions of code can be separated from other portions of code, with code from other modules or other functions located in between.

[0037] The disclosed technology can be used to improve the quality of polygenic risk scoring for rare variants. The disclosed technology can be used to identify patients at the extremes of the phenotypic spectrum who are at highest risk for severe early-onset disease and who will benefit most from clinical intervention. The researchers presenting these disclosures show that outlier phenotypes associated with severe early-onset disease are better explained by rare penetrant variants than by the collective action of many common variants with small effect sizes. Unlike common variants, whose mostly deleterious consequences have been removed by natural selection over many generations, rare variants remain largely unfiltered and retain the potential to exert highly penetrant effects in complex traits and diseases. The disclosed technology includes a novel weighted sum model for quantifying the contribution of rare pathogenic variants to phenotypic response.

[0038] Introduction The disclosed technology develops a complementary rare variant polygenic risk score model based on a weighted sum of rare deleterious variants from multiple phenotype-associated genes. Genome-wide association study (GWAS) is a tool for linking genetic mutations or variants with phenotypes of genetic diseases and complex traits, such as delayed development disorder (DDD). GWAS tools perform better with common or semi-common variants than with very rare variants. For example, the UK Biobank cohort contains 200,643 exomes and associated patient data. Nearly half of the rare variants considered deleterious appeared in only one individual in the dataset. Traditional statistical analysis and summation tests break down when there is only one instance of a rare variant in a large cohort.

[0039] Although individual genome-wide association study (GWAS) variants confer effects that tend to be too mild for clinical actionability, polygenic risk scores that combine signals from hundreds to millions of common variants have shown significant effectiveness for predicting phenotypic extremes of patients at high risk of disease. However, existing common variant polygenic risk score models largely exclude rare variants due to challenges in interpreting variants of unknown significance and imprecision in estimating their effect size. Compared to existing common variant polygenic risk score models, the disclosed technology includes a rare variant polygenic risk score model that is configured to aggregate risk across genes based on whether an individual carries a rare deleterious variant in each gene.

[0040] An individual's genotype contributes to the individual's phenotype (phenotype is defined as one or more observable physical characteristics of an individual), and thus there is a statistical correlation between genotype and phenotype that can be used to aid in predicting the presence of physical traits, abnormalities, or diseases in individuals with a particular variant or group of variants. While certain genetic disorders can be caused by variants in a single gene (i.e., monogenic genetic disorders), genetic disorders caused by variants in multiple genes (i.e., polygenic genetic disorders) are more common. Thus, there is a need for a robust method to determine polygenic risk scores for rare deleterious variants that are likely to cause severe disease.

[0041] Recent large-scale genome and exome sequencing studies of healthy individuals from the general population have revealed that the average human harbors dozens of potentially deleterious rare variants that have arisen through recent mutations. Unlike common variants, whose mostly deleterious consequences have been removed by natural selection over many generations, rare variants remain largely unfiltered and retain the potential to exert highly penetrant effects in complex traits and diseases. The public release of 200,643 UK Biobank exomes, together with rapid advances in the accuracy of variant pathogenicity prediction, creates an opportunity to investigate the impact of rare penetrant variants on a comprehensive set of common human diseases and complex traits, and provides insight into the potential utility of personal genome sequencing for the general population.

[0042] The recent exponential growth of the human population has created a large number of rare variants through randomly occurring mutations, without providing sufficient time for natural selection to weed out those with deleterious consequences, in contrast to the common variants identified in GWAS. We set out to test that in each GWAS locus that contains a common variant associated with mild clinical risk, there should also be rare deleterious variants with much greater severity, forming an allelic series in which the severity of the variant is inversely proportional to its frequency in the population.

[0043] Despite the high penetrance of rare deleterious variants, their rarity limits their ability to explain phenotypic variance to a small percentage of the population, since the majority of individuals do not carry rare deleterious variants for a given phenotype.Thus, one embodiment of the rare variant PRS described in the disclosed technology performs only about 1 / 20th of the common variant PRS in terms of explained variance across the entire population.However, when considering individuals with extreme phenotypes, this trend is reversed: individuals with outlier phenotypes (z-score≧3) are 3 times more likely to have a rare variant PRS score at the 1st or 99th percentile than the baseline population, compared to 1.8 times the common variant PRS.

[0044] A notable barrier to the clinical adoption of common variant PRS models has been their limited generalizability across populations with different ancestries. These issues stem from the incorporation of disease-associated, noncausal variants as predictors in common variant PRS models due to linkage disequilibrium. Although the effects of causal variants can be expected to generalize across cohorts, there is no guarantee that correlations between true causal variants and variants used in the models will hold; rather, these correlations are influenced by differences in population ancestry and technical artifacts in genotyping and imputation. In comparison, rare variant PRS models directly use rare deleterious variants as predictors and are not as affected by linkage disequilibrium issues that make it difficult to distinguish causation from correlation for common variants, given their rarity and recent history in the population. Rather, the challenge for predicting rare variant polygenic risk lies in the accurate interpretation of the effects of variants of uncertain significance (VUS), which is critical for both identifying genes with significant rare variant associations and estimating effect sizes of rare deleterious variants for a given clinical phenotype. With nearly half of the rare deleterious variants occurring in rare variant PRS genes found in only one individual across the entire UK Biobank exome cohort of 200,643, this problem may seem intractable for traditional statistical analysis; however, recent advances employing deep learning, high-throughput experimental assays, and variant information from closely related primate species have each demonstrated progress toward solving the VUS problem on a genome-wide scale.The exceptional allelic heterogeneity underlying rare variant PRS genes can be a challenge in terms of variant effect prediction, but paradoxically, it is also key to the robustness and portability of rare variant PRS models across different cohorts; by integrating signal across thousands of unique rare deleterious variants per gene, rare variant PRS achieves the desired property of smoothing the effects of any variant-specific artifacts that may be present in the data.

[0045] The extent to which rare and common genetic variants contribute to the risk of common diseases and complex traits has been debated for decades in the field of human genetics. The disclosed technology helps reconcile these perspectives by showing the presence of rare penetrant variants in a large proportion of GWAS loci and quantifying the relative contributions of rare and common variants to total population variance as well as outlier individuals at highest risk for severe, early-onset disease.

[0046] From the perspective of precision medicine, the UK Biobank cohort is generally a representative cross-section of adults from the UK population, and also presents a unique opportunity to characterize the sum of rare penetrant variants and their effects in the general population.Across 500 genes that were significant for one or more of the 90 clinical and quantitative phenotypes studied in one embodiment of the disclosed technology, 86% of individuals carried at least one rare deleterious variant, with an average of 2.03 rare deleterious variants per individual.Some embodiments of the disclosed technology show that these rare deleterious variants have high penetrant effects, contributing on average 10 times larger effect size than common GWAS variants present at the same locus.On average, 5.2% of individuals carried rare deleterious variants for a given phenotype, and 0.4% of individuals carried rare penetrant variants for a given gene.

[0047] The disclosed technology augments the rare variant summation test approach to maximize power and filter out benign missense variants. In some embodiments of the disclosed technology, summation test power is maximized by performing a first grid search across the allele counts of the observed variants to find a cutoff for the maximum allele count that maximizes the significance of the summation test, and benign missense variants are filtered out by performing a second grid search across multiple cutoffs for the pathogenicity score (e.g., PrimateAI-3D score) threshold to find an optimal threshold for the pathogenicity score that maximizes the significance of the summation test.

[0048] The disclosed technology develops a rare deleterious variant polygenic risk score model based on the weighted sum of rare deleterious variants from multiple phenotype-related genes. Individual genome-wide association study variants tend to give overly mild effects on clinical viability and genetic disorders caused by variants in multiple genes compared to single-gene disorders. Thus, polygenic risk score models that combine signals from multiple variants are often more effective for predicting patients at high risk of disease. However, existing polygenic risk score models are primarily optimized for common variants and exclude rare variants due to the difficulty of interpreting variants of uncertain significance, challenges from low-powered studies, and inaccurate effect size estimates due to small sample sizes of rare variants. The disclosed technology improves prior polygenic risk score models to include rare variants by aggregating risk on a gene resolution basis based on whether an individual carries at least one rare deleterious variant in a particular gene.

[0049] The disclosed rare polygenic risk scores are calculated as summation tests in which the strength of association quantifies the contribution of rare variants to genes associated with a phenotype and the phenotypic response. Data representing variant carrier status and specific measured phenotypic values ​​are obtained for a cohort of individuals. A summation score is calculated for each of the associated genes using the obtained data, and the summation score identifies the resulting non-random associations in the cohort between the carrier status of each of the associated genes and the phenotypic response to the presence of one or more rare pathogenic variants in the associated gene. The carrier status is a Boolean variable determined by the presence or absence of one or more rare pathogenic variants in a particular gene. The pathogenicity of a particular rare variant is determined by the predicted effect of the particular rare variant on the function of the gene when expressed as a protein, and the predicted effect is determined as a pathogenicity score by a convolutional neural network pathogenicity classifier. The rarity of a particular rare variant is determined by the occurrence of the particular rare variant in a population below a predetermined threshold (i.e., maximum allele count).

[0050] The strength of the association is quantified as an effective intensity score for the carrier status of the rare pathogenic variant for each gene, and a respective phenotypic response for the carrier status of the rare pathogenic variant at a gene-by-gene resolution for each gene. A score for the contribution of the carrier status of the subjects in the cohort across the associated genes to the phenotypic response can be determined based on the effective intensity score. In some embodiments of the disclosed technology, the particular phenotype is a quantitative biomarker phenotype, and the effective intensity score for the resulting non-random association for the gene is determined for the quantitative biomarker phenotype using a two-tailed t-test on the linear regression component of the sum score. The two-tailed t-test generates a p-value to determine whether the difference between the mean phenotypic measurements of carriers and non-carriers is significant at a predetermined significance level.

[0051] In another embodiment of the disclosed technology, the specified phenotype is a categorical clinical diagnostic phenotype, and a significance strength score for the resulting non-random association for the gene is determined for the categorical clinical diagnostic phenotype as a beta coefficient for carrier status, the beta coefficient being determined using a logistic regression component of the sum score, the logistic regression component being adapted to predict a clinical diagnostic label from a binary indicator variable corresponding to carrier status and a plurality of covariates.

[0052] In some embodiments of the disclosed technology, each summation test for a particular gene and a particular phenotype is optimized for a particular gene and corrected for a particular phenotype value. Some embodiments of the disclosed technology may involve either optimization for a particular gene or correction for a particular phenotype value, but not the other.

[0053] Further embodiments of the disclosed technology optimize rare variant summation tests by using nested t-tests that maximize separation between carriers and non-carriers, where summation test parameters are optimized for a particular gene. The summation test power is maximized by performing a first grid search across the allele counts of the observed variants to find a cutoff for the maximum allele count that maximizes the significance of the summation test, and benign missense variants are filtered out by performing a second grid search across multiple cutoffs for the pathogenicity score (e.g., PrimateAI-3D score) threshold to find an optimal threshold for the pathogenicity score that maximizes the significance of the summation test.

[0054] In some embodiments of the disclosed technology, the grid search procedure involves searching a plurality of allele counts and a plurality of pathogenicity score thresholds, the grid search including generating a plurality of combinations of allele counts and pathogenicity score thresholds from the plurality of allele counts and the plurality of pathogenicity score thresholds, identifying a plurality of groups of rare pathogenic variants corresponding to the plurality of combinations of allele counts and pathogenicity score thresholds, summation testing of the plurality of groups of rare pathogenic variants depending on carrier status separating carriers of the particular group of rare pathogenic variants in the cohort of individuals from non-carriers of the particular group of rare pathogenic variants in the cohort of individuals, and determining a plurality of effect sizes and p-values ​​corresponding to the plurality of combinations of allele counts and pathogenicity score thresholds.

[0055] The particular combination of allele count and pathogenicity score threshold with the most significant p-value is selected to be used as the optimal combination for the particular gene. Each gene may or may not have a similar maximum allele count threshold and / or minimum pathogenicity score threshold that results in optimal separation between carriers and non-carriers of each respective gene. The pathogenicity score threshold in the plurality of pathogenicity score thresholds corresponds to a pathogenicity score quantile of the pathogenicity scores determined for the rare pathogenic variants in the group of rare pathogenic variants. In some embodiments of the disclosed technology, the pathogenicity score is generated by a convolutional neural network pathogenicity classifier such as PrimateAI-3D. In other embodiments of the disclosed technology, a wide range of additional AI, machine learning, and deep learning models are employed to generate the pathogenicity score.

[0056] The grid search procedure above requires correction for multiple testing on the generated false discovery rate corrected p-values. In some embodiments of the disclosed technology, multiple testing correction is performed using the Benjamini-Hochberg procedure. In other embodiments of the disclosed technology, multiple testing correction is performed using an adaptive permutation testing procedure.

[0057] Further embodiments of the disclosed technology apply covariate correction for drug use across time-distributed detection points. The summation test input data is corrected for the measured phenotypic values. The phenotypic response values ​​are covariate-corrected for multiple covariates and drug-use corrected for multiple drugs. In some embodiments of the disclosed technology, phenotypic shifts in response to multiple drugs for multiple phenotypes of a cohort of individuals with multiple confounders are predicted to generate covariate-corrected and drug-use corrected phenotypic measurements.

[0058] In some embodiments of the disclosed technology, phenotypic shift is predicted in response to use of multiple drugs for multiple phenotypes of a cohort of individuals with multiple confounding factors (e.g., age, sex, genetic principal components, diet, and smoking status) for each phenotype for the cohort of individuals, and phenotypic measures for the multiple phenotypes, covariate measures for the multiple confounding factors, and drug use patterns for the multiple drugs for each individual in the cohort at two separate time points.

[0059] The phenotypic measurements are covariate-corrected for the first and second time points based on the covariate measurements by fitting a first regression model and regressing the covariate measurements, thereby generating covariate-corrected phenotypic measurements for the first and second time points. A delta is determined based on the difference between the covariate-corrected phenotypic measurements for the first and second time points. For each of the drug use patterns, a second regression model is fitted that uses the delta to predict a phenotypic shift in response to use of the multiple drugs against the covariate-corrected phenotypic measurements.

[0060] The second regression model is a forward selective stepwise regression model that iteratively predicts delta by successively and cumulatively including phenotypic shifts for each of the drug use patterns. The second regression model has a binary indicator independent variable for each of the drug use patterns, the drug use patterns including not taking the drug at the first and second time points, starting to take the drug between the first and second time points, stopping to take the drug between the first and second time points, and taking the drug at the first and second time points. In some embodiments of the disclosed technology, the covariate correction, delta determination, phenotypic shift prediction, and drug use correction are performed on a per-drug basis for drugs in the plurality of drugs. In other embodiments of the disclosed technology, the covariate correction, delta determination, phenotypic shift prediction, and drug use correction are performed on a per-drug basis for drug categories within the plurality of drug categories (e.g., statins, NSAIDs, opioids, etc.).

[0061] The second regression model further models the phenotypic shift in response to the time elapsed between the first and second time points for each individual in the cohort of individuals, and in response to regression to the mean between the first and second time points. The second regression model is fitted jointly to the set of related drugs in the multiple drugs.

[0062] In some embodiments of the disclosed technology, the drug use correction is implemented by fitting a third regression model and includes drug use correcting the phenotypic measure for the first time point based on a first binary indicator independent variable for a first drug use pattern of starting to take a drug between the first and second time points, a second binary indicator independent variable for a second drug use pattern of not taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the first time point.

[0063] A fourth regression model is fitted, the fourth regression model including drug use correction of the phenotype measure for the second time point based on a third binary indicator independent variable for a third drug use pattern of ceasing taking the drug between the first and second time points, a fourth binary indicator independent variable for a fourth drug use pattern of taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a particular drug at the second time point.

[0064] A rank-based inverse normal transformation may be applied to the drug-use-corrected phenotypic measures for the first and second time points to generate normalized drug-use-corrected phenotypic measures for the first and second time points. The normalized drug-use-corrected phenotypic measures are then covariate-corrected for the first and second time points to generate covariate-corrected normalized drug-use-corrected phenotypic measures for the first and second time points. The covariate-corrected normalized drug-use-corrected phenotypic measures may be used to generate a rare variant polygenic risk score, where the measures for the phenotype of interest are corrected for phenotypic shifts in response to covariates and drug use patterns.

[0065] In some embodiments of the disclosed technology, the constructed regression models can be used to detect drug-phenotype associations, such as undesirable side effects and desired target effects.

[0066] Rare Variant Polygenic Risk Score (PRS) FIG. 1A is a schematic diagram illustrating a method for determining weighted rare variant PRS. Database 102 includes genomic and phenotypic data corresponding to a group of individuals belonging to cohort i. Cohort i 104 includes N individuals who have been genetically sequenced and tested for multiple phenotypic measurements. The sequencing data for the individuals in cohort i 104 includes data corresponding to multiple variant carrier states. A particular gene may carry three known possible rare deleterious variants (e.g., variant A, variant B, and variant C), of which an individual may be a carrier (i.e., an individual carries each variant and therefore is a carrier of that variant) or a non-carrier (i.e., an individual does not carry each variant and therefore is not a carrier of that variant). Those skilled in the art will recognize that variant A, variant B, and variant C are provided as illustrative examples, and that a gene may carry any number of rare deleterious variants. The disclosed technology defines a variant as "rare" if the variant occurrence in a population is below a predetermined threshold, and "detrimental" if the pathogenicity of the variant is predicted to have a measured effect on the effect when the gene is expressed as a protein above the predetermined threshold. Each predetermined threshold for both rarity and pathogenicity is gene-specific (i.e., each individual gene will correspond to an individually determined rarity threshold and pathogenicity threshold that may differ from another gene).

[0067] In addition to the genomic data obtained from genome and exome sequencing for each individual in cohort i104, phenotypic data is also available for each individual (i.e., individual N has x n(having a measured value of p = 0.01). By constructing a model 122 for the relationship between a particular genotype and a particular phenotype, it is possible to measure the effect size of a particular genotype on a particular phenotype. Box plot 124 illustrates a representative genetic association test in which the influence of carrier status of one or more particular rare deleterious variants on phenotype D is measured in the form of both a p-value (i.e., the degree of significance when comparing the mean phenotypic values ​​between carriers and non-carriers) and a beta coefficient (i.e., a weighted estimate obtained from a regression line connecting the mean phenotypic values ​​of carriers and non-carriers, where the underlying data has been standardized so that the variances of the dependent and independent variables are equal to one).

[0068] In some embodiments of the present technology where the particular phenotype is a quantitative biomarker phenotype, the strength of the association is quantified as a p-value determined by a two-tailed t-test. In other embodiments of the disclosed technology where the particular phenotype is a qualitative phenotype (e.g., a categorical clinical diagnostic phenotype), the strength of the association is quantified as the beta coefficient of a logistic regression.

[0069] In the disclosed technology, carrier status is defined on a gene resolution basis, not a variant resolution basis. For a particular gene, an individual is defined as a carrier if the individual carries at least one rare deleterious variant in the particular gene (i.e., an individual who is a carrier for two rare deleterious variants is not distinguished from an individual who is a carrier for one rare deleterious variant, but an individual who is a carrier for zero rare deleterious variants is distinguished as a non-carrier). A determination 142 can be made of an effective strength score of the strength of association in the cohort i104 between carrier status and phenotypic response of multiple rare variants.

[0070] Graph 144 includes a t-test for determining an effective strength score of the strength of the association in cohort i 104 between carrier status and phenotypic response for a plurality of rare variants. The null hypothesis of the t-test states that the absolute value of the difference between the sample mean of carriers in cohort i and the sample mean of non-carriers in cohort i 104 is equal to zero. The alternative hypothesis of the t-test states that the absolute value of the difference between the sample mean of carriers in cohort i 104 and the sample mean of non-carriers in cohort i 104 is not equal to zero. The decision to accept or reject the null hypothesis is driven by a particular significance level α and the resulting p-value corresponding to α. Each t-test can be performed for each of a plurality of genes to obtain a plurality of gene-specific effective strength scores (measured by p-values ​​or β coefficients, as represented in chart 124) for a particular shared phenotypic response.

[0071] A weighted sum test 162 can be generated from the plurality of respective gene-specific effective strength scores. This calculation, described as the rare variant PRS, is shown in equation 164, where the rare variant PRS is equal to the sum of the products of the effect size and carrier status for a particular gene for the plurality of genes.

[0072] Here, the discussion provides an example of using a plurality of respective gene-specific effective intensity scores to generate a rare variant PRS for a particular individual, a particular plurality of respective genes sequenced for a particular individual, and a particular phenotype of interest.

[0073] 1B illustrates an exemplary calculation of a rare variant polygenic risk score for a particular plurality of genes and a particular phenotype. Header row 108 includes a group of column labels corresponding to genes, effect sizes, and carrier status for a particular individual. The genes included within the gene column are gene A (row 118), gene B (row 128), gene C (row 138), gene D (row 148), and gene E (row 158), where genes A-E may be hypothesized to be genes responsible for a particular phenotype Y. With respect to the effect size for each gene, effect size A corresponds to gene A (row 118), effect size B corresponds to gene B (row 128), effect size C corresponds to gene C (row 138), effect size D corresponds to gene D (row 148), and effect size E corresponds to gene E (row 158). Carrier status is a Boolean variable that describes the presence of at least one identified rare deleterious variant in a particular gene (e.g., assume variants A, B, and C detected in cohort i104 are specific to gene A (row 118), and therefore each individual listed in the table for cohort i104 would be considered a carrier for gene A). The individuals whose genomes are characterized in FIG. 1B are carriers for at least one rare deleterious variant in gene A (row 118), gene C (row 138), gene D (row 148), and gene E (row 158), respectively. The individuals are not carriers for any rare deleterious variants in gene B (row 128).

[0074] Equation 164 is applied to output weighted sum score 162. Equation 164 for a particular individual may be considered equivalent to equation 182, where the particular individual has a rare variant PRS for a particular phenotype Y that is equivalent to the sum of effect size A corresponding to gene A (row 118), effect size C corresponding to gene C (row 138), effect size D corresponding to gene D (row 148), and effect size E corresponding to gene E (row 158). For each gene for which the particular individual is a carrier for at least one rare deleterious variant, the carrier status is equal to 1 (i.e., effect size *Effect size B corresponding to gene B (row 128) is not shown in equation 182 because the particular individual is not a carrier for at least one rare deleterious variant, and therefore the carrier status is equal to zero (i.e., effect size * 0=0), the effect size B corresponding to gene B (row 128) is not present in equation 182.

[0075] 2 illustrates an exemplary computer system 200 that can be used to implement the disclosed techniques. Computer system 200 includes at least one central processing unit (CPU) 244 that communicates with a number of peripheral devices via a bus subsystem 242. These peripheral devices may include, for example, a storage subsystem 210 including memory devices and file storage subsystem 236, user interface input devices 238, user interface output devices 248, and a network interface subsystem 246. The input and output devices enable user interaction with computer system 200. Network interface subsystem 246 provides an interface to external networks, including interfaces to corresponding interface devices in other computer systems.

[0076] In one embodiment, the rare variant PRS model 240 is communicatively linked to the storage subsystem 210 and the user interface input device 238 .

[0077] The user interface input devices 238 can include pointing devices such as a keyboard, a mouse, a trackball, a touchpad, or a graphics tablet, a scanner, a touch screen integrated into a display, audio input devices such as a voice recognition system and a microphone, as well as other types of input devices. In general, use of the term "input device" is intended to include all possible types of devices and manners for inputting information into the computer system 200.

[0078] The user interface output devices 248 may include a display subsystem, a printer, a fax machine, or a non-visual display such as an audio output device. The display subsystem may include a flat panel device such as an LED display, a cathode ray tube (CRT), a liquid crystal display (LCD), a projection device, or some other mechanism for producing a visible image. The display subsystem may also provide non-visual displays such as an audio output device. In general, use of the term "output device" is intended to include all possible types of devices and manners for outputting information from computer system 200 to a user or to another machine or computer system.

[0079] The storage subsystem 210 stores programming and data constructs that provide the functionality of some or all of the modules and methods described herein. These software modules are generally executed by the processor 249.

[0080] The processor 249 can be a graphics processing unit (GPU), a field-programmable gate array (FPGA), an application-specific integrated circuit (ASIC), and / or a coarse-grained reconfigurable architecture (CGRA). The processor 249 can be hosted by a deep learning cloud platform such as Google Cloud Platform™, Xilinx™, and Cirrascale™. Examples of processors 249 include Google's Tensor Processing Unit (TPU)™, rackmount solutions such as the GX4 Rackmount Series™, GX2 Rackmount Series™, NVIDIA DGX-1™, Microsoft' Stratix V FPGA™, Graphcore's Intelligent Processor Unit (IPU)™, Qualcomm's Zeroth Platform™ with Snapdragon processors™, NVIDIA's Volta™, NVIDIA's DRIVE PX™, NVIDIA's JETSON TX1 / TX2 MODULE™, Intel's Nirvana™, Movidius VPU™, Fujitsu DPI™, ARM's DynamicIQ™, IBM TrueNorth™, Lambda GPU Server with Testa V100s™, and others.

[0081] The memory subsystem 222 used in the storage subsystem 210 may include multiple memories including a main random access memory (RAM) 232 for storing instructions and data during program execution, and a read only memory (ROM) 234 in which fixed instructions are stored. The file storage subsystem 236 may provide persistent storage for program and data files and may include a hard disk drive, associated removable media, a CD-ROM drive, an optical drive, or a removable media cartridge. Modules that implement the functionality of a particular embodiment may be stored by the file storage subsystem 236 in the storage subsystem 210 or in another machine accessible by the processor.

[0082] Bus subsystem 242 provides a mechanism for allowing the various components and subsystems of computer system 200 to communicate with each other as intended. Although bus subsystem 242 is shown generally as a single bus, alternative implementations of the bus subsystem may use multiple buses.

[0083] The computer system 200 itself can be of a variety of types, including a personal computer, a portable computer, a workstation, a computer terminal, a network computer, a television, a mainframe, a server farm, a loosely distributed set of loosely networked computers, or any other data processing system or user device. Due to the ever-changing nature of computers and networks, the description of computer system 200 shown in Figure 2 is intended only as a specific example for purposes of illustrating a preferred embodiment of the present invention. Many other configurations of computer system 200 can have more or fewer components than the computer system shown in Figure 2.

[0084] Alternative rare deleterious polygenic risk score models In one embodiment, the disclosed technology may use a rare deleterious polygenic risk score model that differs from the model discussed above. Instead of using a weighted sum of rare deleterious variants, where the weight assigned to a variant in a gene is equal to the effect size observed across all rare deleterious variants in that gene, i.e., the same for all variants in the gene, the disclosed technology may estimate a weight for each variant in each significant gene. The weight for each variant may be estimated from individuals carrying rare deleterious variants for significant genes via a linear regression of the log10 transformed allele frequency against the variant pathogenicity score (e.g., PrimateAI-3D) and the trait value of the individual. This may improve performance by 18% in terms of explained variance in some embodiments.

[0085] Genotype The method of calculating a polygenic risk score illustrated in Figure 1 involves data associated with an individual's genotype as determined by genotyping arrays (i.e., whole genome sequencing of a particular individual) and exome sequencing (i.e., sequencing of the protein coding regions of a particular individual's genome). Discussion now turns to the details of the genomic data, their associated characteristics, and the relationships between genotype and phenotype for a particular individual.

[0086] FIG. 3 illustrates an example of a genetic variant. An exemplary reference gene sequence A 302 has a length of N, and each sequence position contains a nucleotide base of either thymine, adenine, guanine, or cytosine. The reference gene sequence A 302 can be considered as a natural or normal gene sequence for a particular gene segment. The reference gene sequence A 302 is compared to two exemplary variant sequences 1A 322 and 1B 342, where the variant sequences each have a single nucleotide variant at the fifth base position, but otherwise have the same composition as the reference sequence. For example, a single nucleotide substitution is shown as adenine 324 in variant 1A 322 and thymine 344 in variant 1B 342 compared to cytosine 304 in reference gene sequence A 302. When reference gene sequence A 302, variant 1A 322, and variant 1B 342 are transcribed and translated into amino acid units of a protein, the single nucleotide variants 324 and 344, respectively, may result in at least one similar or different amino acid compared to the amino acid present in the native protein. If variant 1A 322 or variant 1B 342, respectively, results in at least one altered amino acid in the converted protein, then each variant may also result in an altered phenotype for an individual carrying the respective variant.

[0087] FIG. 4 is an example illustrating phenotypic effects in response to genotypic changes. An exemplary reference gene sequence A 402 has a length of N, and each sequence position contains a nucleotide base of either thymine, adenine, guanine, or cytosine. The reference gene sequence A 402 can be considered as the natural or normal gene sequence for a particular gene segment. The reference gene sequence A 402 is compared to three exemplary variant sequences 1A 422, 1B 442, and 1C 462, where the variant sequences each have a single nucleotide variant at the fifth base position, but otherwise have the same composition as the reference sequence. For example, the single nucleotide substitutions are shown as adenine 404 in variant 1A 422, thymine 444 in variant 1B 442, and guanine 464 in variant 1C 462, compared to cytosine 404 in reference gene sequence A 402.

[0088] Reference gene sequence A 402 results in protein A 406, which is the natural protein composition of the particular gene that comprises reference gene sequence A 402. Individuals carrying natural protein A 406 will exhibit phenotype A 408, which is a healthy phenotype. Variant 1A 422 results in mutant protein 1A 426, which contains a missense mutation that results in a protein structure and function that is different from natural protein A 406. Individuals carrying missense protein 1A 406 will exhibit phenotype 1A 428, which is a disease phenotype. Variant 1B 442 results in mutant protein 1B 446, which contains a synonymous mutation that does not result in a change in protein structure and function compared to natural protein A 406 (i.e., no amino acid change occurs in response to the single nucleotide variant at position 5 of variant 1B 442). Individuals carrying the synonymous protein 1B 446 will display phenotype 1B 448, a healthy phenotype similar to phenotype A 408. Variant 1C 462 results in mutant protein 1C 466, which contains a nonsense mutation resulting in a truncated, non-functional protein structure and function compared to native protein A 406. Individuals carrying the nonsense protein 1C 466 will display phenotype 1C 468, likely resulting in non-viable embryos or a significantly shortened life span.

[0089] Those skilled in the art will recognize that variants 1A 422, 1B 442, and 1C 462 are listed as simplified examples, and that the potential phenotypic expressions span a broad spectrum rather than a limited number of distinct expressions. Furthermore, those skilled in the art will also recognize that many phenotypic responses arise due to multiple variants within a single gene, or combined polygenic variant effects. Many embodiments of the disclosed technology specifically address polygenic risk scoring for specific phenotypic responses associated with severe genetic disorders known to cause substantial harm to an individual's quality of life and / or life expectancy.

[0090] Phenotype Considering now the phenotypic representation in more detail: Figure 4 illustrates examples of multiple phenotypes at a broad organism level, while Figure 5 illustrates examples of multiple phenotypes at a targeted measurable level, which can now be described as either quantitative biomarker values ​​or categorical clinical diagnoses.

[0091] Please note that this application uses "quantitative biomarker value" and "quantitative phenotype" interchangeably. Please note that this application also uses "categorical clinical diagnosis" and "qualitative phenotype" interchangeably.

[0092] FIG. 5 illustrates a plurality of phenotypes corresponding to a cardiovascular disease patient X 500. The cardiovascular disease patient X 500 may be described by a plurality of observable phenotypes in response to a plurality of genetic variants in the genome of the patient X. These phenotypes may be qualitative or quantitative. Qualitative phenotypes may include demographic outputs 502 (e.g., ancestry, biological sex, etc.), categorical biomarkers 504 (e.g., binary variables for the presence of a particular protein in blood, binary classification of a patient's metabolite levels in blood determined by a decision threshold, or multi-class variables for the presence of a particular histological morphology), clinical diagnoses 506 (e.g., binary variables for cardiovascular disease diagnoses, binary variables for a particular cancer diagnosis, etc.), or measures of disease severity (e.g., advanced stage of tumor metastasis, classification of disease subtypes such as refractory celiac disease type 1 vs. type 2, etc.).

[0093] Quantitative phenotypes may include accurate blood or urine biomarker values ​​522 (e.g., creatine, cholesterol, low-density lipoprotein (LDL), triglycerides, glucose, etc.), body mass measurements 542 (e.g., height, weight, body fat percentage, body mass index, etc.), or vital signs (e.g., resting heart rate, systolic and diastolic blood pressure, respiratory rate, etc.). One of skill in the art will recognize that the quantitative and qualitative phenotypes enumerated for cardiovascular disease patient X 500 are non-limiting examples and that there are an infinite number of observable values ​​that may be employed to describe the health and composition of an individual.

[0094] Gene-phenotype association tests The previous discussion encompasses definitions and contexts for genomic and phenotypic data used to describe a particular individual. Figure 1 describes at a high level how to determine genetic associations for a particular rare deleterious variant carrier status at a genetic resolution level for a particular phenotype. Here, we revisit this method in more detail through the introduction of genome-wide association studies and their implementation for both common and rare variants. In contrast to traditional genome-wide association studies for rare variants, the discussion is directed to a novel methodology with improved utility in which multiple rare variant summation tests for multiple specific genes, including at least one specific rare deleterious variant, are implemented to determine rare variant polygenic risk scores for multiple specific genes and specific phenotypes. This discussion focuses specifically on rare deleterious variants, due to the statistical possibility that rare deleterious variants are involved in severe early-onset genetic diseases.

[0095] FIG. 6 graphically contrasts genome-wide association studies for common variants versus rare variants. Graphs 602 and 604 illustrate genome-wide associations for cholesterol levels and common variants or rare variants, respectively. For both graphs 602 and 604, alternating colors indicate successive chromosomes. Graph 602 shows individual common variants as data points, whereas graph 604 shows individual genes (i.e., all rare variants for a particular gene condensed to a single carrier state as shown in model 122). The common variant genome-wide association study shown in graph 602 finds many locations in the genome that are strongly associated with cholesterol levels. Due to the fact that common variants close to each other are frequently correlated, the associations occur as peaks with many associated variants. As a result, it may require significant effort to detect any causal variants within a particular associated locus.

[0096] Graph 604 shows rare variant genome-wide association study results for 19,500 protein-coding genes in the human genome, with a small number of rare variants strongly associated. Due to their rarity, rare variant genome-wide association studies are less powerful than the common variant genome-wide association study shown in graph 602. Rare variants typically do not correlate with each other, and therefore the effect will not extend to nearby genes. Thus, significant genes appear isolated rather than clustered in graph 602. Variants that change protein sequences are specifically tested, which means that causal variants can be identified for significant genes identified by rare variant genome-wide association studies.

[0097] While Figure 6 illustrates multiple genetic associations at the genome-wide level, the discussion now shifts to genetic association tests for specific individual variants (i.e., the resulting non-random association between specific individual common variant allele dosages and the phenotypic response to the presence in the associated gene of a specific individual common variant allele dosage value), or for specific aggregations of variants (i.e., the resulting non-random association in a cohort between the carrier status of each of the associated genes and the phenotypic response to the presence in the associated gene of one or more rare pathogenic variants).

[0098] FIG. 7 illustrates a genetic association test for a particular individual common variant or a particular aggregated rare variant at gene resolution, respectively. For a single common variant genetic association test 702 from a genotyping array, each individual variant can be evaluated to determine the relationship between the number of minor alleles for each individual common variant (box plots 704, 706, and 708 correspond to phenotypic value distributions for zero, one, or two minor alleles, respectively) and the phenotype of interest (e.g., cholesterol). A linear regression model can be constructed between the average values ​​for each minor allele count value to obtain a p-value and a beta coefficient that measure the effective strength score of the association for a particular variant and phenotypic response.

[0099] In contrast, for rare variants obtained from exome sequencing, individual variant-by-variant testing is not useful because the incidence rate for each rare variant is too low. To test the genetic association of rare variants, a summation test 722 is performed on aggregated variants within a particular gene, as shown in chart 122. Compared to the single common variant genetic association test 702, the aggregated rare variant summation test 722 measures the strength of association between the carrier status (box plots 724 and 726 correspond to the phenotype value distribution for non-carriers and carriers, respectively) for at least one rare variant within a gene and the phenotype of interest. In addition to filtering for rarity, only deleterious variants are included as determined by a pathogenicity measurement threshold.

[0100] In some embodiments of the disclosed technology, the phenotypic response of interest is a categorical clinical diagnosis. Thus, the effective strength score for the resulting non-random association for a gene is determined for the categorical clinical diagnosis phenotype as a beta coefficient for the carrier status, where the beta coefficient is determined using the logistic regression component of the sum score.

[0101] Optimizing rare variant summation tests Here, we revisit the concept of individually determined rarity thresholds (i.e., maximum allele counts) and pathogenicity score thresholds for a particular gene. The diagram above describes a genetic association test of a particular multiple rare deleterious variant for a particular gene on a particular phenotype, where the maximum rarity threshold and minimum pathogenicity score threshold are specific to a particular gene and may differ from the maximum rarity threshold and minimum pathogenicity score threshold for different particular genes.

[0102] In the embodiment described below with respect to FIG. 8, "optimization" refers to identifying a maximum threshold for rarity and a minimum score for pathogenicity, since the target goal is to identify the combination of parameters that results in the most significant t-test statistic for a particular gene. However, in other embodiments, the target goal for optimization may differ due to the genetic and medical context. Those skilled in the art will recognize that the interpretation of a test statistic (such as a p-value) is context-dependent, and thus the target value or range of values ​​associated with significance is influenced by various parameters, such as sample size, statistical power, variation in effect size per gene or per variant, error variance, base rate of true effect, and payoff calculations associated with possible summation test results.

[0103] Please note that this application uses the terms "rarity threshold", "maximum allele count", "allele count threshold", and "AC" interchangeably. Please note that this application also uses the terms "pathogenicity threshold", "pathogenicity score threshold", "minimum pathogenicity score threshold", and "PST" interchangeably.

[0104] The discussion then turns to the disclosed technique for determining the particular combination of allele counts and pathogenicity score thresholds that has the most significant p-value and using that particular combination as the optimal combination. Each p-value for each combination of allele counts and pathogenicity score thresholds within the multiple combinations of allele counts and pathogenicity score thresholds is corrected to account for multiple testing (i.e., t-tests for genetic association testing of specific multiple rare deleterious variants for a specific gene on a specific phenotype performed within a grid search across multiple combinations of allele counts and pathogenicity score thresholds).

[0105] FIG. 8 is a schematic diagram of a method for optimizing rare variant summation testing. Database 802 includes genomic and phenotypic data corresponding to a group of individuals belonging to cohort i. Cohort i 804 includes N individuals who have been genetically sequenced and tested for multiple phenotypic measurements. The sequencing data for the individuals in cohort i 804 includes data corresponding to multiple variant carrier states. A particular gene may carry three known possible rare deleterious variants (e.g., variant A, variant B, and variant C), of which an individual may be a carrier (i.e., an individual carries each variant and thus they are a carrier of that variant) or a non-carrier (i.e., an individual does not carry each variant and thus they are not a carrier of that variant). One of skill in the art will recognize that variant A, variant B, and variant C are provided as illustrative examples and that a gene may carry any number of rare deleterious variants. The disclosed technology defines a variant as "rare" if the variant occurrence in a population is below a predetermined threshold, and "detrimental" if the pathogenicity of the variant is predicted to have a measured effect on the effect when the gene is expressed as a protein above the predetermined threshold. Each predetermined threshold for both rarity and pathogenicity is gene-specific (i.e., each individual gene will correspond to an individually determined rarity threshold and pathogenicity threshold that may differ from another gene).

[0106] In addition to the genomic data obtained from genome and exome sequencing for each individual in cohort i804, phenotypic data are also available for each individual (i.e., individual N has x n (having a measured value of ).

[0107] To determine optimal thresholds for rarity and pathogenicity, a grid search is performed to search across a space containing all possible pathogenicity score thresholds (PST) and maximum allele count thresholds (AC) 822. Within the grid 822, each individual test (corresponding to a PST value m and an AC value N) is a t-test 824. For each t-test 824, the null hypothesis states that the difference between the phenotypic values ​​of individuals who are carriers of the variant group and the phenotypic values ​​of individuals who are not carriers of the variant group is not statistically significant. The alternative hypothesis states that the difference between the phenotypic values ​​of individuals who are carriers of the variant group and the phenotypic values ​​of individuals who are not carriers of the variant group is statistically significant. The t-tests are performed iteratively across the grid for each possible combination of PST and AC values ​​to obtain a p-value for each t-test at a unequal significance level. Within multiple PST and AC value combinations and associated p-values ​​842, subcohorts (PSTs, AC ... m , A.C. n ) 844 contains genomic and phenotypic data for N individuals who meet a specified PST value m and AC value N (i.e., if variant B does not meet a specified threshold for pathogenicity or allele count for a particular t-test, then the carrier status for variant B is not considered within the aggregated binary carrier status variable).

[0108] The PST and AC combination corresponding to the most significant p-value is output as the optimal PST and AC combination for the particular gene 862. In some embodiments of the disclosed technology, following optimization on the determined optimal PST and AC combination 862, a rare variant summation test is performed. Thus, the method of performing an optimized summation test is performed such that all optimized summation tests for a particular gene start with the optimal combination of maximum allele count and minimum pathogenicity score threshold that maximizes the significance of the summation test effect of rare pathogenic variants in the particular gene on a particular phenotype, and is taken as the parameters for the summation test for the particular gene.

[0109] Determining pathogenicity based on protein structure Thus, the discussion so far has described a rare variant polygenic risk score model that is configured to aggregate risk across genes based on whether an individual carries a rare deleterious variant in each gene. The rare variant polygenic risk score models illustrated in Figures 1, 6, 7, and 8 are implemented with the optimal combination of maximum allele count and minimum pathogenicity score threshold that maximizes the significance of the sum test effect of rare pathogenic variants in a particular gene on a particular phenotype.

[0110] The discussion now turns to a description of a system for pathogenicity determination based on protein structure, where one embodiment of the disclosed technology implements a convolutional neural network pathogenicity classifier configured to determine the pathogenicity of a variant. The pathogenicity score output corresponding to a particular variant for a particular gene may be implemented as a pathogenicity score threshold for the particular gene in a particular t-test 824. A plurality of m pathogenicity score thresholds may be tested within a grid search 822 to determine an optimal minimum pathogenicity score threshold, such as an optimal pathogenicity score threshold 862.

[0111] FIG. 9 is a flow diagram 900 illustrating the process of a system for determining pathogenicity of a variant. In step 902, a sequence accessor 904 of the system accesses a reference amino acid sequence and an alternative amino acid sequence. In 912, a 3D structure generator 914 of the system generates a 3D protein structure of the reference amino acid sequence. In some embodiments, the 3D protein structure is a homology model of a human protein. In one embodiment, the so-called SwissModel homology modeling pipeline provides a public repository of predicted human protein structures. In another embodiment, the so-called HHpred homology modeling uses a tool called Modeller to predict the structure of a target protein from a template structure.

[0112] Proteins are represented by a collection of atoms and their coordinates in 3D space. Amino acids can have various atoms such as carbon atoms, oxygen (O) atoms, nitrogen (N) atoms, and hydrogen (H) atoms. The atoms can be further classified as side chain atoms and backbone atoms. Backbone carbon atoms can include alpha carbon (Cα) atoms and beta carbon (Cβ) atoms.

[0113] In step 922, the coordinate classifier 924 of the system classifies the 3D atomic coordinates of the 3D protein structure on an amino acid basis. In one embodiment, the amino acid-by-amino acid classification includes assigning the 3D atomic coordinates to 21 amino acid categories (including stop or gap amino acid categories). In one example, the amino acid-by-amino acid classification of the alpha carbon atoms can list each of the alpha carbon atoms under each of the 21 amino acid categories. In another example, the amino acid-by-amino acid classification of the beta carbon atoms can list each of the beta carbon atoms under each of the 21 amino acid categories.

[0114] In yet another example, the amino acid classification of oxygen atoms may list each oxygen atom under each of the 21 amino acid categories. In yet another example, the amino acid classification of nitrogen atoms may list each nitrogen atom under each of the 21 amino acid categories. In yet another example, the amino acid classification of hydrogen atoms may list each hydrogen atom under each of the 21 amino acid categories.

[0115] Those skilled in the art will appreciate that in various implementations, the amino acid classification can include a subset of the 21 amino acid categories and a subset of the different atomic elements.

[0116] This application uses voxels and voxelization as non-limiting examples, and one of ordinary skill in the art will appreciate that in various embodiments, different formats for arranging and processing data, such as features of various dimensions, vectors, arrays, etc., may be used alternatively or in combination.

[0117] In step 932, a voxel grid generator 934 of the system instantiates a voxel grid. The voxel grid can have any resolution, e.g., 3x3x3, 5x5x5, 7x7x7, etc. The voxels within the voxel grid can be any size, e.g., 1 Angstrom (Å) on each side, 2 Å on each side, 3 Å on each side, etc. Those skilled in the art will appreciate that since voxels are cubes, these exemplary dimensions refer to cubic dimensions. Those skilled in the art will also appreciate that these exemplary dimensions are non-limiting and voxels can have any cubic dimensions.

[0118] In step 942, the system's voxel grid centerer 944 centers the voxel grid at the amino acid level on the reference amino acid that experiences the target variant. In one embodiment, the voxel grid is centered on the atomic coordinates of a particular atom of the reference amino acid that experiences the target variant, for example, the 3D atomic coordinates of the alpha carbon atom of the reference amino acid that experiences the target variant.

[0119] Distance Channel A voxel in the voxel grid can have multiple channels (or features). In one embodiment, a voxel in the voxel grid has multiple distance channels (e.g., 21 distance channels each for one of the 21 amino acid categories (including a stop or gap amino acid category). In step 952, a distance channel generator 954 of the system generates a distance channel for each amino acid for the voxels in the voxel grid. A distance channel is generated independently for each of the 21 amino acid categories.

[0120] For example, consider the alanine (A) amino acid category. Further, for example, consider a voxel grid that is 3×3×3 in size and has 27 voxels. Then, in one embodiment, the alanine distance channel includes 27 distance values ​​for each of the 27 voxels in the voxel grid. The 27 distance values ​​in the alanine distance channel are measured from the center of each of the 27 voxels in the voxel grid to their respective nearest atoms in the alanine amino acid category.

[0121] In one example, the alanine amino acid category contains only alpha carbon atoms, and therefore the nearest atom is the alanine alpha carbon atom that is closest to each of the 27 voxels in the voxel grid. In another example, the alanine amino acid category contains only beta carbon atoms, and therefore the nearest atom is the alanine beta carbon atom that is closest to each of the 27 voxels in the voxel grid.

[0122] In yet another example, the alanine amino acid category contains only oxygen atoms, and therefore the nearest atom is the alanine oxygen atom that is closest to each of the 27 voxels in the voxel grid. In yet another example, the alanine amino acid category contains only nitrogen atoms, and therefore the nearest atom is the alanine nitrogen atom that is closest to each of the 27 voxels in the voxel grid. In yet another example, the alanine amino acid category contains only hydrogen atoms, and therefore the nearest atom is the alanine hydrogen atom that is closest to each of the 27 voxels in the voxel grid.

[0123] Similar to the alanine distance channel, distance channel generator 954 generates distance channels (i.e., sets of distance values ​​per voxel) for each of the remaining amino acid categories. In other embodiments, distance channel generator 954 generates distance channels for only a subset of the 21 amino acid categories.

[0124] In other embodiments, the selection of the closest atom is not limited to a particular atom type, i.e., the closest atom to a particular voxel within the target amino acid category is selected without regard to the atomic element of the closest atom, and the distance value of the particular voxel is calculated for inclusion in the distance channel of the target amino acid category.

[0125] In yet another embodiment, the distance channels are generated on an atomic element basis. Instead of or in addition to having a distance channel for the amino acid category, distance values ​​can be generated for the atomic element category regardless of the amino acid to which the atom belongs. For example, consider that the atoms of the amino acids in the reference amino acid sequence span the seven atomic elements, carbon, oxygen, nitrogen, hydrogen, calcium, iodine, and sulfur. The voxels in the voxel grid are then configured to have seven distance channels, such that each of the seven distance channels has 27 per-voxel distance values ​​that specify the distance to the nearest atom only in the corresponding atomic element category. In another embodiment, distance channels for only a subset of the seven atomic elements can be generated. In yet another embodiment, the atomic element category and distance channel generation are performed for the same atomic element, e.g., alpha carbon (C α ) atoms and β carbon (C β ) can be further hierarchized into variants of atoms.

[0126] In yet other embodiments, distance channels can be generated on an atom type basis, for example, a distance channel for only side chain atoms and a distance channel for only backbone atoms.

[0127] The closest atom can be searched within a predefined maximum scan radius (e.g., 6 Angstroms (Å)) from the voxel center, and multiple atoms may be closest to the same voxel in the voxel grid.

[0128] The distance is calculated between the 3D coordinates of the voxel center and the 3D atomic coordinates of the atom, and the distance channel is generated using a voxel grid centered on the same location (e.g., centered on the 3D atomic coordinates of the alpha carbon atom of the reference amino acid experiencing the target variant).

[0129] The distance can be a Euclidean distance. The distance can also be parameterized by the atomic size (or influence of the atom) (e.g., by using the Lennard-Jones potential and / or van der Waals atomic radius of the atom in question). The distance value can also be normalized by the maximum scan radius or by the maximum observed distance value of the farthest nearest atom in the target amino acid category or target atomic element category or target atom type category. In some embodiments, the distance between a voxel and an atom is calculated based on the polar coordinates of the voxel and the atom. The polar coordinates are parameterized by the angle between the voxel and the atom. In one embodiment, this angle information is used to generate an angle channel for the voxel (i.e., independent of the distance channel). In some embodiments, the angle between the nearest atom and a neighboring atom (e.g., a backbone atom) can be used as a feature to be encoded with the voxel.

[0130] Reference allele and alternative allele channels A voxel in the voxel grid may also have a reference allele and an alternative allele channel. In step 962, a one-hot encoder 964 of the system generates a reference one-hot encoding of a reference amino acid in the reference amino acid sequence and an alternative one-hot encoding of an alternative amino acid in the alternative amino acid sequence. The reference amino acid experiences a target variant. The alternative amino acid is a target variant. The reference amino acid and the alternative amino acid are located at the same positions in the reference amino acid sequence and the alternative amino acid sequence, respectively. The reference amino acid sequence and the alternative amino acid sequence have the same position-by-position amino acid composition, with one exception. The exception is a position having a reference amino acid in the reference amino acid sequence and an alternative amino acid in the alternative amino acid sequence.

[0131] In step 972, a concatenator 974 of the system concatenates the distance channel for each amino acid with the reference and alternative one-hot encodings. In another embodiment, the concatenator 974 concatenates the distance channel for each atomic element with the reference and alternative one-hot encodings. In yet another embodiment, the concatenator 974 concatenates the distance channel for each atom type with the reference and alternative one-hot encodings.

[0132] In step 982, the system's runtime logic 984 processes the concatenated per-amino acid / per-atom element / per-atom type distance channels and the reference and alternative one-hot encodings through a pathogenicity classifier (pathogenicity determination engine) to determine the pathogenicity of the target variant, which is then inferred as the pathogenicity determination of the underlying nucleotide variant that generates the target variant at the amino acid level. The pathogenicity classifier is trained using a labeled dataset of benign and pathogenic variants, for example, using a backpropagation algorithm. Further details regarding the labeled dataset of benign and pathogenic variants, as well as exemplary architectures and training of the pathogenicity classifier, can be found in commonly owned U.S. Patent Application Nos. 16 / 160,903, 16 / 160,986, 16 / 160,968, and 16 / 407,149.

[0133] Figure 10 illustrates an exemplary processing architecture 1000 of the pathogenicity classifier 900 according to one embodiment of the disclosed technology. The processing architecture 1000 includes a cascade of processing modules 1006, 1010, 1014, 1018, 1022, 1010, 1030, 1034, 1038, and 1042, each of which may include 1D convolution (1x1x1 CONV), 3D convolution (3x3x3 CONV), ReLU nonlinearity, and batch normalization (BN). Other examples of processing modules include a fully-connected (FC) layer, a dropout layer, a flattening layer, and a final softmax layer that generates exponentially normalized scores for target variants belonging to benign and pathogenic classes. In Figure 10, "64" indicates the number of convolution filters applied by a particular processing module. In Figure 10, the size of the input voxel 1002 is 15 x 15 x 15 x 8. Figure 10 also shows the volumetric dimensions of each of the intermediate inputs 1004, 1008, 1012, 1016, 1020, 1024, 1028, 1032, 1036, and 1040 generated by the processing architecture 1000.

[0134] 11 illustrates an exemplary computer system 1100 that can be used to implement the disclosed techniques. The computer system 1100 includes at least one central processing unit (CPU) 1144 that communicates with a number of peripheral devices via a bus subsystem 1142. These peripheral devices can include, for example, a storage subsystem 1110 including memory devices and a file storage subsystem 1136, user interface input devices 1138, user interface output devices 1148, and a network interface subsystem 1146. The input and output devices enable user interaction with the computer system 1100. The network interface subsystem 1146 provides an interface to external networks, including interfaces to corresponding interface devices in other computer systems.

[0135] In one embodiment, the pathogenicity classifier 1000 is communicatively linked to a memory subsystem 1110 and a user interface input device 1138 .

[0136] The user interface input devices 1138 can include pointing devices such as a keyboard, a mouse, a trackball, a touch pad, or a graphics tablet, a scanner, a touch screen integrated into a display, audio input devices such as a voice recognition system and a microphone, as well as other types of input devices. In general, use of the term "input device" is intended to include all possible types of devices and manners for inputting information into the computer system 1100.

[0137] The user interface output devices 1148 may include a display subsystem, a printer, a fax machine, or a non-visual display such as an audio output device. The display subsystem may include a flat panel device such as an LED display, a cathode ray tube (CRT), a liquid crystal display (LCD), a projection device, or some other mechanism for creating a visible image. The display subsystem may also provide non-visual displays such as an audio output device. In general, use of the term "output device" is intended to include all possible types of devices and manners for outputting information from computer system 1100 to a user or to another machine or computer system.

[0138] The storage subsystem 1110 stores programming and data constructs that provide the functionality of some or all of the modules and methods described herein. These software modules are generally executed by the processor 1148.

[0139] The processor 1148 can be a graphics processing unit (GPU), a field programmable gate array (FPGA), an application specific integrated circuit (ASIC), and / or a coarse-grained reconfigurable architecture (CGRA). The processor 1148 can be hosted by a deep learning cloud platform such as Google Cloud Platform™, Xilinx™, and Cirrascale™. Examples of processors 1148 include Google's Tensor Processing Unit (TPU)™, rackmount solutions such as the GX4 Rackmount Series™, GX11 Rackmount Series™, NVIDIA DGX-1™, Microsoft' Stratix V FPGA™, Graphcore's Intelligent Processor Unit (IPU)™, Qualcomm's Zeroth Platform™ with Snapdragon processors™, NVIDIA's Volta™, NVIDIA's DRIVE PX™, NVIDIA's JETSON TX1 / TX11 MODULE™, Intel's Nirvana™, Movidius VPU™, Fujitsu DPI™, ARM's DynamicIQ™, IBM TrueNorth™, Lambda GPU Server with Testa V100s™, and others.

[0140] The memory subsystem 1122 used in the storage subsystem 1110 may include multiple memories including a main random access memory (RAM) 1132 for storing instructions and data during program execution, and a read only memory (ROM) 1134 in which fixed instructions are stored. The file storage subsystem 1136 may provide persistent storage for program and data files and may include a hard disk drive, associated removable media, a CD-ROM drive, an optical drive, or a removable media cartridge. Modules that implement the functionality of a particular embodiment may be stored by the file storage subsystem 1136 in the storage subsystem 1110 or in another machine accessible by the processor.

[0141] Bus subsystem 1140 provides a mechanism for allowing the various components and subsystems of computer system 1100 to communicate with each other as intended. Although bus subsystem 1142 is shown generally as a single bus, alternative implementations of the bus subsystem may use multiple buses.

[0142] The computer system 1100 itself can be of various types, including a personal computer, a portable computer, a workstation, a computer terminal, a network computer, a television, a mainframe, a server farm, a loosely distributed set of loosely networked computers, or any other data processing system or user device. Due to the ever-changing nature of computers and networks, the description of computer system 1100 shown in Figure 11 is intended only as a specific example for purposes of illustrating a preferred embodiment of the present invention. Many other configurations of computer system 1100 can have more or fewer components than the computer system shown in Figure 11.

[0143] 12 illustrates one embodiment 1200 of determining a final pathogenicity score. In action 1202, in one embodiment, the pathogenicity classifier 1000 generates a first pathogenicity score for a first alternative amino acid that is the same as the first reference amino acid. In action 1212, in one embodiment, the pathogenicity classifier 1000 generates a second pathogenicity score for a second alternative amino acid that is different from the first reference amino acid. In action 1222, in one embodiment, the final pathogenicity score for the second alternative amino acid is the second pathogenicity score for the second alternative amino acid.

[0144] In another alternative, the final pathogenicity score for the second alternative amino acid is based on a combination of the first pathogenicity score and the second pathogenicity score. In the first alternative in 1222a, in one embodiment, the final pathogenicity score for the second alternative amino acid is the ratio of the second pathogenicity score to the sum of the first pathogenicity score and the second pathogenicity score. In the second alternative in 1222b, in one embodiment, the final pathogenicity score for the second alternative amino acid is determined by subtracting the first pathogenicity score from the second pathogenicity score.

[0145] Those skilled in the art will appreciate that other current and future artificial intelligence, machine learning, and deep learning models, datasets, and learning techniques can be incorporated into the disclosed variant pathogenicity classifiers without departing from the spirit of the disclosed technology.

[0146] Correction for multiple testing The above discussion describes an embodiment of the disclosed technology including a method for performing an optimized summation test in which the strength of association between genes associated with a phenotype and the contribution of rare variants to the phenotypic response is quantified, the summation test being optimized for a particular combination of allele counts and pathogenicity score thresholds with the most significant p-values. The optimized summation test is based on multiple nested t-tests that maximize the separation between carriers and non-carriers of at least one rare deleterious variant in a particular gene. In this embodiment of the disclosed technology, the multiple t-tests performed within the optimization method ensure the need to correct for multiple testing within each gene.

[0147] The discussion then turns to a description of the correction of the most significant p-value corresponding to the optimal combination of PST and AC values.

[0148] FIG. 13 is a flow diagram 1300 of a process for correcting for multiple testing within each gene. Following false discovery rate (FDR) correction, the p-values ​​are further corrected for multiple testing to account for summation testing optimization for both PST and AC. The FDR corrected p-values ​​1302 are separated into value ranges. In one embodiment of the disclosed technology, the following algorithm for multiple testing correction 1300 is followed: In other embodiments of the disclosed technology, other forms of multiple testing correction may be used. For values ​​in the range (0, 1e-5], the FDR corrected p-values ​​1302 undergo a Benjamini-Hochberg FDR correction 1310.

[0149] For values ​​in the range (0.01, 1], the FDR corrected p-value 1302 is subjected to a permutation testing step 1312. In step 1312, in one embodiment, a count of the total number of permutations generated (N) is set to zero. In other embodiments, N may be set to any predetermined starting value. Following step 1312, step 1322 involves generating 1,000 permutations of phenotypic labels and setting N equal to N+1000. In the next step 1342 Then, a summation test is performed for each permutation of the data and the fraction p of permutations that result in a more significant p-value than the original observed data is counted. If N is greater than 100 / p, proceed to step 1362. At step 1362, stop and output p. The stopping point at step 1362 ensures N<10,000 to maintain computational efficiency; that is, the number of permutations (N) is limited to less than 10,000, thereby preventing reaching computationally infeasible or computationally inefficient or computationally expensive counts.

[0150] For values ​​in the range (1e-5, 0.01], the FDR corrected p-values ​​1302 are subjected to a permutation testing step 1314. In step 1314, 10,000 permutations of phenotype labels are generated. Following step 1314, step 1324 involves fitting a generalized extreme value distribution to the absolute value of the test statistic for each respective sum test from the permuted data. In the next step 1344, a corrected p-value is estimated from the area under the curve of the fitted distribution. That is, in step 1362, an initial p-value is output, whereas in step 1344, a corrected p-value is estimated.

[0151] Correction of phenotypic values ​​for covariates The discussion so far encompasses quantifying the strength of association of genes associated with a phenotype and the contribution of rare variants to the phenotypic response, where the effective strength scores of multiple genes contributing to the rare polygenic risk score are determined by a Boolean carrier state variable, and the rare polygenic risk score model (i.e., summation test) is optimized for the particular optimal combination of allele counts and pathogenicity score thresholds with the most significant p-values. Here, we revisit the concept of phenotypic measures to introduce additional methods of optimization in which the phenotypic measures are covariate-corrected and drug use-corrected.

[0152] Note that the application uses the terms "drug use," "drug use pattern," "drug category use," and "drug category use pattern" interchangeably. Also, note that the application uses the term "confounder" to describe a covariate associated with both one or more other covariates and the outcome of interest (e.g., tobacco use is a confounder associated with drug use patterns and phenotypic measures of cardiovascular disease). One of skill in the art will understand that the disclosed technology may be applied to any disease, such as oncological diseases, diabetes, etc.

[0153] FIG. 14 is a schematic diagram of a method 1400 of correcting phenotypic values ​​for covariates. In step 1402, a quantitative phenotypic dataset is pruned into a non-redundant dataset. In step 1404, a plurality of drugs of interest are grouped into a set of drug categories. In step 1406, each phenotype in the non-redundant quantitative phenotypic dataset is corrected for drug use for each drug category in the set of drug categories (e.g., statins, blood pressure medications, etc.). In step 1408, each corrected phenotype is inverse rank normal transformed. In step 1410, the corrected phenotype is further corrected for a plurality of confounding factors (e.g., age, sex, genetic principal components, diet, smoking status, etc.). With regard to confounding factors, a confounding factor is a factor that confounds the relationship between a test drug and an outcome. Some studies may allow for supplemental medications or supplemental treatments, such as psychotherapy, during the study. These are confounding factors because the outcome may be attributable to the supplemental agent rather than the drug being tested. As an example, pain medication trials often ignore participants' use of non-pharmacological pain medications, such as splints, creams, massages, hot baths, and chiropractor manipulations. Responses to these additional study treatments can significantly confound the results of the study. Common confounders are participant attributes, such as body mass index, smoking status, age at onset of disease, socioeconomic status, educational status, and the extent of support networks. Life events are also potential confounders. They can cause significant changes in participants' moods and symptom levels, and this applies to even the smallest events, especially in the relationship between participants and assessors. All known confounders that are not matched between experimental and control groups can be taken into account in the statistical analysis.

[0154] The discussion now turns to a detailed description of one embodiment of the disclosed technology configured to predict phenotypic shift in response to the use of multiple drugs for multiple phenotypes of a cohort of individuals at a first and second time point, where the individuals in the cohort are grouped based on permuted drug use patterns measured at two time points (e.g., no drug at the first and second time points, start taking drug between the first and second time points, stop taking drug between the first and second time points, take drug at the first and second time points). In contrast to other strategies for covariate correction commonly used in biostatistics and population genetics studies, the use of multiple time points rather than a single time point results in lower error variance and a more accurate approximation of the effect of a particular confounding factor on the phenotypic shift.

[0155] FIG. 15 is a schematic diagram depicting predicting phenotypic shift in response to use of multiple drugs for multiple phenotypes. A database 1502 includes genomic and phenotypic data corresponding to a population belonging to a cohort j. In step 1522, covariate measurements for multiple confounders, phenotypic measurements for multiple phenotypes, and drug use patterns for multiple drug categories are accessed from the database 1502 at two time points t1 and t2. Steps 1542, 1562, 1582, and 1592 are performed for each phenotype in the multiple phenotypes. Those skilled in the art will understand that the disclosed technology can be extended to any number of time points, such as three, four, five, etc.

[0156] In step 1542, the phenotypic measurements are covariate corrected for multiple confounders at t1 and t2 to generate confounder-corrected values ​​by fitting a first regression model. In step 1562, the first regression model is used to determine delta δ for the confounder-corrected values. Thus, covariate correction is implemented by removing and regressing the covariate measurements by fitting the first regression model. Steps 1582 and 1592 are performed for each drug category in the multiple drug categories. In step 1582, δ is used to predict phenotypic shift in response to use of a particular drug category by fitting a second regression model. In step 1592, the confounder-corrected values ​​are further corrected for drug use for a particular drug category to generate confounder-corrected, drug category use-corrected values ​​for a particular phenotype. Thus, phenotypic shift prediction is implemented by fitting a second regression model that models phenotypic shift for each of the drug use patterns. The second regression model has a binary indicator variable for each of the drug use patterns.

[0157] FIG. 16 is a schematic diagram showing an experimental setup for obtaining drug use patterns and phenotypic data for a particular cohort. A database 1602 includes genomic and phenotypic data corresponding to a group of individuals belonging to a cohort j. The database 1602 is segmented into two sub-cohorts, each including an individual 1622 not taking a particular drug category z and an individual 1624 taking a particular drug category z at a first time point t1. The sub-cohort including the individual 1622 not taking a particular drug category z can be further segmented into a smaller sub-cohort including an individual 1642 who did not take a particular drug category z continuously and an individual 1644 who started taking a particular drug category z after the first time point t1 and before the second time point t2, respectively.

[0158] The sub-cohorts including individuals 1646 taking a particular drug category z may be further segmented into smaller sub-cohorts including individuals 1646 who stopped taking the particular drug category z and individuals 1648 who continued to take the particular drug category z after a first time point t1 and before a second time point t2.

[0159] Phenotype value trend graph 1662 illustrates example phenotype values ​​measured at t1 and t2 for sub-cohort 1642 (i.e., patients who were not taking drug category z at t1 or t2), where a particular phenotype value starts high on the y-axis at t1 and increases slightly at t2. Phenotype value trend graph 1664 illustrates example phenotype values ​​measured at t1 and t2 for sub-cohort 1644 (i.e., patients who were not taking drug category z at t1 but started taking drug category z before t2), where a particular phenotype value starts high on the y-axis at t1 and decreases at t2. Phenotype value trend graph 1666 illustrates example phenotype values ​​measured at t1 and t2 for sub-cohort 1646 (i.e., patients who were taking drug category z at t1 but stopped taking drug category z before t2), where a particular phenotype value starts low on the y-axis at t1 and increases at t2. The phenotype value trend graph 1668 illustrates an example of phenotype values ​​measured at t1 and t2 for subcohort 1648 (i.e., patients who were taking drug category z at t1 and t2), where a particular phenotype value begins to decline on the y-axis at t1 and remains low at t2.

[0160] In some embodiments of the disclosed technology, the drug effect on the phenotype value (i.e., the phenotype shift in response to drug category z) can be learned by model 1682, and the delta (δ) for the phenotype value at t1 and t2 can be calculated using the model 1682. y) is modeled as a regression with a β coefficient for drug category z. This model can be tested for significance where the null hypothesis states that the β coefficient for drug category z is equal to zero and the alternative hypothesis states that the β coefficient for drug category z is not equal to zero.

[0161] In some embodiments of the disclosed technology, the regression models from steps 1542 and 1582 that are constructed based on data from database 1602 are constructed via the following protocol.

[0162] Phenotypic shift effects of 34 drug categories (e.g., statins, NSAIDs, opioids, etc.) were estimated from a cohort of participants who had their quantitative phenotype values ​​and corresponding covariates measured at t1 and t2, 5 years apart. For each quantitative phenotype, covariate-adjusted values ​​at t1 and t2 were generated by regression analysis removing all covariates except for drug use (e.g., age, sex, genetic principal components, diet, smoking status, etc.). The difference between the two time points is calculated as follows:

[0163]

number

[0164] For each drug category X, all individuals in the cohort are divided into four groups according to their drug use at t1 and t2 (as described in sub-cohorts 1642, 1644, 1646, and 1648), and a binary indicator variable is introduced for each group that encodes whether the individual belongs to that group or not. 1) Individuals who were not taking drug X at both t1 and t2 (indicator variable: X 00 ) 2) Individuals who started taking drug X between time points t1 and t2 (indicator variable: X 01 ) 3) Individuals who stopped taking drug X between time points t1 and t2 (indicator variable: X 10 ) 4) Individuals who were taking drug X at both time points t1 and t2 (indicator variable: X 11 )

[0165] To determine which drugs have a significant effect on phenotype Y, a forward selection stepwise regression of the following form is fitted, repeating across all drug categories:

[0166]

number

[0167] In the above equation, the term β t t models the effect of time elapsed between t1 and t2 (t=t2-t1 for each individual in the cohort), and the term

[0168]

number

[0169] After the set D of relevant drugs has been determined, their individual effects can be estimated jointly by fitting the following regression:

[0170]

number

[0171] The raw values ​​for phenotype Y at t1 across all individuals in the cohort can be corrected as follows:

[0172]

number

[0173] X (d) is a binary indicator variable that encodes whether an individual was taking drug X at the first visit, t1. After adjusting for drug use, the value

[0174]

number

[0175] In some embodiments of the disclosed technology, the drug use correction is implemented by fitting a third regression model, which includes drug use correction of the phenotypic measure for the first time point based on a first binary indicator independent variable for a first drug use pattern of starting to take a drug between the first and second time points, a second binary indicator independent variable for a second drug use pattern of not taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a particular drug at the first time point. A fourth regression model is fitted, which includes drug use correction of the phenotypic measure for the second time point based on a third binary indicator independent variable for a third drug use pattern of stopping to take a drug between the first and second time points, a fourth binary indicator independent variable for a fourth drug use pattern of taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a particular drug at the second time point.

[0176] A rank-based inverse normal transformation may be applied to the drug-use-corrected phenotypic measures for the first and second time points to generate normalized drug-use-corrected phenotypic measures for the first and second time points. The normalized drug-use-corrected phenotypic measures are then covariate-corrected for the first and second time points to generate covariate-corrected normalized drug-use-corrected phenotypic measures for the first and second time points. The covariate-corrected normalized drug-use-corrected phenotypic measures may be used to generate a rare variant polygenic risk score, where the measures for the phenotype of interest are corrected for phenotypic shifts in response to covariates and drug use patterns.

[0177] FIG. 17 illustrates a number of phenotypes corresponding to a cardiovascular disease patient X 1700. The measured phenotype value for the cardiovascular disease patient X 1700 will be influenced by a number of confounding factors. Examples of these confounding factors include demographics 1702 (e.g., age, sex, ethnicity, etc.), psychosocial factors 1722 (mental health related confounders, socioeconomic class, etc.), tobacco and alcohol consumption 1742, drug use 1704 (e.g., drugs not included in formal drug effect analysis such as illegal drugs), body mass measurements 1724 (e.g., body mass index, body fat percentage, abdominal fat density, etc.), or diet and lifestyle factors 1744 (e.g., dietary restrictions, eating habits, exercise, etc.). One skilled in the art will recognize that these are non-limiting confounding factors and that there are numerous additional confounding factors that may affect the measured phenotype value.

[0178] Performance measurements as objective indicators of nonobviousness and inventive step The preceding discussion encompasses multiple embodiments of the disclosed technology for rare variant summation testing that includes a polygenic risk score model constructed for rare deleterious variant carrier status in a particular gene for a particular phenotype, where model parameters (i.e., pathogenicity score threshold and maximum allele count) are optimized via nested t-tests (i.e., grid search), and where both phenotypic values ​​are covariate-corrected for multiple covariates and drug-use-corrected for multiple drug-use patterns. The discussion now turns to performance results of various embodiments of the disclosed technology.

[0179] Figure 18 is a graph 1800 quantifying the total number of significant gene-phenotype pairs identified for different types of summation tests. The total number of significant gene-phenotype pairs identified across 90 phenotypes for different types of summation tests shows that the combined tests for rare LoF and deleterious missense variants prioritized by PrimateAI-3D outperform other approaches. As a negative control, the number of significant genotype-phenotype pairs for summation tests performed on synonymous variants is also shown.

[0180] Figure 19 is a set of graphs illustrating rare deleterious variants that affect disease severity and age of onset identified by pathogenicity classifier PrimateAI-3D. Graph 1902 corresponds to the LDLR gene and its association with disorders of LDL cholesterol and lipoprotein metabolism. For carriers of rare missense variants in the LDLR gene, LDL cholesterol levels (y-axis) are positively correlated with PrimateAI-3D percentile scores (x-axis). Hereinafter, PrimateAI-3D scores refer to PrimateAI-3D percentile scores normalized from 0 to 1 within each gene to facilitate comparison between genes. Graph 1904 illustrates that PrimateAI-3D scores predict the age of onset of dyslipidemia for carriers of rare missense variants in the LDLR gene.

[0181] Graph 1922 corresponds to the PCSK9 gene, a downregulator of LDLR. LDL cholesterol levels of carriers of missense variants in PCSK9, a downregulator of LDLR, are negatively correlated with Primate AI-3D scores. Graph 1924 shows that LDL cholesterol levels of carriers increase with age at a similar rate to non-carriers, but are on average lower than non-carriers of the same age group. Graph 1942 corresponds to the GCK gene and shows that HbA1c levels correlate with Primate AI-3D for carriers of rare missense variants in the GCK gene. Graph 1944 shows that HbA1c levels of carriers increase with age at a similar rate to non-carriers, but on average, reach the prediabetes threshold earlier in the life of carriers.

[0182] Figure 20 is a graph 2000 of the mean absolute Spearman correlation of different pathogenicity scores with phenotypic values. Graph 2000 shows the mean absolute Spearman correlation of different pathogenicity scores with phenotypic values ​​in a benchmark set of 34 gene-phenotype pairs compared to a theoretical upper bound set by carriers of the same variants.

[0183] Figure 21 is a heatmap 2100 comparison of rare deleterious variants and common genome-wide association study variants. For each pair of phenotypes, the heatmap shows the statistical significance of the overlap between the GWAS genes associated with the phenotype on the y-axis and the rare variant genes associated with the phenotype on the x-axis.

[0184] Figure 22 is a further comparison of rare deleterious variants with common genome-wide association study variants. Graph 2200A shows the relationship between variant effect size and variant allele frequency for LoF, deleterious missense (PrimateAI-3D>0.5), and cryptic splice (SpliceAI>0.2) variants. Synonymous and benign missense (PrimateAI-3D<0.5) variants are shown as negative controls. Dot size is proportional to the square root of the number of variants in each dot. Graph 2200B contains the distribution of the effect size of the rarest LoF variant divided by the effect size of the lead GWAS variant for the same gene. Histograms of both high pLI (top plot) and low pLI (bottom plot) genes are shown. Red vertical lines indicate the mean of each distribution. Graph 2200C shows the percentage of positively mapped GWAS genes with rare variant enrichment (nominal p-value ≦0.05) stratified by genes with different GWAS significance. Results are plotted separately for high pLI genes (pLI>0.5) and low pLI genes (pLI<0.5). High confidence GWAS genes have lead variants with p-values ​​<10-100 and strong LD (r 2 >=0.9). The dashed line indicates the expected percentage of genes that would meet the nominal p-value threshold (p<0.05) by chance.

[0185] Figure 23 illustrates the cholesterol pathway and total cholesterol distribution across all individuals in the UK Biobank cohort. Example 2302 includes the cholesterol pathway, with the genes in the rare variant PRS model overlaid on the picture. For each gene, the numbers and arrows indicate the effect size and direction of the effect. Graph 2304 shows the distribution of total cholesterol across all individuals in the UK Biobank cohort, with the average effect size of common and rare variants shown.

[0186] Figure 24 includes measures of rare variant PRS performance. Graph 2402 shows the average correlation of common variant PRS (AF>1%), rare variant PRS (AF<0.1%), and composite PRS with 72 phenotypes, using 90% of UKBB for PRS training and 10% for testing. Box plots show the distribution of Pearson correlations (R) between traits. Graph 2404 is a comparison of PRS performance (average correlation with phenotype) using different significance cutoffs for determining inclusion of genes or loci in the model. Data point size corresponds to the number of genes or loci used in each PRS model.

[0187] 25 is a graph 2500 of enrichment of PRS outliers. Graph 2500 shows enrichment of PRS outliers in individuals who are phenotypic outliers. The x-axis shows the z-score threshold used to define phenotypic outliers, and the y-axis shows enrichment of PRS outliers in phenotypic outlier individuals relative to the baseline population.

[0188] 26 is a graph 2600 of PRS outliers for quantitative phenotypes. Graph 2600 shows PRS outliers for 56 quantitative phenotypes. Bars indicate statistical significance for either common variant PRS or rare variant PRS to describe individuals at the 99th and 99.9th percentile phenotype outlier thresholds.

[0189] FIG. 27 includes graphs 2702 and 2704 of normalized total cholesterol distributions from two separate cohorts. Graph 2702 shows normalized total cholesterol distributions from UK Biobank individuals in the bottom 0.5% (low), 0.5-99.5% (medium), and top 0.5% (high) groups for the cholesterol PRS, with 50% of the cohort used to train the PRS and the remainder used for validation. Graph 2704 shows normalized total cholesterol distributions from MGB individuals in the bottom 0.5% (low), 0.5-99.5% (medium), and top 0.5% (high) groups for the cholesterol PRS, with 50% of the cohort used to train the PRS and the remainder used for validation.

[0190] Figure 28 includes graphs 2802 and 2804 comparing rare variant PRS outliers and phenotypes between two separate cohorts. Graph 2802 shows rare variant PRS outlier tests for 17 quantitative phenotypes measured in both cohorts. PRS and phenotype outliers were defined as individuals in the top and bottom 0.5% of the population. Graph 2804 shows the correlation between rare variant PRS and phenotype for the 17 quantitative phenotypes.

[0191] Figure 29 includes graphs 2902 and 2904 illustrating performance results by ethnicity. Graph 2902 shows normalized cholesterol distribution for low and high PRS groups in the MGB cohort, with results shown for both white and non-white individuals. Graph 2904 shows a comparison of the mean z-score distance between the low PRS (<0.5%) and high PRS (>99.5%) groups for each of the 11 phenotypes in the MGB, for both white and non-white individuals.

[0192] FIG. 30 is a table 3000 comparing effect sizes and frequencies for common PRS variants and rare PRS genes used for normal cholesterol levels.

[0193] Figure 31 is a graph 3100 of the average percentage of phenotypic variance explained by different pathogenicity scoring methods for a set of 34 gene-phenotype pairs, which were selected based on their enrichment for rare missense and LoF variants. The percentage of phenotypic variance explained is calculated as the squared Spearman correlation between each score from each scoring method and the phenotypic value of carriers of the corresponding variant. The method is compared in graph 3100 to two theoretical lower limits: a theoretical upper limit calculated by using the phenotypic value of carriers of the same missense variant, and carriers of the same synonymous variant and a random score.

[0194] Figure 32 is a heat map 3200 of enrichment of rare variants in GWAS genes for all pairwise comparisons between quantitative and clinical phenotypes. For each pair of phenotypes, the heat map 3200 shows the statistical significance of the ranking of the most significant subset of GWAS genes of one phenotype (y-axis) by their enrichment for rare deleterious variants affecting the second phenotype (x-axis). After removing the effect of all significant GWAS variants from each phenotype and performing regression analysis, a summation test for rare deleterious variants is calculated. The heat map 2100 shows the subset of phenotypes enriched on the main diagonal of the heat map 3200.

[0195] Figure 33 includes a graph 3300 of the distribution of the absolute value of the ratio between the average effect size of singleton LoF variants and the effect size of the most significant GWAS variant for the same gene. Graph 3300A shows the distribution of the absolute value of the ratio between the average effect size of singleton LoF variants and the effect size of the most significant GWAS variant for the same gene for high pLI genes from graph 2200B. Graph 3300B shows the distribution of the absolute value of the ratio between the average effect size of singleton LoF variants and the effect size of the most significant GWAS variant for the same gene for low pLI genes from graph 2200B.

[0196] Figure 34 shows the number of variants per individual in any of the genes from the gene-phenotype pairs that were significant at a 5% false discovery rate. Graph 3400 includes 500 genes from 1031 gene-phenotype pairs that passed the significance threshold. The average number of variants per individual is 2.03 with a standard deviation of 1.51.

[0197] Figure 35 shows a comparison 3500 of effect sizes in training vs. test data splits. Rare variant testing was performed using 50% of the cohort as a training set and the other 50% as a test set. In graph 3500, p-values ​​<10 in the training set -7 For genes with , the effect sizes from the training set versus the test set are compared.

[0198] Figure 36 includes a graph 3600 comparing common variant PRS subsets with rare variant PRS subsets. Graph 3600A shows common and rare PRS correlations when a locus was required to be significant in both common and rare tests. PRS was trained on 90% and tested on the remainder. Graph 3600B shows common and rare PRS correlations to maximum by number of training samples. Graph 3600C shows effect sizes and ratios of variance explained from rare variant PRS models using PrimateAI-3D versus rare variant PRS models built without PrimateAI-3D.

[0199] Terms The disclosed technology, particularly the provisions disclosed in this section, may be implemented as a system, method, or product. One or more features of the embodiments may be combined with the base embodiment. Non-mutually exclusive embodiments are taught as combinable. One or more features of the embodiments may be combined with other embodiments. The present disclosure will periodically inform users of these options. The omission from some embodiments of the enumeration of repeating these options should not be interpreted as limiting the combinations taught in the preceding section. These descriptions are incorporated herein by reference in each of the following implementations.

[0200] One or more embodiments and provisions of the disclosed technology, or elements thereof, can be implemented in the form of a computer product including a non-transitory computer-readable storage medium with computer usable program code for performing the illustrated method steps. Furthermore, one or more embodiments and provisions of the disclosed technology, or elements thereof, can be implemented in the form of an apparatus including a memory and at least one processor coupled to the memory and operative to perform the illustrated method steps. Furthermore, in another aspect, one or more embodiments and provisions of the disclosed technology, or elements thereof, can be implemented in the form of a means for performing one or more of the method steps described herein, which means can include (i) a hardware module, (ii) a software module running on one or more hardware processors, or (iii) a combination of hardware and software modules, any of (i)-(iii) implementing the particular technology described herein, and the software module is stored in a computer-readable storage medium (or multiple such media).

[0201] The clauses described in this section can be combined as features. For the sake of brevity, combinations of features are not listed separately and are not repeated for each base set of features. The reader will understand how features specified in the clauses described in this section can be easily combined with sets of basic features specified as embodiments in other sections of this application. These clauses are not meant to be mutually exclusive, exhaustive, or restrictive, and the disclosed technology is not limited to these clauses, but rather encompasses all possible combinations, modifications, and variations within the scope of the claimed technology and its equivalents.

[0202] Other implementations of the provisions described in this section may include a non-transitory computer-readable storage medium storing instructions executable by a processor to perform any of the provisions described in this section. Yet another implementation of the provisions described in this section may include a system including a memory and one or more processors operable to execute instructions stored in the memory to perform any of the provisions described in this section.

[0203] The inventors disclose the following provisions: 1. A computer-implemented method for predicting phenotypic shift in response to use of multiple drugs for multiple phenotypes in a cohort of individuals with multiple confounding factors, comprising: For a cohort of individuals and for a first and second time point, accessing phenotypic measurements for a plurality of phenotypes; Access to covariate measurements for multiple confounding factors; Accessing drug use patterns for multiple drugs; and For each phenotype, covariate-correcting the phenotypic measures for the first and second time points based on the covariate measures, thereby generating covariate-corrected phenotypic measures for the first and second time points; determining a delta based on the difference between the covariate-adjusted phenotypic measurements for the first and second time points; For each of the drug use patterns, fitting a second regression model that uses delta to predict phenotypic shifts in response to use of the multiple drugs on the covariate-adjusted phenotypic measures; and drug-use correcting the phenotypic measures for the first and second time points based on the phenotypic shift, thereby generating a drug-use corrected phenotypic measure for the first time point. 2. The computer-implemented method of clause 1, wherein the plurality of confounding factors includes age, sex, genetic principal components, diet, and smoking status. 3. The computer-implemented method of clause 1, wherein the covariate correction is implemented by removing the covariate measurements by fitting a first regression model and performing a regression analysis. 4. The computer-implemented method of clause 1, wherein the phenotypic shift prediction is implemented by fitting a second regression model that models the phenotypic shift for each of the drug use patterns. 5. The computer-implemented method of clause 4, wherein the second regression model is a forward selective stepwise regression model that iteratively predicts delta by successively and cumulatively including phenotypic shifts for each of the drug use patterns. 6. Drug use patterns: Being drug-free at time points 1 and 2; commencing taking the medication between a first time point and a second time point; ceasing to take the drug between the first time point and the second time point; taking the medication at the first and second time points. 7. The computer-implemented method of claim 6, wherein the second regression model has a binary indicator independent variable for each of the drug use patterns. 8. The computer-implemented method of clause 1, wherein covariate correction, delta determination, phenotype shift prediction, and drug use correction are performed on a drug-by-drug basis for drugs in the plurality of drugs. 9. The computer-implemented method of claim 8, wherein a second regression model is iteratively fitted for each of the drugs. 10. The computer-implemented method of claim 8, further comprising grouping the drugs into a set of drug categories. 11. The computer-implemented method of claim 10, wherein covariate correction, delta determination, phenotype shift prediction, and drug use correction are performed for each drug category for drug categories in the set of drug categories. 12. The computer-implemented method of claim 11, wherein a second regression model is iteratively fitted to each of the drug categories. 13. The computer-implemented method of clause 1, wherein the second regression model further models phenotypic shift in response to the time elapsed between the first and second time points for each individual in the cohort of individuals. 14. The computer-implemented method of clause 1, wherein the second regression model further models phenotypic shifts in response to regression to the mean between the first and second time points. 15. The computer-implemented method of claim 8, wherein the second regression model is fitted jointly to the set of related drugs in the multiple drugs. 16. The computer-implemented method of clause 1, wherein the drug use correction is implemented by fitting a third regression model. 17. The computer-implemented method of clause 16, further comprising drug-use correcting the phenotype measure for the first time point based on a first binary indicator independent variable for a first drug use pattern of starting to take a drug between the first and second time points, a second binary indicator independent variable for a second drug use pattern of not taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the first time point. 18. The computer-implemented method of clause 1, wherein the drug use correction is implemented by fitting a fourth regression model. 19. The computer-implemented method of clause 18, further comprising drug-use correcting the phenotype measure for the second time point based on a third binary indicator independent variable for a third drug use pattern of ceasing taking the drug between the first and second time points, a fourth binary indicator independent variable for a fourth drug use pattern of taking the drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the second time point. 20. The computer-implemented method of clause 1, further comprising applying a rank-based inverse normal transformation to the drug-use corrected phenotype measures for the first and second time points, and generating normalized drug-use corrected phenotype measures for the first and second time points. 21. The computer-implemented method of clause 12, further comprising covariate-correcting the normalized drug-use corrected phenotypic measures for the first and second time points, and generating covariate-corrected normalized drug-use corrected phenotypic measures for the first and second time points. 22. The computer-implemented method of clause 12, further comprising using the covariate-corrected and normalized drug use-corrected phenotypic measures to generate a rare variant polygenic risk score. 23. The computer-implemented method of clause 1, wherein the plurality of phenotypes corresponds to a plurality of quantitative phenotypes. 24. The computer-implemented method of clause 23, wherein the quantitative phenotype in the plurality of quantitative phenotypes is a quantitative biomarker measurement. 25. The computer-implemented method of claim 23, further comprising pruning the plurality of quantitative phenotypes into a non-redundant set for use in covariate correction, delta determination, phenotype shift, and drug use correction. 26. The computer-implemented method of clause 25, wherein each pair of quantitative phenotypes in the non-redundant set has an absolute pairwise Pearson correlation below an upper threshold. 27. The computer-implemented method of claim 26, wherein the upper threshold is 0.95. 28. The computer-implemented method of clause 26, further comprising selecting, from among each group of redundant quantitative phenotypes in the plurality of quantitative phenotypes, the phenotype with the most samples for inclusion in the non-redundant set. 29. The computer-implemented method of clause 1, wherein the plurality of phenotypes corresponds to a plurality of categorical phenotypes. 30. The computer-implemented method of clause 29, wherein a categorical phenotype in the plurality of categorical phenotypes is a clinical diagnosis. 31. The computer-implemented method of clause 4, further comprising using a second regression model to detect drug-phenotype associations. 32. The computer-implemented method of claim 31, wherein the drug-phenotype associations include potential undesirable side effects and desired target effects. 33. A system comprising one or more processors coupled to a memory, the memory being loaded with computer instructions for predicting phenotypic shift in response to use of multiple drugs for multiple phenotypes in a cohort of individuals having multiple confounding factors, the instructions, when executed on the processor, For a cohort of individuals and for a first and second time point, accessing phenotypic measurements for a plurality of phenotypes; Access to covariate measurements for multiple confounding factors; Accessing drug use patterns for multiple drugs; and For each phenotype, covariate-correcting the phenotypic measures for the first and second time points based on the covariate measures, thereby generating covariate-corrected phenotypic measures for the first and second time points; determining a delta based on the difference between the covariate-adjusted phenotypic measurements for the first and second time points; The system implements actions including: for each of the drug use patterns, using the delta to predict a phenotypic shift in response to use of a plurality of drugs relative to the covariate-corrected phenotypic measures; and drug-use-correcting the phenotypic measures for the first time point based on the phenotypic shift, thereby generating drug-use-corrected phenotypic measures for the first and second time points. 34. The system of clause 33, wherein the multiple confounding factors include age, sex, genetic components, diet, and smoking status. 35. The system of clause 33, wherein the covariate correction is implemented by removing the covariate measurements by fitting a first regression model and performing regression analysis. 36. The system of clause 33, wherein the phenotypic shift prediction is implemented by fitting a second regression model that models phenotypic shift for each of the drug use patterns. 37. The system of clause 36, wherein the second regression model is a forward selective stepwise regression model that iteratively predicts delta by successively and cumulatively including phenotypic shifts for each of the drug use patterns. 38. Your drug use patterns are: Being drug-free at time points 1 and 2; commencing taking the medication between a first time point and a second time point; ceasing to take the drug between the first time point and the second time point; and taking the medication at the first and second time points. 39. The system of clause 38, wherein the second regression model has a binary indicator independent variable for each of the drug use patterns. 40. The system of clause 33, wherein covariate correction, delta determination, phenotype shift prediction, and drug use correction are performed on a drug-by-drug basis for drugs in the plurality of drugs. 41. The system of claim 40, wherein a second regression model is iteratively fitted for each of the drugs. 42. The system of clause 40, further implementing an action including grouping the drugs into a set of drug categories. 43. The system of clause 42, wherein covariate correction, delta determination, phenotype shift prediction, and drug use correction are performed for each drug category for drug categories in the set of drug categories. 44. The system of claim 43, wherein a second regression model is iteratively fitted to each of the drug categories. 45. The system of clause 33, wherein the second regression model further models phenotypic shift in response to the time elapsed between the first and second time points for each individual in the cohort of individuals. 46. ​​The system of clause 33, wherein the second regression model further models phenotypic shifts in response to regression to the mean between the first and second time points. 47. The system of claim 40, wherein the second regression model is fitted jointly to the set of related drugs in the plurality of drugs. 48. The system of clause 33, wherein the drug use correction is implemented by fitting a third regression model. 49. The system of clause 48, further implementing an action including drug use correcting the phenotype measure for the first time point based on a first binary indicator independent variable for a first drug use pattern of starting to take a drug between the first and second time points, a second binary indicator independent variable for a second drug use pattern of not taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the first time point. 50. The system of clause 33, wherein the drug use correction is implemented by fitting a fourth regression model. 51. The system of clause 50, further implementing an action including drug use correcting the phenotype measure for the second time point based on a third binary indicator independent variable for a third drug use pattern of ceasing taking the drug between the first and second time points, a fourth binary indicator independent variable for a fourth drug use pattern of taking the drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the second time point. 52. The system of clause 33, further implementing the action of applying a rank-based inverse normal transformation to the drug use corrected phenotype measures for the first and second time points, and generating normalized drug use corrected phenotype measures for the first and second time points. 53. The system of clause 44, further implementing the actions of covariate-correcting the normalized drug use corrected phenotype measures for the first and second time points, and generating covariate-corrected normalized drug use corrected phenotype measures for the first and second time points. 54. The system of clause 44, further implementing an action including generating a rare variant polygenic risk score using the covariate-corrected and normalized drug use-corrected phenotypic measures. 55. The system of clause 33, wherein the plurality of phenotypes corresponds to a plurality of quantitative phenotypes. 56. The system of clause 55, wherein the quantitative phenotype in the plurality of quantitative phenotypes is a quantitative biomarker measurement value. 57. The system of clause 55, further implementing an action including pruning the plurality of quantitative phenotypes into a non-redundant set for use in covariate correction, delta determination, phenotypic shift prediction, and drug use correction. 58. The system of clause 57, wherein each pair of quantitative phenotypes in the non-redundant set has an absolute pairwise Pearson correlation below an upper threshold. 59. The system of claim 58, wherein the upper threshold is 0.95. 60. The system of clause 58, further implementing an action including selecting, from among each group of redundant quantitative phenotypes in the plurality of quantitative phenotypes, the phenotype with the most samples for inclusion in a non-redundant set. 61. The system of clause 33, wherein the plurality of phenotypes corresponds to a plurality of categorical phenotypes. 62. The system of clause 61, wherein a categorical phenotype in the plurality of categorical phenotypes is a clinical diagnosis. 63. The system of clause 36, further implementing an action including using a second regression model to detect drug-phenotype associations. 64. The system of clause 63, wherein the drug-phenotype associations include potential undesirable side effects and desired target effects. 65. A non-transitory computer-readable storage medium storing computer program instructions for predicting phenotypic shift in response to use of multiple drugs for multiple phenotypes in a cohort of individuals with multiple confounding factors, the instructions, when executed on a processor, comprising: For a cohort of individuals and for a first and second time point, accessing phenotypic measurements for a plurality of phenotypes; Access to covariate measurements for multiple confounding factors; Accessing drug use patterns for multiple drugs; and For each phenotype, covariate-correcting the phenotypic measures for the first and second time points based on the covariate measures, thereby generating covariate-corrected phenotypic measures for the first and second time points; determining a delta based on the difference between the covariate-adjusted phenotypic measurements for the first and second time points; A non-transitory computer-readable storage medium implementing a method including: using delta to predict, for each of the drug use patterns, a phenotypic shift in response to use of a plurality of drugs relative to the covariate-corrected phenotypic measures; and drug-use-correcting the phenotypic measures for the first and second time points based on the phenotypic shift prediction, thereby generating drug-use-corrected phenotypic measures for the first and second time points. 66. The non-transitory computer-readable storage medium of clause 65, wherein the plurality of confounding factors includes age, sex, genetic basis, diet, and smoking status. 67. The non-transitory computer-readable storage medium of clause 65, wherein the covariate correction is implemented by removing the covariate measurements by fitting a first regression model and performing regression analysis. 68. The non-transitory computer-readable storage medium of clause 65, wherein the phenotype shift prediction is implemented by fitting a second regression model that models the phenotype shift for each of the drug use patterns. 69. The non-transitory computer-readable storage medium of clause 68, wherein the second regression model is a forward selective stepwise regression model that iteratively predicts delta by successively and cumulatively including phenotypic shifts for each of the drug use patterns. 70. Drug use patterns: Being drug-free at time points 1 and 2; commencing taking the medication between a first time point and a second time point; ceasing to take the drug between the first time point and the second time point; and taking the medication at the first and second times. 71. The non-transitory computer-readable storage medium of clause 70, wherein the second regression model has a binary indicator independent variable for each of the drug use patterns. 72. The non-transitory computer-readable storage medium of clause 65, wherein covariate correction, delta determination, phenotype shift prediction, and drug use correction are performed on a drug-by-drug basis for drugs in the plurality of drugs. 73. The non-transitory computer-readable storage medium of clause 72, wherein a second regression model is iteratively fitted for each of the drugs. 74. The non-transitory computer-readable storage medium of clause 72, implementing a method further comprising grouping drugs into a set of drug categories. 75. The non-transitory computer-readable storage medium of clause 74, wherein covariate correction, delta determination, phenotype shift prediction, and drug use correction are performed on a drug category-by-drug category basis for drug categories within the set of drug categories. 76. The non-transitory computer-readable storage medium of clause 75, wherein a second regression model is iteratively fitted to each of the drug categories. 77. The non-transitory computer-readable storage medium of clause 65, wherein the second regression model further models phenotypic shift in response to the time elapsed between the first and second time points for each individual in the cohort of individuals. 78. The non-transitory computer-readable storage medium of clause 65, wherein the second regression model further models phenotypic shifts in response to regression to the mean between the first time point and the second time point. 79. The non-transitory computer-readable storage medium of clause 72, wherein the second regression model is fitted jointly to the set of related drugs in the plurality of drugs. 80. The non-transitory computer-readable storage medium of clause 65, wherein the drug use correction is implemented by fitting a third regression model. 81. The non-transitory computer-readable storage medium of clause 80, implementing a method further comprising drug-use correcting the phenotype measure for the first time point based on a first binary indicator independent variable for a first drug use pattern of starting to take a drug between the first and second time points, a second binary indicator independent variable for a second drug use pattern of not taking a drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the first time point. 82. The non-transitory computer-readable storage medium of clause 65, wherein the drug use correction is implemented by fitting a fourth regression model. 83. The non-transitory computer readable storage medium of clause 82, implementing the method further comprising drug use correcting the phenotype measure for the second time point based on a third binary indicator independent variable for a third drug use pattern of ceasing taking the drug between the first and second time points, a fourth binary indicator independent variable for a fourth drug use pattern of taking the drug at the first and second time points, and a drug-specific binary indicator independent variable encoding whether the individual was taking a specific drug at the second time point. 84. The non-transitory computer-readable storage medium of clause 65, implementing a method further comprising applying a rank-based inverse normal transformation to the drug-use corrected phenotype measures for the first and second time points, and generating normalized drug-use corrected phenotype measures for the first and second time points. 85. The non-transitory computer-readable storage medium of clause 76, implementing a method further comprising covariate-correcting the normalized drug-use corrected phenotypic measures for the first and second time points, and generating covariate-corrected normalized drug-use corrected phenotypic measures for the first and second time points. 86. The non-transitory computer-readable storage medium of clause 76, implementing a method further comprising generating a rare variant polygenic risk score using the covariate-corrected and normalized drug use-corrected phenotypic measures. 87. The non-transitory computer-readable storage medium of clause 65, wherein the plurality of phenotypes corresponds to a plurality of quantitative phenotypes. 88. The non-transitory computer-readable storage medium of clause 87, wherein the quantitative phenotype in the plurality of quantitative phenotypes is a quantitative biomarker measurement value. 89. The non-transitory computer-readable storage medium of clause 87, implementing a method further comprising pruning the plurality of quantitative phenotypes into a non-redundant set for use in covariate correction, delta determination, phenotypic shift prediction, and drug use correction. 90. The non-transitory computer-readable storage medium of clause 89, wherein each pair of quantitative phenotypes in the non-redundant set has an absolute pairwise Pearson correlation below an upper threshold. 91. The non-transitory computer-readable storage medium of clause 90, wherein the upper threshold is 0.95. 92. The non-transitory computer-readable storage medium of clause 90, implementing a method further comprising selecting, from among each group of redundant quantitative phenotypes in the plurality of quantitative phenotypes, a phenotype with the most samples for inclusion in a non-redundant set. 93. The non-transitory computer-readable storage medium of clause 65, wherein the plurality of phenotypes corresponds to a plurality of categorical phenotypes. 94. The non-transitory computer-readable storage medium of clause 93, wherein a categorical phenotype in the plurality of categorical phenotypes is a clinical diagnosis. 95. The non-transitory computer-readable storage medium of clause 68, implementing a method further comprising using a second regression model to detect drug-phenotype associations. 96. The non-transitory computer-readable storage medium of clause 95, wherein the drug-phenotype associations include potential undesirable side effects and desired target effects.

Claims

1. 1. A computer-implemented method for correcting phenotypic measures for use in generating a polygenic risk score for rare variants, comprising: For a cohort of individuals and for the first and second time points: accessing phenotypic measurements for a plurality of phenotypes; Access to covariate measures for multiple confounders; Accessing drug use patterns for multiple drugs; and For each phenotype, covariate-correcting the phenotypic measures for the first and second time points based on the covariate measures, thereby generating covariate-corrected phenotypic measures for the first and second time points; determining a delta based on the difference between the covariate-corrected phenotypic measurements for the first and second time points; for each of the drug use patterns, fitting a second regression model using the delta to predict a phenotypic shift in response to use of the plurality of drugs relative to the covariate-adjusted phenotypic measure; drug-use correcting the phenotypic measures for the first and second time points based on the phenotypic shift, thereby generating a drug-use corrected phenotypic measure for the first time point; normalizing the drug-use-corrected phenotypic measures for the first and second time points to generate normalized drug-use-corrected phenotypic measures for the first and second time points; covariate-correcting the normalized drug-use-corrected phenotypic measures for the first and second time points to generate covariate-corrected normalized drug-use-corrected phenotypic measures for the first and second time points; and generating a polygenic risk score for rare variants using the covariate-corrected normalized drug-use-corrected phenotypic measures.

2. The computer-implemented method of claim 1 , wherein the plurality of confounding factors comprises age, sex, genetic principal components, diet, and smoking status.

3. The computer-implemented method of claim 1 , wherein the covariate correction is implemented by performing a regression analysis removing the covariate measurements by fitting a first regression model.

4. 2. The computer-implemented method of claim 1, wherein the phenotype shift prediction is implemented by fitting a second regression model that models phenotype shift for each of the drug use patterns.

5. The drug use pattern is: being drug-free at the first and second time points; commencing taking the medication between the first time point and the second time point; ceasing to take the medication between the first time point and the second time point; and taking the medication at the first and second time points.

6. 2. The computer-implemented method of claim 1, wherein the covariate correction, the delta determination, the phenotype shift prediction, and the drug use correction are performed on a drug-by-drug basis for drugs in the plurality of drugs.

7. 7. The computer-implemented method of claim 6, wherein the second regression model is iteratively fitted to each of the drugs.

8. The computer-implemented method of claim 6 , further comprising grouping the drugs into a set of drug categories.

9. 2. The computer-implemented method of claim 1, wherein the second regression model further models phenotypic shift in response to time elapsed between the first and second time points for each individual in the cohort of individuals.

10. 2. The computer-implemented method of claim 1, wherein the second regression model further models a phenotypic shift in response to regression to the mean between the first time point and the second time point.

11. 7. The computer-implemented method of claim 6, wherein the second regression model is fitted jointly to a set of related drugs in the plurality of drugs.

12. 2. The computer-implemented method of claim 1, wherein normalizing the drug-use-corrected phenotypic measures comprises applying a rank-based inverse normal transformation to the drug-use-corrected phenotypic measures for the first and second time points.

13. The computer-implemented method of claim 1 , wherein the plurality of phenotypes corresponds to a plurality of quantitative phenotypes.

14. 14. The computer-implemented method of claim 13, wherein a quantitative phenotype in the plurality of quantitative phenotypes is a quantitative biomarker measurement.

15. 14. The computer-implemented method of claim 13, further comprising pruning the plurality of quantitative phenotypes to a non-redundant set for use in the covariate correction, the delta determination, the phenotype shift, and the drug use correction.

16. The computer-implemented method of claim 1 , wherein the plurality of phenotypes corresponds to a plurality of categorical phenotypes.

17. 1. A system comprising one or more processors coupled to a memory, the memory loaded with computer instructions for correcting phenotypic measurements for use in generating a polygenic risk score for rare variants, the instructions, when executed on the processor, performing: For a cohort of individuals and for the first and second time points: accessing phenotypic measurements for a plurality of phenotypes; Access to covariate measures for multiple confounders; Accessing drug use patterns for multiple drugs; and For each phenotype, covariate-correcting the phenotypic measures for the first and second time points based on the covariate measures, thereby generating covariate-corrected phenotypic measures for the first and second time points; determining a delta based on the difference between the covariate-corrected phenotypic measurements for the first and second time points; using the delta to predict, for each of the drug use patterns, a phenotypic shift in response to use of the plurality of drugs relative to the covariate-adjusted phenotypic measure; drug-use correcting the phenotypic measure for the first time point based on the phenotypic shift, thereby generating drug-use corrected phenotypic measures for the first and second time points; normalizing the drug-use-corrected phenotypic measures for the first and second time points to generate normalized drug-use-corrected phenotypic measures for the first and second time points; covariate-correcting the normalized drug-use-corrected phenotypic measures for the first and second time points to generate covariate-corrected normalized drug-use-corrected phenotypic measures for the first and second time points; generating a polygenic risk score for rare variants using the covariate-corrected normalized drug use-corrected phenotypic measure.

18. 1. A non-transitory computer-readable storage medium storing computer program instructions for correcting phenotypic measurements for use in generating a polygenic risk score for rare variants, the instructions, when executed on a processor, comprising: For a cohort of individuals and for the first and second time points: accessing phenotypic measurements for a plurality of phenotypes; Access to covariate measures for multiple confounders; Accessing drug use patterns for multiple drugs; and For each phenotype, covariate-correcting the phenotypic measures for the first and second time points based on the covariate measures, thereby generating covariate-corrected phenotypic measures for the first and second time points; determining a delta based on the difference between the covariate-corrected phenotypic measurements for the first and second time points; using the delta to predict, for each of the drug use patterns, a phenotypic shift in response to use of the plurality of drugs relative to the covariate-adjusted phenotypic measure; drug-use correcting the phenotypic measures for the first and second time points based on the phenotypic shift prediction, thereby generating drug-use corrected phenotypic measures for the first and second time points; normalizing the drug-use-corrected phenotypic measures for the first and second time points to generate normalized drug-use-corrected phenotypic measures for the first and second time points; covariate-correcting the normalized drug-use-corrected phenotypic measures for the first and second time points to generate covariate-corrected normalized drug-use-corrected phenotypic measures for the first and second time points; and generating a polygenic risk score for rare variants using the covariate-corrected normalized drug use-corrected phenotypic measures.

Images