Methods and systems for identifying functional non-coding variants

Abstract

The presently disclosed techniques are generally directed to facilitating the identification of functional rare variants (236), such as variants within promoter regions. In certain embodiments, promoter regions to be assessed for functional variants may be identified or localized using a TSS predictor model (530) integrating multiple data sources to accurately identify transcription start site (TSS) positions, which are used to localize promoter regions for a gene. In further aspects as discussed herein, allele specific expression is evaluated for genes of interest, as opposed to total gene expression at the gene of interest. This may be done by selecting a suitable allele-expression single nucleotide polymorphism (aeSNP) for the gene. In certain embodiments discussed herein, the aeSNP closest to the TSS determined for the gene is selected as the best or most suitable SNP for assessing allele-specific expression.

Description

METHODS AND SYSTEMS FOR IDENTIFYING FUNCTIONAL NON-CODING VARIANTSFIELD OF THE TECHNOLOGY DISCLOSED

[0001] The technology disclosed relates to processor-implemented approaches to identify functional, non-coding variants (e g., rare variants) within a limited or expansive genome. Such approaches may be implemented using general- or special-purpose computers and / or other digital data processing systems and corresponding data processing methods and techniques. In certain examples, the technology disclosed relates to processor-driven analytics for identifying variants within promoter regions that may impact gene expression.INCORPORATIONS

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

[0003] U.S. Provisional Patent Application No.: 63 / 395,774, titled “ARTIFICIAL INTELLIGENCE-BASED EPIGENETICS AT BASE RESOLUTION,” filed August 5, 2022;

[0004] U.S. Provisional Patent Application No.: 63 / 395,778, titled “ARTIFICIAL INTELLIGENCE-BASED DETECTION OF GENE CONSERVATION AND EXPRESSION PRESERVATION AT BASE RESOLUION,” filed August 5, 2022;

[0005] U.S. Provisional Patent Application No.: 62 / 903,700, titled “ARTIFICIAL INTELLIGENCE-BASED EPIGENETICS,” filed September 20, 2019;

[0006] Sundaram, L. et al. Predicting the clinical impact of human mutation with deep neural networks. Nat. Genet. 50, 1161-1170 (2018);

[0007] Jaganathan, K. et al. Predicting splicing from primary sequence with deep learning. Cell 176, 535-548 (2019);

[0008] US Provisional Patent Application No. 62 / 573,144, titled “TRAINING A DEEP PATHOGENICITY CLASSIFIER USING LARGE-SCALE BENIGN TRAINING DATA,” filed October 16, 2017;

[0009] US Provisional Patent Application No. 62 / 573, 149, titled “PATHOGENICITY CLASSIFIER BASED ON DEEP CONVOLUTIONAL NEURAL NETWORKS (CNNs),” filed October 16, 2017;

[0010] US Provisional Patent Application No. 62 / 573,153, titled “DEEP SEMISUPERVISED LEARNING THAT GENERATES LARGE-SCALE PATHOGENIC TRAINING DATA,” filed October 16, 2017;

[0011] US Provisional Patent Application No. 62 / 582,898, titled “PATHOGENICITY CLASSIFICATION OF GENOMIC DATA USING DEEP CONVOLUTIONAL NEURAL NETWORKS (CNNs),” filed November 7, 2017;

[0012] US Patent Application No. 16 / 160,903, titled “DEEP LEARNING-BASED TECHNIQUES FOR TRAINING DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed on October 15, 2018;

[0013] US Patent Application No. 16 / 160,986, titled “DEEP CONVOLUTIONAL NEURAL NETWORKS FOR VARIANT CLASSIFICATION,” filed on October 15, 2018;

[0014] US Patent Application No. 16 / 160,968, titled “SEMI- SUPERVISED LEARNINGFOR TRAINING AN ENSEMBLE OF DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed on October 15, 2018;

[0015] US Patent Application No. 16 / 407,149, titled “DEEP LEARNING-BASED TECHNIQUES FOR PRE-TRAINING DEEP CONVOLUTIONAL NEURAL NETWORKS,” filed May 8, 2019;

[0016] US Patent Application No. 17 / 232,056, titled “DEEP CONVOLUTIONAL NEURALNETWORKS TO PREDICT VARIANT PATHOGENICITY USING THREE-DIMENSIONAL (3D) PROTEIN STRUCTURES,” filed on April 15, 2021;

[0017] US Provisional Patent Application No. 63 / 175,495, titled “MULTI-CHANNEL PROTEIN VOXELIZATION TO PREDICT VARIANT PATHOGENICITY USING DEEPCONVOLUTIONAL NEURAL NETWORKS,” filed on April 15, 2021;

[0018] US Provisional Patent Application No. 63 / 175,767, titled “EFFICIENT VOXELIZATION FOR DEEP LEARNING,” filed on April 16, 2021); and

[0019] US Patent Application No. 17 / 468,411 , titled “ARTIFICIAL INTELLIGENCE¬BASED ANALYSIS OF PROTEIN THREE-DIMENSIONAL (3D) STRUCTURES,” filed on September 7, 2021.BACKGROUND

[0020] This section is intended to introduce the reader to various aspects of art that may be related to various aspects of the present disclosure, which are described and / or claimed below. This discussion is believed to be helpful in providing the reader with background information to facilitate a better understanding of the various aspects of the present disclosure. Accordingly, it should be understood that these statements are to be read in this light, and not as admissions of prior art.

[0021] An aspect of studying the human genome, as well as other genomes, relates to the characterization of the function or functions of the different genomic elements of a respective organism or population of organisms. To this end, one goal of studying and characterizing certain elements of a genome is finding associations between the genotype (i.e., underlying code for a given genomic element) and phenotype (i.e., how such a genomic element presents itself functionally, structurally, anatomically, physiologically, and so forth).

[0022] However, certain tasks within genomics may remain challenging. In the context of humans, each human has a unique genetic code, a large portion of which is shared. Despite the extensive amount of shared code between humans, the different respective genetic codes of humans contain nucleic acid sequence “variants” or “genetic variants” which may be present within limited number of individuals or within small groups of a larger population and which may, in some instances, lead to “outlier” effects on the molecular products stemming from the genome. In practice, such outlier effects may correspond to different versions of a molecule(e g., a protein) encoded by the variant sequence relative to what is statistically the common or normal version of the molecule. For example, an outlier effect may be attributable to a variant sequence that encodes a protein that differs by one amino acid in an otherwise same specific sequence. However, in other contexts the outlier effect may originate from a variant in the noncoding portion of the genome, such as a promoter or enhancer region which, impacts the expression of a nearby gene without itself encoding for the expressed product.

[0023] Such genetic variants will often be effectively neutral, with either no discernible difference in the expression of the encoded product or with a discernible difference in the expression of the encoded product, but where the difference has little or no phenotypic effect. Conversely, in other instances such genetic variations may be considered pathogenic, and may be associated with a negative phenotypic effect, such as a genetic-based disease or disorder.Though such pathogenic genetic variants have often been depleted from genomes by natural selection, rare variants especially have largely arisen in the human population recently without time for selection to act or the variant effect may not impact reproductive fitness. Thus, an ability to identify which genetic variants are likely to be pathogenic can facilitate the investigation of such genetic variants so as to gain an understanding of corresponding disease states, associated diagnostics, treatments, and / or cures.

[0024] In practice the majority of variants identified by genome-wide association studies of complex phenotypes are non-coding, and instead of determining the sequence of a downstream biomolecule, may perform a regulatory function with respect to the expression of a nearby genomic element. In such a non-coding context, it may be challenging to estimate the effect of such variations on phenotype.SUMMARY

[0025] The present techniques provide for a software platform (e.g., a software application implemented locally (e.g., on-premise) or in a distributed (e.g., cloud implementation) manner and that provides tools for analyzing genetic data to identify functional rare variants, including rare variants in non-coding regions, such as promoter regions. In certain aspects, such a software platform may include one or more tools that allow a user to use a transcription start site predictor model or tool that integrates multiple types of data (e.g., existing or canonical TSS annotations,Cap Analysis of Gene Expression (CAGE or CAGE-seq), histone modification data, transcription factor (TF) data, chromatin immunoprecipitation sequencing (ChlP-seq) data, assay for transposase-accessible chromatin with sequencing (ATAC-Seq) data, methylation data, sequence conservation data, and so forth) to derive a consensus TSS location or window within which the TSS is located for a respective gene. In certain embodiments, an output of the software platform may be an updated or customized annotation file that includes the consensus TSS location for one or more genes. In practice, the software platform may be generic with respect to a sequencing device used to generate sequence data. Technical effects associated with the creation and use of such a customized TSS annotation file include, but are not limited to, improved identification of functional variants within the promoter region as well as generation of improved datasets of rare promoter variants that may be used in training or using machine learning tools or other models directed toward predicting gene expression effects of variants within the promoter region of a gene.

[0026] In a further aspect, the TSS annotation data in such a file may be used, via the software platform or a downstream application or platform, to select an allelic expression single nucleotide polymorphism (aeSNP) that may in turn be used to identify variants in the promoter region that might impact gene expression, leading to either under- or over-expression of the gene product. Such an aeSNP may be useful in the context of analyzing allele specific gene expression (as opposed to total gene expression) for a given gene. As may be appreciated, allele specific gene expression, in this context, refers to measuring or observing either the absolute or quantitative gene expression of two different alleles at a given gene locus or at the relative expression (e.g., the ratio) of the differing alleles at the locus. By measuring allele specific expression, useful data may be obtained for individuals that are heterozygous at the gene locus of interest as the allele specific expression data allows comparison of the expression of one allele (e g., variant) versus the other allele (e g., a reference allele) within the same individual. In particular, analysis of allele specific expression in this context may provide additional validation of expression effects, may facilitate the discovery of such expression effects that may not be evident when analyzing only the total gene expression data, and / or may provide improved resolution relative to total expression analysis.

[0027] By way of example, in certain approaches an allelic expression single nucleotide polymorphism (aeSNP) may be selected for a gene for which an individual is heterozygous. A comparison may then be made across reads for each variant to assess which reads are associated with each allele. As discussed herein, an aeSNP may be selected that is the closest to the TSS, and thus may benefit from improved localization of the TSS. Such an allele specific expression approach represents a technical improvement over prior techniques based on used of total-gene expression data alone and allows identification of rare variants within the promoter region that impact gene expression and which would be missed based on analysis of total expression alone.

[0028] In addition, a cis-eQTL correction may be employed on allele specific expression measurements that further improves rare variant discovery. Such a correction controls for common variants that impact allelic imbalance to further highlight or emphasize rare variant effects in the promoter region. With this in mind, a further technical improvement includes incorporation of a cis-eQTL correction step by which the signal or effects of one or more common variants within the region are predicted and removed so as to improve the signal relative to background noise and to facilitate identification and / or selection of rare variants within the promoter region.

[0029] With the preceding in mind, in one embodiment a method is provided for identifying a transcription start site (TSS). In accordance with this embodiment, a genetic sequence annotation file containing prior TSS location annotations is accessed. A TSS for a gene is identified based on estimated TSS positions for the gene across a plurality of different tissue types or different cell types. Identifying the consensus TSS for the gene utilizes Cap Analysis of Gene Expression (CAGE) data, Cap Analysis of Gene Expression (CAGE) data with functional measurements to identify the consensus TSS. The consensus TSS for the gene is validated based on enrichment of promoter variants leading to outlier expression effects nearby that TSS as compared to other potential TSSes. The genetic sequence annotation file is updated with the consensus TSS to generate a customized TSS annotation file.

[0030] In a further embodiment, a method is provided for using allele-specific expression (ASE) for variant assessment. In accordance with this embodiment, an allele expression single nucleotide polymorphism (aeSNP) is selected that is present at a gene of an individual where theindividual is heterozygous at the gene. The aeSNP is selected based on location relative to a transcription start site (TSS). The gene expression associated with the different alleles of the aeSNP is compared and one or more variant-level biochemical or phenotype effects are identified based on the expression of the different alleles.BRIEF DESCRIPTION OF THE DRAWINGS

[0031] These and other features, aspects, and advantages of the present invention will become better understood when the following detailed description is read with reference to the accompanying drawings, in which like characters represent like parts throughout the drawings, wherein:

[0032] FIG. 1 shows one implementation of a processing pipeline that identifies rare variants that cause extreme levels of gene expression, in accordance with aspects of the present disclosure;

[0033] FIG. 2 illustrates one implementation of a processing system that generates causality scores for outlier variants that cause extreme levels of gene expression, in accordance with aspects of the present disclosure;

[0034] FIG. 3 depicts one implementation of a fitted causality model that determines causal relationships between the rare variants and the extreme levels of gene expression while controlling for a plurality of confounders, in accordance with aspects of the present disclosure;

[0035] FIG. 4 provides a graphical illustration of the impact of TSS position on functional variant discovery, in accordance with aspects of the present disclosure;

[0036] FIG. 5 illustrates that the CAGE-based approach yields a different TSS than existing annotations files for a subset of genes, in accordance with aspects of the present disclosure;

[0037] FIG. 6 depicts signal-to-nose (y-axis) for prior canonical annotation TSS locations and for TSS location refined using CAGE data, in accordance with aspects of the present disclosure;

[0038] FIG. 7 depicts an aspect of a TSS predictor tool as a visual representation that incorporates multiple TSS localization data sets, in accordance with aspects of the present disclosure;

[0039] FIG. 8 depicts precision-recall (PR) curves for both the training data and test data derived using the TSS predictor model, in accordance with aspects of the present disclosure;

[0040] FIG. 9 depicts bar graphs are comparing the predicted TSS based on the TSS predictor model to the CAGE derived TSS, in accordance with aspects of the present disclosure;

[0041] FIG. 10 illustrates that, when matching window sizes for promoter region analysis, the TSS predictor model outperforms CAGE-based TSS localization, in accordance with aspects of the present disclosure;

[0042] FIG. 11 depicts a bar plot illustrating the number of variants that would be discovered using the CAGE-derived TSS locations versus the TSS predictor model, in accordance with aspects of the present disclosure;

[0043] FIGS. 12A, 12B, and 12C depict analyses of outliers from CAGE-derived TSS, in accordance with aspects of the present disclosure;

[0044] FIGS. 13A, 13B, and 13C depict analyses of outliers from a TSS predictor model, in accordance with aspects of the present disclosure;

[0045] FIG. 14 depicts a bar graph demonstrating enrichment versus p-value where the p- values were derived from T-tests between normalized adjusted allele reference ratios for variant carriers relative to all other individuals, in accordance with aspects of the present disclosure;

[0046] FIG. 15 depicts a bar graph in which the outlier count is illustrated against different techniques for deriving a reference ratio, in accordance with aspects of the present disclosure;

[0047] FIG. 16 depicts a bar graph illustrating that genes having significant allele specific expression have a substantial overlap with total expression effects, in accordance with aspects of the present disclosure;

[0048] FIG. 17 depicts a bar graph in which data is binned based in the ratio of the counts of the closest aeSNP to the top aeSNP and present in the context of enrichment versus outlier threshold, in accordance with aspects of the present disclosure;

[0049] FIG. 18 depicts a bar graph in which data is binned based on the distance between the closest aeSNP to the TSS and the top aeSNP and presented in the context of enrichment versus outlier threshold, in accordance with aspects of the present disclosure;

[0050] FIG. 19 illustrates that observing higher coverage at the closest aeSNP to the TSS increases enrichments of rare promoter variants with significant allelic expression imbalance, in accordance with aspects of the present disclosure;

[0051] FIG. 20 depicts in bar graph form total gene expression (TE), allele specific expression (ASE), and intersection binned data in the context of enrichment and p-value, in accordance with aspects of the present disclosure;

[0052] FIG. 21 depicts intersection data in bar graph form and in the context of counts and p- value, in accordance with aspects of the present disclosure;

[0053] FIG. 22 depicts a bar graph illustrating that improvement in the outlier count was obtained after correcting for the effects of the top aeQTL, in accordance with aspects of the present disclosure;

[0054] FIG. 23 depicts a bar graph illustrating that improvement in the outlier count was obtained after correcting for the effects of multiple aeQTL, in accordance with aspects of the present disclosure;

[0055] FIG. 24 depicts a bar graph in data is binned based upon expression and in the context of counts of the number of unique genes with any significant variants versus p-values, in accordance with aspects of the present disclosure; and

[0056] FIG. 25 depicts an example computer system that can be used to implement various aspects of the technology disclosed.DETAILED DESCRIPTION

[0057] One or more specific embodiments will be described below. In an effort to provide a concise description of these embodiments, not all features of an actual implementation are described in the specification. It should be appreciated that in the development of any such actual implementation, as in any engineering or design project, numerous implementationspecific decisions must be made to achieve the developers’ specific goals, such as compliance with system-related and business-related constraints, which may vary from one implementation to another. Moreover, it should be appreciated that such a development effort might be complex and time consuming, but would nevertheless be a routine undertaking of design, fabrication, and manufacture for those of ordinary skill having the benefit of this disclosure.

[0058] As used herein, the terms “application”, “pipeline”, “module”, “model”, “program”, and so forth may refer to one or more sets of computer software instructions (e.g., computer programs and / or scripts) executable by one or more processors of a computing system to provide particular functionality. Computer software instructions can be written in any suitable programming languages, such as C, C++, C#, Pascal, Fortran, Perl, Python, MATLAB, SAS, SPSS, JavaScript, AJAX, and JAVA. Such computer software instructions can comprise an independent or integrated application with data input and data display modules. Alternatively, the disclosed computer software instructions can be classes that are instantiated as distributed objects. The disclosed computer software instructions can also be component software. Additionally, the disclosed applications or engines can be implemented in computer software, computer hardware, or a combination thereof.

[0059] The detailed description of various implementations will be better understood when read in conjunction with the accompanying figures. To the extent that the figures illustrate concepts and functional blocks, modules, or processes, the depicted blocks or processes are not necessarily indicative of an implementation solely in hardware, software, or firmware. That is, such functionality may or may not be embodied in a discrete software module or modules, hardware circuit or circuits, and so forth, or indeed may be embodied or implemented in a combination of hardware and executable software. Further, applications, programs, modules, and / or computational pipelines as used herein may be implemented as a single executable unit(e g., a stand-alone program) or may be implemented as a groups or series of separate or discrete such executable units, portions of which may operate in sequence or parallel or may operate when called upon as subroutines or callable modules. It should be understood that the various implementations are not limited to the arrangements and instrumentality shown in the figures, but may instead be implemented in other arrangements while providing the same or similar described functionality. In practice the functional units or blocks described in the figures may be implemented on different processors, computers, servers, or virtual machines and / or may be spread among local and / or remote execution platforms (e.g., a workstation, a local server, a cloud-based data center, and so forth).

[0060] As used herein, the term “computing system” or “processor-based system” refers to an electronic computing device such as, but not limited to, a single computer, virtual machine, virtual container, host, server, laptop, and / or mobile device, or to a plurality of electronic computing devices working together to perform the function described as being performed on or by the computing system. As used herein, the term “medium” refers to one or more non- transitory, computer-readable physical media that together store the contents described as being stored thereon. Embodiments may include non-volatile secondary storage, read-only memory (ROM), and / or random-access memory (RAM).

[0061] As discussed herein, the computational analysis of genomic data (e.g., sequence data) may be affected by confounding variation that is not directly related to the coding region or sequence in question. For example, identification of non-coding variants (e.g., variants in promoter regions) that cause extreme levels of gene expression, either high or low, is relevant to the diagnosis of the pathogenicity of genetic diseases. However, there are numerous confounding factors that can interfere with the identification of pathogenic variants. Isolating variants by examining rare variants that can be associated with specific pathologies can simplify the problem. Further, removing noise introduced by confounders can increase the signal-to-noise ratio and improve sensitivity of the analysis, making it easier to identify such impactful variants. With this in mind, it may be useful to apply computational techniques, including techniques incorporating or based on artificial intelligence, to epigenetic studies so as to improve the sensitivity of such techniques in determining genetic associations between variable genetic loci and the expression levels of individual genes.

[0062] With the preceding in mind, the presently disclosed techniques are generally directed to facilitating the identification of functional rare variants, such as variants within promoter regions. Such identified rare variants may be utilized in various ways including, but not limited to, their use as training data for downstream models or analytics that may be used to predict the impact of a variant or variants on gene expression. Factors that need to be addressed in performing such analyses include potentially small sample sizes, other factors that may separately be impacting gene expression, and the ability to associate a given variant in a noncoding region, such as a promoter region, with a respective gene expression level or product. In certain embodiments, promoter regions to be assessed for functional variants may be identified or localized using techniques discussed herein, which utilize multiple data sources to accurately identify transcription start site (TSS) position(s), which are used to localize promoter regions for a gene. In further aspects as discussed herein, allele specific expression is evaluated for genes of interest, as opposed to total gene expression at the gene of interest. This may be done by selecting a suitable allele-expression single nucleotide polymorphism (aeSNP) for the gene. In certain embodiments discussed herein, the aeSNP closest to the TSS determined for the gene is selected as the best or most suitable SNP for assessing allele-specific expression. Computational pipelines incorporating such allele specific expression data may be used both to identify variant effects in a promoter region and / or to facilitate discovery of additional variants that may be pathogenic in nature.Technological Implementation Examples of the Technology Disclosed

[0063] Gene expression is the process by which the instructions in DNA are converted into a functional product, such as an RNA molecule or a protein. The technology disclosed facilitates identification of rare variants in non-coding regions that cause extreme levels of gene expression (e.g., under-expression or over-expression or the gene product). The rare variants are identified by association with nearby genes that have the extreme levels of gene expression. In one implementation, the technology disclosed identifies those individuals who have a particular variant in the promoter region of a gene and also have significantly different gene expression for that gene compared to individuals who do not have the particular variant. Based on identifying such individuals, the technology disclosed classifies the particular variant as a gene-expressionaltering variant. In certain embodiments, the technology disclosed further uses artificial intelligence or other computational techniques to train a one or more models using the identified rare variants as training data and their quantitative or qualitative phenotype of under-expression or over-expression as ground truth labels. In this manner, functional rare variants may be identified with high-confidence and, in some aspects, may be used as training data for downstream artificial-intelligence (Al) or other computational models that may in turn be used to predict the impact of variants on gene expression.

[0064] As discussed herein, accurate localization of the TSS (or TSSs in some instances) is useful in identifying a respective promoter region of a gene and, correspondingly, for identifying functional and rare variants within such promoter regions that may give rise to over- or underexpression. The presently described approaches utilize multiple sources of TSS localization data along with a TSS predictor model having improved accuracy relative to conventional approaches of identifying a TSS. In certain implementations, the TSS localization data may be integrated across multiple tissues (e.g., tens, hundreds, or thousands of tissues), such as in instances where a rare variant is present in a single individual or small number of individuals, so as to counter the lack of repeat examples across individuals. Among other uses, and as discussed herein, the TSS predictor model may be utilized in determining a heterozygous single nucleotide polymorphism (SNP) that may be used in assessing allele-specific expression at the gene, as opposed to total expression.

[0065] As used herein, the terms “chromatin sequence”, “chromatin input sequence”, and “output sequence”, may be used, such as in relation to expression measurement. Chromatin refers to a region of DNA (i.e., a DNA sequence) with bound proteins and / or RNA. The DNA sequence in a section of chromatin may be protected by protein and RNA and later sequenced as in DNA footprinting. The chromatin sequence may also be chemically modified. For example, DNA sequences often have methyl groups attached to the nucleotides in the sequence and thus are methylated. Enzymes that bind directly to chromatin components and catalyze chemical modifications of the chromatin components can alter chromatin structure, and changes in chromatin structure can also alter the ability of chromatin-interacting enzymes to access their target ligands and function. The structure of chromatin and the enzymes that alter structure directly influence the accessibility of a gene for transcription and expression. As discussedherein, DNA variants can cause changes in chromatin structure which subsequently can lead to epigenetic effects related to such processes as transcription factor binding and enzymatic reactions necessary for the proper regulation of gene expression and gene suppression.

[0066] Conversely, epigenetic effects on chromatin such as methylation and protein binding events can affect mutation rate, potentially introducing variants that may be silent or pathogenic. The study of evolutionary constraint on a gene and pathogenicity of variants related to that gene is significantly more comprehensive and accurate when augmented by the inclusion of epigenetic features, the study of which may be facilitated by use of the presently described techniques. Overall, the technology disclosed encompasses various embodiment and implementations which are amenable to a range of training and learning strategies to generate outputs which can be applied to the prediction of gene expression and gene pathogenicity for a target genetic sequence. The presently described tools facilitate or otherwise provide for the study of epigenetic signals affecting biological replication and transcription processes. With this in mind, the presently described tools facilitate or otherwise provide for the study of epigenetic signals affecting biological replication and transcription processes.Rare Variant IdentificationProcessing Pipeline

[0067] Turning to FIG. 1, one implementation of a processing pipeline 100 is depicted that identifies rare variants that cause extreme levels of gene expression. In this example, gene expression levels are accessed (step 110) for a group of individuals. In one implementation, the gene expression levels are accessed from Genotype-Tissue Expression (GTEx) data stores. The gene expression levels in the depicted example are normalized (step 120), for example, by calculating a mean and a plurality of standard deviations from the mean. Outlier individuals from the group of individuals that have extreme levels of gene expression are identified (step 130). The extreme levels of gene expression are, in certain embodiments, determined from tail quantiles 124 of the normalized gene expression levels 122. Examples of the tail quantiles 124 include one or more standard deviations from the mean, both in the positive and the negative directions. For example, the outlier individuals have gene expression levels with z-scores of atleast 11.2| from the mean. Rare variants from gene sequences of the outlier individuals are selected (step 140). The rare variants are selected based on an allele frequency cutoff. For example, a rare variant may be defined or otherwise characterized as having a minor allele frequency (MAF) of less than 1.0%, 0.5%, 0.1%, 0.05%, and so forth.

[0068] A causality model may be fitted (step 150) to characterize causal relationships between the rare variants and the extreme levels of gene expression in the outlier individuals while controlling for a plurality of confounders. In one implementation, the fitted causality model characterizes or models the causal relationship(s), allowing a prediction or estimate of a gene expression level of a particular gene in a particular chromosome to be derived based upon a variant-driven gene expression level caused by a particular rare variant. In one implementation, the fitted causality model measures a contribution of the variant-driven gene expression level as a variant effect covariate.

[0069] By way of example, such causality models may be implemented as a logistic regression model, a linear regression model, an analysis of covariance (ANCOVA) model, or a multivariate analysis of covariance (MANCOVA) model, and / or another suitable statistical model. Examples of the plurality of confounders include distal / / / / AS-expression quantitative trait loci (eQTLs) effects, local cA-eQTLs effects, genotype-based principal components (gPCs), expression residuals (PEER (Probabilistic Estimation of Expression Residuals)) effects, environmental effects, population structure and ancestry effects, and gender effects, batch effects, genotyping platform effects, and library construction protocol effects.

[0070] Causality scores for the rare variants are generated (step 160) based on the determined causal relationships. A particular causality score of a particular rare variant indicates a likelihood of the particular rare variant causing an extreme level of gene expression in those outlier individuals whose gene sequences contain the particular rare variant. In one implementation, the causality scores are probability values (p-values), such as may be determined using a Pearson correlation coefficient or other statistical approach.Processing System

[0071] Turning to FIG. 2, one implementation of a processing system 200 is illustrated that generates causality scores for outlier variants that cause extreme levels of gene expression. In the depicted example, datastore(s) 202 accessed by the processing system 200 stores gene expression data, which may be sourced from GTEx, RNA-Seq, whole genome sequencing (WGS), or other suitable sources. In some implementations, the gene expression levels are further stratified into tissue-specific gene expression levels for a plurality of tissues. In the embodiment of FIG. 2, a normalizer 204 (e.g., an executable normalization application or routine or suitably programmed or configured normalization firmware or dedicated circuitry) normalizes the gene expression data and stores the normalized gene expression data in datastore(s) 206. In one implementation, the gene expression levels for each gene in each tissue are normalized using quantile normalization. The normalized gene expression data may be characterized by a mean and one or more standard deviations from the mean. In one implementation, extreme levels of gene expression are determined from tail quantiles of the normalized gene expression data.

[0072] In the depicted example, a causality model 212 is fitted to determine causal relationships between variants 216 and differing levels of gene expression while controlling for a plurality of confounders 226. In practice, the causality model may be implemented as or by an executable application or routine or suitably programmed or configured firmware or dedicated circuitry. In one implementation, the fitted causality model 212 determines the causal relationships by predicting a particular gene expression level of a particular gene in a particular chromosome in dependence upon a variant-driven gene expression level caused by a particular variant, which is referred to herein as gene expression caused by variants 214. In some implementations, the causality model 212 is fitted separately for different tissues or cell types or for each tissue or cell type of interest. Further, in some instances the causality model is fitted using stratification. In one implementation, the fitted causality model 212 measures a contribution of the variant-driven gene expression level 214 (caused by the variants 216) as a variant effect covariate (See FIG. 3, gene expression caused by rare variants 308). In one implementation, the fitted causality model 212 measures a contribution of the confounders- driven gene expression level 224 (caused by the confounders 226) as a plurality of confoundereffect covariates (See FIG. 3, gene expression caused by distal trans-eQTLs effects 304, gene expression caused by local cis-eQTLs effects 306, and gene expression caused by other effects 310).

[0073] Examples of suitable causality model 212 include, but are not limited to, a logistic regression model, a linear regression model, an analysis of covariance (ANCOVA) model, and / or a multivariate analysis of covariance (MANCOVA) model. Examples of the confounders 226 include distal / / ' s-expression quantitative trait loci (eQTLs) effects, local c / .s-eQTLs effects, genotype-based principal components (gPCs), expression residuals (PEER) effects, environmental effects, population structure and ancestry effects, and gender effects, batch effects, genotyping platform effects, and library construction protocol effects.

[0074] Turning back to FIG. 2, in the depicted example, the causality model 212 generates, as an output, confounder-corrected-normalized gene expression 232 caused by the variants 216. In one such example, the effects of the confounders 226 are regressed out.

[0075] A rare variant identifier 234 (e.g., an executable identifier application or routine or suitably programmed or configured identifier firmware or dedicated circuitry) may be used to identify rare variants 236 from among the variants 216. The rare variants 216 can be selected based on an allele frequency cutoff. For example, the rare variants 236 can have a minor allele frequency (MAF) of less than 1.0%, 0.5%, 0.1%, 0.05%, and so forth. In some implementations, the rare variants may be or may include singleton variants, such as singleton variants that occur in only one outlier individual from the outlier individuals. In practice, the rare variants may be non-coding variants, such as, but not limited to, 5' untranslated region (UTR) variants, 3' UTR variants, enhancer variants, and promoter variants.

[0076] A causality score generator 242 (e.g., an executable score generation application or routine or suitably programmed or configured score generation firmware or dedicated circuitry) generates causality scores 246 for the rare variants 236 based on the confounder-corrected normalized gene expression 232. A particular causality score of a particular rare variant indicates a likelihood of the particular rare variant causing an extreme level of gene expression in those outlier individuals whose gene sequences contain the particular rare variant. In one implementation, the causality scores 246 are probability values (p-values). In oneimplementation, the p-values are determined by a Pearson correlation coefficient or other suitable statistical techniques. In the depicted example, a rare variant ranker 256 (e.g., an executable ranking application or routine or suitably programmed or configured ranking firmware or dedicated circuitry) ranks the rare variants 236 based at least on the causality scores 246 and stores the ranked rare variants in datastore 252. In practice other variables or factors may go into the ranking process implemented by the variant ranker 256, including but not limited to statistical confidence, degree of rarity (or lack thereof) of the respective variant, the distribution of expression associated with the respective gene, sequence similarity to other variants.Causality Model

[0077] With the preceding in mind, and turning to FIG. 3, a conceptual example of a fitted causality model 212 is illustrated that determines causal relationships between rare variants 236 and levels of gene expression while controlling for a plurality of confounders. The illustrated example of a fitted causality model 212 has a dependent variable 302 titled “G” (gene expression) and a plurality of independent variables corresponding to: gene expression caused by distal Zra / 7.s-eQTLs effects 304 (denoted Gl), gene expression caused by local cz.s-eQTLs effects 306 (denoted G2), gene expression caused by the rare variants 308 (denoted G3), and gene expression caused by other effects 310 (denoted G4). With respect to the gene expression caused by other effects 310, such effects may include, but are not limited to, genotype-based principal components (gPCs) 320, expression residuals (PEER) effects 330, environmental effects 340, population structure and ancestry effects 350, and gender effects, batch effects, genotyping platform effects, and library construction protocol effects 360.

[0078] The techniques discussed herein may facilitate the estimation or calculation of the gene expression caused by rare variants, such as by removing or otherwise correcting for the effects of gene expression caused by distal / ra / z.s-eQTLs effects 304, gene expression caused by local czv-eQTLs effects 306, and / or gene expression caused by other effects 310 from gene expression 302. That is, in accordance with the high-level overview provided by FIG. 3:(1) G3= G - Gl - G2 - G4so as to eliminate (or correct for) confounding effects and to accurately determine the gene expression caused by the rare variants 308. By way of example, in certain embodiments, the fitted causality model 212 may control for the distal Zraws-eQTLs effects by predicting the particular gene expression level (G) 302 in dependence upon a trans gene expression level (Gl) 304 caused by other genes in other chromosomes. In one implementation, the fitted causality model 212 measures a contribution of the trans gene expression level as a trans effect covariate.

[0079] In certain implementations, and as discussed in greater detail herein, the fitted causality model 212 may control for the local c / .s-eQTLs effects by predicting the particular gene expression level (G) 302 in dependence upon a cis gene expression level (G2) 306 caused by a presence of a plurality of common variants in a neighborhood of the particular gene. In some implementations, the neighborhood is defined by an offset from a transcription start site (TSS) in the particular gene. In one implementation, and as discussed herein, the fitted causality model 212 measures a contribution of the cis gene expression level as a cis effect covariate.

[0080] Further, in certain implementations the fitted causality model 212 may control for population structure and ancestry effects by predicting the particular gene expression level (G) 302 in dependence upon a gPC gene expression level (G4) 310 caused by the gPCs. In one implementation, the fitted causality model 212 measures a contribution of the gPC gene expression level (G4) 310 as a population structure and ancestry effect covariate.

[0081] Similarly, the fitted causality model 212 may control for the PEER effects by predicting the particular gene expression level (G) 302 in dependence upon a PEER gene expression level (G4) 310 caused by the PEER. In one implementation, the fitted causality model 212 measures a contribution of the PEER gene expression level (G4) 310 as a PEER effect covariate.

[0082] In addition, the fitted causality model 212 may control for environmental effects by predicting the particular gene expression level (G) 302 in dependence upon an environment gene expression level (G4) 310 caused by the environmental effects. In one implementation, the fitted causality model 212 measures a contribution of the environment gene expression level (G4) 310 as an environmental effect covariate.

[0083] As used herein, extreme levels of gene expression may include one or both of gene over-expression and gene under-expression.

[0084] In practice, and as discussed herein, the fitted causality model 212 may be used to determine the causal relationships between the rare variants and gene over- or under-expression while controlling for the plurality of confounders. In some implementations, the fitted causality model 212 generates “over causality scores” and / or “under causality scores” for the rare variants. A particular causality score of a respective rare variant indicates a likelihood of the rare variant causing under- or over-gene expression, respectively, in those outlier individuals whose gene sequences contain the respective rare variant.

[0085] In one implementation, the causality scores are probability values (p-values), such as may be determined by a Pearson correlation coefficient or other suitable statistical technique. In some implementations, the p-values specify statistically unconfounded likelihoods of the rare variants increasing or decreasing gene expression with respect to the gene in question. With respect to the manner in which such causality score are reported or otherwise provided to a reviewer, such causality scores may be reported along with a particular chromosome on which the respective rare variants is located, along with a location or position on such a chromosome for the rare variant, an identification of reference nucleotides (or bases) corresponding to the rare variants, and / or an identification of alternate nucleotides (or bases) (relative to a reference sequence) representing the rare variants.

[0086] With the preceding in mind, certain further aspects discussed herein are directed to identifying functional rare variants, including variants within a promoter region of a gene, with a high degree of confidence. As noted herein, such identified rare variants may be utilized in various ways including, but not limited to, their use as training data for downstream models or analytics that may be used to predict the impact of a variant or variants on gene expression. Particular, factors addressed by the present techniques include, but are not limited to, the ability to associate a given variant in a non-coding region (e g., a promotor, enhancer, or UTR) with a respective gene. By way of example, in certain embodiments promoter regions to be assessed for functional variants may be identified or localized using the presently disclosed techniques, which utilize multiple data sources to accurately identify transcription start site (TSS)position(s), which are used to localize promoter regions for a gene. In further aspects as discussed herein, allele specific expression is evaluated for genes of interest, as opposed to total gene expression at the gene of interest. This may be done by selecting a suitable alleleexpression single nucleotide polymorphism (aeSNP) for the gene. In certain embodiments discussed herein, the aeSNP closest to the TSS determined for the gene is selected as the best or most suitable SNP for assessing allele-specific expression.Transcription Start Site (TSS) Identification and / or Prediction

[0087] With the preceding computational pipeline and causality modeling in mind, one aspect of the presently described techniques relates to the identification, with high-confidence, of functional rare variants that can be employed as positive training data for downstream models, such as promoter identification and / or evaluation models. In particular, functional rare variants identified within promoter regions may be employed to model and / or predict the impact of such variants on gene expression.

[0088] In certain approaches, it may be desirable to identify a promoter region associated with a gene. In such approaches, identification of the promoter region may initially involve identifying a transcription start site (TSS) or sites associated with the gene of interest. A “window” (e.g., ± 100 bp, ± 500 bp, or ± 1 kb around an identified or selected center point or consensus location) extending upstream and downstream from the TSS may also be defined that is likely to have a regulatory role for the gene of interest. Conventionally, the TSS may be identified for a respective gene based on an annotation file or similar resource that may be accessed, such as from a datastore, and that, among other annotations potentially, defines one or more TSSs for a given gene based on prior work or studies.

[0089] In practice, however, TSS information derived from such annotation files may be inaccurate and / or may vary (such as due to the use of different sources or techniques to identify TSS locations). This may be problematic for approaches, such as promoter region localization, where the location of the TSS has significance and, correspondingly, may impact efforts to localize a promoter region and / or to identify promoter region variants, which may be a source of under- or over-expression of the respective gene.

[0090] With the preceding in mind, in certain embodiments discussed herein a TSS predictor model or tool is described that employs or otherwise integrates TSS localization data from multiple sources to predict or provide accurate TSS positions for a given gene. By way of example, in one such implementation data from one or more annotation files may be combined with experimental assay results across multiple different tissue and / or cell types as well as with sequence conservation data so as to identify the most frequently used TSS per gene and to validate that selection via enrichment of rare variants.

[0091] In a practical implementation, the output of such a process may be embodied in a customized or modified annotation file that may subsequently be used to provide accurate TSS location data for a given gene. By way of example, such a customized or modified annotation file may integrate available annotation file data, conservation data, and experimental assay data (such as may be related to functional measurements) to identify consensus TSS locations or regions. Such a custom annotation file may be used to improve functional variant identification within the promoter region and to improve model predictions related to variants in such promoter regions.

[0092] Prior to further discussion of the TSS predictor model, a brief discussion of the outlier techniques and methodology employed herein is provided so as to introduce certain terminology. As used herein, an outlier test may be employed for the purpose of determining the presence of over- or under-expression for a given gene in an individual. In the present context, to identify functional variants, all singleton single nucleotide variants (SNVs) present in promoter regions across a representative set of individuals within a dataset (e.g., GTEx) are tested for significant effects on gene expression. For each promoter variant, a vector is created that contains all expression values for the proximal gene (i.e., the gene associated with that promoter region) in that individual across tissues. A T-test is performed to assess if individuals with the variant exhibit gene expression that is significantly different from the population. This may be compared to a threshold p-value to identify outlier variants for training of the model. To assess the signal-to-noise ratio (SNR) in the set of identified significant variants, a shuffled background may be created by swapping genomes between individuals, thereby providing an estimate of how many significant hits would be expected just due to chance. Using the same methodology, thenumber of significant variants using the shuffled background may be determined, providing an assessment of signal-to-noise.

[0093] In particular, and in accordance with aspects discussed herein, in assessing the gene expression effects of rare variants in a promoter region, a determination may initially be made as to a TSS and promotor region (e g., a window selected based on the TSS) to be evaluated using a suitable dataset, such as GTEx. Based on this promoter window or windows, individuals are identified in the dataset that carry a gene variant in the promoter region defined by the window or windows, which may be a single individual in rare variant contexts. In this manner, a candidate list of variants that may have a functional impact on proximal gene expression may be developed and these candidate variants may be tested or otherwise assessed for such gene expression impacts.

[0094] With the preceding in mind, and in support of the preceding discussion, FIG. 4 provides a graphical illustration of the impact of TSS position on functional variant discovery. In particular, FIG. 4 graphically illustrates the impact on variant discovery of the annotation file employed to localize TSS positions. To generate the illustrated data, existing annotation files were used to identify the TSS for each gene investigated and to then determine which variants in the vicinity of the respective TSS to test based on a window around the reported TSS position. As shown, the y-axis illustrates a measure of how accurate the functional variant test set is (i.e., how many identified variants within the window actually do impact gene expression, with higher number denoting higher accuracy). In particular, the y-axis denotes the outlier count compared to what is seen in a shuffled context and is a measure of how clean (i.e., noise-free) the data set is, thereby corresponding to a signal-to-noise ratio or measure. As shown, the annotation file employed does impact the ability to detect functional variants in the promoter region. In particular, the more accurate the TSS location estimate is, the more accurate the localization for the promoter region also is.

[0095] In certain implementations discussed herein, data acquired using Cap Analysis of Gene Expression (CAGE or CAGE-seq) may be used to generate data that may be used to generate or update a TSS predictor as discussed herein. In particular, CAGE may be employed to generate data related to transcription initiation frequency of transcription start sites at singlebase pair resolution. In practice, CAGE techniques trap (via oxidation and biotinylation) the RNA 5 '-cap structure enabling the isolation of cDNAs from the 5' ends of mRNAs. This further enables the enrichment of and selective sequencing of complete transcripts, thereby allowing mapping of initiation sites of both capped coding and noncoding RNAs. This technique allows the precise position (e.g., location or position on a chromosome or nucleic acid sequence) of every read to be mapped, allowing the start locations of every read to be determined across tissues and, correspondingly, allowing TSS positions to be more accurately identified.

[0096] With the preceding in mind, CAGE data can be acquired for different transcripts for the same gene and can be used to map, down to single-nucleotide resolution, the specific start location for each respective transcript. In this manner, for each gene being analyzed, a determination can be made as to where the most reads originate in a given tissue, i.e., which TSS is used most frequently in that tissue.

[0097] With this in mind, in certain implementations CAGE reads may be determined for a given gene which provide all of the different transcripts possible for a respective gene. The different read counts for each transcript may be summed across tissues to select the transcript having the greatest read count support as the TSS across the analyzed tissues. In this manner, data may be integrated to inform TSS selection based on read counts summarized across many different tissues.

[0098] While TSSs determined in this manner often map to known annotations for the respective TSS, for a subset of genes this is not the case. For this subset of genes the CAGE- based analysis yields a different presumptive TSS than that identified in existing annotations. This is illustrated graphically in FIG. 5, which illustrates that the CAGE-based approach yields a different TSS than existing annotations files for a subset of genes having a difference greater than 10° between the canonical TSS and CAGE-determined TSS. With respect to this subset of genes, it is likely that the genes in question have more than one TSS position and that the CAGE- determined TSS differs from the canonical annotated TSS location. As may be appreciated, such results further support the utility of a TSS predictor model configured to utilize multiple sources of information to facilitate accurate localization of a respective TSS.

[0099] With this in mind, and turning to FIG. 6, signal-to-noise (y-axis) for prior canonical annotation TSS locations is illustrated by the leftmost and rightmost bars, while an assessment of TSS location refined using CAGE data is illustrated by the middle bar. As in the prior illustrations, the greater the height of the bar, the better the signal-to-noise of the measure. As illustrated, supplementing the canonical annotation data with CAGE-based TSS location data yields improved signal-to-noise in instances where the CAGE-based results yield a different location than the canonical annotation TSS data, which in turn equates to improved signal in selecting variants of interest (e.g., rare variants of interest) in the promoter region that have effects on gene expression.

[0100] While the preceding illustrates the efficacy of CAGE-based techniques in identifying a TSS, in practice useful information is still present in other indicators. Correspondingly, and as discussed herein, a suitable TSS predictor tool may incorporate CAGE data with non-CAGE data to further improve localization of a TSS, which may in turn facilitate identification of rare variants in a promoter region which may be of interest for model training or for identification of variants in the promoter region that alter or otherwise affect gene expression. By way of example, such a TSS predictor tool may incorporate existing or canonical TSS annotations, CAGE-based TSS location data, histone modification data, transcription factor (TF) data, chromatin immunoprecipitation sequencing (ChlP-seq) data, assay for transposase-accessible chromatin with sequencing (ATAC-Seq) data, methylation data, sequence conservation data, and so forth. In accordance with such an approach, multi-track data in addition to CAGE-data may be employed to improve TSS localization relative to the use of CAGE-data or canonical annotation data alone. By way of example, in one embodiment approximately 6,000 tracks for every 5 nucleotides (NT) may be used to predict the probability of a genetic loci being a TSS.

[0101] By way of example, and turning to FIG. 7, a representation of an analysis or output of a TSS predictor tool is visually represented incorporating CAGE-based TSS estimation, canonical annotation data, histone modification data, TF data, ChlP-seq data, ATAC-Seq data, and conservation data. In this example, the TSS predictor tool incorporates 6,065 tracks in total and may be used to identify a consensus TSS region based on the multiple data sources. In one implementation, the estimated TSS position associated with each track or technique may be annotated on the representation, so as to allow a reviewer to review the different TSSannotations, though a consensus TSS estimate may still be provided as well. Signals are averaged for every ± 5 nt, 10 nt, and 25 nt. In the depicted example, the predictor line 400 illustrates the consensus TSS location. For further reference, track 404 corresponds to the CAGE-derived TSS location, track 408 corresponds to the TSS predictor tool consensus (illustrating the probability of the TSS prediction tool’s consensus), and track 412 comprises the CAGE data (e.g., read count data), which illustrates a strong correlation with the TSS predictor of track 408.

[0102] In certain embodiments a hierarchy or rule set may be used in evaluating the different data sets and applying them to the estimation of the consensus TSS position. By way of example, TSS data may have priority to ATAC-Seq or ChlP-Seq data, which in turn may be prioritized over sequence conservation data. In practice this may be implemented as a set of rules prioritizing or ranking the different data sources. In addition or in the alternative, such a scheme may be implemented as a weighting scheme, with different data sets being given different numeric weighting in estimating the TSS location. Such a weighting scheme may be fixed or pre-defined such that weights are not changed or may be dynamic in nature (e.g., determining or assigning weights based on a measure of data quality and / or statistical confidence). Similarly, benchmarks or thresholds may be employed based on data quality either for the inclusion or exclusion of certain data sources or to establish weights in a dynamic weighting scheme.

[0103] As may be appreciated, in addition to indicating possible additional locations of a TSS a given gene, the TSS predictor model may also be used to define a window or region within which the TSS is present. By way of example, the region of the TSS predictor track 408 exhibiting peaks may be used to define such a window 416, with the respective TSS likely present within the span defined by or encompassing the peaks. Such a window 416 around a tight or limited range around the TSS may be useful in identifying a region in which variants may likely impact gene expression. By way of example, such a window may be defined as ± 100 bp, ± 500 bp, or ± 1 kb around an identified or selected center point or consensus location. In practice, such a window may be a fixed window across genes and / or tissues. Alternatively, in other embodiments the window may be defined on a gene-by-gene basis, such as on a gene-by- gene basis across multiple tissues, or on a tissue-by-tissue basis.

[0104] Performance of the TSS predictor model is illustrated in FIG. 8 in the form of precision-recall (PR) curves for both training data and test data derived using the TSS predictor model. Further, and turning to FIG. 9, bar graphs are illustrated demonstrating the proximity of the predicted TSS based on the TSS predictor model (as represented by track 408 of FIG. 7) relative to the CAGE derived TSS (track 404 of FIG. 7). As illustrated, for a probability cutoff of 0.07 (P > 0.7), approximately 71% of the predicted TSS are within 10 nt distance to the CAGE TSS location (distance to the CAGE-derived TSS or d CAGE TSS) and about 91% are within 100 nt. This illustrate that the TSS predictor model is close to the results obtained via CAGE-derived TSS techniques but, further, that there are still additional TSS locations that can be found by the TSS predictor model in addition to the ones found using CAGE data alone. In practice this may be due to the CAGE-derived techniques selecting a single TSS per gene, when in some instances there may be multiple active TSSs, such as in multiple or differing tissues. In such circumstances, the TSS predictor model is able to identify additional TSSs, as opposed to identifying only one per gene.

[0105] In a further comparison between the results of the TSS predictor model and CAGE- derived singlets annotations, FIG. 10 illustrates that, when matching window sizes for promoter region analysis, the TSS predictor model as described herein outperforms CAGE-based TSS localization. In particular, and turning to FIG. 10, a plot of enrichment (y-axis) versus distance around peak (x-axis) is provided comparing the TSS predictor model TSS location results against CAGE-derived TSS positions. In this context, enrichment provides useful signal-to-noise context and is based on how many more tests were statistically significant relative to background noise or measurements. At 0 distance around the peak (i.e., at the peak) the variant in question falls exactly at the CAGE-annotated TSS or falls in the peak window 416 that had highly predictive values in the TSS predictor model. Moving rightward along the x-axis, greater distances correspond to greater distances away from that peak. As illustrated, enrichment decreases as distance from peak increases for both the TSS predictor and the CAGE-derived TSS locations. However, a consistent increase in the enrichment signal for the TSS predictor model locations relative to the CAGE-derived TSS locations is observed. This result is indicative that significant new variants of interest (e.g., rare variants in the promoter region of a gene of interest) would be missed using the CAGE-derived TSS alone but identified using the TSSpredictor model. This is further illustrated in FIG. 11, which is a bar plot illustrating the number of variants that would be discovered using the CAGE-derived TSS locations and the presently described TSS predictor model locations respectively. As will be appreciated, to the extent that it is desirable to improve the number of significant variants discovered over background, the TSS predictor model outperforms CAGE-based TSS location techniques.

[0106] Turning to FIGS. 12A, 12B, 12C, 13A, 13B, and 13C, respective analyses of outliers from CAGE-derived TSS (FIGS. 12A, 12B and 12C) and from the present TSS predictor model (FIGS. 13A, 13B, and 13C) are illustrated. In particular, for each technique, model predictions based on GTEx data were performed and outlier results plotted. In particular, a plot is provided of sensitivity versus specificity (for under- versus over-expression), of precision versus recall (for under- versus no expression), and of precision versus recall (for over- versus no expression). Per this analysis the TSS predictor model was able to identify more outliers than the CAGE- based model (3769 versus 2925) without a loss in performance.

[0107] As discussed herein, in practice the TSS predictor model may be used to output a new or custom TSS annotation file (or to modify an existing file) that may be used to more accurately localize TSS locations or regions associated with the genes having annotations within the file. Such an annotation may include, for example, a chromosome and a TSS location or window, though other information such as a proximal gene and locations of other regulatory features may also be present. Unlike prior efforts, such a custom file integrates existing annotation data, CAGE-derived data, other functional measurement data or experimental assays, and / or sequence conservation data to more accurately locate one or more TSS locations associated with the genes referenced in the file. Technical effects associated with the creation and use of such a customized TSS annotation file include, but are not limited to, improved identification of functional variants within the promoter region as well as generation of improved datasets of rare promoter variants that may be used in training or using machine learning tools or other models directed toward predicting gene expression effects of variants within the promoter region of a gene.Allele Specific Expression

[0108] As noted herein, one use to which improved TSS localization may be applied is identifying a promoter or promoter region associated with a gene and, correspondingly, identifying variants in the promoter region that might impact gene expression, leading to either under- or over-expression of the gene product. In certain embodiments discussed herein, and in contrast to prior approaches, allele specific gene expression (as opposed to total gene expression) for a given gene may be measured. As may be appreciated, allele specific gene expression, in this context, refers to measuring or observing either the absolute or quantitative gene expression of two different alleles at a given gene locus or at the relative expression (e.g., the ratio) of the differing alleles at the locus. This may be contrasted with approaches in which total gene expression is measured as such measures collapse the gene expression of both alleles into a single measurement for a given gene. By measuring allele specific expression, useful data may be obtained for individuals that are heterozygous at the gene locus of interest as the allele specific expression data allows comparison of the expression of one allele (e.g., variant) versus the other allele (e.g., a reference allele) within the same individual. In particular, analysis of allele specific expression in this context may provide additional validation of expression effects and may facilitate the discovery of such expression effects that may not be evident when analyzing only the total gene expression data, and / or may provide improved resolution relative to total expression analysis (e.g., 1.5 x greater resolution).

[0109] By way of example, in certain approaches an allelic expression single nucleotide polymorphism (aeSNP) may be selected for a gene for which an individual is heterozygous. A comparison may then be made across reads for each variant to assess which reads are associated with each allele. In practice, it is useful to select an SNP that is contained in a high number of reads, such as above a threshold amount corresponding to pre-defined or user-specified degree of coverage. Reference to the aeSNP selected in this manner effectively associates or “tags” which reads are derived from or associated with each allele of the heterozygous gene. In certain implementations, one or more variant-level biochemical or phenotype effects may be identified based on the expression of the different alleles.

[0110] As discussed herein, selection of a suitable aeSNP from the available aeSNPs for a given gene impacts analysis of and determination of variant effects within the promoter region on gene expression. In particular, as described below, selecting an aeSNP based on distance tothe TSS location provides more accurate data than taking the top-covered aeSNP or combined read counts across all possible phased aeSNPs. With the preceding in mind, selection of an aeSNP based on distance from the TSS for a gene may benefit from improved localization of the TSS, as discussed herein. In addition, a czs-eQTL correction may be employed on allele specific expression measurements, as discussed herein, that further improves rare variant discovery. Such a correction controls for common variants that impact allelic imbalance to further highlight or emphasize those rare variant effects in the promoter region that may be of greater interest.

[0111] With the preceding in mind and with respect to the selection of an aeSNP for use in assessing allele specific expression, enrichments of rare promoter variants with significant allelic imbalance can be improved by using a reference ratio derived for the closest aeSNP to the TSS for a gene. By way of example, and turning to FIG. 14 a bar graph is illustrated demonstrating enrichment versus p-value where the p-values were derived from T-tests between normalized adjusted allele reference ratios for variant carriers relative to all other individuals. The data illustrated considers singletons within 100 base pairs (bp) of the TSS. In particular, the y-axis shows enrichment (which may be considered as corresponding to signal-to-noise in the present context). As noted herein, the T-tests in question relate to allelic ratios, as opposed to total gene expression.

[0112] Turning to FIG. 15, a further bar graph is illustrated in which the outlier count is illustrated as the y-axis (again corresponding to a signal-to-noise metric) while the three different techniques for deriving the reference ratio (i.e., how the SNP that represents the different alleles is chosen) are provided along the x-axis. On the left, a bar corresponding to phased expression (which integrates over many different heterozygous SNPs contained in the gene to map each read) is illustrated. In the middle, a bar corresponding to the top aeSNP is illustrated, where the top aeSNP is the SNP having the most reads that contain that SNP over any other potential heterozygous SNP for that individual (i.e., the top SNP by read). The top aeSNP may, therefore, be characterized as the SNP having the highest coverage in this context. On the right, a bar corresponding to the aeSNP being selected that is closest to the TSS is illustrated. As shown, the aeSNP closest to the TSS has superior signal relative to both the phased expression and to aeSNP. With this in mind, and as discussed herein, performance is improved by defining theallelic reference ratio based on the ratio of reads containing the closest heterozygous SNP (i.e., aeSNP) to the TSS for a gene in a given individual.

[0113] Further support is illustrated in FIG. 16, where a further bar graph illustrates that genes having significant allele specific expression (e.g., p < 0.01) have a substantial overlap with total expression effects (e.g., p < 0.01, as shown on the y-axis). As illustrated in the depicted graph, the x-axis is configured so as to differentiate allele specific expression by outlier thresholds, with a bar shown for closest aeSNP and top aeSNP for each outlier threshold. In this manner, expression effects can be compared for top aeSNP relative to the aeSNP closest to the TSS. In this example, similar results are seen when looking at both total and allele specific expression effects

[0114] With respect to enrichment, in accordance with the presently described studies it was determined that the enrichment of rare promoter variants within promoter regions was strongest when the closest aeSNP to the TSS and the top aeSNP are within 1 kb on average. By way of example, and with reference to FIGS. 17 and 18, various factors that influence how robust the allelic ratios are were evaluated. Turning to FIG. 17, data is binned based on the ratio of the counts of the closest aeSNP to the top aeSNP and then presented in the context of enrichment versus outlier threshold. Bars at each break point correspond to low-, mid-, and high-count bins (i.e., differing levels of coverage). With respect to the information conveyed, the coverage of the closest aeSNP to the TSS corresponds to how many reads were actually mapped that contained that aeSNP, which may be understood to convey how robust the reference estimate is. With respect to FIG. 18, data is binned based on the distance between the closest aeSNP to the TSS and the top aeSNP and then presented in the context of enrichment versus outlier threshold. Of particular note, the bin corresponding to the shortest distance between the closest aeSNP and the top aeSNP estimate (distance < 1 kb) consistently outperforms. Similarly, and turning to FIG.19, it is illustrated that observing higher coverage at the closest aeSNP to the TSS increases enrichments of rare promoter variants with significant allelic expression imbalance. That is, the greater the coverage (i.e., reads), the greater the observed enrichment. Further, it was observed that there was minimal impact from imposing a higher read count threshold (e.g., a minimum of thirty reads) when selecting the closest aeSNP to the TSS.

[0115] In further embodiments, combining signals associated with total gene expression and allele specific expression, as described herein, may allow for the identification of additional significant promoter region variants. By way of example, turning to FIGS. 20 and 21, combined total gene expression (TE) and allele specific expression (ASE) were analyzed along with the intersection of the two categories. With respect to the data employed, all variants having an allele frequency < 0.5% within 500 bp were considered where expression had a p < 0.01. As may be observed in FIG. 20, additional variants may be identified for analysis by including both total gene expression and allele specific gene expression signal as there is non-intersection signal attributable to each of these categories that would be lost absent consideration of the respective category. Thus, additional variants may be identified using analyses based on allelic expression in addition to total gene expression.

[0116] In further embodiments, the reference ratio may be improved for use in variant identification and selection by correcting for the signal attributable to common variants (i.e., non-rare variants) since this portion of the signal is not of interest. By way of example, cis- eQTLs may also contribute to differences in allelic expression, but their effects may be of little or no interest. In certain aspects the signal attributable to czs-eQTLs variants may be corrected for by calculating a reference ratio that is a proportion of reads coming from aeSNP relative to total expression. An expected reference ratio is calculated based on the c / .s-eQTL or c / .s-eQTLs. The expected reference ratio can then be subtracted from the calculated reference ratio, which may then be used in performing outlier tests. In such an approach, the effects attributable to a common variant, which are not of interest in the present context, are corrected and / or controlled for. The remaining effects are, thereby, more likely due to rare or non-common variants, such as singleton variants carried by the individual in question.

[0117] By way of example, and turning to FIG. 22, a bar graph is provided illustrating that improvement in the outlier count was obtained after correcting for the effects of the top aeQTL with consistent directionality. In effect, such a correction, by removing the effects of the largest contributing common variant, effectively increases the signal contributed by the other variants, including the rare or non-common variants of interest, essentially increasing the signal -to-noise ratio. Turning to FIG. 23, it may be observed that this effect can be even further improved bycorrecting for the effects of multiple aeQTLs, as opposed to correcting for only the top contributing aeQTL.

[0118] With the preceding discussion in mind regarding the use of allelic imbalance variants, and the correction for common variants as appropriate, it was determined that inclusion of significant allelic imbalance variants, alone or in conjunction with analysis of total gene expression, identified over 100 (e.g., between 100 and 1,500) additional promoter region variants for testing or training purposes. This is illustrated in FIG. 24, in which data is binned based upon expression (i.e., allelic imbalance, total gene expression, or the intersection of the two) and in the context of counts of the number of unique genes with any significant variants versus p-values. As will be noted, hundreds to over a thousand of variants were identified based on the allelic imbalance analyses described herein.

[0011] With the preceding in mind, it may be appreciated that various technological improvements and solutions are described herein. In one instance, the technological improvement of more accurate localization of a transcription start site (TSS) is described in which a computationally implemented TSS predictor model is generated and executed and which incorporates data from multiple sources. The data incorporated into the TSS predictor model includes some or all of canonical or existing TSS annotation data, CAGE-based TSS location data, histone modification data, transcription factor (TF) data, chromatin immunoprecipitation sequencing (ChlP-seq) data, assay for transposase-accessible chromatin with sequencing (ATAC-Seq) data, conservation data, and so forth. A TSS predictor model implemented in this manner allows identification of rare variants within a promoter region that would be missed using conventional methodology and / or allows more accurate localization of a respective TSS.

[0120] In other aspects, such an improved TSS localization or placement may be used in selecting an aeSNP for use in analyzing allelic expression imbalances associated with heterozygous variants within a promoter region. In particular, an aeSNP is selected that is closest to the localized TSS, as opposed to selection based on top coverage. Such an approach represents a technical improvement over prior techniques based on total-gene expression data alone and allows identification of rare variants within the promoter region that impact gene expression and which would be missed based solely on analysis of total gene expression. Afurther technical improvement includes incorporation of a c / .s-eQTL correction step by which the signal or effects of one or more common variants within the region are predicted and removed so as to improve the signal relative to background noise and to facilitate identification and / or selection of rare variants within the promoter region.

[0121] With the preceding in mind, and by way of technical summary, the techniques discussed and described herein facilitate localization of transcription start sites (TSS) and, thereby, the identification of variants proximate to the TSS, including but not limited to promoter regions. By way of example, an output or tangible product of the present TSS localization techniques may include a list (such as may be provided as a modified or curated annotation file) of consensus TSS locations that maybe used in the identification of functional rare variants within a promoter region and which may have associated or discernible effects on associated gene expression. Such a modified annotation file, as described herein solves a technical problem with respect to prior annotation files in that such prior annotation files were typically designed to be complete and list all possible transcripts for a given gene, inadvertently leading to noise in the present context, while the present TSS predictor generated annotation files have improved signal-to-noise characteristics in this regard.

[0122] Rare functional variants identified in this manner may be used for training and / or testing processor-implemented models (e.g., deep learning or machine learning models) designed to localize TSSs, localize promoter regions or windows, identify or predict variants within promoter regions, and / or identify or predict the impact on gene expression of such variants. In terms of technical improvements, the list of functional rare variants may be quantitatively superior to those based on prior techniques, may be qualitatively superior to those based on prior techniques, or both. Further, the TSS predictor model as described herein, by identifying where transcription starts for a gene, may be useful for determining or evaluating what non-coding regions are likely to impact expression of the respective gene.

[0123] In certain implementations rare variants identified within an individual may be used in disease risk prediction, burden testing, and so forth. For example, the techniques described herein may be used to facilitate identifying which variant within a promoter region (in a large promoter region, many variants may be present) is most likely contributing to over- or under-expression of an associated gene. Further, in certain embodiments rare variants identified within an individual may be identified and localized for treatment using advanced, targeted gene therapeutic techniques, such as CRISPR or similar gene modification technologies. In this manner, rare promoter variants associated with gene expression levels of an associated gene that are considered outside the bounds of clinically normal or otherwise in need of treatment may be addressed. Use of such treatment techniques may be facilitated or improved by use of the TSS locations based on the techniques disclosed herein.Computer System

[0124] With the preceding in mind, FIG. 25 is an example computer system 500 that can be used to implement various aspects of the technology disclosed. Computer system 500 includes at least one central processing unit (CPU) 504 that communicates with a number of peripheral devices via bus subsystem 508. These peripheral devices can include a storage subsystem 512 including, for example, memory devices and a file storage subsystem 516, user interface input devices 520, user interface output devices 524, and a network interface subsystem 528. The input and output devices allow user interaction with computer system 500. Network interface subsystem 528 provides an interface to outside networks, including an interface to corresponding interface devices in other computer systems.

[0125] In one implementation, a TSS predictor model 530 is communicably linked to the storage subsystem 512 and the user interface input devices 520. In other aspects, an allele specific expression routine or module and / or a c / .s-eQTL correction routine or module may also be present or may be present in place of the TSS predictor model 530.

[0126] User interface input devices 520 can include a keyboard; pointing devices such as a mouse, trackball, touchpad, or graphics tablet; a scanner; a touch screen incorporated into the display; audio input devices such as voice recognition systems and microphones; and other types of input devices. In general, use of the term “input device” is intended to include all possible types of devices and ways to input information into computer system 500.

[0127] User interface output devices 524 can include a display subsystem, a printer, a fax machine, or non-visual displays such as audio output devices. The display subsystem can includean LED display, a cathode ray tube (CRT), a flat-panel device such as a liquid crystal display (LCD), a projection device, or some other mechanism for creating a visible image. The display subsystem can also provide a non-visual display such as audio output devices. In general, use of the term “output device” is intended to include all possible types of devices and ways to output information from computer system 500 to the user or to another machine or computer system.

[0128] Storage subsystem 512 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 deep learning processors 532 or CPU 504.

[0129] Deep learning processors 532 can be graphics processing units (GPUs), field- programmable gate arrays (FPGAs), application-specific integrated circuits (ASICs), and / or coarse-grained reconfigurable architectures (CGRAs). Processors 532 can be hosted by a deep learning cloud platform such as Google Cloud Platform™, Xilinx™, and Cirrascale™. Examples of processors 2078 include Google’s Tensor Processing Unit (TPU)™, rackmount solutions like GX4 Rackmount Series™, GX20 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 VlOOs™, and others.

[0130] Memory subsystem 536 used in the storage subsystem 512 can include a number of memories including a main random access memory (RAM) 540 for storage of instructions and data during program execution and a read only memory (ROM) 544 in which fixed instructions are stored. A file storage subsystem 516 can provide persistent storage for program and data files, and can include a hard disk drive, a floppy disk drive along with associated removable media, a CD-ROM drive, an optical drive, or removable media cartridges. The modules implementing the functionality of certain implementations described herein can be stored by file storage subsystem 516 in the storage subsystem 512, or in other machines accessible by the processor.

[0131] Bus subsystem 508 provides a mechanism for letting the various components and subsystems of computer system 500 communicate with each other as intended. Although bussubsystem 508 is shown schematically as a single bus, alternative implementations of the bus subsystem can use multiple busses.

[0132] Computer system 500 itself can be of varying types including a personal computer, a portable computer, a workstation, a computer terminal, a network computer, a television, a mainframe, a server farm, a widely-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 500 depicted in FIG. 24 is intended only as a specific example for purposes of illustrating the preferred implementations of the present invention. Many other configurations of computer system 500 are possible having more or less components than the computer system depicted in FIG. 24.

[0133] The technology disclosed can be practiced as a system, method, or article of manufacture. One or more features of an implementation can be combined with the base implementation. Implementations that are not mutually exclusive are taught to be combinable. One or more features of an implementation can be combined with other implementations. This disclosure periodically reminds the user of these options. Omission from some implementations of recitations that repeat these options should not be taken as limiting the combinations taught in the preceding sections.

[0134] One or more implementations of the technology disclosed 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 actions indicated. Furthermore, one or more implementations and clauses of the technology disclosed or elements thereof can be implemented in the form of an apparatus including a memory and at least one processor that is coupled to the memory and operative to perform routine encoded within the memory as executable code. Yet further, in another aspect, one or more implementations and clauses of the technology disclosed or elements thereof can be implemented in the form of means for carrying out one or more of the method steps described herein; the means can include (i) hardware module(s), (ii) software module(s) executing on one or more hardware processors, or (iii) a combination of hardware and software modules; any of (i)-(iii) implement the specifictechniques set forth herein, and the software modules are stored in a computer readable storage medium (or multiple such media).

[0135] Other implementations can include a non-transitory computer readable storage medium storing instructions executable by a processor to perform any of the actions or implements any of the models or routines described herein. Yet another implementation of the implementations described herein can include a system including memory and one or more processors operable to execute instructions, stored in the memory.

[0136] This written description uses examples to disclose the invention, including the best mode, and also to enable any person skilled in the art to practice the invention, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insubstantial differences from the literal languages of the claims.

Claims

What is claimed is:

1. A method for identifying a transcription start site (TSS), comprising: accessing a genetic sequence annotation file containing prior TSS location annotations; identifying a consensus TSS for a gene based on estimated TSS positions for the gene across a plurality of different tissue types or different cell types, wherein identifying the consensus TSS for the gene utilizes Cap Analysis of Gene Expression (CAGE) data, Cap Analysis of Gene Expression (CAGE) data with functional measurements to identify the consensus TSS; validating the consensus TSS for the gene based on enrichment of promoter variants; and updating the genetic sequence annotation file with the consensus TSS to generate a customized TSS annotation file.

2. The method of claim 1, wherein the identifying the consensus TSS for the gene utilizing CAGE data comprises using read counts summarized across the plurality of different tissue types or different cell types with CAGE sequence data.

3. The method of claim 2, wherein the consensus TSS is identified by selecting the transcript with the most read count support or the highest CAGE count ratio as the TSS across the plurality of different tissue types or different cell types.

4. The method of claim 1, wherein variants within a window around the consensus TSS are associated with changes in an expression level of the gene.

5. The method of claim 4, wherein the window comprises a region ± 100 bp, ± 500 bp, or ± 1 kb about the TSS.

6. The method of claim 4, wherein the window is a fixed window.

7. The method of claim 4, wherein the window is defined gene-by-gene or tissue-by-tissue.

8. The method of claim 1, wherein identifying the consensus TSS comprises applying a hierarchy or weighting to different data sources used in identifying the consensus TSS.

9. The method of claim 1, wherein identifying the TSS for the gene further comprises analysis of sequence conservation data.

10. The method of claim 1, wherein the customized TSS annotation file has increased signal- to-noise relative to the genetic sequence annotation files.

11. The method of claim 1, further comprising: identifying one or more functional promoter variants using the customized TSS annotation file.

12. The method of claim 1, wherein the functional measurements comprise one or more of histone modification data, transcription factor (TF) data, chromatin immunoprecipitation sequencing (ChlP-seq) data, assay for transposase-accessible chromatin with sequencing (ATAC-Seq) data, or methylation data.

13. The method of claim 1, further comprising: mapping promoter variants to affected genes using the customized TSS annotation file.

14. The method of claim 1, comprising identifying one or more functional variants using the customized TSS annotation file; and generating a treatment or therapy, a disease risk prediction, or a burden estimate for the one or more functional variants using the TSS.

15. The method of claim 14, further comprising: prescribing or administering the treatment or therapy to a subject having the one or more functional variants, wherein the treatment or therapy comprises a targeted gene therapy treatment.

16. A method for using allele-specific expression (ASE) for variant assessment, comprising: selecting an allele expression single nucleotide polymorphism (aeSNP) present at a gene of an individual, wherein the individual is heterozygous at the gene, and wherein the aeSNP is selected based on location relative to a transcription start site (TSS); comparing the gene expression associated with the different alleles of the aeSNP; and identifying one or more variant-level biochemical or phenotype effects based on the expression of the different alleles.

17. The method of claim 16, further comprising: performing a c / .s-expression quantitative trait loci (eQTL) correction on the expression of the different alleles at the gene.

18. The method of claim 17, wherein performing the cZs-eQTL correction comprises: calculating a reference ratio that is a proportion of reads coming from the aeSNP relative to total gene expression; predicting an expected ratio based the / .s-eQTL; and subtracting the predicted expected ratio from the reference ratio.

19. The method of claim 16, further comprising: quantifying ASE at the gene to generate a quantified measure of ASE for the aeSNP.

20. The method of claim 19, wherein the quantified measure of ASE is a corrected ratio of sequence reads attributable to each allele.

21. The method of claim 19, further comprising: verifying or calibrating one or more variant effects based on the quantified measure of ASE.

22. The method of claim 21, wherein the one or more variant effects comprise the one or more variant-level biochemical or phenotype effects.

23. The method of claim 19, further comprising: discovering one or more promoter variants using the quantified measure of ASE.

24. The method of claim 16, further comprising: identifying one or more promoter variants based on an allelic imbalance.

Images