Prioritizing disease-related rare variants by integrating gene expression data

Prioritizing disease-related rare variants by integrating gene expression data | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Method Article Prioritizing disease-related rare variants by integrating gene expression data Hanmin Guo, Alexander Eckehart Urban, Wing Hung Wong This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4355589/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Rare variants, comprising a vast majority of human genetic variations, are likely to have more deleterious impact on human diseases compared to common variants. Here we present carrier statistic, a statistical framework to prioritize disease-related rare variants by integrating gene expression data. By quantifying the impact of rare variants on gene expression, carrier statistic can prioritize those rare variants that have large functional consequence in the diseased patients. Through simulation studies and analyzing real multi-omics dataset, we demonstrated that carrier statistic is applicable in studies with limited sample size (a few hundreds) and achieves substantially higher sensitivity than existing rare variants association methods. Application to Alzheimer's disease reveals 16 rare variants within 15 genes with extreme carrier statistics. We also found strong excess of rare variants among the top prioritized genes in diseased patients compared to that in healthy individuals. The carrier statistic method can be applied to various rare variant types and is adaptable to other omics data modalities, offering a powerful tool for investigating the molecular mechanisms underlying complex diseases. Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Introduction Rare variants (minor allele frequency (MAF) < 1%) constitute a vast majority of human genetic variations 1 , 2 . They are on average more deleterious compared with common variants, and thus undergoes stronger selection and remains at low frequency in the general population. By analyzing large cohorts of whole genome sequencing/whole exome sequencing (WGS/WES) data, researchers have identified a handful of rare variants-trait association 3 – 6 and shown that rare variants contribute to a large proportion of missing heritability that cannot be explained by common variants 7 . Rare variants on average confer larger effects on gene expression and complex diseases and are easier to map to causal genes than common variants 8 . However, statistical power to identify disease-associated rare variants, especially for ultra-rare variants or even singletons, is limited, given that the sample size is not too high or the effect size is not too large. Variants collapsing methods (burden test, variance component test, omnibus test) are proposed to circumvent this obstacle 6 , 9 – 11 , which evaluate association for multiple variants in a biologically relevant region, such as a gene, instead of testing the effect of single variant. These methods work well under the assumption that multiple variants in a gene cumulatively contribute to the disease risk with each individual allele explaining only a small fraction of the cases, but resolution to pinpoint the risk variants may be diluted when the assumption is violated. Recent studies based on large scale WES have revealed rare protein truncating variants associated with a wide range of phenotypes 4 , 5 , 12 , 13 . These variants exert extreme effects on the function of genes and their encoded proteins, underscoring the importance of considering how rare variants are related to gene expression. Complementary to genome sequencing assay, RNA-seq can quantitatively measure gene expression level and provide molecular cause of complex diseases, especially rare diseases. Previous studies have shown that rare variants are enriched near genes with aberrant gene expression 14 – 16 . We posit that those rare variants will be more prone to disease pathology. In this work, we propose carrier statistic, a statistical framework for prioritizing those rare variants with large functional consequence in diseased patients by integrating gene expression data. We demonstrate superior performance of our method through extensive simulations and application study to the Alzheimer’s disease, where given a limited sample size, existing rare variants association methods without functional gene expression data cannot provide positive findings. Results Method Overview Our method stems from the expectation that diseased population shows enrichment in rare variants that have large impact on expression for disease-related genes. Suppose we have genotypes (e.g. variants call from WGS data) and gene expression measurements (e.g. reads count from RNA-seq data) on a disease relevant cell type for both diseased patients and healthy controls. For each rare variant-gene pair, we calculate the expression association z-score for rare variant carriers by using the gene expression from individuals without the variant as the null distribution (Fig. 1 a). The expression association z-score, which we term as carrier statistic, is calculated separately within the case and the control group. We only consider rare variant-gene pair wherein the variant is located within the exon of the gene throughout this study. The carrier statistic quantifies the degree to which the rare variant impacts gene expression level. We assume that most rare variants do not have a large impact on the gene function, so the distribution of carrier statistic will be centered around 0. Now, if a gene is relevant for the disease, then conditioning on having the disease will bias the sampling towards people carrying rare variants with large functional impact on the gene expression. Thus, the carrier statistics for disease-related rare variant-gene pairs will tend to be more extreme compared to those for non-related pairs in the case group (Fig. 1 b). We prioritize those rare variants and genes with outlier carrier statistic in the case group. False discovery rate (FDR) can be computed as the ratio of tail probability for carrier statistic between two groups ( Methods ). Simulation Results We first carried out simulations to assess whether the carrier statistic-based method would produce false positive findings. We simulated genotypes based on whole exome sequencing (WES) data from Genome Aggregation Database (gnomAD) 2 and simulated gene expression profiles based on RNA-seq data in whole blood tissue from the Genotype-Tissue Expression (GTEx) project ( Methods ). We perturbed the expression level of the causal genes for causal variants carriers by assuming that causal rare variants have large functional impact on disease-related genes. FDR for carrier statistic was well-calibrated in all simulation settings with varying penetrance of causal variant, prevalence in causal variant noncarriers, and number of causal variants per causal gene ( Supplementary Fig. 1 ). We also checked if using gene expression from noncarriers in two groups together rather than separately as null distribution will produce false positive findings ( Methods ). In this case, we observed substantial inflation in FDR, particularly when diseased patients have systematic change in gene expression profile from healthy controls. In contrast, using gene expression from noncarriers in the same group as null distribution consistently gave well-calibrated error rates ( Supplementary Fig. 2 ). Next, based on the simulated data, we benchmarked the performance of carrier statistic with three existing rare variants association methods: burden test 9 , Sequence Kernel Association Test (SKAT) 10 , and SKAT-O 11 . Burden test counts the number of rare variants within a gene followed by the association test with the disease. SKAT computes a gene-level variance component score statistic which allows bidirectional effect of different variants. The unified test SKAT-O implements a linear combination of burden test statistic and variance component test statistic, which is preferred when the underlying genetic architecture of the disease is not known. While all four methods successfully controlled FDR at the nominal level ( Supplementary Fig. 1 ), carrier statistic achieved higher sensitivity than the three variants collapsing methods under all simulation settings (Fig. 2 ). We believe the low sensitivity of burden-like statistics is due to the small number of case samples that can be attributed to the causal variants in any gene region, which makes it difficult to attain statistical significance of enrichment of rare variants burden for any causal gene region ( Supplementary Tables 1–3 ). We also performed empirical power analysis, which provides guidance for designing new disease association study with genome sequencing and RNA-seq data. Simulations were repeated 50 times to determine the sample size required for achieving 80% sensitivity, given that the penetrance of causal variant is 70%, the prevalence of causal variant noncarrier is 1%, and there are 5 causal variants per causal gene. When the effect size of causal variant on gene expression is large (e.g. Z = 5), 80% causal genes can be identified based on a cohort of 500 cases and 500 controls (Fig. 3 ). The necessary sample size gradually increases as causal variants become less deleterious and have smaller functional consequence on gene expression, e.g. 2,000 samples will be needed to achieve 80% sensitivity for Z = 4.5. In contrast, given the same sample sizes standard GWAS will not have sufficient power to detect rare causal variants regardless of their effect sizes. Application to Alzheimer’s disease Alzheimer’s disease is highly heritable, with heritability estimated to be as high as 60%-80% based on twin studies 17 . Large scale GWASs have identified multiple loci contributing to Alzheimer’s disease, but the genetic variance explained by these loci is far below the level suggested by the disease heritability 18 . Additionally, there is limited understanding regarding the molecular mechanism through which these GWAS variants affect the disease, with the exception of the well-known APOE locus. To investigate whether rare variants (single nucleotide variants [SNVs] and short indels) confer functional consequences in Alzheimer’s disease, we applied carrier statistic to a harmonized multi-omics dataset ( \({n}_{case}=444\) , \({n}_{control}=234\) ) consisting of WGS and RNA-seq from prefrontal cortex in four aging cohort studies: the Religious Orders Study (ROS) and Memory and Aging Project (MAP), the Mount Sinai Brain Bank (MSBB), and the Mayo Clinic ( Methods ). We found significant excess of large carrier statistic in the diseased patients (Fig. 4 ). Controlling FDR with a cutoff of 0.2, we prioritized 16 rare variants within 15 genes with large carrier statistic in the case group (Table 1 ), implicating them as candidate variants that may contribute to Alzheimer’s disease through up-regulating gene expression in the brain. Table 1 16 rare variants within 15 genes with large carrier statistic in the Alzheimer’s disease patients. To see if existing methods can also detect these variants, we applied burden test, SKAT, SKAT-O to the same Alzheimer’s disease dataset. These three variants collapsing methods did not identify any significant genes (FDR < 0.2), possibly due to insufficient sample size. To further evaluate the performance of carrier statistic, we assessed the enrichment of rare variants burden within the top prioritized genes in case group compared to that in controls. Among the top 100 genes with largest carrier statistic, 67 genes have fold enrichment larger than 1, 32 genes have fold enrichment larger than 4/3, while only 6 genes have fold enrichment smaller than 3/4 ( Supplementary Fig. 3 ). Consistent with results from the simulations, enrichment of rare variants burden within each of those genes was moderate and did not pass significance threshold by the variants collapsing methods. The significant genes prioritized by the carrier statistics may shed light on the genetic etiology of Alzheimer’s disease (Fig. 5 ). COCH has the largest carrier statistic of 5.78. Missense mutations within this gene were found to cause the late-onset DFNA9 deafness disorder 19 , 20 . Furthermore, deposits of Cochlin encoded by COCH is associated with age-related glaucomatous trabecular meshwork but absent in healthy controls. Additionally, SNPs inside COCH are associated with cortical thickness 21 , changes in which through neuroimaging techniques is commonly used in early detection and monitoring of Alzheimer's disease progression 22 , 23 . Gene ARHGAP11A has a carrier statistic of 4.62. Transcribed mRNAs of the gene subcellularly localize and are locally translated in radial glia cells of human cerebral cortex and further regulate cortical development 24 . More importantly, ARHGAP11A may contribute to Alzheimer’s disease pathology by mediating Amyloid-β generation and Amyloid-β oligomer neurotoxicity 25 . PPP1R17 (carrier statistic = 4.61) functions as a suppressor of phosphatase complexes 1 (PP1) and 2A (PP2A). A recent study suggests that a subpopulation of neurons in the dorsomedial hypothalamus regulate aging and lifespan in mice through hypothalamic-adipose inter-tissue communication and this regulation depends on Ppp1r17 expression 26 . Interestingly, PPP1R17 is also involved in human-specific cortical neurodevelopment regulated by enhancers in human accelerated regions 27 . ZIC4 (carrier statistic = 5.02) plays an important role in the embryonal development of the cerebellum. Heterozygous deletions encompassing the ZIC4 locus are associated with a rare congenital cerebellar malformation known as the Dandy–Walker malformation 28 . Notably, mutations in proximity to the ZIC4 loci are implicated in multiple system atrophy, a rare neurodegenerative disease 29 . Large scale GWAS study of brain morphology has also identified associations with ZIC4 , underscoring its significance in diverse neurological processes 30 . MDGA1 (carrier statistic = 4.71) encodes a glycosylphosphatidylinositol (GPI)-anchored cell surface glycoprotein. It has been reported that MDGA1 can contribute to cognitive deficits through altering inhibitory synapse development and transmission in the hippocampus 31 . Of note, MDGA1 is one of the 96 genes from the Olink neurology panel with established links to neurobiological processes and neurological diseases. Discussion We presented carrier statistic, a statistical framework to perform multi-omics data analysis, for prioritization of disease-related rare variants and their regulated genes. Through simulations and analyses of real multi-omics dataset, we demonstrated that carrier statistic overcomes sample size limitation and achieves substantial gain in statistical power compared to existing variants collapsing methods. The superior performance of carrier statistic can be attributed to incorporation of functional gene expression data, which allows quantitatively measuring the impact of rare variants that cannot be determined by looking at the variants alone. We applied carrier statistic to Alzheimer’s disease and highlighted several novel risk genes, providing insights into the molecular etiology of the complex disease. Carrier statistic serves as a general approach to study how rare variants affect complex disease through mediating gene expression. There exist several methods such as transcriptome-wide association study (TWAS) 32 , 33 or colocalisation 34 , 35 that can also perform integrative analysis across multiple data modalities (genotype, gene expression, and phenotype). However, those methods focus exclusively on effects of common SNPs and will have limited power for rare variants. In addition to SNVs and short indels that we included in this study, the statistical framework can be also applied to other types of rare variants (simple structural variants [SVs], complex SVs, mobile element insertions, tandem repeat expansions), which in general have larger effect size than SNVs 36 . Finally, carrier statistic can be adapted to other omics data, such as epigenomics and proteomics. Over the last fifteen years, abundant disease-associated loci have been identified based on genome sequences of biobank-scale sample size (e.g. hundreds of thousand) 37 , 38 . However, even with such large sample sizes it is still difficult for GWAS analysis to detect rare causal variants. Another important objective is to understand the biological roles of the detected loci, which remains challenging. We showed here that by integrating RNA-seq data, the carrier statistic approach offers a study design that may overcome the sample size limitation and may help to associate the functional rare variants and their target genes. This approach will be especially useful for the study of diseases for which biospecimen of disease relevant tissue are easy to obtain and RNA-seq can be performed, such as autoimmune disease (relevant to blood) or skin-related disease. Furthermore, when the tissue sample is available, adding RNA-seq to a WGS-based GWAS will not increase the cost of the study significantly. As multi-omics data accumulates alongside genome sequencing data, we anticipate that carrier statistics will become an effective approach to dissect the molecular mechanism of complex diseases. Methods Carrier statistic For each rare variant-gene pair (the variant is located within the exon of the gene), we used the expression of that gene in the rare variant noncarriers as the null distribution and computed a z-score for each rare variant carrier, then average over carriers of that variant. The carrier statistic was computed separately within the case group and the control group. Rare variants were defined as SNVs and short indels whose allele count was no larger than 5 within the case group or within the control group. Therefore, the rare variants and thus the number of carrier statistics are not the same between two groups. The carrier statistic quantifies the degree to which the rare variant impacts gene expression level. We assume that diseased population shows enrichment in rare variants that have large impact on expression for disease-related genes, thus there will also be enrichment of extreme carrier statistic for disease-related rare variant-gene pairs in the case group. We prioritize those rare variants and genes with outlier carrier statistic in the case group. For rare variant-gene pairs with positive carrier statistic, false discovery rate at a given threshold of carrier statistic, denoted by \({z}_{0}\) , can be computed as \(\frac{\text{Pr}\left[{z}_{ctrl} \ge {z}_{0}\right]}{\text{Pr}\left[{z}_{case} \ge {z}_{0}\right]}\) , where \({z}_{case}\) and \({z}_{ctrl}\) denote the carrier statistic in the case group and control group, respectively. Duplicative carrier statistics were removed (i.e. multiple rare variants occurring in the same individuals are counted as the same rare variant). Similarly, for rare variant-gene pairs with negative carrier statistic, false discovery rate at threshold of \({z}_{0}\) can be computed as \(\frac{\text{Pr}\left[{z}_{ctrl} \le {z}_{0}\right]}{\text{Pr}\left[{z}_{case} \le {z}_{0}\right]}\) . Simulations We simulated genotypes for a large population consisting of 125,748 individuals based on the alternative allele count from 125,748 exomes in the gnomAD v2.1.1 dataset 2 . Only exonic variants that passed all variant filters in the gnomAD dataset were retained. Then we simulated gene expression data for the large population as follows. We first simulated background gene expression profile for these 125,748 individuals while matching the mean and standard deviation of normalized gene expression in the reference expression dataset. In this study, we used \({\text{log}}_{2}\left(\text{r}\text{e}\text{a}\text{d}\text{s} \text{c}\text{o}\text{u}\text{n}\text{t} + 1\right)\) as normalized gene expression and RNA-seq in the whole blood tissue from GTEx project v8 as the reference expression dataset 39 . Genes whose median number of reads count in the large population < 10 were removed. We randomly selected \(m\) causal genes and \(l\) causal variants for each causal gene, where \(m\) was set as 50 and \(l\) was set to vary from 1 to 10. For each causal gene, we perturbed the gene expression for causal variant carriers by \(z*sd\) fold, where \(z\) was set as 4 and \(sd\) was the standard deviation of normalized gene expression in the reference dataset. Next, we simulated the disease status for the large population by assuming penetrance of causal variant as \({p}_{carrier}\) and prevalence in causal variant noncarrier as \({p}_{noncarrier}\) . Here \({p}_{carrier}\) varied from 0.5 to 0.9 and \({p}_{noncarrier}\) varied from 0.005 to 0.02. Finally, we randomly sampled 500 cases and 500 controls from the affected and nonaffected population respectively to mimic the sample recruitment procedure in the disease study. Each simulation setting was repeated for 100 times. We evaluated the performance of different methods using two metrics: FDR and sensitivity. FDR was defined as the proportion of falsely identified genes among all identified ones. If no gene was identified then FDR was set as 0. Sensitivity was defined as the proportion of truly identified genes among all underlying causal genes. Note that carrier statistic was computed by using gene expression from rare variant noncarriers in the same group (i.e. case or control) as the carriers as null distribution. We also checked if using gene expression from noncarriers in both case group and control group as the null distribution will produce false positive findings. We perturbed expression level for all genes in the case group. In this case, FDR showed substantial inflation for using gene expression from all individuals in two groups as null distribution, especially when there is large systematic difference in the transcriptome between two groups ( Supplementary Fig. 2b ). On the contrary, using gene expression from rare variant noncarriers in the same group of carriers as null distribution consistently controlled FDR at the nominal level ( Supplementary Fig. 2a ). Implementation of different methods Variants collapsing methods were performed using the R package SKAT v.2.2.5. All the parameters were set as default value. Both common and rare variants were included in the analysis. Multi-omics data analysis for Alzheimer’s disease WGS data for the four aging cohorts (ROS/MAP, MSBB, and the Mayo Clinic) were obtained from the Whole Genome Sequence Harmonization Study (Synapse ID: syn22264775). RNA-seq data for the same four cohorts were downloaded from the RNAseq Harmonization Study (syn21241740). Only white people with both WGS and RNA-seq data were included in the analysis. We determined disease status following description in previous publications 40 , 41 . For the ROS/MAP cohorts, individuals with a Braak neurofibrillary tangle score \(\ge\) 4, a CERAD neuritic and cortical plaque score \(\le\) 2, and a cognitive diagnosis of probable Alzheimer’s disease with no other causes (cogdx = 4) were classified as cases, while individuals with a Braak score \(\le\) 3, a CERAD score \(\ge\) 3, and a cognitive diagnosis of no cognitive impairment (cogdx = 1) were classified as controls. For MSBB, individuals with a Braak score \(\ge\) 4, a CERAD score \(\ge\) 2, and a Clinical Dementia Rating (CDR) score \(\ge\) 1 were classified as cases, while individuals with a Braak score \(\le\) 3, a CERAD score \(\le\) 1, and a CDR score \(\le\) 0.5 were classified as controls. For the Mayo Clinic cohort, individuals with a Braak score \(\ge\) 4 and a CERAD score \(\ge\) 2 were classified as cases, while individuals with a Braak score \(\le\) 3 and a CERAD score \(\le\) 1 were classified as controls. Of note, definition of CERAD score in the ROS/MAP cohort is different from that in the MSBB and the Mayo Clinic cohorts. After harmonization across cohorts, 444 cases and 234 controls in total were identified and used for downstream analysis. Next, we performed quality control on the WGS and RNA-seq data. For WGS data, only exonic variants with missing genotypes \(<\) 10% were retained. For RNA-seq data, we selected prefrontal cortex as the target brain tissue. If a donor does not have RNA-seq in the prefrontal cortex, then RNA-seq in other tissues will be used based on the following order: dorsolateral prefrontal cortex > posterior cingulate cortex > head of caudate nucleus in the ROS/MAP cohort, prefrontal cortex > frontal pole > superior temporal gyrus > inferior frontal gyrus > parahippocampal gyrus in the MSBB cohort, and temporal cortex > cerebellum in the Mayo Clinic cohort. Genes with zero reads count in more than 10% of samples or with median reads count \(<\) 10 across samples were excluded. \({\text{log}}_{2}\left(\text{r}\text{e}\text{a}\text{d}\text{s} \text{c}\text{o}\text{u}\text{n}\text{t} + 1\right)\) was used as normalized gene expression. Then we applied carrier statistic to perform downstream analysis. Declarations Declaration of interests The authors declare no competing interests. Author Contribution W.H.W. and A.E.U. supervised the study. H.G. and W.H.W. designed the idea. H.G. performed data analysis. H.G., A.E.U., and W.H.W. interpreted the results. H.G. and W.H.W. wrote the manuscript. Acknowledgement We thank Dr. Bo Zhou and Dr. Hua Tang for discussion and providing feedbacks. This work was partially supported by NIH grants R01 HG010359, R01 MH116529, and P50 HG007735. Data Availability The WGS data in the Whole Genome Sequence Harmonization Study (https://www.synapse.org/#!Synapse:syn22264775) and the RNA-seq data in the RNAseq Harmonization Study (https://www.synapse.org/#!Synapse:syn21241740) are publicly available on the AD Knowledge Portal platform through completion of a data use certificate. The gnomAD v2.1.1 data consisting of 125,748 exomes were downloaded from https://gnomad.broadinstitute.org/. The gene expression data (gene reads count) from GTEx project version 8 were downloaded from the GTEx Portal, https://gtexportal.org/home/downloads/adult-gtex/bulk_tissue_expression. References Lek M, et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature. 2016;536:285–91. Karczewski KJ, et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature. 2020;581:434–43. Jurgens SJ, et al. Analysis of rare genetic variation underlying cardiometabolic diseases and traits among 200,000 individuals in the UK Biobank. Nat Genet. 2022;54:240–50. Wang Q, et al. Rare variant contribution to human disease in 281,104 UK Biobank exomes. Nature. 2021;597:527–32. Flannick J, et al. Exome sequencing of 20,791 cases of type 2 diabetes and 24,440 controls. Nature. 2019;570:71–6. Li X, et al. Powerful, scalable and resource-efficient meta-analysis of rare variant associations in large whole genome sequencing studies. Nat Genet. 2023;55:154–64. Wainschtein P, et al. Assessing the contribution of rare variants to complex trait heritability from whole-genome sequence data. Nat Genet. 2022;54:263–73. Claussnitzer M, et al. A brief history of human disease genetics. Nature. 2020;577:179–89. Li B, Leal SM. Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. Am J Hum Genet. 2008;83:311–21. Wu MC, et al. Rare-variant association testing for sequencing data with the sequence kernel association test. Am J Hum Genet. 2011;89:82–93. Lee S, et al. Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies. Am J Hum Genet. 2012;91:224–37. Liu D, et al. Schizophrenia risk conferred by rare protein-truncating variants is conserved across diverse human populations. Nat Genet. 2023;55:369–76. DeBoever C, et al. Medical relevance of protein-truncating variants across 337,205 individuals in the UK Biobank study. Nat Commun. 2018;9:1612. Li X, et al. The impact of rare variation on gene expression across tissues. Nature. 2017;550:239–43. Zeng Y, et al. Aberrant gene expression in humans. PLoS Genet. 2015;11:e1004942. Brechtmann F, et al. OUTRIDER: a statistical method for detecting aberrantly expressed genes in RNA sequencing data. Am J Hum Genet. 2018;103:907–17. Gatz M, et al. Role of genes and environments for explaining Alzheimer disease. Arch Gen Psychiatry. 2006;63:168–74. Wightman DP, et al. A genome-wide association study with 1,126,563 individuals identifies new risk loci for Alzheimer’s disease. Nat Genet. 2021;53:1276–82. Bhattacharya SK. Focus on molecules: cochlin. Exp Eye Res. 2006;82:355. Robertson NG, et al. Cochlin immunostaining of inner ear pathologic deposits and proteomic analysis in DFNA9 deafness and vestibular dysfunction. Hum Mol Genet. 2006;15:1071–85. Van Der Meer D, et al. The genetic architecture of human cortical folding. Sci Adv. 2021;7:eabj9446. Querbes O, et al. Early diagnosis of Alzheimer's disease using cortical thickness: impact of cognitive reserve. Brain. 2009;132:2036–47. Schwarz CG, et al. A large-scale comparison of cortical thickness and volume methods for measuring Alzheimer's disease severity. NeuroImage: Clin. 2016;11:802–12. Pilaz L-J et al. Subcellular mRNA localization and local translation of Arhgap11a in radial glial cells regulates cortical development. bioRxiv , 2020.07. 30.229724 (2020). Huang Y-r et al. ArhGAP11A mediates amyloid-β generation and neuropathology in an Alzheimer’s disease-like mouse model. Cell Rep 42(2023). Tokizane K, Brace CS, Imai S. -i. DMHPpp1r17 neurons regulate aging and lifespan in mice through hypothalamic-adipose inter-tissue communication. Cell Metabol. 2024;36:377–92. e11. Girskis KM, et al. Rewiring of human neurodevelopmental gene regulatory programs by human accelerated regions. Neuron. 2021;109:3239–51. Grinberg I, et al. Heterozygous deletion of the linked genes ZIC1 and ZIC4 is involved in Dandy-Walker malformation. Nat Genet. 2004;36:1053–5. Hopfner F, et al. Common Variants Near ZIC1 and ZIC4 in Autopsy-Confirmed Multiple System Atrophy. Mov Disord. 2022;37:2110–21. Zhao B, et al. Genome-wide association analysis of 19,629 individuals identifies variants influencing regional brain volumes and refines their genetic co-architecture with cognitive and mental health traits. Nat Genet. 2019;51:1637–44. Connor SA, et al. Loss of synapse repressor MDGA1 enhances perisomatic inhibition, confers resistance to network excitation, and impairs cognitive function. Cell Rep. 2017;21:3637–45. Gamazon ER, et al. A gene-based association method for mapping traits using reference transcriptome data. Nat Genet. 2015;47:1091–8. Wainberg M, et al. Opportunities and challenges for transcriptome-wide association studies. Nat Genet. 2019;51:592–9. Giambartolomei C, et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet. 2014;10:e1004383. Hormozdiari F, et al. Colocalization of GWAS and eQTL signals detects target genes. Am J Hum Genet. 2016;99:1245–60. Chiang C, et al. The impact of structural variation on human gene expression. Nat Genet. 2017;49:692–9. Tam V, et al. Benefits and limitations of genome-wide association studies. Nat Rev Genet. 2019;20:467–84. Abdellaoui A, Yengo L, Verweij KJ, Visscher PM. 15 years of GWAS discovery: realizing the promise. Am J Hum Genet. 2023;110:179–94. Consortium G. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369:1318–30. Wan Y-W et al. Meta-analysis of the Alzheimer’s disease human brain transcriptome and functional dissection in mouse models. Cell Rep 32(2020). Vialle RA, de Paiva Lopes K, Bennett DA, Crary JF, Raj T. Integrating whole-genome sequencing with multi-omic data reveals the impact of structural variants on gene regulation in the human brain. Nat Neurosci. 2022;25:504–14. Additional Declarations No competing interests reported. Supplementary Files SupplementaryFigures0430.pdf SupplementaryTables0430.xlsx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4355589","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Method Article","associatedPublications":[],"authors":[{"id":299763742,"identity":"aa7fc35b-7616-4a08-8eeb-f7f2f4bc7131","order_by":0,"name":"Hanmin Guo","email":"","orcid":"","institution":"Stanford University","correspondingAuthor":false,"prefix":"","firstName":"Hanmin","middleName":"","lastName":"Guo","suffix":""},{"id":299763743,"identity":"586f70e6-70c1-4eed-93cb-3d4c94fa20ce","order_by":1,"name":"Alexander Eckehart Urban","email":"","orcid":"","institution":"Stanford University School of Medicine","correspondingAuthor":false,"prefix":"","firstName":"Alexander","middleName":"Eckehart","lastName":"Urban","suffix":""},{"id":299763744,"identity":"c3643b3e-43cd-49f3-8049-c7aad6be8684","order_by":2,"name":"Wing Hung Wong","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA10lEQVRIiWNgGAWjYDACCTBpAyY/MDAwA6kEorSkMfAwMDDOIEXLYRK0GNzuMXxc8Ou8vT372YMNH/dYM/Cz5xjg13LnjLHxzL7biT08eYmNM56lM0j2vMGvxexG7jZp3p7bCTwMOeaPeQ4cZjC4QcAWqJZz9jz8bwybQVrsidLC8+MAY49EDkSLgQQBLfY38j8b8zYkJ/bceGPYOONAOo/EmWcFeLVIzkhLfMzzx86evT/HsOHDAWs5/vbkDXi1gAFjG4LNQ1g5GPwhUt0oGAWjYBSMTAAA+B5JReiQ494AAAAASUVORK5CYII=","orcid":"","institution":"Stanford University","correspondingAuthor":true,"prefix":"","firstName":"Wing","middleName":"Hung","lastName":"Wong","suffix":""}],"badges":[],"createdAt":"2024-05-01 18:24:07","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4355589/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4355589/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":56282638,"identity":"8d815ba3-ae31-4cfd-a340-6fb7291cb269","added_by":"auto","created_at":"2024-05-10 21:34:42","extension":"jpeg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":912366,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCarrier statistics for disease-related rare variant-gene pairs will tend to be more extreme compared to those for non-related ones in the case group\u003c/strong\u003e\u003c/p\u003e","description":"","filename":"floatimage1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/097b984ed725225f7afa2d25.jpeg"},{"id":56282700,"identity":"0b4fa8d8-b771-4aba-9a10-e6d8aa15e9e2","added_by":"auto","created_at":"2024-05-10 21:35:16","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":90976,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCarrier statistic achieves higher sensitivity than variants collapsing methods in simulations with varying (a)\u003c/strong\u003e \u003cstrong\u003epenetrance of causal variant, (b) prevalence in causal variant noncarriers, and (c) number of causal variants per causal gene\u003c/strong\u003e\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/2ae11eb651ed5bda07aeba32.png"},{"id":56282732,"identity":"7dfa1c34-78b9-49e7-9bb8-7cf4181e5c68","added_by":"auto","created_at":"2024-05-10 21:35:38","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":36328,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eRequired sample size for carrier statistic to attain 80% sensitivity based on simulations.\u003c/strong\u003e Y-axis denotes total sample size with 50% case-control ratio and x-axis denotes effect size of causal variant on gene expression. Log-transformed sample size is regressed on Z and the fitted curve is shown in red.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/f4778ba0539ab96aa2d44fda.png"},{"id":56282631,"identity":"c6a474f4-3310-4bf2-b503-20a435b608ee","added_by":"auto","created_at":"2024-05-10 21:34:40","extension":"jpeg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":460145,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAlzheimer’s disease patients show significant excess of large carrier statistic\u003c/strong\u003e\u003c/p\u003e","description":"","filename":"floatimage5.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/3748f28e695298d4dfa897bc.jpeg"},{"id":56282694,"identity":"3c6c04da-208c-456b-9776-7105ce3f23fa","added_by":"auto","created_at":"2024-05-10 21:35:11","extension":"jpeg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":946179,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eRare variant carriers show outlier expression level for (a) \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eCOCH\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e, (b) \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eARHGAP11A\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e, and (c) \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003ePPP1R17\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e. \u003c/strong\u003ePseudo count 1 was added to the RNA reads count for visualization purpose. Y-axis is on log scale\u003c/p\u003e","description":"","filename":"floatimage6.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/33932ffcadd837e05156a98a.jpeg"},{"id":56592367,"identity":"75b1f258-9866-41c4-82d9-0431f29ec079","added_by":"auto","created_at":"2024-05-16 09:30:57","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3014686,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/4cebb8e2-6332-4c8b-af5a-7a7797141d9c.pdf"},{"id":56282630,"identity":"8ff0189f-bdb2-43f1-b6fa-404764c262ac","added_by":"auto","created_at":"2024-05-10 21:34:40","extension":"pdf","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":497841,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFigures0430.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/5ccc0049c726e39188dba50a.pdf"},{"id":56282729,"identity":"9f487dd1-593a-447a-9a42-1ae6c0d42894","added_by":"auto","created_at":"2024-05-10 21:35:28","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":14124,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTables0430.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4355589/v1/738de0fae8ab6beac9a3c7fd.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Prioritizing disease-related rare variants by integrating gene expression data","fulltext":[{"header":"Introduction","content":"\u003cp\u003eRare variants (minor allele frequency (MAF)\u0026thinsp;\u0026lt;\u0026thinsp;1%) constitute a vast majority of human genetic variations\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. They are on average more deleterious compared with common variants, and thus undergoes stronger selection and remains at low frequency in the general population. By analyzing large cohorts of whole genome sequencing/whole exome sequencing (WGS/WES) data, researchers have identified a handful of rare variants-trait association\u003csup\u003e\u003cspan additionalcitationids=\"CR4 CR5\" citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e and shown that rare variants contribute to a large proportion of missing heritability that cannot be explained by common variants\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eRare variants on average confer larger effects on gene expression and complex diseases and are easier to map to causal genes than common variants\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. However, statistical power to identify disease-associated rare variants, especially for ultra-rare variants or even singletons, is limited, given that the sample size is not too high or the effect size is not too large. Variants collapsing methods (burden test, variance component test, omnibus test) are proposed to circumvent this obstacle\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e,\u003cspan additionalcitationids=\"CR10\" citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e, which evaluate association for multiple variants in a biologically relevant region, such as a gene, instead of testing the effect of single variant. These methods work well under the assumption that multiple variants in a gene cumulatively contribute to the disease risk with each individual allele explaining only a small fraction of the cases, but resolution to pinpoint the risk variants may be diluted when the assumption is violated. Recent studies based on large scale WES have revealed rare protein truncating variants associated with a wide range of phenotypes\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e,\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e. These variants exert extreme effects on the function of genes and their encoded proteins, underscoring the importance of considering how rare variants are related to gene expression.\u003c/p\u003e \u003cp\u003eComplementary to genome sequencing assay, RNA-seq can quantitatively measure gene expression level and provide molecular cause of complex diseases, especially rare diseases. Previous studies have shown that rare variants are enriched near genes with aberrant gene expression\u003csup\u003e\u003cspan additionalcitationids=\"CR15\" citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e. We posit that those rare variants will be more prone to disease pathology. In this work, we propose carrier statistic, a statistical framework for prioritizing those rare variants with large functional consequence in diseased patients by integrating gene expression data. We demonstrate superior performance of our method through extensive simulations and application study to the Alzheimer\u0026rsquo;s disease, where given a limited sample size, existing rare variants association methods without functional gene expression data cannot provide positive findings.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eMethod Overview\u003c/h2\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eOur method stems from the expectation that diseased population shows enrichment in rare variants that have large impact on expression for disease-related genes. Suppose we have genotypes (e.g. variants call from WGS data) and gene expression measurements (e.g. reads count from RNA-seq data) on a disease relevant cell type for both diseased patients and healthy controls. For each rare variant-gene pair, we calculate the expression association z-score for rare variant carriers by using the gene expression from individuals without the variant as the null distribution (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea). The expression association z-score, which we term as carrier statistic, is calculated separately within the case and the control group. We only consider rare variant-gene pair wherein the variant is located within the exon of the gene throughout this study. The carrier statistic quantifies the degree to which the rare variant impacts gene expression level. We assume that most rare variants do not have a large impact on the gene function, so the distribution of carrier statistic will be centered around 0. Now, if a gene is relevant for the disease, then conditioning on having the disease will bias the sampling towards people carrying rare variants with large functional impact on the gene expression. Thus, the carrier statistics for disease-related rare variant-gene pairs will tend to be more extreme compared to those for non-related pairs in the case group (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb). We prioritize those rare variants and genes with outlier carrier statistic in the case group. False discovery rate (FDR) can be computed as the ratio of tail probability for carrier statistic between two groups (\u003cb\u003eMethods\u003c/b\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eSimulation Results\u003c/h2\u003e \u003cp\u003eWe first carried out simulations to assess whether the carrier statistic-based method would produce false positive findings. We simulated genotypes based on whole exome sequencing (WES) data from Genome Aggregation Database (gnomAD)\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e and simulated gene expression profiles based on RNA-seq data in whole blood tissue from the Genotype-Tissue Expression (GTEx) project (\u003cb\u003eMethods\u003c/b\u003e). We perturbed the expression level of the causal genes for causal variants carriers by assuming that causal rare variants have large functional impact on disease-related genes. FDR for carrier statistic was well-calibrated in all simulation settings with varying penetrance of causal variant, prevalence in causal variant noncarriers, and number of causal variants per causal gene (\u003cb\u003eSupplementary Fig.\u0026nbsp;1\u003c/b\u003e). We also checked if using gene expression from noncarriers in two groups together rather than separately as null distribution will produce false positive findings (\u003cb\u003eMethods\u003c/b\u003e). In this case, we observed substantial inflation in FDR, particularly when diseased patients have systematic change in gene expression profile from healthy controls. In contrast, using gene expression from noncarriers in the same group as null distribution consistently gave well-calibrated error rates (\u003cb\u003eSupplementary Fig.\u0026nbsp;2\u003c/b\u003e).\u003c/p\u003e \u003cp\u003eNext, based on the simulated data, we benchmarked the performance of carrier statistic with three existing rare variants association methods: burden test\u003csup\u003e\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e, Sequence Kernel Association Test (SKAT)\u003csup\u003e\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e, and SKAT-O\u003csup\u003e\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. Burden test counts the number of rare variants within a gene followed by the association test with the disease. SKAT computes a gene-level variance component score statistic which allows bidirectional effect of different variants. The unified test SKAT-O implements a linear combination of burden test statistic and variance component test statistic, which is preferred when the underlying genetic architecture of the disease is not known. While all four methods successfully controlled FDR at the nominal level (\u003cb\u003eSupplementary Fig.\u0026nbsp;1\u003c/b\u003e), carrier statistic achieved higher sensitivity than the three variants collapsing methods under all simulation settings (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). We believe the low sensitivity of burden-like statistics is due to the small number of case samples that can be attributed to the causal variants in any gene region, which makes it difficult to attain statistical significance of enrichment of rare variants burden for any causal gene region (\u003cb\u003eSupplementary Tables\u0026nbsp;1\u0026ndash;3\u003c/b\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe also performed empirical power analysis, which provides guidance for designing new disease association study with genome sequencing and RNA-seq data. Simulations were repeated 50 times to determine the sample size required for achieving 80% sensitivity, given that the penetrance of causal variant is 70%, the prevalence of causal variant noncarrier is 1%, and there are 5 causal variants per causal gene. When the effect size of causal variant on gene expression is large (e.g. Z\u0026thinsp;=\u0026thinsp;5), 80% causal genes can be identified based on a cohort of 500 cases and 500 controls (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). The necessary sample size gradually increases as causal variants become less deleterious and have smaller functional consequence on gene expression, e.g. 2,000 samples will be needed to achieve 80% sensitivity for Z\u0026thinsp;=\u0026thinsp;4.5. In contrast, given the same sample sizes standard GWAS will not have sufficient power to detect rare causal variants regardless of their effect sizes.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eApplication to Alzheimer\u0026rsquo;s disease\u003c/h2\u003e \u003cp\u003eAlzheimer\u0026rsquo;s disease is highly heritable, with heritability estimated to be as high as 60%-80% based on twin studies\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e. Large scale GWASs have identified multiple loci contributing to Alzheimer\u0026rsquo;s disease, but the genetic variance explained by these loci is far below the level suggested by the disease heritability\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e. Additionally, there is limited understanding regarding the molecular mechanism through which these GWAS variants affect the disease, with the exception of the well-known APOE locus. To investigate whether rare variants (single nucleotide variants [SNVs] and short indels) confer functional consequences in Alzheimer\u0026rsquo;s disease, we applied carrier statistic to a harmonized multi-omics dataset (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({n}_{case}=444\\)\u003c/span\u003e\u003c/span\u003e, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({n}_{control}=234\\)\u003c/span\u003e\u003c/span\u003e) consisting of WGS and RNA-seq from prefrontal cortex in four aging cohort studies: the Religious Orders Study (ROS) and Memory and Aging Project (MAP), the Mount Sinai Brain Bank (MSBB), and the Mayo Clinic (\u003cb\u003eMethods\u003c/b\u003e). We found significant excess of large carrier statistic in the diseased patients (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). Controlling FDR with a cutoff of 0.2, we prioritized 16 rare variants within 15 genes with large carrier statistic in the case group (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e), implicating them as candidate variants that may contribute to Alzheimer\u0026rsquo;s disease through up-regulating gene expression in the brain.\u003c/p\u003e \u003cp\u003e\u003cb\u003eTable 1\u003c/b\u003e 16 rare variants within 15 genes with large carrier statistic in the Alzheimer\u0026rsquo;s disease patients.\u003c/p\u003e \u003cp\u003e\u003cimg src=\"https://myfiles.space/user_files/127393_c7e80a1c9bb65875/127393_custom_files/img1715180965.png\"\u003e\u003cbr\u003e\u003c/p\u003e \u003cp\u003eTo see if existing methods can also detect these variants, we applied burden test, SKAT, SKAT-O to the same Alzheimer\u0026rsquo;s disease dataset. These three variants collapsing methods did not identify any significant genes (FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.2), possibly due to insufficient sample size. To further evaluate the performance of carrier statistic, we assessed the enrichment of rare variants burden within the top prioritized genes in case group compared to that in controls. Among the top 100 genes with largest carrier statistic, 67 genes have fold enrichment larger than 1, 32 genes have fold enrichment larger than 4/3, while only 6 genes have fold enrichment smaller than 3/4 (\u003cb\u003eSupplementary Fig.\u0026nbsp;3\u003c/b\u003e). Consistent with results from the simulations, enrichment of rare variants burden within each of those genes was moderate and did not pass significance threshold by the variants collapsing methods.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe significant genes prioritized by the carrier statistics may shed light on the genetic etiology of Alzheimer\u0026rsquo;s disease (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e). \u003cem\u003eCOCH\u003c/em\u003e has the largest carrier statistic of 5.78. Missense mutations within this gene were found to cause the late-onset DFNA9 deafness disorder\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e,\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e. Furthermore, deposits of Cochlin encoded by \u003cem\u003eCOCH\u003c/em\u003e is associated with age-related glaucomatous trabecular meshwork but absent in healthy controls. Additionally, SNPs inside \u003cem\u003eCOCH\u003c/em\u003e are associated with cortical thickness\u003csup\u003e\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e, changes in which through neuroimaging techniques is commonly used in early detection and monitoring of Alzheimer's disease progression\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e,\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e. Gene \u003cem\u003eARHGAP11A\u003c/em\u003e has a carrier statistic of 4.62. Transcribed mRNAs of the gene subcellularly localize and are locally translated in radial glia cells of human cerebral cortex and further regulate cortical development\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u003c/sup\u003e. More importantly, \u003cem\u003eARHGAP11A\u003c/em\u003e may contribute to Alzheimer\u0026rsquo;s disease pathology by mediating Amyloid-β generation and Amyloid-β oligomer neurotoxicity\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. \u003cem\u003ePPP1R17\u003c/em\u003e (carrier statistic\u0026thinsp;=\u0026thinsp;4.61) functions as a suppressor of phosphatase complexes 1 (PP1) and 2A (PP2A). A recent study suggests that a subpopulation of neurons in the dorsomedial hypothalamus regulate aging and lifespan in mice through hypothalamic-adipose inter-tissue communication and this regulation depends on Ppp1r17 expression\u003csup\u003e\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e\u003c/sup\u003e. Interestingly, \u003cem\u003ePPP1R17\u003c/em\u003e is also involved in human-specific cortical neurodevelopment regulated by enhancers in human accelerated regions\u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e. \u003cem\u003eZIC4\u003c/em\u003e (carrier statistic\u0026thinsp;=\u0026thinsp;5.02) plays an important role in the embryonal development of the cerebellum. Heterozygous deletions encompassing the \u003cem\u003eZIC4\u003c/em\u003e locus are associated with a rare congenital cerebellar malformation known as the Dandy\u0026ndash;Walker malformation\u003csup\u003e\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. Notably, mutations in proximity to the \u003cem\u003eZIC4\u003c/em\u003e loci are implicated in multiple system atrophy, a rare neurodegenerative disease\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e. Large scale GWAS study of brain morphology has also identified associations with \u003cem\u003eZIC4\u003c/em\u003e, underscoring its significance in diverse neurological processes\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. \u003cem\u003eMDGA1\u003c/em\u003e (carrier statistic\u0026thinsp;=\u0026thinsp;4.71) encodes a glycosylphosphatidylinositol (GPI)-anchored cell surface glycoprotein. It has been reported that \u003cem\u003eMDGA1\u003c/em\u003e can contribute to cognitive deficits through altering inhibitory synapse development and transmission in the hippocampus\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. Of note, \u003cem\u003eMDGA1\u003c/em\u003e is one of the 96 genes from the Olink neurology panel with established links to neurobiological processes and neurological diseases.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eWe presented carrier statistic, a statistical framework to perform multi-omics data analysis, for prioritization of disease-related rare variants and their regulated genes. Through simulations and analyses of real multi-omics dataset, we demonstrated that carrier statistic overcomes sample size limitation and achieves substantial gain in statistical power compared to existing variants collapsing methods. The superior performance of carrier statistic can be attributed to incorporation of functional gene expression data, which allows quantitatively measuring the impact of rare variants that cannot be determined by looking at the variants alone. We applied carrier statistic to Alzheimer\u0026rsquo;s disease and highlighted several novel risk genes, providing insights into the molecular etiology of the complex disease.\u003c/p\u003e \u003cp\u003eCarrier statistic serves as a general approach to study how rare variants affect complex disease through mediating gene expression. There exist several methods such as transcriptome-wide association study (TWAS)\u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e,\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e or colocalisation\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e,\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e that can also perform integrative analysis across multiple data modalities (genotype, gene expression, and phenotype). However, those methods focus exclusively on effects of common SNPs and will have limited power for rare variants. In addition to SNVs and short indels that we included in this study, the statistical framework can be also applied to other types of rare variants (simple structural variants [SVs], complex SVs, mobile element insertions, tandem repeat expansions), which in general have larger effect size than SNVs\u003csup\u003e\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e. Finally, carrier statistic can be adapted to other omics data, such as epigenomics and proteomics.\u003c/p\u003e \u003cp\u003eOver the last fifteen years, abundant disease-associated loci have been identified based on genome sequences of biobank-scale sample size (e.g. hundreds of thousand)\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e,\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e. However, even with such large sample sizes it is still difficult for GWAS analysis to detect rare causal variants. Another important objective is to understand the biological roles of the detected loci, which remains challenging. We showed here that by integrating RNA-seq data, the carrier statistic approach offers a study design that may overcome the sample size limitation and may help to associate the functional rare variants and their target genes. This approach will be especially useful for the study of diseases for which biospecimen of disease relevant tissue are easy to obtain and RNA-seq can be performed, such as autoimmune disease (relevant to blood) or skin-related disease. Furthermore, when the tissue sample is available, adding RNA-seq to a WGS-based GWAS will not increase the cost of the study significantly. As multi-omics data accumulates alongside genome sequencing data, we anticipate that carrier statistics will become an effective approach to dissect the molecular mechanism of complex diseases.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eCarrier statistic\u003c/h2\u003e \u003cp\u003eFor each rare variant-gene pair (the variant is located within the exon of the gene), we used the expression of that gene in the rare variant noncarriers as the null distribution and computed a z-score for each rare variant carrier, then average over carriers of that variant. The carrier statistic was computed separately within the case group and the control group. Rare variants were defined as SNVs and short indels whose allele count was no larger than 5 within the case group or within the control group. Therefore, the rare variants and thus the number of carrier statistics are not the same between two groups. The carrier statistic quantifies the degree to which the rare variant impacts gene expression level. We assume that diseased population shows enrichment in rare variants that have large impact on expression for disease-related genes, thus there will also be enrichment of extreme carrier statistic for disease-related rare variant-gene pairs in the case group. We prioritize those rare variants and genes with outlier carrier statistic in the case group. For rare variant-gene pairs with positive carrier statistic, false discovery rate at a given threshold of carrier statistic, denoted by \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({z}_{0}\\)\u003c/span\u003e\u003c/span\u003e, can be computed as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\frac{\\text{Pr}\\left[{z}_{ctrl} \\ge {z}_{0}\\right]}{\\text{Pr}\\left[{z}_{case} \\ge {z}_{0}\\right]}\\)\u003c/span\u003e\u003c/span\u003e, where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({z}_{case}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({z}_{ctrl}\\)\u003c/span\u003e\u003c/span\u003e denote the carrier statistic in the case group and control group, respectively. Duplicative carrier statistics were removed (i.e. multiple rare variants occurring in the same individuals are counted as the same rare variant). Similarly, for rare variant-gene pairs with negative carrier statistic, false discovery rate at threshold of \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({z}_{0}\\)\u003c/span\u003e\u003c/span\u003e can be computed as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\frac{\\text{Pr}\\left[{z}_{ctrl} \\le {z}_{0}\\right]}{\\text{Pr}\\left[{z}_{case} \\le {z}_{0}\\right]}\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eSimulations\u003c/h2\u003e \u003cp\u003eWe simulated genotypes for a large population consisting of 125,748 individuals based on the alternative allele count from 125,748 exomes in the gnomAD v2.1.1 dataset\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. Only exonic variants that passed all variant filters in the gnomAD dataset were retained. Then we simulated gene expression data for the large population as follows. We first simulated background gene expression profile for these 125,748 individuals while matching the mean and standard deviation of normalized gene expression in the reference expression dataset. In this study, we used \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{log}}_{2}\\left(\\text{r}\\text{e}\\text{a}\\text{d}\\text{s} \\text{c}\\text{o}\\text{u}\\text{n}\\text{t} + 1\\right)\\)\u003c/span\u003e\u003c/span\u003e as normalized gene expression and RNA-seq in the whole blood tissue from GTEx project v8 as the reference expression dataset\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u003c/sup\u003e. Genes whose median number of reads count in the large population \u0026lt; 10 were removed. We randomly selected \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(m\\)\u003c/span\u003e\u003c/span\u003e causal genes and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(l\\)\u003c/span\u003e\u003c/span\u003e causal variants for each causal gene, where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(m\\)\u003c/span\u003e\u003c/span\u003e was set as 50 and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(l\\)\u003c/span\u003e\u003c/span\u003e was set to vary from 1 to 10. For each causal gene, we perturbed the gene expression for causal variant carriers by \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(z*sd\\)\u003c/span\u003e\u003c/span\u003e fold, where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(z\\)\u003c/span\u003e\u003c/span\u003e was set as 4 and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(sd\\)\u003c/span\u003e\u003c/span\u003e was the standard deviation of normalized gene expression in the reference dataset. Next, we simulated the disease status for the large population by assuming penetrance of causal variant as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({p}_{carrier}\\)\u003c/span\u003e\u003c/span\u003e and prevalence in causal variant noncarrier as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({p}_{noncarrier}\\)\u003c/span\u003e\u003c/span\u003e. Here \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({p}_{carrier}\\)\u003c/span\u003e\u003c/span\u003e varied from 0.5 to 0.9 and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({p}_{noncarrier}\\)\u003c/span\u003e\u003c/span\u003e varied from 0.005 to 0.02. Finally, we randomly sampled 500 cases and 500 controls from the affected and nonaffected population respectively to mimic the sample recruitment procedure in the disease study. Each simulation setting was repeated for 100 times.\u003c/p\u003e \u003cp\u003eWe evaluated the performance of different methods using two metrics: FDR and sensitivity. FDR was defined as the proportion of falsely identified genes among all identified ones. If no gene was identified then FDR was set as 0. Sensitivity was defined as the proportion of truly identified genes among all underlying causal genes.\u003c/p\u003e \u003cp\u003eNote that carrier statistic was computed by using gene expression from rare variant noncarriers in the same group (i.e. case or control) as the carriers as null distribution. We also checked if using gene expression from noncarriers in both case group and control group as the null distribution will produce false positive findings. We perturbed expression level for all genes in the case group. In this case, FDR showed substantial inflation for using gene expression from all individuals in two groups as null distribution, especially when there is large systematic difference in the transcriptome between two groups (\u003cb\u003eSupplementary Fig.\u0026nbsp;2b\u003c/b\u003e). On the contrary, using gene expression from rare variant noncarriers in the same group of carriers as null distribution consistently controlled FDR at the nominal level (\u003cb\u003eSupplementary Fig.\u0026nbsp;2a\u003c/b\u003e).\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eImplementation of different methods\u003c/h3\u003e\n\u003cp\u003eVariants collapsing methods were performed using the R package SKAT v.2.2.5. All the parameters were set as default value. Both common and rare variants were included in the analysis.\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eMulti-omics data analysis for Alzheimer\u0026rsquo;s disease\u003c/h2\u003e \u003cp\u003eWGS data for the four aging cohorts (ROS/MAP, MSBB, and the Mayo Clinic) were obtained from the Whole Genome Sequence Harmonization Study (Synapse ID: syn22264775). RNA-seq data for the same four cohorts were downloaded from the RNAseq Harmonization Study (syn21241740). Only white people with both WGS and RNA-seq data were included in the analysis.\u003c/p\u003e \u003cp\u003eWe determined disease status following description in previous publications\u003csup\u003e\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e,\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e. For the ROS/MAP cohorts, individuals with a Braak neurofibrillary tangle score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 4, a CERAD neuritic and cortical plaque score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 2, and a cognitive diagnosis of probable Alzheimer\u0026rsquo;s disease with no other causes (cogdx = 4) were classified as cases, while individuals with a Braak score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 3, a CERAD score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 3, and a cognitive diagnosis of no cognitive impairment (cogdx = 1) were classified as controls. For MSBB, individuals with a Braak score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 4, a CERAD score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 2, and a Clinical Dementia Rating (CDR) score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 1 were classified as cases, while individuals with a Braak score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 3, a CERAD score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 1, and a CDR score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 0.5 were classified as controls. For the Mayo Clinic cohort, individuals with a Braak score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 4 and a CERAD score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\ge\\)\u003c/span\u003e\u003c/span\u003e 2 were classified as cases, while individuals with a Braak score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 3 and a CERAD score \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\le\\)\u003c/span\u003e\u003c/span\u003e 1 were classified as controls. Of note, definition of CERAD score in the ROS/MAP cohort is different from that in the MSBB and the Mayo Clinic cohorts. After harmonization across cohorts, 444 cases and 234 controls in total were identified and used for downstream analysis.\u003c/p\u003e \u003cp\u003eNext, we performed quality control on the WGS and RNA-seq data. For WGS data, only exonic variants with missing genotypes \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\u0026lt;\\)\u003c/span\u003e\u003c/span\u003e 10% were retained. For RNA-seq data, we selected prefrontal cortex as the target brain tissue. If a donor does not have RNA-seq in the prefrontal cortex, then RNA-seq in other tissues will be used based on the following order: dorsolateral prefrontal cortex \u0026gt; posterior cingulate cortex \u0026gt; head of caudate nucleus in the ROS/MAP cohort, prefrontal cortex \u0026gt; frontal pole \u0026gt; superior temporal gyrus \u0026gt; inferior frontal gyrus \u0026gt; parahippocampal gyrus in the MSBB cohort, and temporal cortex \u0026gt; cerebellum in the Mayo Clinic cohort. Genes with zero reads count in more than 10% of samples or with median reads count \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\u0026lt;\\)\u003c/span\u003e\u003c/span\u003e 10 across samples were excluded. \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\text{log}}_{2}\\left(\\text{r}\\text{e}\\text{a}\\text{d}\\text{s} \\text{c}\\text{o}\\text{u}\\text{n}\\text{t} + 1\\right)\\)\u003c/span\u003e\u003c/span\u003e was used as normalized gene expression. Then we applied carrier statistic to perform downstream analysis.\u003c/p\u003e \u003c/div\u003e "},{"header":"Declarations","content":" \u003ch2\u003eDeclaration of interests\u003c/h2\u003e \u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eW.H.W. and A.E.U. supervised the study. H.G. and W.H.W. designed the idea. H.G. performed data analysis. H.G., A.E.U., and W.H.W. interpreted the results. H.G. and W.H.W. wrote the manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe thank Dr. Bo Zhou and Dr. Hua Tang for discussion and providing feedbacks. This work was partially supported by NIH grants R01 HG010359, R01 MH116529, and P50 HG007735.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe WGS data in the Whole Genome Sequence Harmonization Study (https://www.synapse.org/#!Synapse:syn22264775) and the RNA-seq data in the RNAseq Harmonization Study (https://www.synapse.org/#!Synapse:syn21241740) are publicly available on the AD Knowledge Portal platform through completion of a data use certificate. The gnomAD v2.1.1 data consisting of 125,748 exomes were downloaded from https://gnomad.broadinstitute.org/. The gene expression data (gene reads count) from GTEx project version 8 were downloaded from the GTEx Portal, https://gtexportal.org/home/downloads/adult-gtex/bulk_tissue_expression.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eLek M, et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature. 2016;536:285\u0026ndash;91.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKarczewski KJ, et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature. 2020;581:434\u0026ndash;43.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJurgens SJ, et al. Analysis of rare genetic variation underlying cardiometabolic diseases and traits among 200,000 individuals in the UK Biobank. Nat Genet. 2022;54:240\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang Q, et al. Rare variant contribution to human disease in 281,104 UK Biobank exomes. Nature. 2021;597:527\u0026ndash;32.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFlannick J, et al. Exome sequencing of 20,791 cases of type 2 diabetes and 24,440 controls. Nature. 2019;570:71\u0026ndash;6.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi X, et al. Powerful, scalable and resource-efficient meta-analysis of rare variant associations in large whole genome sequencing studies. Nat Genet. 2023;55:154\u0026ndash;64.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWainschtein P, et al. Assessing the contribution of rare variants to complex trait heritability from whole-genome sequence data. Nat Genet. 2022;54:263\u0026ndash;73.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eClaussnitzer M, et al. A brief history of human disease genetics. Nature. 2020;577:179\u0026ndash;89.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi B, Leal SM. Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. Am J Hum Genet. 2008;83:311\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWu MC, et al. Rare-variant association testing for sequencing data with the sequence kernel association test. Am J Hum Genet. 2011;89:82\u0026ndash;93.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLee S, et al. Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies. Am J Hum Genet. 2012;91:224\u0026ndash;37.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu D, et al. Schizophrenia risk conferred by rare protein-truncating variants is conserved across diverse human populations. Nat Genet. 2023;55:369\u0026ndash;76.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDeBoever C, et al. Medical relevance of protein-truncating variants across 337,205 individuals in the UK Biobank study. Nat Commun. 2018;9:1612.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi X, et al. The impact of rare variation on gene expression across tissues. Nature. 2017;550:239\u0026ndash;43.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZeng Y, et al. Aberrant gene expression in humans. PLoS Genet. 2015;11:e1004942.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBrechtmann F, et al. OUTRIDER: a statistical method for detecting aberrantly expressed genes in RNA sequencing data. Am J Hum Genet. 2018;103:907\u0026ndash;17.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGatz M, et al. Role of genes and environments for explaining Alzheimer disease. Arch Gen Psychiatry. 2006;63:168\u0026ndash;74.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWightman DP, et al. A genome-wide association study with 1,126,563 individuals identifies new risk loci for Alzheimer\u0026rsquo;s disease. Nat Genet. 2021;53:1276\u0026ndash;82.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBhattacharya SK. Focus on molecules: cochlin. Exp Eye Res. 2006;82:355.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRobertson NG, et al. Cochlin immunostaining of inner ear pathologic deposits and proteomic analysis in DFNA9 deafness and vestibular dysfunction. Hum Mol Genet. 2006;15:1071\u0026ndash;85.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVan Der Meer D, et al. The genetic architecture of human cortical folding. Sci Adv. 2021;7:eabj9446.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eQuerbes O, et al. Early diagnosis of Alzheimer's disease using cortical thickness: impact of cognitive reserve. Brain. 2009;132:2036\u0026ndash;47.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSchwarz CG, et al. A large-scale comparison of cortical thickness and volume methods for measuring Alzheimer's disease severity. NeuroImage: Clin. 2016;11:802\u0026ndash;12.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePilaz L-J et al. Subcellular mRNA localization and local translation of Arhgap11a in radial glial cells regulates cortical development. \u003cem\u003ebioRxiv\u003c/em\u003e, 2020.07. 30.229724 (2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHuang Y-r et al. ArhGAP11A mediates amyloid-β generation and neuropathology in an Alzheimer\u0026rsquo;s disease-like mouse model. Cell Rep 42(2023).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTokizane K, Brace CS, Imai S. -i. DMHPpp1r17 neurons regulate aging and lifespan in mice through hypothalamic-adipose inter-tissue communication. Cell Metabol. 2024;36:377\u0026ndash;92. e11.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGirskis KM, et al. Rewiring of human neurodevelopmental gene regulatory programs by human accelerated regions. Neuron. 2021;109:3239\u0026ndash;51.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGrinberg I, et al. Heterozygous deletion of the linked genes ZIC1 and ZIC4 is involved in Dandy-Walker malformation. Nat Genet. 2004;36:1053\u0026ndash;5.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHopfner F, et al. Common Variants Near ZIC1 and ZIC4 in Autopsy-Confirmed Multiple System Atrophy. Mov Disord. 2022;37:2110\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhao B, et al. Genome-wide association analysis of 19,629 individuals identifies variants influencing regional brain volumes and refines their genetic co-architecture with cognitive and mental health traits. Nat Genet. 2019;51:1637\u0026ndash;44.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eConnor SA, et al. Loss of synapse repressor MDGA1 enhances perisomatic inhibition, confers resistance to network excitation, and impairs cognitive function. Cell Rep. 2017;21:3637\u0026ndash;45.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGamazon ER, et al. A gene-based association method for mapping traits using reference transcriptome data. Nat Genet. 2015;47:1091\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWainberg M, et al. Opportunities and challenges for transcriptome-wide association studies. Nat Genet. 2019;51:592\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGiambartolomei C, et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet. 2014;10:e1004383.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHormozdiari F, et al. Colocalization of GWAS and eQTL signals detects target genes. Am J Hum Genet. 2016;99:1245\u0026ndash;60.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChiang C, et al. The impact of structural variation on human gene expression. Nat Genet. 2017;49:692\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTam V, et al. Benefits and limitations of genome-wide association studies. Nat Rev Genet. 2019;20:467\u0026ndash;84.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAbdellaoui A, Yengo L, Verweij KJ, Visscher PM. 15 years of GWAS discovery: realizing the promise. Am J Hum Genet. 2023;110:179\u0026ndash;94.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eConsortium G. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369:1318\u0026ndash;30.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWan Y-W et al. Meta-analysis of the Alzheimer\u0026rsquo;s disease human brain transcriptome and functional dissection in mouse models. Cell Rep 32(2020).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVialle RA, de Paiva Lopes K, Bennett DA, Crary JF, Raj T. Integrating whole-genome sequencing with multi-omic data reveals the impact of structural variants on gene regulation in the human brain. Nat Neurosci. 2022;25:504\u0026ndash;14.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-4355589/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4355589/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eRare variants, comprising a vast majority of human genetic variations, are likely to have more deleterious impact on human diseases compared to common variants. Here we present carrier statistic, a statistical framework to prioritize disease-related rare variants by integrating gene expression data. By quantifying the impact of rare variants on gene expression, carrier statistic can prioritize those rare variants that have large functional consequence in the diseased patients. Through simulation studies and analyzing real multi-omics dataset, we demonstrated that carrier statistic is applicable in studies with limited sample size (a few hundreds) and achieves substantially higher sensitivity than existing rare variants association methods. Application to Alzheimer's disease reveals 16 rare variants within 15 genes with extreme carrier statistics. We also found strong excess of rare variants among the top prioritized genes in diseased patients compared to that in healthy individuals. The carrier statistic method can be applied to various rare variant types and is adaptable to other omics data modalities, offering a powerful tool for investigating the molecular mechanisms underlying complex diseases.\u003c/p\u003e","manuscriptTitle":"Prioritizing disease-related rare variants by integrating gene expression data","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-05-10 21:19:09","doi":"10.21203/rs.3.rs-4355589/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"9cf1bc79-1125-4c5d-9744-334cb8d66967","owner":[],"postedDate":"May 10th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-05-16T09:22:49+00:00","versionOfRecord":[],"versionCreatedAt":"2024-05-10 21:19:09","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4355589","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4355589","identity":"rs-4355589","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}