Impact of Rare Non-coding Variants on Human Diseases through Alternative Polyadenylation Outliers

Impact of Rare Non-coding Variants on Human Diseases through Alternative Polyadenylation Outliers | 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 Article Impact of Rare Non-coding Variants on Human Diseases through Alternative Polyadenylation Outliers Lei Li, Xudong Zou, Zhaozhao Zhao, Yu Chen, Kewei Xiong, Zeyang Wang, and 6 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3907149/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 16 Jan, 2025 Read the published version in Nature Communications → Version 1 posted You are reading this latest preprint version Abstract Although rare non-coding variants (RVs) play crucial roles in human complex traits and diseases, understanding their functional mechanisms and identifying those most closely associated with diseases continue to be major challenges. Here, we constructed the first comprehensive atlas of alternative polyadenylation (APA) outliers (aOutliers) from 15,201 samples across 49 human tissues. Strikingly, these aOutliers exhibit unique characteristics markedly distinct from those of outliers based on transcriptional abundance or splicing. This is evidenced by a pronounced enrichment of RVs specifically within aOutliers. Mechanistically, aOutlier RVs frequently alter poly(A) signals and splicing sites, and experimental perturbation of these RVs indeed triggers APA events. Furthermore, we developed a Bayesian-based APA RV prediction model, which successfully pinpointed a specific set of RVs with significantly large effect sizes on complex traits or diseases. A particularly intriguing discovery was the observed convergence effect on APA between rare and common cancer variants, exemplified by the combinatorial regulation of APA in the DDX18 gene. Together, this study introduces a novel APA-enhanced framework for individual genome annotation and underscores the importance of APA in uncovering previously unrecognized functional non-coding RVs linked to human complex traits and diseases. Biological sciences/Genetics/Gene regulation Biological sciences/Computational biology and bioinformatics/Data mining Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction The human genome harbors numerous rare genetic variants 1 , each with a minor allele frequency (MAF) of less than 1%. Many of these rare variants strongly contribute to human diseases 2–5 . While exome sequencing of large population cohorts has identified numerous rare protein-coding variants associated with both common and rare diseases 6 , the vast majority of rare variants (RVs) are located in non-coding regions. These non-coding RVs do not function through altering the protein sequences, thereby posing a significant challenge in interpreting their functions. To address this challenge, analysis of population-scale transcriptomic data has been used to uncover functional rare non-coding variants affecting gene expression or splicing outliers 7–10 . Despite these efforts, a significant portion of disease-associated RVs remain uncharacterized. Alternative polyadenylation (APA) of mRNA is a widespread post-transcriptional regulatory mechanism observed across various species. By employing different polyadenylation sites within 3′untranslated regions (3′ UTRs), genes can produce various mRNA isoforms with either shortened or extended 3′ UTRs. These 3′ UTRs contain many regulatory elements that modulate the abundance or localization of the mRNA and protein 11–15 . Moreover, APA can also occur in intronic regions, leading to truncated mRNA or proteins 16,17 . Accordingly, disruptions in APA events have been increasingly implicated in many human diseases 17–19 . For example, altered APA leading to 3′ UTR shortening of competing-endogenous RNAs for tumor suppressor genes can result in the release of microRNAs, inhibiting tumor suppressor genes and potentially leading to tumorigenesis 20 . Moreover, recent studies have reported the ubiquitous genetic regulation of APA, highlighting its importance in the functional interpretation of disease-associated non-coding variants 21–23 . A notable example is a single-nucleotide polymorphism (SNP; rs10954213) within the 3′ UTR of interferon regulatory factor 5 ( IRF5) , which can alter the length and stability of its 3′ UTR, thereby contributing to systemic lupus erythematosus susceptibility 24 . In our previous study, we built an atlas of human 3′ UTR APA quantitative trait loci (3′aQTLs) across human tissues, identifying approximately 0.4 million common SNPs associated with interindividual APA changes, which colocalize with 16.1% of trait-associated genetic variants 25 . Yet, these studies mainly focus on assessing the APA regulation of common variants. To our knowledge, the effect of RVs on APA has not been explored. Here, to better understand the impact of RVs on APA, we systematically analyzed aberrant APA events across 49 human tissues from the Genotype-Tissue Expression Project (GTEx). We identified 1,534 multi-tissue APA outliers (aOutliers) from European individuals. Intriguingly, 74.2% of these aOutliers are associated with genes not previously identified in outlier analysis of other molecular phenotypes (e.g., expression or splicing). These aOutliers exhibit distinct characteristics, such as unique 3′ UTR length and GC-contents, setting them apart from other types of molecular outliers. Moreover, a significant enrichment of deleterious RVs was observed in regions proximal to these aOutliers. To prioritize functional RVs impacting APA, we developed a Bayesian hierarchical model and identified a distinct set of RVs with large effect sizes on human complex traits and disease phenotypes. Intriguingly, we observed and demonstrated strong convergence effects between prioritized RVs and common variants in regulating 3′ UTR APA, exemplified by the combinatorial regulation of APA in DDX18 . Lastly, to facilitate broad access to aOutliers-associated RVs, we have constructed a user-friendly portal at http://bioinfo.szbl.ac.cn/rareAPA/index.php . Collectively, our findings indicate that APA highlights a specific set of RVs with significant impacts on human traits and diseases, providing a new avenue for interpreting rare human non-coding genetic variants. Results The landscape of APA outliers across 49 human tissues We first conducted a comprehensive identification of 3′ UTR and intronic APA events in 15,201 GTEx RNA-seq samples from 49 human tissues of 838 individuals (Fig. 1 a) using our Dapars2 25,26 and IPAFinder 18 algorithms, respectively (see Materials and Methods) (Supplementary Fig. 1). Considering the potential influence of many known and unknown technical confounders on APA usage among samples, we regressed out these confounders, such as age, sex, sequencing platform, and other hidden confounders inferred by using probabilistic estimation of expression residuals (PEER) factors (Supplementary Fig. 2). We then calculated Z-scores for the PEER-adjusted 3′ UTR and intronic APA usage in each tissue to identify individuals with aberrant APA usage for a specific gene, which we refer to as APA outliers (aOutliers) with an absolute Z-score > 3. The individuals and genes were designated as “aOutlier individuals” and “aOutlier genes”, respectively. Importantly, a single gene could be associated with multiple outlier individuals, and conversely, one individual could be an aOutlier individual for multiple genes. Our analysis of these aOutliers revealed that, on average, 68.5% of all transcripts per tissue were present in at least one outlier individual (Supplementary Fig. 3a). The number of aOutlier genes strongly correlated (Spearman’s correlation rho = 0.91, P < 2.2 × 10 ‒16 ) with sample size across tissues (Supplementary Fig. 3b), suggesting that additional aOutlier genes might be discovered as more RNA-seq samples become available. This strong sample size correlation was further confirmed by down-sampling analyses in representative tissues (Supplementary Fig. 3c). Moreover, we noticed that the incidence of an aOutlier identified in one tissue being replicated in another was as low as 14.3% (Supplementary Fig. 4), indicating a significant degree of tissue-specificity among these single-tissue aOutliers. We further defined multi-tissue aOutliers based on aberrant APA usage across five or more tissues (see Materials and Methods). From this analysis, we identified a total of 2,147 multi-tissue aOutliers, comprising 1,930 3′ UTR aOutliers and 217 intronic aOutliers based on the genomic location of the APA event. Focusing specifically on the 715 European individuals, in whom we detected 1,534 multi-tissue aOutliers, including 1,334 3′ UTR and 200 intronic aOutliers (Fig. 1 b and Supplementary Figs. 5 and 6). In our further investigation into the distribution of multi-tissue aOutliers across different tissues, we found that intronic aOutliers exhibited a broader replication pattern than 3′ UTR aOutliers (one-sided Wilcoxon rank–sum test P = 3.35 × 10 ‒14 ; Fig. 1 c, d). Notably, among these aOutliers, several significant genes were identified (Fig. 1 e–g and Supplementary Fig. 7a–f), including SUGP1 , known for its crucial role in mRNA splicing regulation in cancer 27,28 . In certain outlier individual(s), EIF2A , FLYWCH , TP53RK , and SUGP1 exhibited increased usage of distal poly(A) sites, whereas genes such as UNC5A, RAB31 , and LSS preferentially use proximal poly(A) sites. Additionally, genes like COL4A2 (Fig. 1 g), ADCY4 , and HMGCL (Supplementary Figs. 7g, h) were found to utilize intronic poly(A) sites in outlier individuals. Altogether, the single and multi-tissue aOutliers we identified represent the first comprehensive atlas of aberrant APA events across 49 human tissues. aOutliers represent a unique gene set with characteristics distinct from other molecular outliers To determine the extent of sharing between aOutliers genes and those identified as expression outlier or splicing outlier genes (i.e., eOutliers and sOutliers, respectively), we conducted a comparative analysis using the same datasets. Remarkably, we found that 74.2% of multi-tissue aOutlier genes were not detected by analysis of multi-tissue eOutliers or sOutliers (Fig. 2 a and Supplementary Fig. 8a). For example, TRIT1 , a human tRNA isopentenyl transferase 1 gene, is an aOutlier-only gene that preferentially utilizes a distal poly(A) site in outlier individuals across multiple tissues (median Z-score > 11) (Fig. 2 b). This finding suggests that multi-tissue aOutliers represent a novel set of aberrant genes not detectable by traditional eOutlier and sOutlier analyses. Further comparisons between the genomic lengths of multi-tissue aOutliers and eOutliers disclosed that aOutlier genes have significantly longer 3′ UTRs than eOutlier genes (one-sided Wilcoxon rank–sum test, P = 1.4 × 10 ‒16 ) (Fig. 2 c and Supplementary Fig. 8b). In contrast, aOutlier genes have only slightly longer 5′ UTRs than eOutliers (one-sided Wilcoxon rank–sum test, P = 0.004; Supplementary Fig. 8c), and no significant difference was observed in coding sequence length (two-sided Wilcoxon rank–sum test, P = 0.19). Furthermore, aOutlier genes have a lower GC-content (Fig. 2 d) in their 3′ UTR regions (one-sided Wilcoxon rank–sum test, P = 6.8 × 10 –6 ) than eOutlier genes. Gene ontology enrichment analysis 29 on multi-tissue aOutliers further highlighted specific biological processes and signaling pathways unique to these genes (Supplementary Fig. 9). Collectively, these data indicate that aOutliers comprise a distinct gene set with unique molecular and functional characteristics, thereby significantly distinguishing them from other types of molecular outliers. RVs are significantly enriched among APA outliers To assess the impact of RVs (MAF < 0.01) on aberrant APA usage, we computed odds ratios (ORs) for RVs located within varying proximity of the gene body (window size: 1 kb, 2 kb, or 10 kb) to multi-tissue aOutlier genes in outlier individuals compared to those in nonoutlier individuals. Our analysis revealed strong enrichment of nearby RVs in multi-tissue aOutliers (Supplementary Fig. 10a). Interestingly, we observed higher ORs for the enrichment of insertion and deletions (indels) than for single-nucleotide variants (SNVs) (Supplementary Fig. 10a, b). Furthermore, the degree of enrichment became more pronounced when we considered RVs located in closer proximity to the aOutlier genes or employed increased Z-score thresholds (Supplementary Fig. 10b, c). To gain further functional insights into aOutliers-associated RVs, we first determined the proportions of these RVs with functional category using Variant Effect Predictor (VEP) 30 . A higher proportion of aOutliers-associated RVs had function annotation than nonoutliers, increasing with higher Z-score thresholds (Fig. 2 e). The functional categories of aOutliers-associated RVs were largely distinct from those associated with eOutliers and sOutliers. For example, aOutliers-associated RVs are strongly enriched in the 3′ UTR region (OR = 4.6 and 10.1, respectively; Fig. 2 f and Supplementary Fig. 10d). To examine whether aOutliers-associated RVs are more likely to be deleterious and potentially pathogenic, we further employed Combined Annotation-Dependent Depletion (CADD) scores 31 to stratify RVs into three groups: ( 1 ) lowly deleterious, CADD score 0–15; ( 2 ) moderately deleterious, CADD score ≥ 15 but < 25; and ( 3 ) highly deleterious, CADD score ≥ 25. Highly deleterious RVs showed significantly higher enrichment (20-fold increase for singletons and 11-fold increase for RVs with MAF < 1%; Fig. 2 g) in aOutliers compared to moderately deleterious RVs (10-fold increase for singletons and 6-fold increase for RVs with MAF < 1%) and lowly deleterious RVs (2-fold increase for singletons and RVs with MAF < 1%). In total, we identified 179 rare SNVs with CADD scores ≥ 15 near 155 aOutlier genes (two-sided Fisher’s exact test, P = 5.2 × 10 ‒107 ; Supplementary Table 1). In two examples, the rare SNV rs557639120 in SUGP1 (CADD score = 18.4, MAF in GTEx = 0.0056, and gnomAD = 0.0033) leads to an increase in distal poly(A) site usage in its 3′ UTR. Similarly, the rare SNV rs759305120 in COL4A2 (CADD score = 34, MAF in GTEx = 0.0007 and gnomAD = 0.000031) leads to preferential use of its intronic poly(A) site (Supplementary Table 1). We also identified 211 indels near 186 aOutlier genes (two-sided Fisher’s exact test, P = 1.9 × 10 ‒16 ; Supplementary Table 2), including 49 located in 3′ UTR. For example, an indel variant (C > CAAAT, rs112906978) at the 3′ UTR of ACSF3 introduces a canonical "AAUAAA" motif near a poly(A) site, leading to three aOutliers (Supplementary Fig. 10e, f). Enrichment of RVs was also observed in single-tissue aOutliers across nearly all individual tissues (including SNVs and Indels) (Fig. 2 h and Supplementary Fig. 11). Considered collectively, our analyses reveal that a distinct class of RVs is significantly associated with aOutlier genes. Rare APA variants frequently alter the 3′ UTR PAS, 5′ splice sites, and RNA binding proteins (RBPs) binding sites We next investigated the potential regulatory mechanisms of aOutliers-associated RVs on aberrant APA usage. We first focused on 3′ UTR aOutliers-associated RVs and performed motif enrichment analysis to determine the prevalence of RVs altering 3′end processing. Our results show that 3′ UTR aOutliers-associated RVs frequently alter polyadenylation signals (PAS) and AU-rich motifs, such as "AWUAAA" and "AAUAAA" (Fig. 3 a). Additionally, by using saturation mutagenesis data 32 , we found that RVs associated with aOutliers have a more significant impact on poly(A) site usage than RVs associated with nonoutliers (one-sided Wilcoxon rank–sum test P = 1.32 × 10 ‒23 ; Supplementary Fig. 12a). Notably, we observed a significant proportion of large-effect RVs (fold change, LFC > 1) associated with aOutliers compared to nonoutliers (50.3% vs. 6.6%; one-sided Wilcoxon rank–sum test P = 6.1 × 10 ‒44 ; Supplementary Fig. 12b), indicating their pronounced effects on 3′ UTR APA. To further experimentally validate these findings, we selected four top-ranked 3′ UTR aOutlier genes by median Z-score and utilized a minigene reporter system containing reference allele and alternative allele of four rare variants in selected genes, including MKKS (Fig. 3 b), SUGP1 , TP53RK , and ATP5F1E . In all four cases, we could detect significant changes in the poly(A) site usage, which agreed well with the predicted effects of these RVs (Figs. 3 c, d and Supplementary Fig. 13a, b). Further investigation into multi-tissue intronic aOutliers revealed a higher incidence of RVs at 5′ splice donor sites than at acceptor sites (Fig. 3 e). Compared to nonoutlier RVs, aOutlier RVs are 19 to 441 times more prevalent at donor sites, and up to 47 times more prevalent at acceptor sites. Specifically, aOutlier RVs are 441 times more prevalent in the "D + 1" site and “D + 4” site and 302 times more prevalent in the "D + 2" site relative to the nonoutlier RVs. For example, RVs that alter the first nucleotide of the "GT" sequence in the intron of COL4A2 (Fig. 1 h) and the intron of TXNRD2 lead to significant intronic APA events in these genes (Fig. 3 f and Supplementary Fig. 13c). We also found that RVs altering the last base of exon 11 in ADCY4 and exon 4 in HMGCL resulted in intronic APA events (Fig. 3 g and Supplementary Fig. 7e, f). Based on these findings, we hypothesized that RVs affecting canonical donor sites drive intronic aOutliers. This hypothesis is also supported by our recent finding that mutations near the donor sites can promote IPA usage, potentially by blocking U1 small-nuclear RNP binding 33 . Predicting the strength of donor sites with MAXENT 34 showed a reduced strength of mutant donor sites compared to wild type (Fig. 3 h, i). We then performed intronic APA minigene reporter assays for TXNRD2 and COL4A2 with RVs at the conserved donor sites, as well as HMGCL and ADCY4 with RVs at the last base of the exons. For these assays, we cloned fragments containing full-length intronic sequences, including the donor sites, and upstream and downstream exons into the pcDNA3.1 vector. Results from 3′ Rapid Amplification of cDNA Ends (3′ RACE) assays indicate that all four RVs significantly increase alter IPA regulation relative to the wild-type sequence (Fig. 3 j, k and Supplementary Fig. 13d, e). Lastly, we investigated whether aOutlier-associated RVs impact other transcriptional and posttranscriptional regulation of target genes. DeepBind 35 analysis of 927 binding motifs revealed 11 significantly enriched motifs in aOutlier-associated RVs (Supplementary Fig. 14a) using randomly shuffled RVs as control, including known APA regulator PABPN1 36 . Furthermore, we analyzed 166 publicly accessible RBPs cross-linking immunoprecipitation sequencing (CLIP-seq) datasets from the Encyclopedia of DNA Elements (ENCODE) project 37 . We found seven RBPs's CLIP-seq data are strongly enriched with multi-tissue aOutlier RVs compared to nonoutlier RVs (Fig. 3 l and Supplementary Fig. 14b), including LARP4 , an APA regulator identified in our previous study 25 , and a known APA regulator CSTF2T . Knockdown of the two RBPs resulted in widespread APA dysregulation (Supplementary Fig. 14c, d), affecting two aOutlier genes, SREBF2 (Supplementary Fig. 14e) and TOLLIP (Fig. 3 m), in which the associated RVs were inside binding peaks of LARP4 (Supplementary Fig. 14f) and CSTF2T (Fig. 3 n), respectively. Beyond these known APA regulators, other RBPs such as TIA1 , UPF1 , and SAFB2 were also identified as potential new APA regulators (Supplementary Figs. 14g-i). Collectively, these results suggested that aOutlier-associated RVs trigger aberrant APA usage through altering PAS, splice sites, or RBP binding sites. Inclusion of APA significantly improves functional RV effect prediction To prioritize potentially impactful RVs for the interpretation of individual genomes, we repurposed the traditional Watershed 7 method into an APA-included version (aWatershed). This revised aWatershed model is an unsupervised probabilistic Bayesian hierarchical graphical model incorporating three RNA outlier signals, including aOutliers, eOutliers, and sOutliers, and annotations of a matched individual genome (Supplementary Table 3). The aWatershed model can allow us to quantify the posterior probability of an RV leading to a functional effect on APA usage (Supplementary Figs. 15a, b; Materials and Methods). To evaluate the aWatershed performance on the GTEx v8 data, we used held-out individual pairs with the same RVs as the evaluation dataset. By applying aWatershed prediction on the first individual of each pair and evaluating this prediction using the outlier status of the second individual as a label, we observed that our model significantly outperforms both the RIVER (RNA-informed variant effect on regulation) model 8 , a simplification of the Watershed model which integrates genomic features with aOutlier signals alone, and the GAM (genomic annotation model), a generalized logistic regression model based on genomic features alone (Fig. 4 a and Supplementary Fig. 15c). 93% of aWatershed prioritized RVs have low posterior probabilities in the GAM (Fig. 4 b), highlighting the importance of transcriptomic aOutlier signals in functional RVs prioritization. Moreover, aWatershed successfully captures the regulatory mechanisms underlying the effect of RVs on aOutlier signal (Fig. 4 c). Strikingly, the integrated aWatershed model can prioritize RVs associated with 73.8% of aOutliers, in contrast to only 12.4% when relying on the genomic features alone (Fig. 4 d). Next, we used the saturation mutagenesis data 32 to further evaluate the efficacy of aWatershed in prioritizing RVs with significant effects on APA regulation. In this analysis, we stratified RVs into two groups based on aWatershed APA posterior probabilities and compared poly(A) usage between them. We found that RVs in the group with high posterior probability had significantly larger effects on APA than those in the low posterior probabilities group (Fig. 4 e, f), suggesting our aWatershed model is effective in identifying RVs with substantial APA effects. Furthermore, our analysis revealed that aWatershed successfully identified many functional RVs overlooked by the previous variant prediction model 38 , as exemplified by two RVs in RPL13A and PAAF1 , respectively (Supplementary Fig. 15d). Overall, aWatershed prioritized 1,799 RVs predicted to impact 278 APA genes (Supplementary Table 4). Interestingly, there was minimal overlap between RVs impacting APA and those affecting gene expression or splicing, as only 60 of these 1,799 RVs were common to those categories. For example, the RV rs191575428 within the 3′ UTR of MTHFD2 , which exhibited a high aWatershed APA posterior probability of 0.997 based on aOutliers, showed considerably lower posterior probabilities for expression and splicing (0.055 and 0.008, respectively). This variant is associated with 3′ UTR lengthening in outlier individuals without changing gene expression levels (Supplementary Fig. 15e, f). Further extending the aWatershed model to prioritize tissue-specific functional RVs by integrating genomic features with single-tissue aOutliers signals, we observed that the tissue-aWatershed model outperforms both the tissue-RIVER model and tissue-GAM model (Supplementary Figs. 16–17). In summary, by leveraging these aOutliers, we have implemented a robust Bayesian hierarchical variant effect prediction model aWatershed that effectively prioritizes rare functional variants with significant effects on APA regulation. Analysis of aOutliers prioritizes RVs impacting complex traits and diseases To test the hypothesize that aWatershed RVs could be used to interpret the complex traits and diseases, we first examined the 278 genes prioritized by aWatershed and cross-referenced with genes annotated in the Online Mendelian Inheritance in Man (OMIM) database 39 . We identified 21.2% of the prioritized genes were well-known disease genes (Supplementary Fig. 18a). For example, we identified a prioritized RV, rs79940214, associated with MKKS (Supplementary Fig. 18b), which encoded a centrosome-shuttling protein and was associated with many genetic diseases, including McKusick-Kaufman syndrome (OMIM id: 236770) 40,41 and Bardet-Biedl syndrome 6 (OMIM id: 605231) 42,43 . Another example is one prioritized intronic RV, rs76984877, that is associated with gene EXT2 (Supplementary Fig. 18b), which was associated with hereditary multiple exostosis, type 2 44,45 . We also identified five prioritized RVs associated with gene BCR (Supplementary Fig. 18b), which has been frequently reported to be associated with chronic myeloid leukemia 46,47 . We further cross-referenced aWatershed-prioritized RVs with trait variants from 1,234 well-powered GWAS summary statistics from UK Biobank (UKBB) and literature (Supplementary Table 5), resulting in 1,385 RVs associated with 171 aOutlier genes in 1,186 traits (Supplementary Table 6). We focused on the subset of 623 traits, which also have evidence of colocalization with 3′aQTLs (Supplementary Table 7). Notably, aOutlier prioritized RVs fell in or nearby genes had evidence of colocalization with 3′aQTLs having larger trait effect size than the non-colocalized RVs ( P = 0.0014, one-sided Wilcoxon rank–sum test; Fig. 5 a). We also conducted a permutation test to determine whether these prioritized RVs exhibit larger effect sizes on these complex traits and diseases. We found that the mean odds of aOutlier-prioritized RVs had a more significant effect size than non-prioritized RVs ( P = 2.5 × 10 ‒15 , one-sided and paired Wilcoxon rank–sum test; Supplementary Fig. 19a). To exemplify the larger effect size in aOutlier-prioritized RVs, we focused on two traits: height related traits (UKBB trait ID: 50_irnt and 20015_irnt) and high blood pressure (UKBB trait ID: 6150_4). This analysis revealed a significant shift in the odds favoring RVs with higher aWatershed posterior probabilities over those with lower ones ( P = 1.6 × 10 ‒9 and P = 2.3 × 10 ‒54 , respectively; one-sided Wilcoxon rank–sum test; Fig. 5 b, c; Supplementary Fig. 19b). In the case of height related traits and high blood pressure, these aOutlier prioritized RVs had larger effect sizes on the trait than other variants within a 1Mb of the RV, including RVs prioritized by eOutliers or sOutliers. Notably, for the height related traits, the RV (rs112567314), located in the intron of CUL3 , had a greater effect size than other variants within 1 Mb and RVs prioritized by eOutlier or sOutlier (Fig. 5 d and Supplementary Table 6). Similarly, for high blood pressure, the RV (rs893929), located in the intron of USP38 , also had a greater effect size than 99.6% of variants within 1 Mb, including the nearest trait-associated significant variants as well as eOutlier or sOutlier RVs (Fig. 5 e). Collectively, our results demonstrate the capability of aWatershed in prioritizing RVs with large effect sizes on APA, significantly impacting complex traits and diseases. Strong convergence between rare and common variants on DDX18 links APA regulation with cancer susceptibility Emerging evidence suggests potential interactions between rare and common variants in affecting the same disease genes 48–50 . As expected, we also observed the strong convergence effect on 3' UTR APA regulation between RVs and common variants (Supplementary Fig. 20). To further mechanistically examine their convergence effects on disease, we focused on aWatershed prioritized RVs and their associated genes. We found 126 out of the 278 aOutlier RV associated genes were also identified as susceptibility to disease risks, including cancer risks through 3′aQTLs in our gene-based association studies 51,52 (Fig. 6 a). Among the top-ranked APA genes that were prioritized by both RV and 3′aQTLs analyses (Fig. 6 b), we noticed several highly constrained genes (pLI score > 0.9), and we particularly focused on the gene DDX18 , a member of the DEAD-box RNA helicase family, that was identified as an APA-mediated susceptibility gene across many cancer types 53,54 . Moreover, CRISPR-Cas9 based gene essentiality screens also demonstrated that DDX18 has an essential role in cancer cell proliferation 55,56 (Fig. 6 c). Examining our 3′aQTLs data revealed significant associations between common variants and 3′ UTR APA of DDX18 across tissues, with the most significant one was found near the 3′ end (Fig. 6 d and Supplementary Fig. 21a-c). Intriguingly, an aWatershed prioritized RV, rs1680042046, located near the distal poly(A) site of DDX18 , was identified in the outlier individual (Fig. 6 e, f). This RV alters the hexamer motif "AUUAAA" to "AUUAAG" (Supplementary Fig. 21b) and has a highly deleterious effect (CADD = 17.5) (Supplementary Table 1). To further experimentally validate the convergence effect of RVs and common variants on DDX18 , we designed minigenes introducing the APA variants by PCR-based site-directed mutagenesis in HEK293T and MCF7 cells (Fig. 6 g). We then performed 3′ RACE to quantitatively evaluate the effect of the common variant (rs1052628; A > G) alone, the RV (rs1680042046) alone, or their joint effect on APA. We first mutate the reference A allele to the alternative G allele for either RV or the common variants. In HEK293T cells, this mutation decreased the use of the distal poly(A) site (dPAS) for both the common variant and RV (two-sided Student’s t-test P = 1.5 × 10 ‒6 and 2.3 × 10 ‒7 ; Figs. 6 h, i), indicating that both variants indeed trigger DDX18 APA regulation. A similar APA effect was also observed in MCF7 cells (Fig. 6 j). To further assess the functional roles of DDX18 APA regulation in breast cancer cells, we measured DDX18 protein level using luciferase reporter assays and assessed the effect of gene silencing on the proliferation of MCF7 cells proliferation. We observed lower luciferase activities in the short 3′ UTR isoform of DDX18 and the reporter containing RV or both RV and common variant (Supplementary Fig. 21d, e). Knockdown of DDX18 in MCF7 results in inhibition of cell proliferation (Supplementary Fig. 21f, g). Collectively, these findings highlight the critical role of rare variants in understanding the risk of common diseases and offer a novel approach to linking functional rare variants to complex diseases. Discussion The human genome contains a plethora of rare genetic variants whose functional effects and underlying molecular mechanisms are challenging to interpret. In this study, we introduce the aOutlier as an emerging molecular phenotype reflecting aberrant 3′ UTR or intronic APA usage across multiple samples. aOutlier can be used to identify functional rare APA variants. By analyzing population-scale transcriptomics data using our DaPars2 25,26 and IPAfinder algorithms 18 , we identified 1,534 multi-tissue aOutliers based on European individuals. These aOutlier genes exhibit unique molecular features, such as genomic lengths and GC-content, setting them apart from other molecular outliers, such as eOutliers and sOutliers. Importantly, aOutliers can aid in identifying a distinct class of rare functional variants. We observed that aOutliers-associated RVs are more likely to be deleterious and are highly enriched in outlier individuals. Mechanistically, these aOutlier-associated RVs can modulate APA usage by either altering PAS, AU-rich elements, or splice donor sites, as confirmed by saturation mutagenesis data and 3′ RACE experiments. To further enhance the utility of our aOutlier atlas, we adapted a Bayesian hierarchical prediction model (aWatershed) by incorporating genomic features with multiple functional signals, including aOutliers, eOutliers, and sOutliers. This integration aims to predict the probability of RV leading to aberrant APA usage. Notably, our aWatershed model outperformed models trained only on genomic features or those combined with aOutlier signals alone. Moreover, aWatershed-prioritized RVs exhibited more significant effects on APA regulation than non-prioritized RVs. The predictive power of aWatershed was validated using GWAS summary data from the UKBB, showing that aWatershed-prioritized RVs had larger trait effect sizes than non-prioritized RVs, as exemplified by RVs near POLR2L and ATP5F1D associated with height and BMI, respectively. Interestingly, we observed a significant proportion of intersection between aOutlier transcripts and 3′aQTL associated transcripts in matched tissue, suggesting the potential interplay of common variants and RV in APA regulation, similar to previous findings in gene expression studies 48,49,57,58 . Additionally, a rare deletion 16p11.2 and common variants in chromosome 16p modulate downstream gene expression and affect the risk for autism 48 . Moreover, using minigene reporters and 3′ RACE assays, we demonstrated the potential additive effect of rare and common APA variants on DDX18 3′ UTR regulation. We further demonstrated that the regulation of DDX18 3′ UTR contributes to DDX18 protein expression level, which is tightly linked to breast cancer cell proliferation. In summary, our study identifies a novel set of rare functional variants that influence APA and connects these RVs to human trait phenotypes, providing valuable information for the identification of novel genes associated with increased disease risk. Materials and Methods GTEx data collection and processing We downloaded both RNA-seq raw sequencing data and whole-genome genotype data of the v8 release of the GTEx project from dbGAP (accession: phs000424.v7.p2). Expression outlier (eOutlier) and splicing outlier (sOutlier) data, and the metadata of samples (filename: GTEx_Analysis_v8_Annotations_SampleAttributesDD.xlsx) and subjects (filename: GTEx_Analysis_v8_Annotations_SubjectPhenotypesDD.xlsx) were downloaded from GTEx Portal ( https://gtexportal.org/home/ ). The GTEx RNA-seq dataset contains 17,832 samples representing 54 biological tissues collected from 838 donors. For this study, we included data from 49 tissues, each with at least 70 samples. Original GTEx RNA-seq reads were aligned with the human genome (hg38/GRCh38) using STAR v.2.7.3a 59 , with the following alignment parameters: outSAMtype, BAM; SortedByCoordinate; outSAMstrandField, intronMotif; outFilterMultimapNmax, 10; outFilterMultimapScoreRange, 1; alignSJDBoverhangMin, 1; sjdbScore, 2; alignIntronMin, 20; and alignSJoverhangMin, 8. The aligned BAM files were sorted and further converted to bedGraph format using BEDTools v.2.27.1 60 . The genotype data in VCF format (filename: GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.vcf.gz) was processed with vcftools v.0.1.13 to calculate MAF across all subjects and extract allele information for each variant. 3′ UTR APA and intronic APA quantification To quantify the 3′ UTR APA, we analyzed alignment files in BAM format using the DaPars2 algorithm. We followed the workflow implemented in our 3′aQTL analysis 25,61 . Briefly, the BAM files were firstly transformed to bedGraph format with a bin size of 1, which records the read coverage of each position in the genome. Before analyzing APA, we downloaded the gene annotation file containing all transcripts of genome build hg38 in RefSeq database through the UCSC Genome Browser, from which we extracted 3′ UTR region of each transcript using script "DaPars_Extract_Anno.py". The DaPars2 algorithm then detects the proximal poly(A) site in the 3′ UTR region of each transcript and calculates the relative usage of the distal poly(A) site by the script “Dapars2_Multi_Sample.py” for all samples in each of the 49 tissues. This is indicated as the Percent of Distal Poly (A) site Usage Index (PDUI) only if a proximal poly(A) site is detected. For intronic APA detection and quantification, we used IPAfinder 18 , which is a python-based tool that enables de novo identification and quantification of intronic APA (IPA) events using RNA-seq data. IPAfinder can identify potential IPA sites and calculate the Intronic poly(A) site Usage Index (IPUI), which represents the proportion of total transcripts that are intronic-polyadenylated for each intronic APA event 18,33,62 . BAM files were analyzed together by IPAfinder and separated by tissues. Covariate correction and normalization To avoid batch effects and unobserved confounders in each tissue, we adjusted the sample genotype and APA usages with known covariates, such as population structure, sex, and sequencing platform. Briefly, for genotype data, we first removed sites marked as "wasSplit" from the GTEx analysis freeze variant call format (VCF) using BCFtools v.1.10.2. We further applied the PEER model 63 with sex, age, sequencing platform, and the top five genotype principal components as known covariates to estimate a set of latent covariates for PDUI or IPUI values in each tissue. The number of PEER factors was optimized based on tissue sample size, as suggested by the GTEx Consortium; 15 PEER factors were chosen for sample sizes 250. Before running the PEER model for inferring hidden covariates, PDUI/IPUI values in each tissue were quantile normalized to remove batch effects. APA outlier calling After inferring the hidden covariates for each tissue, we calculated PDUI/IPUI residuals by regressing out inferred PEER factors and known covariates, including population structure, sex, and sequencing platform, using the function "PEER_getResiduals". In each individual tissue, we obtained normal-distributed \(Z\left(g,t\right)\) score for each gene ( \(g\) ) in the tissue ( \(t\) ) by scaling the PDUI/IPUI residuals across samples with the following equation, \({X}_{r}\left(g,t\right)\) denotes the residuals of PDUI/IPUI values, \(\stackrel{-}{{X}_{r}\left(g,t\right)}\) and \(sd\left({X}_{r}\left(g,t\right)\right)\) represent the mean and standard deviation of the residuals across samples, respectively: $$Z\left(g,t\right)=\frac{{X}_{r}\left(g,t\right)-\stackrel{-}{{X}_{r}\left(g,t\right)}}{sd\left({X}_{r}\left(g,t\right)\right)}$$ We defined two types of aOutliers in the current study. One is single-tissue aOutlier, which is called from a single tissue based on the Z-score of each gene in that tissue. When the absolute Z-score of an individual exceeds a threshold of three for a gene, then the individual is called a single-tissue aOutlier for that gene. The other is multi-tissue aOutlier, for which we calculated a median Z-score 7,8 for each APA event across all tissues when data were available and restricted our analysis to individuals with APA measurements in at least five tissues. Multi-tissue aOutliers were defined as those with an absolute median Z-score > 3. The same threshold was used for eOutlier and sOutlier calling 7,8 . Our method allowed that one gene could have multiple aOutlier individuals, and one individual could also be aOutliers of multiple genes. To account for situations in which widespread aberrant APA might occur in an individual due to non-genetic influences, we removed 11 individuals in which the proportion of tested genes identified as multi-tissue aOutliers exceeded 1.5 times the interquartile range of the distribution for aOutlier gene proportion across all individuals. These 11 individuals were marked as global outliers. Estimation of replication rates of aOutliers To estimate the replication rate of aOutliers between different tissues, we selected one of the 49 GTEx tissues each as discovery tissue, and compared aOutliers detected in it with those of the other tissues. For replication rate calculation, we only considered the shared aOutlier genes in both tissues and an aOutlier to be replicated only when the gene and individual of the aOutlier matched between the compared tissue pairs. For multi-tissue aOutliers replication, we used the cross-validation method described in a previous study 8 to estimate their replication rate. In brief, the 49 human tissues were separated into two groups, one group with 39 tissues as the discovery group, the other group has the remaining ten tissues as the replication group. Each time we randomly sampled t ( t = 10, 15, 20, 25, 30) tissues from the discovery group and called multi-tissue aOutliers in the discovery group using a Z-score threshold of 3 in at least five tissues as described above. Then we estimated the replication rate as the proportion of multi-tissue aOutliers in the discovery group with an absolute median Z-score 2 or 3 in the replication group. We also computed the expected replication rate by randomly selecting individuals in the discovery group with at least five tissues that have APA usage for the gene and determined the replication status in the replication group. For each discovery group size ( t ), we repeated this process 10 times. RV annotation We defined RVs as those with MAF < 1% within the GTEx European individuals and with MAF < 1% in non-Finnish Europeans within gnomAD 64 . Singletons were defined as RVs with minor allele only presents once in GTEx European individuals and were extracted using vcftools. The annotation of RVs was performed by Ensembl VEP (release 104), which assigned 36 different annotation terms to each RV, including protein-coding gene position (e.g., "splice_donor" "splice_acceptor," "frameshift”) and regulatory regions (e.g., "TFBS_ablation", "TF_binding_site"). Annotation terms were grouped into one of the four classes based on predicted impact: "High", "MODERATE", "MODIFIER" and "LOW". The high-impact one was used for variants assigned with two or more annotations. In addition to 36 VEP annotations, we added two other annotations to each RV; "PAS region" describes variants located within 50 bp upstream of the annotated PAS, and "PAS signal" refers to variants located at the PAS motif "AAUAAA" and its additional 14 variants ("AUUAAA", "UAUAAA", "AGUAAA", "AAAAAA", "AACAAA", "AAGAAA", "AAUAUA", "AAUACA", "CAUAAA", "UUUAAA", "ACUAAA", "AAUAGA" and "GAUAAA"). We also used genomic annotations of variants extracted from CADD v.1.5 release ( http://cadd.gs.washington.edu/download ). RV enrichment analysis We examine the enrichment of Rare Variants (RVs), including single-nucleotide variants (SNV) and small insertion and deletion (indel) near aOutlier genes. Only genes with at least one aOutlier individual were considered, and the remaining individuals for the same genes were treated as nonoutlier controls. We first counted the RVs present within 1kb, 2kb, or 10kb of the outlier genes in both outlier and control individuals and built a 2 \(\times\) 2 contingency table for each of the flanking region, containing the number of aOutliers with RVs, the number of nonoutlier controls with RVs, the number of aOutliers without RVs, and the number of controls without RVs. We then calculated Odds Ratios (ORs), P value, and 95% confidence interval (CI) using Fisher’s exact test in R base package. We grouped variants into four groups based on their MAF (0–1%, 1–5%, 5–10%, and 10–25%), and performed enrichment analysis for each group. We also conducted enrichment for RVs that stratified by VEP annotations and CADD scores as described above. Enrichment analysis for RVs that influence PAS and AU-rich motifs To identify potential regulatory variants associated with aberrant APA events, we defined RVs located within the gene body or in the 10-kb region surrounding outlier genes in outlier individuals as aOutlier RVs. Those in nonoutlier individuals in the same region were classified as nonoutlier RVs. For each aOutlier RV and nonoutlier RV located in the 50-bp region upstream (PAS region) of the poly(A) sites annotated in PolyA_DB V.3.2 65,66 , we extracted its upstream and downstream 5 base pairs sequences and examined whether it matched with one of the 15 known PAS motifs by using script "dna-pattern" in RSAT ( https://github.com/rsa-tools ). We then summarized all tested RVs and conducted PAS motif enrichment analysis using Fisher's exact test, which determined the odds ratios (ORs) and 95% confidence intervals (CIs) for each PAS motif. To perform enrichment analysis at the 12 known AU-rich motifs, we restricted RVs to those within the 100 bp flanking the annotated poly(A) sites. We then counted RV enrichment analysis for each of the AU-rich motifs using the same method as for PAS motifs. Identification of aOutlier RVs enriched RNA motifs We focused on multi-tissue aOutlier associated RVs located within the gene body region, which spans from 3 kb downstream of the transcription start site (TSS) to the end of the gene. We extracted the 3 base pairs of sequences flanking each RV from both sides. Next, we used DeepBind v0.11 35 to score these 7-mer sequences using 617 pre-built models, including 515 transcription factors and 102 RNA-binding proteins from Homo sapiens. For each 7-mer sequence, we selected the top three motifs with a DeepBind score of at least 0.1. To validate the enrichment of RVs in predicted binding motifs, we created a control set of RVs by randomly shuffling the genomic locations of multi-tissue aOutlier associated RVs within the same gene body regions. We used Fisher's exact test to estimate the level of enrichment. Identification of aOutlier RVs enriched RBPs We obtained CLIP-seq data for 166 RNA-binding proteins (RBPs) from the Encyclopedia of DNA Elements (ENCODE) data portal for HepG2 and K562 cells. We only considered significant binding peaks with P-values < 0.01, shared by two biological replicates for each RBP. To assess the enrichment of aOutlier RVs in RBP binding peaks, we selected RVs associated with multi-tissue aOutliers within gene body regions representing the region of 3 kb downstream of the transcription start site (TSS) to the end of the gene. We created a control RV set by randomly shuffling the genomic locations of multi-tissue aOutlier associated RV set within the same gene body regions. We then counted the RVs in binding peaks of each RBP using bedtools. Finally, we compared the RVs between the two sets using Fisher's exact test to determine the enrichment. Development of a Bayesian prediction model that integrates APA outlier signals To prioritize rare functional variants with significant impact, we improved the Watershed Bayesian hierarchical model by incorporating APA outlier signals with other layers of transcriptomic outlier signals and genomic annotations. The improved model called aWatershed, includes a layer of genomic annotation features ( G ) which denotes the 40 observed genomic features aggregated over all RVs in the outlier individual that are within 10 kb region of the gene, a fully connected conditional random field (CRF) layer ( Z ) represent the unobserved regulatory variables for each of the three transcriptomic outlier signals (APA, mRNA expression, and splicing), and a layer of variables ( E ) representing the observed outlier status of each transcriptomic data type. The three layers were linked by the following conditional distributions: Where \(K\) represents the three outlier signals (APA, Expression, and Splice), \({\beta }_{k}\) are parameters defining the contribution of the 40 genomic features to the CRF of the three outlier signals, \(\alpha\) defines the intercept of the CRF for each outlier signal, \(\theta\) represent parameters defining the edge weights between pairs of the three outlier signals, \({\phi }_{k}\) are the paramters denoting the categorical distributions of each of the three outlier signal, and \(C\) and \(\lambda\) are hyper-parameters. To train and evaluate aWatershed, we utilized all gene-individual pairs that have at least one of the three multi-tissue outlier signals, which are defined as the absolute value of Z-score greater than 3 or P-value less than 0.0027 for splicing outliers, measured in GTEx v8 data. We also used a set of 38 binary and continuous genomic annotation features aggregated across all rare variants within the 10-kb region, flanking each gene. We then trained aWatershed to learn edge weights connecting the three transcriptomic outlier signals, weights representing the contribution of each genomic annotation for each type of outlier signal, and other parameters, as described previously 7,8 . To evaluate aWatershed, we selected pairs of individuals with the same set of rare variants associated with the same gene (known as "N2pair") from the training dataset. We estimated the posterior probability of a functional rare variant in the first individual of the pair and used the outlier status of the second individual as a label for evaluation. We also trained and evaluated the genomic annotation model on each layer of the three transcriptomic signals to determine whether the integration of transcriptomic outlier signals contributes to the prediction of rare functional variants. We compared the results to those obtained from the aWatershed model. After evaluation, we utilized the aWatershed prediction model to calculate posterior probabilities. 3′aQTL mapping across 49 GTEx tissues Genetic associations between GTEx common variants within 1 Mb of each gene and PEER-corrected APA usage were mapped by Matrix eQTL 67 , as described in our previous study 25 . Known covariates, including sex, RNA integrity number, platform, top five genotype principal components, and unobserved covariates inferred from PEER, were used during 3′aQTL mapping with Matrix eQTL. The number of PEER covariates for each tissue was used as suggested by the GTEx Consortium. We performed 1,000 rounds of permutation to obtain empirical P -values for each gene, which were then adjusted using the R package qvalue. Colocalization analysis between GWAS summary statistics and 3′aQTL We conducted colocalization analysis comparing GWAS summary statistics from the UKBB and literature and 3′aQTL summary data from 49 human tissues using the coloc v.5.1.0.1 package in R 68 . Only GWAS summary data for traits with at least 10,000 cases (binary traits) or 50,000 participates (continuous trait) and with at least 10 SNPs overlapped with aOutlier-associated RVs were kept, which resulted in 1,186 well-powered traits. We extracted the sentinel SNPs for each GWAS trait, defined as GWAS SNPs with P < 5 × 10 ‒8 , located at least 1 Mb away from more significant variants. We then searched for colocalizing signals within the 1-Mb region surrounding each sentinel SNP. The coordinates from 3′aQTL summary data were converted from human genome build 38 (hg38) to build 37 (hg19) by CrossMap software 69 to match the version used in all GWAS summary statistics. As defined by the coloc method, five posterior probabilities under five different null hypotheses were calculated. In detail, PP0 represents the null model of no association. PP1 and PP2 represent the probability that causal genetic variants are associated with disease signals or 3′aQTL only. PP3 represents the probability that the genetic effects of trait signals and 3′aQTL are independent, and PP4 represents the probability that trait signals and 3′aQTL data share causal SNPs. The current study classified colocalized events as those with PP4 > 0.75. 3 ′ UTR APA transcriptome-wide association study (3 ′ aTWAS) analysis We used APA quantitative data that was previously used for 3′aQTL mapping 25,52,61 and genotype data of individual genomes from whole genome sequencing (WGS) of GTEx consortia to construct 3′aTWAS model using FUSION 70 for each of the 49 human tissues. To avoid the effects of confounders, well-established factors used in 3′aQTL mapping, including gender, sequencing platform, and other covariates, were incorporated to adjust APA usages. To build the TWAS model, four different models embedded in FUSION were used for weight calculation, including best linear unbiased predictor (blup), elastic-net regression (enet), lasso regression (lasso), and single best eQTL (top1). Subsequently, the cross-validation approach was employed to choose the optimal 3′aTWAS model for each gene. Of note, only genes exhibiting significant heritability estimates ( cis-h 2 ) (Bonferroni-corrected P < 0.05) were retained for subsequent analysis. The built models were then applied to GWAS summary statistics for gene-based association analysis, and a significant association was defined by the FDR < 0.05. The disease risk genes identified by 3′aTWAS in two or more tissues were used for further analysis. Prioritization of trait-associated RVs To determine the frequency with which randomly selected aWatershed-prioritized RVs exhibit larger GWAS effect sizes than matched non-prioritized RVs, we conducted a random sampling test (n = 1000) on all RVs using posterior probabilities obtained from the aWatershed prediction model and effect sizes from UKBB GWAS summary statistics. We used aWatershed-prioritized RVs based on aOutlier signals, matched non-prioritized RVs, as well as GWAS effect sizes, gene IDs, and prioritized scores as input data. We defined matched non-prioritized RVs as those with a posterior probability of < 0.1 and MAF within ± 0.001 of the selected prioritized RVs in the UKBB cohort. For each gene in each trait, we randomly selected one prioritized RV and one matched non-prioritized RV and then identified the one with the largest absolute GWAS effect size in the pair. By summarizing all genes in the trait, we computed the odds of observing a prioritized RV with a larger absolute effect size than a non-prioritized RV across all genes. To generate a null distribution of odds, we repeated this process for matched non-prioritized variants only and randomly selected and compared two non-prioritized RVs for each gene. Cell culture HEK293T and MCF7 cells were purchased from the Cell Bank of the Type Culture Collection at the Shanghai Institute of Biochemistry & Cell Biology, Chinese Academy of Science. Cells were maintained in Dulbecco's modified Eagle medium (DMEM; Invitrogen, #11960044) supplemented with 10% fetal bovine serum (Gibco), 100-µg/ml streptomycin, and 100-units/ml penicillin at 37°C in a humidified incubator with 5% CO 2 . Plasmid construction All primers used in this study are listed in Supplementary Table 8. For intronic APA (IPA) minigenes, the candidate intron and its flanking exons were amplified from genomic DNA as wild-type fragments. For 3′ UTR APA minigenes, the 3′ UTR of each gene was amplified from genomic DNA as wild-type fragments, and mutations were introduced by PCR-based site-directed mutagenesis. In short, genomic DNA from HEK293T and MCF7 cells was amplified by PCR using primers to generate two 20–25 bp overlapping fragments containing a mutant site. The IPA wild-type and mutant fragments were subcloned into the EcoRI and BamHI sites of the pcDNA3.1 vector, while 3′ UTR APA wild-type and mutant fragments were subcloned into the XhoI and PmeI sites of the mpCHECK2 vector by the One Step Cloning Kit (Vazyme). Two sets of predesigned shRNAs from Sigma against DDX18 were used to clone into pLKO.1-puro vector. Transient transfection For transient transfection, HEK293T and MCF7 cells were plated in a 2-ml culture medium at 6 × 10 5 cells/well in six-well plates. After 24 h of culture, cells were transfected with 2 µg of wild-type or mutant minigene plasmid using Lipofectamine 2000 (Invitrogen), according to the manufacturer's instructions. The culture medium was replaced at 6 h post-transfection, and cells were harvested for RNA extraction at 48 h post-transfection. Total RNA was extracted using TRIzol reagent (Invitrogen), according to the manufacturer's instructions, and cDNA was synthesized using the FastKing RT Kit (Tiangen, KR116) with the S-CDS primer. All cDNA was diluted 4-fold in nuclease-free double-distilled H 2 O for further use. 3′ RACE The total length of 3′ UTR was identified and amplified from the total RNA of NCI-H1299 cells by 3′ RACE using the HiScript-TS 5′/3′ RACE Kit (Vazyme, RA101) following the manufacturer’s protocol. 3′ RACE was performed using the S-PCR primer and pcDNA3.1-F or mpCHECK2-F primer to distinguish minigene RNA from endogenous RNA, respectively. The 3′ RACE PCR products were separated by gel electrophoresis, and excised bands were purified for Sanger sequencing using the Zymoclean Gel DNA Extraction kit. Cleaned DNA fragments were cloned into the PCE2 vector using the 5 min TA/Blunt-Zero Cloning Kit (Vazyme, C601) and bidirectionally sequenced with M13 forward and reverse primers. At least five colonies were sequenced for every gel product that was purified. Primer sequences are listed in Supplementary Table 8. Dual-luciferase reporter assay MCF-7 cells were seeded 1 day prior to transfection. The Renilla luciferase in the mpCHECK-2 vector was transfected into cells using Lipofectamine 3000 Transfection Reagent (Invitrogen, cat#: L3000015) according to the manufacturer′s instructions. Forty-eight hours post-transfection, firefly, and renilla luciferase activities were measured by Dual-Luciferase Assay System (Promega, #E1980) on a BioTek Synergy H1 plate reader with full waveband. Each assay was measured in three independent replicates. Cell viability and proliferation assays for shRNA-mediated knockdown shRNA-expressing lentivirus was produced with the third-generation packaging system in human embryonic kidney (HEK) 293T cells. For lentivirus infection, target cells (MCF7) were seeded in a 6-well plate 16–18 h before infection and were grown to 70–80% confluency upon transduction. The culture medium was removed, and cells were incubated with virus supernatant along with 8 µg/ml polybrene. Puromycin was applied to kill non-infected cells 2days after infection. After two days of selection, when non-infected control cells were all dead, surviving cells were split and maintained with the same concentration of puromycin. Cells were trypsinized, resuspended at 1 × 10 4 cells/ml, and seeded in 96-well plates, with each well containing 100ul medium of 1 × 10 3 cells. Cell viability and proliferation were determined using CCK8 assays (Yeasen, cat#: 40203ES76) at designated time points (day 1, day 3, day 5, and day 7) by measuring the absorbance at 450 nm, following the manufacturer′s instructions. Values were obtained from three replicate wells for each treatment and time point. Results were representative of three independent experiments. The comprehensive data portal for aOutliers We have established a database along with a web interface called rareAPA ( http://bioinfo.szbl.ac.cn/rareAPA/index.php ) on a standard LAMP (Linux + Apache + MySQL + PHP) system, which serves as a comprehensive resource presenting detailed and comprehensive information on rare APA events and their associated RVs. All these data in the rareAPA were stored in MySQL ( www.mysql.com ). The interactive web pages were implemented using HTML, CSS, JavaScript, and PHP languages ( www.php.net ), with several JavaScript libraries (JQuery.js, DataTable.js, and IGV.js) and Bootstrap framework (a popular framework for developing interactive websites) on Red Hat Linux powered by an Apache server ( www.apache.org ). This data portal is valuable for exploring aOutliers and their associated functional rare variants. With rareAPA, users can search, browse, and visualize important information on aOutliers in 49 human tissues. Users can search by gene or tissue name and scrutinize rare APA events among individuals in each tissue. Additionally, users can also visualize aOutliers using a scatter plot or explore them through a genome browser. Furthermore, rareAPA provides a curated list of prioritized RVs using the aWatershed algorithm, allowing users to examine rare variants and their aWatershed posterior scores. Additionally, rareAPA offers batch downloading of all single-tissue aOutliers and multi-tissue aOutliers. The rareAPA is freely available online without registration or login requirements. Declarations Code availability DaPars2 is available at https://github.com/3UTR/DaPars2, and IPAFinder can be accessed through https://github.com/ZhaozzReal/IPAFinder. The codes for mapping 3′aQTL are available at https://github.com/3UTR/3aQTL-pipe. The custom scripts and source codes for data analysis relevant to this study are available, under the MIT license, at Github repository: https://github.com/Xu-Dong/rareAPA and Zenodo: https://doi.org/10.5281/zenodo.10576656. Data availability The raw data of whole transcriptome and genome sequencing data from the GTEx project V8 are available at the database of Genotypes and Phenotypes (dbGaP) under the accession number: phs000424.v7.p2 [https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000424.v7.p2] 71 . All processed GTEx data, including gene expression outlier (eOutlier) and splicing outlier (sOutlier), are available via the GTEx portal (http://gtexportal.org). GWAS summary statistics used in this study were obtained from UK Biobank GWAS (https://www.nealelab.is/uk-biobank), Finn Gen (https://www.finngen.fi/en), and JENGER (http://jenger.riken.jp). The details about the GWAS summary statistics are listed in Supplementary Table 5. Genomic and functional annotations of rare variants are available via the Combined Annotation Dependent Depletion (CADD v1.5, https://cadd.gs.washington.edu/), and gnomAD v3.1(https://gnomad.broadinstitute.org/). The crosslinking and immunoprecipitation (CLIP) assay data for RNA binding proteins used in this study are available at The Encyclopedia of DNA Elements (ENCODE, https://www.encodeproject.org/). The data described in this study are freely available for querying, visualizing, and downloading at http://bioinfo.szbl.ac.cn/rareAPA/index.php, a website portal dedicated to rare APA. Author Contributions L.L. T.N., and W.L. conceived and supervised the project. X.Z., and Z.Z. performed the bioinformatics analysis with the help from K.X., and H.C. X.Z. constructed the website. Y.C. performed the experiments with the help from Z.W. and S.C., X.Z., T.N., W.L., and L.L. interpreted the data and wrote the manuscript. G.W. and S.X. reviewed and revised the manuscript. Competing Interests The authors declare no competing interests. Acknowledgments We thank Dr. Jian Yang from Westlake University for providing feedback on the manuscript. We also thank members of the Li laboratory for helpful discussions. This work was supported by the National Natural Science Foundation of China (no. 32100533, 32370721, 32288101, 32030020) and startup funds from Shenzhen Bay Laboratory to L.L. We also thank Qin Wang at the Shenzhen Bay Laboratory Supercomputing Center for high-level computing support and the Medical Science Data Center of Fudan University. References Taliun, D. et al. Sequencing of 53,831 diverse genomes from the NHLBI TOPMed Program. Nature 590 , 290-299 (2021). Keinan, A. & Clark, A.G. Recent explosive human population growth has resulted in an excess of rare genetic variants. Science 336 , 740-3 (2012). Consortium, U.K. et al. The UK10K project identifies rare variants in health and disease. Nature 526 , 82-90 (2015). Nelson, M.R. et al. An abundance of rare functional variants in 202 drug target genes sequenced in 14,002 people. Science 337 , 100-4 (2012). Tennessen, J.A. et al. Evolution and functional impact of rare coding variation from deep sequencing of human exomes. Science 337 , 64-9 (2012). Wang, Q. et al. Rare variant contribution to human disease in 281,104 UK Biobank exomes. Nature 597 , 527-532 (2021). Ferraro, N.M. et al. Transcriptomic signatures across human tissues identify functional rare genetic variation. Science 369 (2020). Li, X. et al. The impact of rare variation on gene expression across tissues. Nature 550 , 239-243 (2017). Hernandez, R.D. et al. Ultrarare variants drive substantial cis heritability of human gene expression. Nat Genet 51 , 1349-1355 (2019). Fresard, L. et al. Identification of rare-disease genes using blood transcriptome sequencing and large control cohorts. Nat Med 25 , 911-919 (2019). Mayr, C. What Are 3' UTRs Doing? Cold Spring Harb Perspect Biol 11 (2019). Tian, B. & Manley, J.L. Alternative polyadenylation of mRNA precursors. Nat Rev Mol Cell Biol 18 , 18-30 (2017). Mayr, C. Regulation by 3'-Untranslated Regions. Annu Rev Genet 51 , 171-194 (2017). Berkovits, B.D. & Mayr, C. Alternative 3' UTRs act as scaffolds to regulate membrane protein localization. Nature 522 , 363-7 (2015). Di Giammartino, D.C., Nishida, K. & Manley, J.L. Mechanisms and consequences of alternative polyadenylation. Mol Cell 43 , 853-66 (2011). Mitschka, S. & Mayr, C. Context-specific regulation and function of mRNA alternative polyadenylation. Nat Rev Mol Cell Biol 23 , 779-796 (2022). Singh, I. et al. Widespread intronic polyadenylation diversifies immune cell transcriptomes. Nat Commun 9 , 1716 (2018). Zhao, Z. et al. Cancer-associated dynamics and potential regulators of intronic polyadenylation revealed by IPAFinder using standard RNA-seq data. Genome Res 31 , 2095-2106 (2021). Masamha, C.P. et al. CFIm25 links alternative polyadenylation to glioblastoma tumour suppression. Nature 510 , 412-6 (2014). Park, H.J. et al. 3' UTR shortening represses tumor-suppressor genes in trans by disrupting ceRNA crosstalk. Nat Genet 50 , 783-789 (2018). Mittleman, B.E. et al. Alternative polyadenylation mediates genetic regulation of gene expression. Elife 9 (2020). Mariella, E., Marotta, F., Grassi, E., Gilotto, S. & Provero, P. The Length of the Expressed 3' UTR Is an Intermediate Molecular Phenotype Linking Genetic Variants to Complex Diseases. Front Genet 10 , 714 (2019). Li, L., Li, Y., Zou, X., Peng, F., Cui, Y., Wagner, E.J., Li, W. Population-scale genetic control of alternative polyadenylation and its association with human diseases. Quantitative Biology 10 , 44-54 (2022). Graham, R.R. et al. Three functional variants of IFN regulatory factor 5 (IRF5) define risk and protective haplotypes for human lupus. Proc Natl Acad Sci U S A 104 , 6758-63 (2007). Li, L. et al. An atlas of alternative polyadenylation quantitative trait loci contributing to complex trait and disease heritability. Nat Genet 53 , 994-1005 (2021). Feng, X., Li, L., Wagner, E.J. & Li, W. TC3A: The Cancer 3' UTR Atlas. Nucleic Acids Res 46 , D1027-D1030 (2018). Liu, Z. et al. Pan-cancer analysis identifies mutations in SUGP1 that recapitulate mutant SF3B1 splicing dysregulation. Proc Natl Acad Sci U S A 117 , 10305-10312 (2020). Alsafadi, S. et al. Genetic alterations of SUGP1 mimic mutant-SF3B1 splice pattern in lung adenocarcinoma and other cancers. Oncogene 40 , 85-96 (2021). Chen, E.Y. et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics 14 , 128 (2013). McLaren, W. et al. The Ensembl Variant Effect Predictor. Genome Biol 17 , 122 (2016). Rentzsch, P., Witten, D., Cooper, G.M., Shendure, J. & Kircher, M. CADD: predicting the deleteriousness of variants throughout the human genome. Nucleic Acids Res 47 , D886-D894 (2019). Bogard, N., Linder, J., Rosenberg, A.B. & Seelig, G. A Deep Neural Network for Predicting and Engineering Alternative Polyadenylation. Cell 178 , 91-106 e23 (2019). Zhao, Z. et al. Comprehensive characterization of somatic variants associated with intronic polyadenylation in human cancers. Nucleic Acids Res 49 , 10369-10381 (2021). Yeo, G. & Burge, C.B. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. J Comput Biol 11 , 377-94 (2004). Alipanahi, B., Delong, A., Weirauch, M.T. & Frey, B.J. Predicting the sequence specificities of DNA- and RNA-binding proteins by deep learning. Nat Biotechnol 33 , 831-8 (2015). Jenal, M. et al. The poly(A)-binding protein nuclear 1 suppresses alternative cleavage and polyadenylation sites. Cell 149 , 538-53 (2012). Dominguez, D. et al. Sequence, Structure, and Context Preferences of Human RNA Binding Proteins. Mol Cell 70 , 854-867 e9 (2018). Linder, J., Koplik, S.E., Kundaje, A. & Seelig, G. Deciphering the impact of genetic variation on human polyadenylation using APARENT2. Genome Biol 23 , 232 (2022). Hamosh, A., Scott, A.F., Amberger, J.S., Bocchini, C.A. & McKusick, V.A. Online Mendelian Inheritance in Man (OMIM), a knowledgebase of human genes and genetic disorders. Nucleic Acids Res 33 , D514-7 (2005). Slavotinek, A.M. et al. Mutation analysis of the MKKS gene in McKusick-Kaufman syndrome and selected Bardet-Biedl syndrome patients. Hum Genet 110 , 561-7 (2002). Stone, D.L. et al. Mutation of a gene encoding a putative chaperonin causes McKusick-Kaufman syndrome. Nat Genet 25 , 79-82 (2000). Slavotinek, A.M. et al. Mutations in MKKS cause Bardet-Biedl syndrome. Nat Genet 26 , 15-6 (2000). Katsanis, N. et al. Mutations in MKKS cause obesity, retinal dystrophy and renal malformations associated with Bardet-Biedl syndrome. Nat Genet 26 , 67-70 (2000). Wuyts, W. et al. Mutations in the EXT1 and EXT2 genes in hereditary multiple exostoses. Am J Hum Genet 62 , 346-54 (1998). Stickens, D. et al. The EXT2 multiple exostoses gene defines a family of putative tumour suppressor genes. Nat Genet 14 , 25-32 (1996). Quintas-Cardama, A. & Cortes, J. Molecular biology of bcr-abl1-positive chronic myeloid leukemia. Blood 113 , 1619-30 (2009). Salesse, S. & Verfaillie, C.M. BCR/ABL: from molecular mechanisms of leukemia induction to treatment of chronic myelogenous leukemia. Oncogene 21 , 8547-59 (2002). Weiner, D.J. et al. Statistical and functional convergence of common and rare genetic influences on autism at chromosome 16p. Nat Genet 54 , 1630-1639 (2022). Schrode, N. et al. Synergistic effects of common schizophrenia risk variants. Nat Genet 51 , 1475-1485 (2019). Singh, T. et al. Rare coding variants in ten genes confer substantial risk for schizophrenia. Nature 604 , 509-516 (2022). Cui, Y. et al. Alternative polyadenylation transcriptome-wide association study identifies APA-linked susceptibility genes in brain disorders. Nat Commun 14 , 583 (2023). Chen, H. et al. A distinct class of pan-cancer susceptibility genes revealed by alternative polyadenylation transcriptome-wide association study. medRxiv , 2023.02.28.23286554 (2023). Dong, G. et al. DDX18 drives tumor immune escape through transcription-activated STAT1 expression in pancreatic cancer. Oncogene 42 , 3000-3014 (2023). Redmond, A.M. et al. Genomic interaction between ER and HMGB2 identifies DDX18 as a novel driver of endocrine resistance in breast cancer cells. Oncogene 34 , 3871-80 (2015). McFarland, J.M. et al. Improved estimation of cancer dependencies from large-scale RNAi screens using model-based normalization and data integration. Nat Commun 9 , 4610 (2018). Tsherniak, A. et al. Defining a Cancer Dependency Map. Cell 170 , 564-576 e16 (2017). Demontis, D. et al. Genome-wide analyses of ADHD identify 27 risk loci, refine the genetic architecture and implicate several cognitive domains. Nat Genet 55 , 198-208 (2023). Wu, N. et al. TBX6 null variants and a common hypomorphic allele in congenital scoliosis. N Engl J Med 372 , 341-50 (2015). Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29 , 15-21 (2013). Quinlan, A.R. & Hall, I.M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26 , 841-2 (2010). Zou, X. et al. Using population-scale transcriptomic and genomic data to map 3' UTR alternative polyadenylation quantitative trait loci. STAR Protoc 3 , 101566 (2022). Ma, X. et al. ipaQTL-atlas: an atlas of intronic polyadenylation quantitative trait loci across human tissues. Nucleic Acids Res 51 , D1046-D1052 (2023). Stegle, O., Parts, L., Piipari, M., Winn, J. & Durbin, R. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses. Nat Protoc 7 , 500-7 (2012). Gudmundsson, S. et al. Variant interpretation using population databases: Lessons from gnomAD. Hum Mutat 43 , 1012-1030 (2022). Wang, R., Zheng, D., Yehia, G. & Tian, B. A compendium of conserved cleavage and polyadenylation events in mammalian genes. Genome Res 28 , 1427-1441 (2018). Wang, R., Nambiar, R., Zheng, D. & Tian, B. PolyA_DB 3 catalogs cleavage and polyadenylation sites identified by deep sequencing in multiple genomes. Nucleic Acids Res 46 , D315-D319 (2018). Shabalin, A.A. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28 , 1353-8 (2012). Giambartolomei, C. et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet 10 , e1004383 (2014). Zhao, H. et al. CrossMap: a versatile tool for coordinate conversion between genome assemblies. Bioinformatics 30 , 1006-7 (2014). Grishin, D. & Gusev, A. Allelic imbalance of chromatin accessibility in cancer identifies candidate causal risk variants and their mechanisms. Nat Genet 54 , 837-849 (2022). Consortium, G.T. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 369 , 1318-1330 (2020). Additional Declarations There is NO Competing Interest. Supplementary Files SupplementaryInformation.docx SupplementalTables.xlsx Supplementary Table 1. Rare single-nucleotide variants (SNVs) with Combined Annotation-Dependent Depletion (CADD) scores ≥ 15 nearby aOutlier genes. Supplementary Table 2. aOutlier genes with nearby rare indels in the outlier individuals. Supplementary Table 3. Genomic features used in aWatershed model. Supplementary Table 4. aWatershed prioritized functional RVs impacting APA. Supplementary Table 5. The metadata of the 1,234 GWAS summary statistics from the UK Biobank (UKBB) and literature. Supplementary Table 6. aWatershed prioritized APA-related rare variants that overlapped with variants in GWAS summary data of 1,186 UKBB traits. Supplementary Table 7. Genes with evidence of colocalization between GWAS of UKBB traits and 3'aQTL in any GTEx tissues. Supplementary Table 8. Oligos used for plasmid construction and primers used for 3' RACE. Cite Share Download PDF Status: Published Journal Publication published 16 Jan, 2025 Read the published version in Nature Communications → 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-3907149","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":276614910,"identity":"2c52ba69-780f-456a-8c7b-f24293b0b821","order_by":0,"name":"Lei Li","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAu0lEQVRIiWNgGAWjYBACPgbGBoYHBjYMDBLMDcRpYQNpSTBIA2phJFoLECQwHCZFC3sz0JaC8/bysxsbGH7UMMibE9TCcxDksNuJjXMONjD2HGMw3EnIMjaJRLCWBGYQgxfEPkBIi/xDkLJz9iC9jH+J0iIBDrEDjD1ALczE2cKT2HAgwSA5cYbMwYbDMsckDDcQ0sLPfvzhgw9/7IAh1nzw4ZsaG3mCtoDAASSGBBHqR8EoGAWjYBQQBAAbGTtvWXOnZwAAAABJRU5ErkJggg==","orcid":"https://orcid.org/0000-0003-3924-2544","institution":"Shenzhen Bay Laboratory","correspondingAuthor":true,"submittingAuthor":false,"prefix":"","firstName":"Lei","middleName":"","lastName":"Li","suffix":""},{"id":276614911,"identity":"19786b9f-30d4-4fac-aed0-1daaedc80177","order_by":1,"name":"Xudong Zou","email":"","orcid":"https://orcid.org/0000-0002-2958-0438","institution":"Institute of Systems and Physical Biology, Shenzhen Bay Laboratory","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Xudong","middleName":"","lastName":"Zou","suffix":""},{"id":276614912,"identity":"007b7399-1949-46ed-b225-30710f21340d","order_by":2,"name":"Zhaozhao Zhao","email":"","orcid":"","institution":"Fudan university","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Zhaozhao","middleName":"","lastName":"Zhao","suffix":""},{"id":276614913,"identity":"0e23dab7-ff16-4119-97be-b507e7af2f16","order_by":3,"name":"Yu Chen","email":"","orcid":"","institution":"Fudan university","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Yu","middleName":"","lastName":"Chen","suffix":""},{"id":276614914,"identity":"6352281f-c3a6-4fa6-a946-142721f24333","order_by":4,"name":"Kewei Xiong","email":"","orcid":"","institution":"Shenzhen Bay Laboratory","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Kewei","middleName":"","lastName":"Xiong","suffix":""},{"id":276614915,"identity":"82c6ead1-112b-4679-b7a1-3711585aade0","order_by":5,"name":"Zeyang Wang","email":"","orcid":"https://orcid.org/0000-0001-5735-0675","institution":"Shenzhen Bay Laboratory","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Zeyang","middleName":"","lastName":"Wang","suffix":""},{"id":276614916,"identity":"490e0787-259d-40fd-97ea-16e7ce8d2f23","order_by":6,"name":"Shuxin Chen","email":"","orcid":"","institution":"Shenzhen Bay Laboratory","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Shuxin","middleName":"","lastName":"Chen","suffix":""},{"id":276614917,"identity":"08b8b36b-3c5c-4b95-b710-d9deaacb8207","order_by":7,"name":"Hui Chen","email":"","orcid":"","institution":"Institute of Systems and Physical Biology, Shenzhen Bay Laboratory","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Hui","middleName":"","lastName":"Chen","suffix":""},{"id":276614918,"identity":"638d3f28-f970-4537-9a72-21b3b5f79ff4","order_by":8,"name":"Gong-Hong Wei","email":"","orcid":"https://orcid.org/0000-0001-6546-9334","institution":"Fudan University Shanghai Cancer Center \u0026 MOE Key Laboratory of Metabolism and Molecular Medicine and Department of Biochemistry and Molecular Biology of School Basic Medical Sciences, Shanghai Medi","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Gong-Hong","middleName":"","lastName":"Wei","suffix":""},{"id":276614919,"identity":"eddbb36f-d9ec-40f9-977c-45a6b6acd82a","order_by":9,"name":"Shuhua Xu","email":"","orcid":"","institution":"School of Life Sciences, Fudan University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Shuhua","middleName":"","lastName":"Xu","suffix":""},{"id":276614920,"identity":"2bae91f1-4abf-407d-affb-02723f2153c9","order_by":10,"name":"Wei Li","email":"","orcid":"https://orcid.org/0000-0001-9931-5990","institution":"University of California, Irvine","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Wei","middleName":"","lastName":"Li","suffix":""},{"id":276614921,"identity":"9a5845b8-4b9d-4af8-bd59-333c7839fd55","order_by":11,"name":"Ting Ni","email":"","orcid":"https://orcid.org/0000-0001-7007-1072","institution":"Collaborative Innovation Center of Genetics and Development, Human Phenome Institute, School of Life Sciences, Fudan University","correspondingAuthor":false,"submittingAuthor":false,"prefix":"","firstName":"Ting","middleName":"","lastName":"Ni","suffix":""}],"badges":[],"createdAt":"2024-01-28 23:20:37","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3907149/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3907149/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1038/s41467-024-55407-3","type":"published","date":"2025-01-16T05:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":52189063,"identity":"3d5effd3-0763-43ae-8fd1-a2de389c77de","added_by":"auto","created_at":"2024-03-07 19:02:18","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":458389,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAtlas of human alternative polyadenylation (APA) outliers (aOutliers).\u003c/strong\u003e \u003cstrong\u003ea.\u003c/strong\u003e Schematic illustrating the overall design of this study. \u003cstrong\u003eb.\u003c/strong\u003e Distribution of 3′ untranslated region (UTR; blue) and intronic (red) aOutliers across the human genome. Genes with the highest (for positive median Z-scores) or lowest (for negative median Z-scores) Z-score at each chromosome region were labeled. \u003cstrong\u003ec.\u003c/strong\u003eDistribution of the number of tissues in which 3′ UTR aOutliers (deep blue) and intronic aOutliers (red) were detected. \u003cstrong\u003ed.\u003c/strong\u003e Comparison of average tissue counts of 3′ UTR aOutliers (deep blue; n=603 genes) and intronic aOutliers (red; n=100 genes). Box plots show the median and first and third quartiles, and whiskers extend up to 1.5 times the interquartile range. \u003cstrong\u003ee.\u003c/strong\u003e RNA sequencing (RNA-seq) read coverage of the \u003cem\u003eSUGP1\u003c/em\u003e gene 3′ UTR in outlier individuals (red) and nonoutlier individuals (gray) in the Lung and Brain hippocampus. \u003cstrong\u003ef. \u003c/strong\u003eMedian\u003cstrong\u003e \u003c/strong\u003eZ-score distribution of the \u003cem\u003eSUGP1\u003c/em\u003e and \u003cem\u003eCOL4A2\u003c/em\u003e genes across individuals. Outliers are highlighted with red dots. \u003cstrong\u003eg.\u003c/strong\u003e RNA-seq read coverage of the \u003cem\u003eCOL4A2\u003c/em\u003e gene at the region of “exon5-intron-exon6” in outlier individuals (red) and nonoutlier individuals (gray). For the data shown in this figure, significance was calculated using the single-tailed Wilcoxon rank–sum test.\u003c/p\u003e","description":"","filename":"image1.png","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/d505c26cdc61a668ed0c98cd.png"},{"id":52187552,"identity":"be46c1dd-ff90-48bf-aeae-35ab478d0446","added_by":"auto","created_at":"2024-03-07 18:54:17","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":520511,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eOutliers are distinct from other molecular outliers.\u003c/strong\u003e \u003cstrong\u003ea.\u003c/strong\u003e Number of aOutlier genes also detected by analysis of expression outliers (eOutliers) and splicing outliers (sOutlier) in the same dataset. \u003cstrong\u003eb. \u003c/strong\u003eExample of an aOutlier-only gene, not detected by eOutlier and sOutlier analysis. \u003cstrong\u003ec.\u003c/strong\u003e Analysis of 5¢ UTR length, coding sequence (CDS) length, and 3¢ UTR length in aOutlier genes (n=562 genes) compared to eOutlier genes (n=1,833 genes). \u003cem\u003eP\u003c/em\u003e-values were calculated using the one-sided Wilcoxon rank–sum test. Box plots show the median and first and third quartiles, and whiskers extend up to 1.5 times the interquartile range. \u003cstrong\u003ed\u003c/strong\u003e. Analysis of GC-content in 3¢ UTR regions of aOutlier (n=562) and eOutlier genes (n=1,833). \u003cem\u003eP\u003c/em\u003e-values were calculated using the one-sided Wilcoxon rank–sum test. Box plots show the median and first and third quartiles, and whiskers extend up to 1.5 times the interquartile range. \u003cstrong\u003ee.\u003c/strong\u003e The proportion of aOutliers with nearby RVs of different categories. aOutliers were stratified by absolute median Z-score thresholds: Z \u0026lt; 1 (nonoutlier), Z \u0026gt; 3, Z \u0026gt; 4, Z \u0026gt; 5, and Z \u0026gt; 10. RV categories were assigned by VEP (v.104), and some were manually merged. Terms are defined as follows: \"Splice\" includes RVs at the splice donor site, splice acceptor site, and splice region; \"Stop\" includes RVs resulting in stop gained, stop lost, and start lost; \"Coding\" includes missense variant, stop retained variant; \"otherCoding\" includes CDS and synonymous variant; and \"other noncoding\" includes downstream gene variant, upstream gene variant, and non-coding transcript exon variant. \u003cstrong\u003ef.\u003c/strong\u003e Enrichment of RVs of different categories in aOutliers (red), eOutliers (green), and sOutliers (blue). \"PAS50bp\" represents the 50 base pairs upstream of the annotated poly(A) site. \"Conserved\" RVs are defined by mammalian phaseCons score \u0026gt; 0.9. Data are presented as log2 odds ratios (OR) and 95% confidence intervals (CI). \u003cstrong\u003eg.\u003c/strong\u003e Enrichment of deleterious single-nucleotide variants (SNVs) in aOutliers; variants within 1 kb of aOutlier genes were counted. Data are presented as ORs and 95% CIs. \u003cstrong\u003eh.\u003c/strong\u003e Enrichment of RVs in single-tissue aOutliers. Data are presented as odds ratio and 95% CI (y-axis) for each tissue (x-axis).\u003c/p\u003e","description":"","filename":"image2.png","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/a1b6b86c512f78a5db6cad4f.png"},{"id":52189059,"identity":"9052cc6c-77a7-4a7c-afa7-5b844f0a7d5f","added_by":"auto","created_at":"2024-03-07 19:02:17","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":567605,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunctional rare variants (RVs) and RBPs associated with aOutliers. a.\u003c/strong\u003e Enrichment of aOutlier-associated RVs in poly(A) signal (PAS) and AU-rich motifs. Central dots show the log transformed odds ratio, and lines show 95% confidence intervals\u003cstrong\u003e. b.\u003c/strong\u003e The aOutlier of gene \u003cem\u003eMKKS\u003c/em\u003e. Data are presented as median Z score (y-axis) and samples (x-axis) ranked by median Z score. The outlier individual was represented as a red dot. \u003cstrong\u003ec-d.\u003c/strong\u003e The minigenes and 3′ RACE assays for the 3′ UTRs of \u003cem\u003eMKKS\u003c/em\u003e (c) and \u003cem\u003eSUGP1\u003c/em\u003e (d) in HEK293 and HeLa cells. The structures of each minigene reporter are shown at the top, and PCR priming data for both long and short isoforms are presented below. Tested RVs that alter PAS motifs are indicated with red asterisks. \u003cem\u003eGAPDH\u003c/em\u003e was used as a reference in all assays. \u003cstrong\u003ee.\u003c/strong\u003e Enrichment of intronic aOutlier-associated RVs that disrupt splice sites compared with RVs associated with nonoutliers. The splice site was defined as nine bp (indicated as \"D-3\" to \"D+6\") for the donor site and six bp (indicated as \"A-3\" to \"A+3\") for the acceptor site. Enrichments are presented as ORs and 95% CIs. \u003cstrong\u003ef-g.\u003c/strong\u003e aOutliers in gene \u003cem\u003eTXNRD2\u003c/em\u003e (f) and \u003cem\u003eHMGCL\u003c/em\u003e (g). \u003cstrong\u003eh.\u003c/strong\u003e Consensus donor site sequences in outlier and nonoutlier individuals. \u003cstrong\u003ei. \u003c/strong\u003eStrength of donor splice sites in intronic aOutlier individuals (MUT) and controls (WT). The center horizontal lines represent the median values; boxes span from the 25th to 75th percentile, and whiskers extend to 1.5 × interquartile range. Significance was calculated using the single-tailed Wilcoxon rank–sum test. \u003cstrong\u003ej-k.\u003c/strong\u003e Minigenes and 3′ RACE assays for the intronic APA of \u003cem\u003eTXNRD2\u003c/em\u003e (j) and \u003cem\u003eHMGCL\u003c/em\u003e (k) in HEK293 and HeLa cells. The structures of each minigene reporter are shown at the top, and PCR priming data for both long and short isoforms are presented below. Tested RVs that alter PAS motifs are indicated with red asterisks. \u003cem\u003eGAPDH\u003c/em\u003e was used as a reference in all assays. \u003cstrong\u003el\u003c/strong\u003e. Enrichment of RNA-binding protein (RBP) binding regions in aOutlier-associated RVs compared to nonoutlier RVs. Data are presented as -log\u003csub\u003e10\u003c/sub\u003e(\u003cem\u003eP\u003c/em\u003e) (y-axis) and odds ratio (dot size). \u003cstrong\u003em. \u003c/strong\u003eThe scatter plot shows aOutlier in gene TOLLIP. \u003cstrong\u003en.\u003c/strong\u003e One rare variant in aOutlier individual of gene TOLLIP was identified located at the binding regions of CSTF2T. The RNA-seq reads coverage of the aOutlier and three nonoutlier individuals in the 3′ UTR region were presented, and the binding peaks of the CSTF2T regulator by eCLIP assay were presented below. The rare variant was presented in red.\u003c/p\u003e","description":"","filename":"image3.png","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/0ffdb0be076e8f2f634da8a5.png"},{"id":52187554,"identity":"e97345fa-92b1-4c3e-b527-84ab1f03b16e","added_by":"auto","created_at":"2024-03-07 18:54:17","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":427426,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDevelopment and evaluation of the APA-based Watershed (aWatershed) model.\u003c/strong\u003e \u003cstrong\u003ea.\u003c/strong\u003e Performance of the aWatershed model (red) compared to the RIVER model (blue) and the genomic annotation model (orange). Data are presented as the area under the precision–recall curve (AUC-PR). \u003cstrong\u003eb.\u003c/strong\u003e The correlation between aWatershed predicted posteriors and GAM predicted posteriors. RVs with posteriors \u0026gt; 0.5 in either group were filled with red. \u003cstrong\u003ec. \u003c/strong\u003eEdge weights connecting top genomic annotation features to latent regulatory variables in aOutlier signal (red), eOutlier signal (green), and sOutlier signal (blue), ranked by weight in decreasing order. The top three most influential genomic features are highlighted in bold font. \u003cstrong\u003ed. \u003c/strong\u003eThe proportion of RVs leading to aOutliers. RVs were stratified based on aWatershed (red) and GAM (orange) posterior probability for APA signal. \u003cstrong\u003ee, f.\u003c/strong\u003e Evaluation of aWatershed-prioritized RVs using the data estimated from a published massively parallel reporter assay. RVs were stratified into two groups based on aWatershed APA posterior probabilities (i.e., probabilities \u0026gt; 0.5 (red) and £ 0.5 (blue)), and poly(A) site usage change was compared for reference and alternative alleles in each group (e), and proportion of large-effect (absolute log Fold change \u0026gt; 1) was compared between the two groups (f). Box plots show the median and first and third quartiles, and whiskers extend up to 1.5 times the interquartile range. \u003cem\u003eP\u003c/em\u003e-values were calculated using the one-sided Wilcoxon rank–sum test.\u003c/p\u003e","description":"","filename":"image4.png","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/8c52e9232ea257fb833cec96.png"},{"id":52187553,"identity":"51a6b360-7c47-4ec6-b112-604e87d39909","added_by":"auto","created_at":"2024-03-07 18:54:17","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":379727,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTrait effect sizes for aOutlier RVs prioritized by aWatershed. a. \u003c/strong\u003eComparison of the trait effect size of aOutlier prioritized RV (red) nearby genes with evidence of colocalization to non-prioritized RVs (gray) nearby the same genes. Box plots show the median and first and third quartiles, and whiskers extend up to 1.5 times the interquartile range. \u003cem\u003eP\u003c/em\u003evalue was calculated using the one-sided Wilcoxon rank–sum test (n = 77,388). \u003cstrong\u003eb,c.\u003c/strong\u003eDistribution (red) of odds estimated from permutation test assessing how often randomly drawn aWatershed-prioritized RVs have larger effect sizes in GWAS of height (b) or high blood pressure (c) than matched non-prioritized RVs across genes. The null distribution (in gray) of odds was obtained from a permutation test by randomly drawning two RVs from a non-prioritized RV set only. \u003cem\u003eP\u003c/em\u003e value was calculated from the one-sided Wilcoxon rank‒sum test. \u003cstrong\u003ed. \u003c/strong\u003eManhattan plot (left) across 20 Mb in chromosome 2 for GWAS signals of height (50_irnt) in the UKBB. The aOutlier prioritized RV rs112567314 in the colocalized region was highlighted by the red triangle, and the GWAS lead SNP is indicated by a blue diamond. The blue square denotes RVs prioritized by eOutliers or sOutliers in the same region. UKBB MAF \u003cem\u003evs\u003c/em\u003e. effect size for all variants within 1Mb of the aOutlier prioritized RV was shown on the right. \u003cstrong\u003ee.\u003c/strong\u003e Manhattan plot (left) across 20 Mb in chromosome 4 for GWAS signals of high blood pressure (6150_4) in the UKBB, and the scatter plot (right) shows the UKBB MAF \u003cem\u003evs.\u003c/em\u003e effect size for all variants in a 2Mb region cross the aOutleir prioritized RV (rs149094812). The red triangle highlights the aOutlier prioritized RV, and the blue diamond highlights the GWAS lead SNP.\u003c/p\u003e","description":"","filename":"image5.png","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/c9cf51448a36797753e050ab.png"},{"id":52187557,"identity":"7021a4b8-dfdf-47db-a6e6-21ea7afe3fe6","added_by":"auto","created_at":"2024-03-07 18:54:17","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":470626,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe convergence effect of rare and common variants links APA to human diseases. a\u003c/strong\u003e. Intersection of disease risk genes identified by 3′aQTLs using gene-based methods, including colocalization and transcriptome-wide association study and genes associated with aWatershed-prioritized functional APA RVs. 278 genes associated with RVs having aWatershed posterior score \u0026gt; 0.5 were involved. \u003cstrong\u003eb.\u003c/strong\u003e The top 20 intersected genes are ranked by aWatershed posterior score. \u003cstrong\u003ec.\u003c/strong\u003e Distribution of dependency scores estimated from CRISPR-Cas9 essentiality screening assays in cancer cells for genes associated with RVs prioritizing by aWatershed and also identified as cancer risk genes by 3′aQTLs analysis. \u003cstrong\u003ed.\u003c/strong\u003e Boxplot shows the association between the common variant (rs1052628) and 3′ UTR APA of \u003cem\u003eDDX18\u003c/em\u003e. The outlier individual with the rare variant was labeled. \u003cstrong\u003ee. \u003c/strong\u003eaOutliers (left) and eOutliers (right) of gene \u003cem\u003eDDX18\u003c/em\u003e. Data are presented as median\u003cstrong\u003e \u003c/strong\u003eZ-score (y-axis), and individuals (x-axis) are ranked by median Z-score. Outliers are highlighted with red dots. \u003cstrong\u003ef.\u003c/strong\u003e RNA-seq reads coverage of \u003cem\u003eDDX18\u003c/em\u003e 3′ UTR region and the last second exon in the outlier individual (red) and non-outlier (control) individual (gray). \u003cstrong\u003eg.\u003c/strong\u003e Minigenes of \u003cem\u003eDDX18\u003c/em\u003e 3′ UTR containing the common APA variant and the rare APA variant. \u003cstrong\u003eh.\u003c/strong\u003e 3′ RACE assays with \u003cem\u003eDDX18\u003c/em\u003e minigenes containing only the RV or only the common variant, or both the RV and the common variant. Assays were performed in HEK293T cells. \u003cstrong\u003ei.\u003c/strong\u003e The bar plot shows the usage of dPAS in each minigene measured through image J. Error bar represents the standard deviation and a two-sided student t-test was used to test the difference (n=6 independent experimental replicates). \u003cstrong\u003ej.\u003c/strong\u003e The 3′ RACE assays of \u003cem\u003eDDX18\u003c/em\u003e minigenes in MCF7 cells. \u003cem\u003eGAPDH\u003c/em\u003e was used as the loading control.\u003c/p\u003e","description":"","filename":"image6.png","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/249444ef7c298dc8de2e53bb.png"},{"id":74038330,"identity":"1691922c-989a-4d4f-bfae-141f5acda9dd","added_by":"auto","created_at":"2025-01-17 08:05:35","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":4279126,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/bebb4a88-bc33-46d3-a49b-2f67a17fa1ed.pdf"},{"id":52187559,"identity":"76310181-28d6-4f6c-a40a-d01bbad8998a","added_by":"auto","created_at":"2024-03-07 18:54:18","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":8565324,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryInformation.docx","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/8ff1a21b9812e0b5f9b06a4f.docx"},{"id":52187558,"identity":"ddc3bde5-83c5-49a4-9f4c-8a515e8629aa","added_by":"auto","created_at":"2024-03-07 18:54:18","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":31634897,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 1. Rare single-nucleotide variants (SNVs) with Combined Annotation-Dependent Depletion (CADD) scores ≥ 15 nearby aOutlier genes.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 2. aOutlier genes with nearby rare indels in the outlier individuals.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 3. Genomic features used in aWatershed model.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 4. aWatershed prioritized functional RVs impacting APA.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 5. The metadata of the 1,234 GWAS summary statistics from the UK Biobank (UKBB) and literature.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 6. aWatershed prioritized APA-related rare variants that overlapped with variants in GWAS summary data of 1,186 UKBB traits.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 7. Genes with evidence of colocalization between GWAS of UKBB traits and 3'aQTL in any GTEx tissues.\u003c/p\u003e\n\u003cp\u003eSupplementary Table 8. Oligos used for plasmid construction and primers used for 3' RACE.\u003c/p\u003e","description":"","filename":"SupplementalTables.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-3907149/v1/da3b7442576674d68fb2be8f.xlsx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Impact of Rare Non-coding Variants on Human Diseases through Alternative Polyadenylation Outliers","fulltext":[{"header":"Introduction","content":"\u003cp\u003eThe human genome harbors numerous rare genetic variants\u003csup\u003e1\u003c/sup\u003e, each with a minor allele frequency (MAF) of less than 1%. Many of these rare variants strongly contribute to human diseases\u003csup\u003e2\u0026ndash;5\u003c/sup\u003e. While exome sequencing of large population cohorts has identified numerous rare protein-coding variants associated with both common and rare diseases\u003csup\u003e6\u003c/sup\u003e, the vast majority of rare variants (RVs) are located in non-coding regions. These non-coding RVs do not function through altering the protein sequences, thereby posing a significant challenge in interpreting their functions. To address this challenge, analysis of population-scale transcriptomic data has been used to uncover functional rare non-coding variants affecting gene expression or splicing outliers\u003csup\u003e7\u0026ndash;10\u003c/sup\u003e. Despite these efforts, a significant portion of disease-associated RVs remain uncharacterized.\u003c/p\u003e \u003cp\u003eAlternative polyadenylation (APA) of mRNA is a widespread post-transcriptional regulatory mechanism observed across various species. By employing different polyadenylation sites within 3\u0026prime;untranslated regions (3\u0026prime; UTRs), genes can produce various mRNA isoforms with either shortened or extended 3\u0026prime; UTRs. These 3\u0026prime; UTRs contain many regulatory elements that modulate the abundance or localization of the mRNA and protein\u003csup\u003e11\u0026ndash;15\u003c/sup\u003e. Moreover, APA can also occur in intronic regions, leading to truncated mRNA or proteins \u003csup\u003e16,17\u003c/sup\u003e. Accordingly, disruptions in APA events have been increasingly implicated in many human diseases\u003csup\u003e17\u0026ndash;19\u003c/sup\u003e. For example, altered APA leading to 3\u0026prime; UTR shortening of competing-endogenous RNAs for tumor suppressor genes can result in the release of microRNAs, inhibiting tumor suppressor genes and potentially leading to tumorigenesis\u003csup\u003e20\u003c/sup\u003e. Moreover, recent studies have reported the ubiquitous genetic regulation of APA, highlighting its importance in the functional interpretation of disease-associated non-coding variants\u003csup\u003e21\u0026ndash;23\u003c/sup\u003e. A notable example is a single-nucleotide polymorphism (SNP; rs10954213) within the 3\u0026prime; UTR of interferon regulatory factor 5 (\u003cem\u003eIRF5)\u003c/em\u003e, which can alter the length and stability of its 3\u0026prime; UTR, thereby contributing to systemic lupus erythematosus susceptibility\u003csup\u003e24\u003c/sup\u003e. In our previous study, we built an atlas of human 3\u0026prime; UTR APA quantitative trait loci (3\u0026prime;aQTLs) across human tissues, identifying approximately 0.4\u0026nbsp;million common SNPs associated with interindividual APA changes, which colocalize with 16.1% of trait-associated genetic variants\u003csup\u003e25\u003c/sup\u003e. Yet, these studies mainly focus on assessing the APA regulation of common variants. To our knowledge, the effect of RVs on APA has not been explored.\u003c/p\u003e \u003cp\u003eHere, to better understand the impact of RVs on APA, we systematically analyzed aberrant APA events across 49 human tissues from the Genotype-Tissue Expression Project (GTEx). We identified 1,534 multi-tissue APA outliers (aOutliers) from European individuals. Intriguingly, 74.2% of these aOutliers are associated with genes not previously identified in outlier analysis of other molecular phenotypes (e.g., expression or splicing). These aOutliers exhibit distinct characteristics, such as unique 3\u0026prime; UTR length and GC-contents, setting them apart from other types of molecular outliers. Moreover, a significant enrichment of deleterious RVs was observed in regions proximal to these aOutliers. To prioritize functional RVs impacting APA, we developed a Bayesian hierarchical model and identified a distinct set of RVs with large effect sizes on human complex traits and disease phenotypes. Intriguingly, we observed and demonstrated strong convergence effects between prioritized RVs and common variants in regulating 3\u0026prime; UTR APA, exemplified by the combinatorial regulation of APA in \u003cem\u003eDDX18\u003c/em\u003e. Lastly, to facilitate broad access to aOutliers-associated RVs, we have constructed a user-friendly portal at \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://bioinfo.szbl.ac.cn/rareAPA/index.php\u003c/span\u003e\u003cspan address=\"http://bioinfo.szbl.ac.cn/rareAPA/index.php\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. Collectively, our findings indicate that APA highlights a specific set of RVs with significant impacts on human traits and diseases, providing a new avenue for interpreting rare human non-coding genetic variants.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eThe landscape of APA outliers across 49 human tissues\u003c/h2\u003e \u003cp\u003eWe first conducted a comprehensive identification of 3\u0026prime; UTR and intronic APA events in 15,201 GTEx RNA-seq samples from 49 human tissues of 838 individuals (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea) using our Dapars2\u003csup\u003e25,26\u003c/sup\u003e and IPAFinder\u003csup\u003e18\u003c/sup\u003e algorithms, respectively (see Materials and Methods) (Supplementary Fig.\u0026nbsp;1). Considering the potential influence of many known and unknown technical confounders on APA usage among samples, we regressed out these confounders, such as age, sex, sequencing platform, and other hidden confounders inferred by using probabilistic estimation of expression residuals (PEER) factors (Supplementary Fig.\u0026nbsp;2). We then calculated Z-scores for the PEER-adjusted 3\u0026prime; UTR and intronic APA usage in each tissue to identify individuals with aberrant APA usage for a specific gene, which we refer to as APA outliers (aOutliers) with an absolute Z-score\u0026thinsp;\u0026gt;\u0026thinsp;3. The individuals and genes were designated as \u0026ldquo;aOutlier individuals\u0026rdquo; and \u0026ldquo;aOutlier genes\u0026rdquo;, respectively. Importantly, a single gene could be associated with multiple outlier individuals, and conversely, one individual could be an aOutlier individual for multiple genes. Our analysis of these aOutliers revealed that, on average, 68.5% of all transcripts per tissue were present in at least one outlier individual (Supplementary Fig.\u0026nbsp;3a). The number of aOutlier genes strongly correlated (Spearman\u0026rsquo;s correlation rho\u0026thinsp;=\u0026thinsp;0.91, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;2.2 \u0026times; 10\u003csup\u003e‒16\u003c/sup\u003e) with sample size across tissues (Supplementary Fig.\u0026nbsp;3b), suggesting that additional aOutlier genes might be discovered as more RNA-seq samples become available. This strong sample size correlation was further confirmed by down-sampling analyses in representative tissues (Supplementary Fig.\u0026nbsp;3c). Moreover, we noticed that the incidence of an aOutlier identified in one tissue being replicated in another was as low as 14.3% (Supplementary Fig.\u0026nbsp;4), indicating a significant degree of tissue-specificity among these single-tissue aOutliers.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe further defined multi-tissue aOutliers based on aberrant APA usage across five or more tissues (see Materials and Methods). From this analysis, we identified a total of 2,147 multi-tissue aOutliers, comprising 1,930 3\u0026prime; UTR aOutliers and 217 intronic aOutliers based on the genomic location of the APA event. Focusing specifically on the 715 European individuals, in whom we detected 1,534 multi-tissue aOutliers, including 1,334 3\u0026prime; UTR and 200 intronic aOutliers (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb and Supplementary Figs.\u0026nbsp;5 and 6). In our further investigation into the distribution of multi-tissue aOutliers across different tissues, we found that intronic aOutliers exhibited a broader replication pattern than 3\u0026prime; UTR aOutliers (one-sided Wilcoxon rank\u0026ndash;sum test \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;3.35 \u0026times; 10\u003csup\u003e‒14\u003c/sup\u003e; Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ec, d). Notably, among these aOutliers, several significant genes were identified (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ee\u0026ndash;g and Supplementary Fig.\u0026nbsp;7a\u0026ndash;f), including \u003cem\u003eSUGP1\u003c/em\u003e, known for its crucial role in mRNA splicing regulation in cancer\u003csup\u003e27,28\u003c/sup\u003e. In certain outlier individual(s), \u003cem\u003eEIF2A\u003c/em\u003e, \u003cem\u003eFLYWCH\u003c/em\u003e, \u003cem\u003eTP53RK\u003c/em\u003e, and \u003cem\u003eSUGP1\u003c/em\u003e exhibited increased usage of distal poly(A) sites, whereas genes such as \u003cem\u003eUNC5A, RAB31\u003c/em\u003e, and \u003cem\u003eLSS\u003c/em\u003e preferentially use proximal poly(A) sites. Additionally, genes like \u003cem\u003eCOL4A2\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eg), \u003cem\u003eADCY4\u003c/em\u003e, and \u003cem\u003eHMGCL\u003c/em\u003e (Supplementary Figs.\u0026nbsp;7g, h) were found to utilize intronic poly(A) sites in outlier individuals. Altogether, the single and multi-tissue aOutliers we identified represent the first comprehensive atlas of aberrant APA events across 49 human tissues.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eaOutliers represent a unique gene set with characteristics distinct from other molecular outliers\u003c/h2\u003e \u003cp\u003eTo determine the extent of sharing between aOutliers genes and those identified as expression outlier or splicing outlier genes (i.e., eOutliers and sOutliers, respectively), we conducted a comparative analysis using the same datasets. Remarkably, we found that 74.2% of multi-tissue aOutlier genes were not detected by analysis of multi-tissue eOutliers or sOutliers (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea and Supplementary Fig.\u0026nbsp;8a). For example, \u003cem\u003eTRIT1\u003c/em\u003e, a human tRNA isopentenyl transferase 1 gene, is an aOutlier-only gene that preferentially utilizes a distal poly(A) site in outlier individuals across multiple tissues (median Z-score\u0026thinsp;\u0026gt;\u0026thinsp;11) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb). This finding suggests that multi-tissue aOutliers represent a novel set of aberrant genes not detectable by traditional eOutlier and sOutlier analyses.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFurther comparisons between the genomic lengths of multi-tissue aOutliers and eOutliers disclosed that aOutlier genes have significantly longer 3\u0026prime; UTRs than eOutlier genes (one-sided Wilcoxon rank\u0026ndash;sum test, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1.4 \u0026times; 10\u003csup\u003e‒16\u003c/sup\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec and Supplementary Fig.\u0026nbsp;8b). In contrast, aOutlier genes have only slightly longer 5\u0026prime; UTRs than eOutliers (one-sided Wilcoxon rank\u0026ndash;sum test, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.004; Supplementary Fig.\u0026nbsp;8c), and no significant difference was observed in coding sequence length (two-sided Wilcoxon rank\u0026ndash;sum test, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.19). Furthermore, aOutlier genes have a lower GC-content (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed) in their 3\u0026prime; UTR regions (one-sided Wilcoxon rank\u0026ndash;sum test, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;6.8 \u0026times; 10\u003csup\u003e\u0026ndash;6\u003c/sup\u003e) than eOutlier genes. Gene ontology enrichment analysis\u003csup\u003e29\u003c/sup\u003e on multi-tissue aOutliers further highlighted specific biological processes and signaling pathways unique to these genes (Supplementary Fig.\u0026nbsp;9). Collectively, these data indicate that aOutliers comprise a distinct gene set with unique molecular and functional characteristics, thereby significantly distinguishing them from other types of molecular outliers.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eRVs are significantly enriched among APA outliers\u003c/h2\u003e \u003cp\u003eTo assess the impact of RVs (MAF\u0026thinsp;\u0026lt;\u0026thinsp;0.01) on aberrant APA usage, we computed odds ratios (ORs) for RVs located within varying proximity of the gene body (window size: 1 kb, 2 kb, or 10 kb) to multi-tissue aOutlier genes in outlier individuals compared to those in nonoutlier individuals. Our analysis revealed strong enrichment of nearby RVs in multi-tissue aOutliers (Supplementary Fig.\u0026nbsp;10a). Interestingly, we observed higher ORs for the enrichment of insertion and deletions (indels) than for single-nucleotide variants (SNVs) (Supplementary Fig.\u0026nbsp;10a, b). Furthermore, the degree of enrichment became more pronounced when we considered RVs located in closer proximity to the aOutlier genes or employed increased Z-score thresholds (Supplementary Fig.\u0026nbsp;10b, c).\u003c/p\u003e \u003cp\u003eTo gain further functional insights into aOutliers-associated RVs, we first determined the proportions of these RVs with functional category using Variant Effect Predictor (VEP)\u003csup\u003e30\u003c/sup\u003e. A higher proportion of aOutliers-associated RVs had function annotation than nonoutliers, increasing with higher Z-score thresholds (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee). The functional categories of aOutliers-associated RVs were largely distinct from those associated with eOutliers and sOutliers. For example, aOutliers-associated RVs are strongly enriched in the 3\u0026prime; UTR region (OR\u0026thinsp;=\u0026thinsp;4.6 and 10.1, respectively; Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ef and Supplementary Fig.\u0026nbsp;10d).\u003c/p\u003e \u003cp\u003eTo examine whether aOutliers-associated RVs are more likely to be deleterious and potentially pathogenic, we further employed Combined Annotation-Dependent Depletion (CADD) scores\u003csup\u003e31\u003c/sup\u003e to stratify RVs into three groups: (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) lowly deleterious, CADD score 0\u0026ndash;15; (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e) moderately deleterious, CADD score\u0026thinsp;\u0026ge;\u0026thinsp;15 but \u0026lt;\u0026thinsp;25; and (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e) highly deleterious, CADD score\u0026thinsp;\u0026ge;\u0026thinsp;25. Highly deleterious RVs showed significantly higher enrichment (20-fold increase for singletons and 11-fold increase for RVs with MAF\u0026thinsp;\u0026lt;\u0026thinsp;1%; Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eg) in aOutliers compared to moderately deleterious RVs (10-fold increase for singletons and 6-fold increase for RVs with MAF\u0026thinsp;\u0026lt;\u0026thinsp;1%) and lowly deleterious RVs (2-fold increase for singletons and RVs with MAF\u0026thinsp;\u0026lt;\u0026thinsp;1%). In total, we identified 179 rare SNVs with CADD scores\u0026thinsp;\u0026ge;\u0026thinsp;15 near 155 aOutlier genes (two-sided Fisher\u0026rsquo;s exact test, \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;5.2 \u0026times; 10\u003csup\u003e‒107\u003c/sup\u003e; Supplementary Table\u0026nbsp;1). In two examples, the rare SNV rs557639120 in \u003cem\u003eSUGP1\u003c/em\u003e (CADD score\u0026thinsp;=\u0026thinsp;18.4, MAF in GTEx\u0026thinsp;=\u0026thinsp;0.0056, and gnomAD\u0026thinsp;=\u0026thinsp;0.0033) leads to an increase in distal poly(A) site usage in its 3\u0026prime; UTR. Similarly, the rare SNV rs759305120 in \u003cem\u003eCOL4A2\u003c/em\u003e (CADD score\u0026thinsp;=\u0026thinsp;34, MAF in GTEx\u0026thinsp;=\u0026thinsp;0.0007 and gnomAD\u0026thinsp;=\u0026thinsp;0.000031) leads to preferential use of its intronic poly(A) site (Supplementary Table\u0026nbsp;1). We also identified 211 indels near 186 aOutlier genes (two-sided Fisher\u0026rsquo;s exact test, \u003cem\u003eP\u0026thinsp;=\u0026thinsp;1.9\u003c/em\u003e \u0026times; 10\u003csup\u003e‒16\u003c/sup\u003e; Supplementary Table\u0026nbsp;2), including 49 located in 3\u0026prime; UTR. For example, an indel variant (C\u0026thinsp;\u0026gt;\u0026thinsp;CAAAT, rs112906978) at the 3\u0026prime; UTR of \u003cem\u003eACSF3\u003c/em\u003e introduces a canonical \"AAUAAA\" motif near a poly(A) site, leading to three aOutliers (Supplementary Fig.\u0026nbsp;10e, f). Enrichment of RVs was also observed in single-tissue aOutliers across nearly all individual tissues (including SNVs and Indels) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eh and Supplementary Fig.\u0026nbsp;11). Considered collectively, our analyses reveal that a distinct class of RVs is significantly associated with aOutlier genes.\u003c/p\u003e \u003cp\u003e \u003cb\u003eRare APA variants frequently alter the 3\u0026prime; UTR PAS, 5\u0026prime; splice sites, and RNA binding proteins (RBPs) binding sites\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe next investigated the potential regulatory mechanisms of aOutliers-associated RVs on aberrant APA usage. We first focused on 3\u0026prime; UTR aOutliers-associated RVs and performed motif enrichment analysis to determine the prevalence of RVs altering 3\u0026prime;end processing. Our results show that 3\u0026prime; UTR aOutliers-associated RVs frequently alter polyadenylation signals (PAS) and AU-rich motifs, such as \"AWUAAA\" and \"AAUAAA\" (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea). Additionally, by using saturation mutagenesis data\u003csup\u003e32\u003c/sup\u003e, we found that RVs associated with aOutliers have a more significant impact on poly(A) site usage than RVs associated with nonoutliers (one-sided Wilcoxon rank\u0026ndash;sum test \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1.32 \u0026times; 10\u003csup\u003e‒23\u003c/sup\u003e; Supplementary Fig.\u0026nbsp;12a). Notably, we observed a significant proportion of large-effect RVs (fold change, LFC\u0026thinsp;\u0026gt;\u0026thinsp;1) associated with aOutliers compared to nonoutliers (50.3% \u003cem\u003evs.\u003c/em\u003e 6.6%; one-sided Wilcoxon rank\u0026ndash;sum test \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;6.1 \u0026times; 10\u003csup\u003e‒44\u003c/sup\u003e; Supplementary Fig.\u0026nbsp;12b), indicating their pronounced effects on 3\u0026prime; UTR APA. To further experimentally validate these findings, we selected four top-ranked 3\u0026prime; UTR aOutlier genes by median Z-score and utilized a minigene reporter system containing reference allele and alternative allele of four rare variants in selected genes, including \u003cem\u003eMKKS\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb), \u003cem\u003eSUGP1\u003c/em\u003e, \u003cem\u003eTP53RK\u003c/em\u003e, and \u003cem\u003eATP5F1E\u003c/em\u003e. In all four cases, we could detect significant changes in the poly(A) site usage, which agreed well with the predicted effects of these RVs (Figs.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec, d and Supplementary Fig.\u0026nbsp;13a, b).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFurther investigation into multi-tissue intronic aOutliers revealed a higher incidence of RVs at 5\u0026prime; splice donor sites than at acceptor sites (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ee). Compared to nonoutlier RVs, aOutlier RVs are 19 to 441 times more prevalent at donor sites, and up to 47 times more prevalent at acceptor sites. Specifically, aOutlier RVs are 441 times more prevalent in the \"D\u0026thinsp;+\u0026thinsp;1\" site and \u0026ldquo;D\u0026thinsp;+\u0026thinsp;4\u0026rdquo; site and 302 times more prevalent in the \"D\u0026thinsp;+\u0026thinsp;2\" site relative to the nonoutlier RVs. For example, RVs that alter the first nucleotide of the \"GT\" sequence in the intron of \u003cem\u003eCOL4A2\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eh) and the intron of \u003cem\u003eTXNRD2\u003c/em\u003e lead to significant intronic APA events in these genes (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ef and Supplementary Fig.\u0026nbsp;13c). We also found that RVs altering the last base of exon 11 in \u003cem\u003eADCY4\u003c/em\u003e and exon 4 in \u003cem\u003eHMGCL\u003c/em\u003e resulted in intronic APA events (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eg and Supplementary Fig.\u0026nbsp;7e, f). Based on these findings, we hypothesized that RVs affecting canonical donor sites drive intronic aOutliers. This hypothesis is also supported by our recent finding that mutations near the donor sites can promote IPA usage, potentially by blocking U1 small-nuclear RNP binding\u003csup\u003e33\u003c/sup\u003e. Predicting the strength of donor sites with MAXENT\u003csup\u003e34\u003c/sup\u003e showed a reduced strength of mutant donor sites compared to wild type (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eh, i). We then performed intronic APA minigene reporter assays for \u003cem\u003eTXNRD2\u003c/em\u003e and \u003cem\u003eCOL4A2\u003c/em\u003e with RVs at the conserved donor sites, as well as \u003cem\u003eHMGCL\u003c/em\u003e and \u003cem\u003eADCY4\u003c/em\u003e with RVs at the last base of the exons. For these assays, we cloned fragments containing full-length intronic sequences, including the donor sites, and upstream and downstream exons into the pcDNA3.1 vector. Results from 3\u0026prime; Rapid Amplification of cDNA Ends (3\u0026prime; RACE) assays indicate that all four RVs significantly increase alter IPA regulation relative to the wild-type sequence (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ej, k and Supplementary Fig.\u0026nbsp;13d, e).\u003c/p\u003e \u003cp\u003eLastly, we investigated whether aOutlier-associated RVs impact other transcriptional and posttranscriptional regulation of target genes. DeepBind\u003csup\u003e35\u003c/sup\u003e analysis of 927 binding motifs revealed 11 significantly enriched motifs in aOutlier-associated RVs (Supplementary Fig.\u0026nbsp;14a) using randomly shuffled RVs as control, including known APA regulator \u003cem\u003ePABPN1\u003c/em\u003e\u003csup\u003e36\u003c/sup\u003e. Furthermore, we analyzed 166 publicly accessible RBPs cross-linking immunoprecipitation sequencing (CLIP-seq) datasets from the Encyclopedia of DNA Elements (ENCODE) project\u003csup\u003e37\u003c/sup\u003e. We found seven RBPs's CLIP-seq data are strongly enriched with multi-tissue aOutlier RVs compared to nonoutlier RVs (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003el and Supplementary Fig.\u0026nbsp;14b), including \u003cem\u003eLARP4\u003c/em\u003e, an APA regulator identified in our previous study\u003csup\u003e25\u003c/sup\u003e, and a known APA regulator \u003cem\u003eCSTF2T\u003c/em\u003e. Knockdown of the two RBPs resulted in widespread APA dysregulation (Supplementary Fig.\u0026nbsp;14c, d), affecting two aOutlier genes, \u003cem\u003eSREBF2\u003c/em\u003e (Supplementary Fig.\u0026nbsp;14e) and \u003cem\u003eTOLLIP\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003em), in which the associated RVs were inside binding peaks of \u003cem\u003eLARP4\u003c/em\u003e (Supplementary Fig.\u0026nbsp;14f) and \u003cem\u003eCSTF2T\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003en), respectively. Beyond these known APA regulators, other RBPs such as \u003cem\u003eTIA1\u003c/em\u003e, \u003cem\u003eUPF1\u003c/em\u003e, and \u003cem\u003eSAFB2\u003c/em\u003e were also identified as potential new APA regulators (Supplementary Figs.\u0026nbsp;14g-i). Collectively, these results suggested that aOutlier-associated RVs trigger aberrant APA usage through altering PAS, splice sites, or RBP binding sites.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eInclusion of APA significantly improves functional RV effect prediction\u003c/h2\u003e \u003cp\u003eTo prioritize potentially impactful RVs for the interpretation of individual genomes, we repurposed the traditional Watershed\u003csup\u003e7\u003c/sup\u003e method into an APA-included version (aWatershed). This revised aWatershed model is an unsupervised probabilistic Bayesian hierarchical graphical model incorporating three RNA outlier signals, including aOutliers, eOutliers, and sOutliers, and annotations of a matched individual genome (Supplementary Table\u0026nbsp;3). The aWatershed model can allow us to quantify the posterior probability of an RV leading to a functional effect on APA usage (Supplementary Figs.\u0026nbsp;15a, b; Materials and Methods). To evaluate the aWatershed performance on the GTEx v8 data, we used held-out individual pairs with the same RVs as the evaluation dataset. By applying aWatershed prediction on the first individual of each pair and evaluating this prediction using the outlier status of the second individual as a label, we observed that our model significantly outperforms both the RIVER (RNA-informed variant effect on regulation) model\u003csup\u003e8\u003c/sup\u003e, a simplification of the Watershed model which integrates genomic features with aOutlier signals alone, and the GAM (genomic annotation model), a generalized logistic regression model based on genomic features alone (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea and Supplementary Fig.\u0026nbsp;15c). 93% of aWatershed prioritized RVs have low posterior probabilities in the GAM (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb), highlighting the importance of transcriptomic aOutlier signals in functional RVs prioritization. Moreover, aWatershed successfully captures the regulatory mechanisms underlying the effect of RVs on aOutlier signal (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ec). Strikingly, the integrated aWatershed model can prioritize RVs associated with 73.8% of aOutliers, in contrast to only 12.4% when relying on the genomic features alone (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ed).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eNext, we used the saturation mutagenesis data\u003csup\u003e32\u003c/sup\u003e to further evaluate the efficacy of aWatershed in prioritizing RVs with significant effects on APA regulation. In this analysis, we stratified RVs into two groups based on aWatershed APA posterior probabilities and compared poly(A) usage between them. We found that RVs in the group with high posterior probability had significantly larger effects on APA than those in the low posterior probabilities group (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ee, f), suggesting our aWatershed model is effective in identifying RVs with substantial APA effects. Furthermore, our analysis revealed that aWatershed successfully identified many functional RVs overlooked by the previous variant prediction model\u003csup\u003e38\u003c/sup\u003e, as exemplified by two RVs in \u003cem\u003eRPL13A\u003c/em\u003e and \u003cem\u003ePAAF1\u003c/em\u003e, respectively (Supplementary Fig.\u0026nbsp;15d). Overall, aWatershed prioritized 1,799 RVs predicted to impact 278 APA genes (Supplementary Table\u0026nbsp;4). Interestingly, there was minimal overlap between RVs impacting APA and those affecting gene expression or splicing, as only 60 of these 1,799 RVs were common to those categories. For example, the RV rs191575428 within the 3\u0026prime; UTR of \u003cem\u003eMTHFD2\u003c/em\u003e, which exhibited a high aWatershed APA posterior probability of 0.997 based on aOutliers, showed considerably lower posterior probabilities for expression and splicing (0.055 and 0.008, respectively). This variant is associated with 3\u0026prime; UTR lengthening in outlier individuals without changing gene expression levels (Supplementary Fig.\u0026nbsp;15e, f). Further extending the aWatershed model to prioritize tissue-specific functional RVs by integrating genomic features with single-tissue aOutliers signals, we observed that the tissue-aWatershed model outperforms both the tissue-RIVER model and tissue-GAM model (Supplementary Figs.\u0026nbsp;16\u0026ndash;17). In summary, by leveraging these aOutliers, we have implemented a robust Bayesian hierarchical variant effect prediction model aWatershed that effectively prioritizes rare functional variants with significant effects on APA regulation.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eAnalysis of aOutliers prioritizes RVs impacting complex traits and diseases\u003c/h2\u003e \u003cp\u003eTo test the hypothesize that aWatershed RVs could be used to interpret the complex traits and diseases, we first examined the 278 genes prioritized by aWatershed and cross-referenced with genes annotated in the Online Mendelian Inheritance in Man (OMIM) database\u003csup\u003e39\u003c/sup\u003e. We identified 21.2% of the prioritized genes were well-known disease genes (Supplementary Fig.\u0026nbsp;18a). For example, we identified a prioritized RV, rs79940214, associated with \u003cem\u003eMKKS\u003c/em\u003e (Supplementary Fig.\u0026nbsp;18b), which encoded a centrosome-shuttling protein and was associated with many genetic diseases, including McKusick-Kaufman syndrome (OMIM id: 236770)\u003csup\u003e40,41\u003c/sup\u003e and Bardet-Biedl syndrome 6 (OMIM id: 605231)\u003csup\u003e42,43\u003c/sup\u003e. Another example is one prioritized intronic RV, rs76984877, that is associated with gene \u003cem\u003eEXT2\u003c/em\u003e (Supplementary Fig.\u0026nbsp;18b), which was associated with hereditary multiple exostosis, type 2\u003csup\u003e44,45\u003c/sup\u003e. We also identified five prioritized RVs associated with gene \u003cem\u003eBCR\u003c/em\u003e (Supplementary Fig.\u0026nbsp;18b), which has been frequently reported to be associated with chronic myeloid leukemia\u003csup\u003e46,47\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eWe further cross-referenced aWatershed-prioritized RVs with trait variants from 1,234 well-powered GWAS summary statistics from UK Biobank (UKBB) and literature (Supplementary Table\u0026nbsp;5), resulting in 1,385 RVs associated with 171 aOutlier genes in 1,186 traits (Supplementary Table\u0026nbsp;6). We focused on the subset of 623 traits, which also have evidence of colocalization with 3\u0026prime;aQTLs (Supplementary Table\u0026nbsp;7). Notably, aOutlier prioritized RVs fell in or nearby genes had evidence of colocalization with 3\u0026prime;aQTLs having larger trait effect size than the non-colocalized RVs (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.0014, one-sided Wilcoxon rank\u0026ndash;sum test; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ea). We also conducted a permutation test to determine whether these prioritized RVs exhibit larger effect sizes on these complex traits and diseases. We found that the mean odds of aOutlier-prioritized RVs had a more significant effect size than non-prioritized RVs (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;2.5 \u0026times; 10\u003csup\u003e‒15\u003c/sup\u003e, one-sided and paired Wilcoxon rank\u0026ndash;sum test; Supplementary Fig.\u0026nbsp;19a). To exemplify the larger effect size in aOutlier-prioritized RVs, we focused on two traits: height related traits (UKBB trait ID: 50_irnt and 20015_irnt) and high blood pressure (UKBB trait ID: 6150_4). This analysis revealed a significant shift in the odds favoring RVs with higher aWatershed posterior probabilities over those with lower ones (\u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1.6 \u0026times; 10\u003csup\u003e‒9\u003c/sup\u003e and \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;2.3 \u0026times; 10\u003csup\u003e‒54\u003c/sup\u003e, respectively; one-sided Wilcoxon rank\u0026ndash;sum test; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eb, c; Supplementary Fig.\u0026nbsp;19b). In the case of height related traits and high blood pressure, these aOutlier prioritized RVs had larger effect sizes on the trait than other variants within a 1Mb of the RV, including RVs prioritized by eOutliers or sOutliers. Notably, for the height related traits, the RV (rs112567314), located in the intron of \u003cem\u003eCUL3\u003c/em\u003e, had a greater effect size than other variants within 1 Mb and RVs prioritized by eOutlier or sOutlier (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ed and Supplementary Table\u0026nbsp;6). Similarly, for high blood pressure, the RV (rs893929), located in the intron of \u003cem\u003eUSP38\u003c/em\u003e, also had a greater effect size than 99.6% of variants within 1 Mb, including the nearest trait-associated significant variants as well as eOutlier or sOutlier RVs (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003ee). Collectively, our results demonstrate the capability of aWatershed in prioritizing RVs with large effect sizes on APA, significantly impacting complex traits and diseases.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eStrong convergence between rare and common variants on\u003c/b\u003e \u003cb\u003eDDX18\u003c/b\u003e \u003cb\u003elinks APA regulation with cancer susceptibility\u003c/b\u003e\u003c/p\u003e \u003cp\u003eEmerging evidence suggests potential interactions between rare and common variants in affecting the same disease genes\u003csup\u003e48\u0026ndash;50\u003c/sup\u003e. As expected, we also observed the strong convergence effect on 3' UTR APA regulation between RVs and common variants (Supplementary Fig.\u0026nbsp;20). To further mechanistically examine their convergence effects on disease, we focused on aWatershed prioritized RVs and their associated genes. We found 126 out of the 278 aOutlier RV associated genes were also identified as susceptibility to disease risks, including cancer risks through 3\u0026prime;aQTLs in our gene-based association studies\u003csup\u003e51,52\u003c/sup\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea). Among the top-ranked APA genes that were prioritized by both RV and 3\u0026prime;aQTLs analyses (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb), we noticed several highly constrained genes (pLI score\u0026thinsp;\u0026gt;\u0026thinsp;0.9), and we particularly focused on the gene \u003cem\u003eDDX18\u003c/em\u003e, a member of the DEAD-box RNA helicase family, that was identified as an APA-mediated susceptibility gene across many cancer types\u003csup\u003e53,54\u003c/sup\u003e. Moreover, CRISPR-Cas9 based gene essentiality screens also demonstrated that \u003cem\u003eDDX18\u003c/em\u003e has an essential role in cancer cell proliferation\u003csup\u003e55,56\u003c/sup\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). Examining our 3\u0026prime;aQTLs data revealed significant associations between common variants and 3\u0026prime; UTR APA of \u003cem\u003eDDX18\u003c/em\u003e across tissues, with the most significant one was found near the 3\u0026prime; end (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed and Supplementary Fig.\u0026nbsp;21a-c). Intriguingly, an aWatershed prioritized RV, rs1680042046, located near the distal poly(A) site of \u003cem\u003eDDX18\u003c/em\u003e, was identified in the outlier individual (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ee, f). This RV alters the hexamer motif \"AUUAAA\" to \"AUUAAG\" (Supplementary Fig.\u0026nbsp;21b) and has a highly deleterious effect (CADD\u0026thinsp;=\u0026thinsp;17.5) (Supplementary Table\u0026nbsp;1).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo further experimentally validate the convergence effect of RVs and common variants on \u003cem\u003eDDX18\u003c/em\u003e, we designed minigenes introducing the APA variants by PCR-based site-directed mutagenesis in HEK293T and MCF7 cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eg). We then performed 3\u0026prime; RACE to quantitatively evaluate the effect of the common variant (rs1052628; A\u0026thinsp;\u0026gt;\u0026thinsp;G) alone, the RV (rs1680042046) alone, or their joint effect on APA. We first mutate the reference A allele to the alternative G allele for either RV or the common variants. In HEK293T cells, this mutation decreased the use of the distal poly(A) site (dPAS) for both the common variant and RV (two-sided Student\u0026rsquo;s t-test \u003cem\u003eP\u003c/em\u003e\u0026thinsp;=\u0026thinsp;1.5 \u0026times; 10\u003csup\u003e‒6\u003c/sup\u003e and 2.3 \u0026times; 10\u003csup\u003e‒7\u003c/sup\u003e; Figs.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eh, i), indicating that both variants indeed trigger \u003cem\u003eDDX18\u003c/em\u003e APA regulation. A similar APA effect was also observed in MCF7 cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ej). To further assess the functional roles of \u003cem\u003eDDX18\u003c/em\u003e APA regulation in breast cancer cells, we measured \u003cem\u003eDDX18\u003c/em\u003e protein level using luciferase reporter assays and assessed the effect of gene silencing on the proliferation of MCF7 cells proliferation. We observed lower luciferase activities in the short 3\u0026prime; UTR isoform of \u003cem\u003eDDX18\u003c/em\u003e and the reporter containing RV or both RV and common variant (Supplementary Fig.\u0026nbsp;21d, e). Knockdown of \u003cem\u003eDDX18\u003c/em\u003e in MCF7 results in inhibition of cell proliferation (Supplementary Fig.\u0026nbsp;21f, g). Collectively, these findings highlight the critical role of rare variants in understanding the risk of common diseases and offer a novel approach to linking functional rare variants to complex diseases.\u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eThe human genome contains a plethora of rare genetic variants whose functional effects and underlying molecular mechanisms are challenging to interpret. In this study, we introduce the aOutlier as an emerging molecular phenotype reflecting aberrant 3\u0026prime; UTR or intronic APA usage across multiple samples. aOutlier can be used to identify functional rare APA variants. By analyzing population-scale transcriptomics data using our DaPars2\u003csup\u003e25,26\u003c/sup\u003e and IPAfinder algorithms\u003csup\u003e18\u003c/sup\u003e, we identified 1,534 multi-tissue aOutliers based on European individuals. These aOutlier genes exhibit unique molecular features, such as genomic lengths and GC-content, setting them apart from other molecular outliers, such as eOutliers and sOutliers. Importantly, aOutliers can aid in identifying a distinct class of rare functional variants. We observed that aOutliers-associated RVs are more likely to be deleterious and are highly enriched in outlier individuals. Mechanistically, these aOutlier-associated RVs can modulate APA usage by either altering PAS, AU-rich elements, or splice donor sites, as confirmed by saturation mutagenesis data and 3\u0026prime; RACE experiments.\u003c/p\u003e \u003cp\u003eTo further enhance the utility of our aOutlier atlas, we adapted a Bayesian hierarchical prediction model (aWatershed) by incorporating genomic features with multiple functional signals, including aOutliers, eOutliers, and sOutliers. This integration aims to predict the probability of RV leading to aberrant APA usage. Notably, our aWatershed model outperformed models trained only on genomic features or those combined with aOutlier signals alone. Moreover, aWatershed-prioritized RVs exhibited more significant effects on APA regulation than non-prioritized RVs. The predictive power of aWatershed was validated using GWAS summary data from the UKBB, showing that aWatershed-prioritized RVs had larger trait effect sizes than non-prioritized RVs, as exemplified by RVs near \u003cem\u003ePOLR2L\u003c/em\u003e and \u003cem\u003eATP5F1D\u003c/em\u003e associated with height and BMI, respectively.\u003c/p\u003e \u003cp\u003eInterestingly, we observed a significant proportion of intersection between aOutlier transcripts and 3\u0026prime;aQTL associated transcripts in matched tissue, suggesting the potential interplay of common variants and RV in APA regulation, similar to previous findings in gene expression studies\u003csup\u003e48,49,57,58\u003c/sup\u003e. Additionally, a rare deletion 16p11.2 and common variants in chromosome 16p modulate downstream gene expression and affect the risk for autism\u003csup\u003e48\u003c/sup\u003e. Moreover, using minigene reporters and 3\u0026prime; RACE assays, we demonstrated the potential additive effect of rare and common APA variants on \u003cem\u003eDDX18\u003c/em\u003e 3\u0026prime; UTR regulation. We further demonstrated that the regulation of \u003cem\u003eDDX18\u003c/em\u003e 3\u0026prime; UTR contributes to \u003cem\u003eDDX18\u003c/em\u003e protein expression level, which is tightly linked to breast cancer cell proliferation. In summary, our study identifies a novel set of rare functional variants that influence APA and connects these RVs to human trait phenotypes, providing valuable information for the identification of novel genes associated with increased disease risk.\u003c/p\u003e"},{"header":"Materials and Methods","content":"\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003eGTEx data collection and processing\u003c/h2\u003e \u003cp\u003eWe downloaded both RNA-seq raw sequencing data and whole-genome genotype data of the v8 release of the GTEx project from dbGAP (accession: phs000424.v7.p2). Expression outlier (eOutlier) and splicing outlier (sOutlier) data, and the metadata of samples (filename: GTEx_Analysis_v8_Annotations_SampleAttributesDD.xlsx) and subjects (filename: GTEx_Analysis_v8_Annotations_SubjectPhenotypesDD.xlsx) were downloaded from GTEx Portal (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://gtexportal.org/home/\u003c/span\u003e\u003cspan address=\"https://gtexportal.org/home/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). The GTEx RNA-seq dataset contains 17,832 samples representing 54 biological tissues collected from 838 donors. For this study, we included data from 49 tissues, each with at least 70 samples. Original GTEx RNA-seq reads were aligned with the human genome (hg38/GRCh38) using STAR v.2.7.3a\u003csup\u003e59\u003c/sup\u003e, with the following alignment parameters: outSAMtype, BAM; SortedByCoordinate; outSAMstrandField, intronMotif; outFilterMultimapNmax, 10; outFilterMultimapScoreRange, 1; alignSJDBoverhangMin, 1; sjdbScore, 2; alignIntronMin, 20; and alignSJoverhangMin, 8. The aligned BAM files were sorted and further converted to bedGraph format using BEDTools v.2.27.1\u003csup\u003e60\u003c/sup\u003e. The genotype data in VCF format (filename: GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.vcf.gz) was processed with vcftools v.0.1.13 to calculate MAF across all subjects and extract allele information for each variant.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003e3\u0026prime; UTR APA and intronic APA quantification\u003c/h2\u003e \u003cp\u003eTo quantify the 3\u0026prime; UTR APA, we analyzed alignment files in BAM format using the DaPars2 algorithm. We followed the workflow implemented in our 3\u0026prime;aQTL analysis\u003csup\u003e25,61\u003c/sup\u003e. Briefly, the BAM files were firstly transformed to bedGraph format with a bin size of 1, which records the read coverage of each position in the genome. Before analyzing APA, we downloaded the gene annotation file containing all transcripts of genome build hg38 in RefSeq database through the UCSC Genome Browser, from which we extracted 3\u0026prime; UTR region of each transcript using script \"DaPars_Extract_Anno.py\". The DaPars2 algorithm then detects the proximal poly(A) site in the 3\u0026prime; UTR region of each transcript and calculates the relative usage of the distal poly(A) site by the script \u0026ldquo;Dapars2_Multi_Sample.py\u0026rdquo; for all samples in each of the 49 tissues. This is indicated as the Percent of Distal Poly (A) site Usage Index (PDUI) only if a proximal poly(A) site is detected. For intronic APA detection and quantification, we used IPAfinder\u003csup\u003e18\u003c/sup\u003e, which is a python-based tool that enables \u003cem\u003ede novo\u003c/em\u003e identification and quantification of intronic APA (IPA) events using RNA-seq data. IPAfinder can identify potential IPA sites and calculate the Intronic poly(A) site Usage Index (IPUI), which represents the proportion of total transcripts that are intronic-polyadenylated for each intronic APA event\u003csup\u003e18,33,62\u003c/sup\u003e. BAM files were analyzed together by IPAfinder and separated by tissues.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eCovariate correction and normalization\u003c/h2\u003e \u003cp\u003eTo avoid batch effects and unobserved confounders in each tissue, we adjusted the sample genotype and APA usages with known covariates, such as population structure, sex, and sequencing platform. Briefly, for genotype data, we first removed sites marked as \"wasSplit\" from the GTEx analysis freeze variant call format (VCF) using BCFtools v.1.10.2. We further applied the PEER model\u003csup\u003e63\u003c/sup\u003e with sex, age, sequencing platform, and the top five genotype principal components as known covariates to estimate a set of latent covariates for PDUI or IPUI values in each tissue. The number of PEER factors was optimized based on tissue sample size, as suggested by the GTEx Consortium; 15 PEER factors were chosen for sample sizes\u0026thinsp;\u0026lt;\u0026thinsp;150, 30 PEER factors were selected for sample sizes ranging from 150 to 250, and 35 peer factors were chosen for sample sizes\u0026thinsp;\u0026gt;\u0026thinsp;250. Before running the PEER model for inferring hidden covariates, PDUI/IPUI values in each tissue were quantile normalized to remove batch effects.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eAPA outlier calling\u003c/h2\u003e \u003cp\u003eAfter inferring the hidden covariates for each tissue, we calculated PDUI/IPUI residuals by regressing out inferred PEER factors and known covariates, including population structure, sex, and sequencing platform, using the function \"PEER_getResiduals\". In each individual tissue, we obtained normal-distributed \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(Z\\left(g,t\\right)\\)\u003c/span\u003e\u003c/span\u003e score for each gene (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(g\\)\u003c/span\u003e\u003c/span\u003e) in the tissue (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(t\\)\u003c/span\u003e\u003c/span\u003e) by scaling the PDUI/IPUI residuals across samples with the following equation, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({X}_{r}\\left(g,t\\right)\\)\u003c/span\u003e\u003c/span\u003e denotes the residuals of PDUI/IPUI values, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\stackrel{-}{{X}_{r}\\left(g,t\\right)}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(sd\\left({X}_{r}\\left(g,t\\right)\\right)\\)\u003c/span\u003e\u003c/span\u003e represent the mean and standard deviation of the residuals across samples, respectively:\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$Z\\left(g,t\\right)=\\frac{{X}_{r}\\left(g,t\\right)-\\stackrel{-}{{X}_{r}\\left(g,t\\right)}}{sd\\left({X}_{r}\\left(g,t\\right)\\right)}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWe defined two types of aOutliers in the current study. One is single-tissue aOutlier, which is called from a single tissue based on the Z-score of each gene in that tissue. When the absolute Z-score of an individual exceeds a threshold of three for a gene, then the individual is called a single-tissue aOutlier for that gene. The other is multi-tissue aOutlier, for which we calculated a median Z-score\u003csup\u003e7,8\u003c/sup\u003e for each APA event across all tissues when data were available and restricted our analysis to individuals with APA measurements in at least five tissues. Multi-tissue aOutliers were defined as those with an absolute median Z-score\u0026thinsp;\u0026gt;\u0026thinsp;3. The same threshold was used for eOutlier and sOutlier calling\u003csup\u003e7,8\u003c/sup\u003e. Our method allowed that one gene could have multiple aOutlier individuals, and one individual could also be aOutliers of multiple genes. To account for situations in which widespread aberrant APA might occur in an individual due to non-genetic influences, we removed 11 individuals in which the proportion of tested genes identified as multi-tissue aOutliers exceeded 1.5 times the interquartile range of the distribution for aOutlier gene proportion across all individuals. These 11 individuals were marked as global outliers.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eEstimation of replication rates of aOutliers\u003c/h2\u003e \u003cp\u003eTo estimate the replication rate of aOutliers between different tissues, we selected one of the 49 GTEx tissues each as discovery tissue, and compared aOutliers detected in it with those of the other tissues. For replication rate calculation, we only considered the shared aOutlier genes in both tissues and an aOutlier to be replicated only when the gene and individual of the aOutlier matched between the compared tissue pairs. For multi-tissue aOutliers replication, we used the cross-validation method described in a previous study\u003csup\u003e8\u003c/sup\u003e to estimate their replication rate. In brief, the 49 human tissues were separated into two groups, one group with 39 tissues as the discovery group, the other group has the remaining ten tissues as the replication group. Each time we randomly sampled \u003cem\u003et\u003c/em\u003e (\u003cem\u003et\u003c/em\u003e\u0026thinsp;=\u0026thinsp;10, 15, 20, 25, 30) tissues from the discovery group and called multi-tissue aOutliers in the discovery group using a Z-score threshold of 3 in at least five tissues as described above. Then we estimated the replication rate as the proportion of multi-tissue aOutliers in the discovery group with an absolute median Z-score 2 or 3 in the replication group. We also computed the expected replication rate by randomly selecting individuals in the discovery group with at least five tissues that have APA usage for the gene and determined the replication status in the replication group. For each discovery group size (\u003cem\u003et\u003c/em\u003e), we repeated this process 10 times.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eRV annotation\u003c/h2\u003e \u003cp\u003eWe defined RVs as those with MAF\u0026thinsp;\u0026lt;\u0026thinsp;1% within the GTEx European individuals and with MAF\u0026thinsp;\u0026lt;\u0026thinsp;1% in non-Finnish Europeans within gnomAD\u003csup\u003e64\u003c/sup\u003e. Singletons were defined as RVs with minor allele only presents once in GTEx European individuals and were extracted using vcftools. The annotation of RVs was performed by Ensembl VEP (release 104), which assigned 36 different annotation terms to each RV, including protein-coding gene position (e.g., \"splice_donor\" \"splice_acceptor,\" \"frameshift\u0026rdquo;) and regulatory regions (e.g., \"TFBS_ablation\", \"TF_binding_site\"). Annotation terms were grouped into one of the four classes based on predicted impact: \"High\", \"MODERATE\", \"MODIFIER\" and \"LOW\". The high-impact one was used for variants assigned with two or more annotations. In addition to 36 VEP annotations, we added two other annotations to each RV; \"PAS region\" describes variants located within 50 bp upstream of the annotated PAS, and \"PAS signal\" refers to variants located at the PAS motif \"AAUAAA\" and its additional 14 variants (\"AUUAAA\", \"UAUAAA\", \"AGUAAA\", \"AAAAAA\", \"AACAAA\", \"AAGAAA\", \"AAUAUA\", \"AAUACA\", \"CAUAAA\", \"UUUAAA\", \"ACUAAA\", \"AAUAGA\" and \"GAUAAA\"). We also used genomic annotations of variants extracted from CADD v.1.5 release (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://cadd.gs.washington.edu/download\u003c/span\u003e\u003cspan address=\"http://cadd.gs.washington.edu/download\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eRV enrichment analysis\u003c/h2\u003e \u003cp\u003eWe examine the enrichment of Rare Variants (RVs), including single-nucleotide variants (SNV) and small insertion and deletion (indel) near aOutlier genes. Only genes with at least one aOutlier individual were considered, and the remaining individuals for the same genes were treated as nonoutlier controls. We first counted the RVs present within 1kb, 2kb, or 10kb of the outlier genes in both outlier and control individuals and built a 2 \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\times\\)\u003c/span\u003e\u003c/span\u003e 2 contingency table for each of the flanking region, containing the number of aOutliers with RVs, the number of nonoutlier controls with RVs, the number of aOutliers without RVs, and the number of controls without RVs. We then calculated Odds Ratios (ORs), \u003cem\u003eP\u003c/em\u003e value, and 95% confidence interval (CI) using Fisher\u0026rsquo;s exact test in R base package. We grouped variants into four groups based on their MAF (0\u0026ndash;1%, 1\u0026ndash;5%, 5\u0026ndash;10%, and 10\u0026ndash;25%), and performed enrichment analysis for each group. We also conducted enrichment for RVs that stratified by VEP annotations and CADD scores as described above.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eEnrichment analysis for RVs that influence PAS and AU-rich motifs\u003c/h2\u003e \u003cp\u003eTo identify potential regulatory variants associated with aberrant APA events, we defined RVs located within the gene body or in the 10-kb region surrounding outlier genes in outlier individuals as aOutlier RVs. Those in nonoutlier individuals in the same region were classified as nonoutlier RVs. For each aOutlier RV and nonoutlier RV located in the 50-bp region upstream (PAS region) of the poly(A) sites annotated in PolyA_DB V.3.2\u003csup\u003e65,66\u003c/sup\u003e, we extracted its upstream and downstream 5 base pairs sequences and examined whether it matched with one of the 15 known PAS motifs by using script \"dna-pattern\" in RSAT (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/rsa-tools\u003c/span\u003e\u003cspan address=\"https://github.com/rsa-tools\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). We then summarized all tested RVs and conducted PAS motif enrichment analysis using Fisher's exact test, which determined the odds ratios (ORs) and 95% confidence intervals (CIs) for each PAS motif. To perform enrichment analysis at the 12 known AU-rich motifs, we restricted RVs to those within the 100 bp flanking the annotated poly(A) sites. We then counted RV enrichment analysis for each of the AU-rich motifs using the same method as for PAS motifs.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eIdentification of aOutlier RVs enriched RNA motifs\u003c/h2\u003e \u003cp\u003eWe focused on multi-tissue aOutlier associated RVs located within the gene body region, which spans from 3 kb downstream of the transcription start site (TSS) to the end of the gene. We extracted the 3 base pairs of sequences flanking each RV from both sides. Next, we used DeepBind v0.11\u003csup\u003e35\u003c/sup\u003e to score these 7-mer sequences using 617 pre-built models, including 515 transcription factors and 102 RNA-binding proteins from Homo sapiens. For each 7-mer sequence, we selected the top three motifs with a DeepBind score of at least 0.1. To validate the enrichment of RVs in predicted binding motifs, we created a control set of RVs by randomly shuffling the genomic locations of multi-tissue aOutlier associated RVs within the same gene body regions. We used Fisher's exact test to estimate the level of enrichment.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eIdentification of aOutlier RVs enriched RBPs\u003c/h2\u003e \u003cp\u003eWe obtained CLIP-seq data for 166 RNA-binding proteins (RBPs) from the Encyclopedia of DNA Elements (ENCODE) data portal for HepG2 and K562 cells. We only considered significant binding peaks with \u003cem\u003eP-values\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.01, shared by two biological replicates for each RBP. To assess the enrichment of aOutlier RVs in RBP binding peaks, we selected RVs associated with multi-tissue aOutliers within gene body regions representing the region of 3 kb downstream of the transcription start site (TSS) to the end of the gene. We created a control RV set by randomly shuffling the genomic locations of multi-tissue aOutlier associated RV set within the same gene body regions. We then counted the RVs in binding peaks of each RBP using bedtools. Finally, we compared the RVs between the two sets using Fisher's exact test to determine the enrichment.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eDevelopment of a Bayesian prediction model that integrates APA outlier signals\u003c/h2\u003e \u003cp\u003eTo prioritize rare functional variants with significant impact, we improved the Watershed Bayesian hierarchical model by incorporating APA outlier signals with other layers of transcriptomic outlier signals and genomic annotations. The improved model called aWatershed, includes a layer of genomic annotation features (\u003cem\u003eG\u003c/em\u003e) which denotes the 40 observed genomic features aggregated over all RVs in the outlier individual that are within 10 kb region of the gene, a fully connected conditional random field (CRF) layer (\u003cem\u003eZ\u003c/em\u003e) represent the unobserved regulatory variables for each of the three transcriptomic outlier signals (APA, mRNA expression, and splicing), and a layer of variables (\u003cem\u003eE\u003c/em\u003e) representing the observed outlier status of each transcriptomic data type. The three layers were linked by the following conditional distributions:\u003c/div\u003e\n\u003cp\u003e\u003cimg src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbgAAADOCAYAAABMxhfbAAAgAElEQVR4Ae2djZUUN9NGcQg4BQgBx2AywClACCYBHIEdASQACZgAcAAmAhMD+527/q7fopC6e/52dmYfnTNH3VKpfp6SqiT1Yj+6SQkCQSAIBIEgcIUIPLpCm2JSEAgCQSAIBIGbJLhMgiAQBIJAELhKBJLgrtKtMSoIBIEgEASS4DIHgkAQCAJB4CoRSIK7SrfGqCAQBIJAEEiCyxwIAkEgCASBq0QgCe4q3RqjgkAQCAJBIAkucyAIBIEgEASuEoEkuKt0a4wKAkEgCASBJLjMgXuBwNevX2/4pQSBIPAtAlkX3+Kxy1sS3C5ohfYkCPz99983L1++vPnll19Owj9Mg8ClIvDx48ebn3766eb9+/fZAO7hxCS4PUDLkOMhwMJ9/PjxzevXr2++fPlyPMbhFASuBAGS3I8//ni7CcxpbjenJsHthleoj4jAp0+fbh49enTz22+/HZFrWAWB60Pgr7/++i/JXZ91p7MoCe6I2GZ3tRuYT548ufn55593G3REavxVf0dkfWesqv6XPv+uyZZTTIB3797dbgi59UjZhkAS3ASnvthG73Uo1whbTyMjXrSdsmyROaOxfaaf/Ut1H/v27dtbvMDtHOWff/65vfLh6ueHH3643R2/efPmHKrsLZNA9+zZs1v9sYHNAqfiSyvMG07xT58+/c8Wvsnio4dQ+rpZshl/szFM2YZAEtwEJxYYCYvAUX+02c61gWVrgiMAwdvACm8mLG2nKn/++efNixcvvpHJQukynz9/fmub9mm3+o0CDoGpYuIYa/q6HOyE57kWKnZgPx/vP3/+fAu7u+NLuS51g/DHH3/c6s/3S2zie+alfct0PbjZYV1hxzlP96dai51v36T8+uuvneSbd+epWH3TmZfvEEiC+w6SfxsIvjUASsZiJGj3QLglwf3++++3Y/lrQQMr41jMpwr2nErQt8ukrQcQAiN60Kd+tBFER/RiQvKkn8VXi5gYhO0jgEG/tpilP3aN3WCujfJfslGa+1CLa8fPzQb9l1LU+cOHD9+ojI/wxzUX1iYbQdeH66LHlooBcxZcRpvGSpfnfxG47hm0p5eZRJywegB018xJpxeDzmxyOvbVq1d96O1fEI7avyPcscGEOuKNDSNdCfz8ejHg1FOrNCa4jhf9LOA+xqB2jkBsEOnJAV0vJcGJdz+pnRNX58KuNXONk2cvzrfefi3vnNyYbyY37cJufkuFTehojS6Neah9SXAbPc/VIpOKxNcDCyyWEpy7Lk6Ed1WUyXeNrcUxBNBeDDijpMSC2+UEKq8Rjl0u7/UbRe2nfddCYiOw9KTLNS7to8S3q4xT0oMZes58xIZiK66n1HMLb05toyDPFTJr7S7XyxZ9j0nDehklMtrWkpcbnNGG8pg6XgOvJLgNXiSQssskeIwCPCyWEpzXmv0aZoPovUmU2a8NlxgacEYnO+wnGPXgacDtVyZLyYfFDa+1Ag9OvmJPsvYvyDyt7GIf8uDVA4jf5GablzU977JfH/WdPzgxP0e+u0v9dpE12mzgc+bS0lrbRcZ9pPU72mjusi6Yo0vFuT+LRUtjH1rfepR5aIgM7HUhLgWPWYIzAexywqkqsOD91falZ2X2QL40hj7t7AvH083oqtOAS4BVT65Gl2RvWcToY6AzmBMQ4OsOf3RVvGYjsj39oK8f+Tkt9FPdGq9z9Bvc1JX/Cgxtl5bcwK6eVvAFtyT4ho3GXW4G79qP2Mg8dr1Y40vm5+hkV3V0DizFo0r/kJ+T4Fa871352nXJLMG5W9vn6osFj1yCF5Oe91EhuRjw6FdmP1WNxtY2T2l1wcGLgEMy6ac3xpoU0dFfTSKVv8+7LGKTWx2LXfuctrqP/LN0EvKlFOYB+FHwDZjjN4LjpZU6D9g8YQu+ZQNzrcXNp2ul1uDBby1W9Hl8rVgdw64kuAUUCRrstAima/fds0nHgmXSjnakJhLrqgryCMAkLhYF/ykrFoMnJWkJBtBV/ZRJANxaXHjo6qLjmeC5xIf+elqDD+89MVU94Lu0S5XH6NTLWH5LOlVZ9dmdL76iIAc94QeuSwUfjYq+q/WIrrZJW9u2PnfsmB9sgvBDTwzK6XWXZf/W9krH2H3KaL2wRlhrzN9eqpyZvn1MfXdM5WP/qM2+Y9cju5Wxdd0u8ZBX6n8RSIKbzAQmPUGDgDJKTn3YbNK5464nLMYabOHPr1//Mdl7EOedAMAVB8mXgAb/TqdMdNpasBE9DC4Ef/iQrGryrPygYQx0tXDa6/bW/tGY2o890IyuYNbGVj79WVx6u4FlZicnZ+T2pK3P6XNTQFJe+gfjYAaN16Rdl6V35XVc0Bsd9J080KnqZvKoiVCe0FW+Jn7aZ/Nohovyl2rnf+et7+ucRl9sYww/Enqfc0uy6OOfyWALa7oW9ahtp3wW1243Mllr6MgcWSr6DN1TlhFIgpvg49UbJ6daZru92aSbBVV4GlhHwa4HU3UgcbBYCV6c3GogkGZJpjS91t7Kz6Q3uzKxf9eFxiJeClAGnR4ExLi3d1tm7wSQHuCgncmTD7oydmQn2NBnkdcswdMP/ZL98uq1wXG04Rphqn/UBdyQ3RMhNpB0+4kZHTtt1WkJl0o3embOj4K5Pq5YQ1vnIPasfTIYycT33R7W0j7fckf8t7Q5P/oc1lejWND5jjDqNHn/F4EkuMFMWPruRoAY7fRnk44gwEIeFRNcn+wj2l3almSSoEdJ2jHVNk9oo6SAPrPFuqYreMx4MlZdOi6z9jV59BMUkdsDHH1LfuD7HAEQ2TXoKpOA1IMScka0nHg4RRGs4bdrQQ68q4/g4Qmu8zQhVjnQdDr46kuTIWOg6z6Q1xou0s1q1lFPqNB6gqv4zfCE3vnsnO7vVT4y6wZO2krjs33WtlvbTr1LEec+Bqyxc4Z3pZdHxaj25/l/CIwj7//6H9zT0nc3rp5YJKPCxBwtRIMnf4lYCwsDXjN+lXbX5yWZ9PGrxUQ2Cji0YdeokKToW7tS6WOXeEKr/nWxE6wN8LW98jbo1DafDZw1wNHH9Rc+IPH04lwggYwSA/SMrYHGa7uRjuCFfANUl7ekP7TgNvKR/7WafuoHr5rM4M/47n/nIHXtm/l2Cy5LtpiQqyzsY4yBviZabirQu16tih2bBW4zqOn3dgOs67x0g+PmgHevcLuvsA/d6OdXMUdHfOgff4FZP1Ev2a7v1QM7/CvlekrVvlEtj6q3MqlT/ofAOHL9r//BPRm0OcXVSeMkrAGjgsNkGyU4FxJ8WTjwpDZYUx+7zGSyaNGxB3l247SjYw8ijul/hGEgZxzPuxSDWA1AdTz6wZfADV4Ebr5RijG+4bkmFsYzZoandlQ/oDfvBLEepOAHL2Xw3H1v0GRuoKf8RtdnYOz40ckKefSbbCoePJsUsJHxyMNXPKP/SCa8xNCgTSKvwRUbwIBCgFU++NredVnDBfolX+hfdBE79cOW/lkAHaEl0UBvwWfOA3RljmCbfgEbC8/aZpt6+E6ND6GDF/MTGRVb5hHvyKUgV7/KZ8l2dMNGsMaH2IM8eM7Wg3ytwR8Z1Y/wo21rkpTXtddJcMXDBgsm4OxnwCvDbh+Z8EywUT99TGB58szi4uqL+hSly0Q2ctnh1oWBfPWSpupjPwHGgMH4OobnrYsT3mAEVgaJKs9nd+Ls3pWLDN7RZfQ/SIVnDzbyIxARIJEND/VHzkgPkiiBhyRCICJR9ICPXsiUF3wNjMqlJnjDy2BOzbiOGbrTPioEWvoIYGKDXGRiU+dlkFc3x3Y6bDAoOgZZtb3qswUX6JE384XBGBmjdVHl+cycgxa+6FcL/PBttQ26uhZJCvxqob/ryHtvcwzzBL74kx9rg7nY1/CS7fCC3jmo/5SxpWYeMp9qYU6gC+11fVeah/g8Xk0PEYkDbXby10V1IMurHW4gNbAew1AWNYGFoNkLgY8+dt9bCvQETJODNTxqgV9PerXfZwKrPKzh1RMrMnsQlodJYWvwIogiA3rsmQVubKgJAx1oQ48euLfisuQL7HGzoW1ba+QTwHsCgl/1u/K7XX1tMq62OS/ruKobmIApPiSZkFQ6rbKrPpXHoc/yn81lsIEm5V8Evl2xQWVvBJLgdoOOnSvB9FiFQMUOnyDYi77pAbvT+U4yIfhVXvDvCQ7915K0Jx6Cp0V9aoLjmaBZ6aSnJnDtghcBsNITcE14lS80NSBCRxLp7YzZisuSL+CDHrMAXXXjmdNzLYyrCQ4fwa9iyWkbLPWfSQEa+TmOBGUb/fByHHLrlT1yu785ydWyZnul3efZjUu1Vz7MnX69a99DrZPgjuR5FwcTPGUdgaWFuj76WwoCEte9NTBVCpPTLHlUWr41EuT6N0cCG+1+A/KakWA6K37P6cGcpAevOpbTwJJ+0M9Od10+QbmfdOANjxqgtaEGcZMB42vZisuaL1wnW044zJGqL/r0kyUJCrvAlEQFX05YdTOjTOyEHzraRoLCP9gtRnwvhRc0ddOEbDY+8OFnMhOnNdulO6RGPpvDXtBlC6Z93LW/J8EdycMsBhZWEtx2QFmodTe+feRulCQ/dvRrhQDnFSK1Rd/aR7vP1PSPCjKlqzS2LY2t/JS/NYBhrzLqGPCm3TkqDfS1kGz5WXbBxTGzGtnIhedagZYEgw1gwHu1h/EkLJKxWEM/ujakn599JDTbajLkBCQuYACdBZ3hTz9YiqP9p67dGGnDqeVdA/8kuGvw4oXaYOAkUHhNdKGmRO0zIdC/o51JjZOL9TaAP2JK2Y5AEtx2rEJ5AgTYjbKTZider8pOICosrwwBrgS9nrwy074xh6tk1kiS2zewbHpJgtsEU4hOiQDXQCzeei12SnnhfR0IsDnyOrFeAV+Hdf9agV2ja9drsvGUtiTBnRLd8A4CQSAIBIGzIZAEdzboIzgIBIEgEAROiUAS3CnRDe8gEASCQBA4GwJJcGeDPoKDQBAIAkHglAgkwZ0S3fAOAkEgCASBsyGQBHc26CM4CASBIBAETolAEtwp0Q3vIBAEgkAQOBsCSXBngz6Cg0AQCAJB4JQIJMGdEt3wDgJBIAgEgbMhkAR3NugjOAgEgSAQBE6JQBLcKdEN7yAQBIJAEDgbAklwZ4M+goNAEAgCQeCUCCTBnRLd8A4CQSAIBIGzIZAEdzboIzgIBIEgEAROiUAS3CnRDe8gEASCQBA4GwJJcGeDPoKDQBAIAkHglAgkwZ0S3fAOAkEgCASBsyGQBHc26CM4CASBIBAETolAEtwp0Q3vIBAEgkAQOBsCSXBngz6Cg0AQCAJB4JQIJMGdEt3wDgJBIAgEgbMhkAR3NugjOAgEgSAQBE6JQBLcKdEN7yAQBIJAEDgbAklwZ4M+goNAEAgCQeCUCCTBnRLd8A4CQSAIBIGzIZAEdzboIzgIBIEgEAROiUAS3CnRDe8gEASCQBA4GwJJcGeDPoKDQBAIAkHglAgkwZ0S3fAOAkEgCASBsyGQBHc26CM4CASBIBAETolAEtwp0Q3vIHBhCHz9+vXONL5LWXdmVATdKwSS4O6VO/5VhoWfxX8PHXOhKjmfluYUfe/evbv56aefbj5+/HhyS5H14sWLmz///PNkc32L3Sc3NALOikAS3IHw10U0e95VxKNHj25+/vnn4bCZjNo+HHhHjVWP+nxH4i9GjNicWmESyQ8//HCzNKf++eef2/lGcvv8+fOpVfqPP4n06dOnNy9fvjx6koO3dj979uw/med+0O/WVR/brGtfnvdDIAluP9xuR9VFRABhQdUfbfy+fPmyk5SlYEQwWJKF/L/++msnecci/vvvv2+DFUFLHJ48eXK7Uz80cLLor6mAy+PHj+/EpA8fPtzOmV9//XUojwRActt1ng6Z7diITJNcH7pLoB/Rajdr5r4U1oVxgWcLscR21zdtKYchkAR3GH43f/zxx+3EHC0iTmEEsl0LE3x2goPfaKdt4vvtt992FXcU+vfv398G7HrFRVJDX+w5pMAHHgSsayjYQ3D75Zdf7sQcA/0IP+YLifbQDcghhhjcmUO1kJDx+9qGjSQ5WjPy5RR7aDGBjupdeGMLuo7s0t66hnbhHdrvETgs8nzP78G1mOBGi4jg8erVq50xGS1WmBCEfvzxx++C0du3b28XzPPnz3eWdYwBnz59upU/OgWAwaF6/f7777f8zxmEj4HTuXjgA+ZUx4/EQHIbbc7uWldOkX0ziL7oPTt5qqNrsJ94tHstQcpnVnOFiw7Mb28mao2crcUEx/fHWpCB/Wx6znGSrrpc03MS3IHedNfVg8chbGcJbsST5EKQIvGda2F4StsVg7obxjbea+GdP0LAPjDp9JXW8Ws0jql0PDve/l6P6CuN/aM2+dMnXW2rY3yudDPaEc2IltsAMOxllhg6ne+dd5Uvzb61uvRkNEp8XQbzb3TjQRLpdu+qs+uLxNN163pseWcjzFzGXgvJDTtfv35tU+ojIZAEdyCQTMzRIjqE7dYEx2JFPrvJvns9RP4uY12wu5wC0PvNmze3O1Z0JzmzC4ZH3Q1z8qMfPOqOmW82tcCPREhAk46AN0u4fitELvTQcgrufkQGvOkTZ+iRQ+CzcPWHTuip/rSx4+++ZDxt/EaFYAcPv2NCh+weXLnOqzqBHSfdUaCHRz8xIJu2mR5VNzGQB+9VRxJMxaOO3frsyUb8HOf8gj9y+w+/Y8Po+hV/qjP8oIGWH77ZUsAT+461eXRDrD+xi3k4ugHaol9olhEYr7LlMektCLBYWEQuPCYsC+uQhAPPUaAqYm8fXSw9KHS6Le/qT71LIciirwt2y1iCMYHeXSyJyFNgX+izwFflgAP8wIFA5LeXUdLVPwQ4EyDjsAEdeoFHDUDqKi2yuIZGLljgN2zQJ8jxGd7QzfzrTh55zh99XAOsV9LKpc9k1W1ewo95is5LhfkgBvoYnZzjyAYL3g8tI1zgD2/68HH/0d43POiBn+ir2IMv49k4VTxnei9hNxuz1g7e2AOubEjATlzXxqZ/dwSS4HbH7L8RBtK66FhU/A4pjF9LcOzgodu6E13Sh4VPgMQOFhy8R4kOe+tO2SBisF+SYZ/JxORmO/Zij0nHdq+uqlz7qP0+VwMZ7SMMDcYExBrg9GNPDjNdbTcJqQ82ELyWvrsiF91IEr2IQQ14yKgBXF37H6ioU98giF/XdYZR16knt1HQV3bVu/PZ8g4uo7mEDvR1/trWbUYWbYzRbuZ432ys6SR//60ea2L0W+Njv743XszslT714QgcFokPl3/RHFzYLiKMYTGOklNdGGtGM/FHPBzHFRuBlJ1+TwjSUCuzto2e2VUaJFjU8MUOxlt4RqeaaAwi0G4pLHD0HgUxeI/aDW41ISlrxs9rK2ypxYClrfaZNGr7jDdjRn6nHb9h30hXZXlNVnGkTx3WsAQn5HS/q1Nv92Sn/FqvzTM3URXHkT88ZdZ1sHXudX3QqReTat8UMF9Gc4bx6oQenHihq/p1GaN3MTUhzWrothR9Dz1zRF/2ubCFV2i2IfD9bNo2LlQ3N/9N0BrQmLyjCc8unMW7FsAAdinwsGBJSNAsLQyTD3Q96FXnQdd1IqCgL3LcvRIw+skEO+E/srfK8FmdRvTwIRj3shTE5Idu4ML1I7xJ0KNrKHGr/kKedtQTgsGoB1XoDfI1YJqgRvTVJvpHPpFn1aGO4xk/MpbA2AttJNdewA+7R2XGS1rG9YQ94megFlf0JBnAHx9tLdDzG5U+D/RPTb51nL5WjxkGdUx/dl5UP3eaXd71vfy0YTTvd+Eb2jkC49k0p09PQYDFP1s4BNxaWPws3tmCrLTQjYIYNC6S/hdXXZ60o6BXZbHIRgkQfUloJAt+XR48dg0AM3oXesdGzHoCVn+xcGfNFRRXd7OgCq4jf9HWcVJXdOuFYNsDMbrTZvDqY3xHFuN7ob3z7DQmUXSrRZx6oDQhzvBD3myewb/3yw/cLcqu16j06dOlhC0Pa+SNsKFffPUHtvbkKx91Eo+egKVbq5XZ5+XauFn/yMdufEdrcMYn7dsRSILbjtU3lCxcFuQoeBgI6gDbXPAkJH+VjuceWOz3ymj03Y3F3hcJurnI4TGTJ/9da5PALKj3pGug6fS2i416GCRrQK88Hdft1lb5WI9wNWnAq5aZbfq94so43uFPcJ2VHngrHWMJgEtFnSoe0NveA7H4zRI+c2ZJZsdLfiaZKrvL4L1uGrbMvS6vYiF2zGnXUsdBen1qP7rAu+so/azW16M1Phuz1D7ysbqp69L49O2OQBLc7pjdjnB3N1o0LIgeMJ3IimOyzxb0qH3puxt/cl+DiTLYDbtwqOHLaacnGOl3rbVJGXW83z1qG7igQ5UPjqOdLePUWXqCZE0snuB6UuFatZ4y1KHjyh8emCS7DV02PJQ/whCsZ6cP5ffAW5P1bHylGeHNtaz4VZxm+KkLtbbXtvrMRgrMwInS8a5/kVrH8VzXgHbDq+PsOGlGfpMGnsxzaKi736XrvjMhMp5SMXXMrOaqG71Zf4cUNwfdPmzAlqW5g7676HyIntc2NgluT4+6Y2eROwGpZ3/VVxc8dIxn8RiUqho9ENNnEPMvHJXpH1T0hOqilj+7Ua4aWfyzwFB12PIMH65YWKB+qyMYmvwNKPIyQJOQ0R86rkHF0jbaKfKBHr7ws49+bCPZYBNjlY2d/TQIvboSrPiBvzrVU0nljW7wVT7yGFOLWHd7K4088S32oB+8LQZl/Ys85JBk9BdykI+vsRf7scGx+IA2v5XWduzrwVV8u+3q5Jxh7sGbml/FuernOGoCNvIpzhP0woZRUZeR36Q3SYDhEtajxI3v0QncGOu6kPesrhjoG3jU32ys7dA6x6vP7VdfNoWjAuassVn/aEza/kUgCW6PmcDCJdAs/XrQYJISYEiI1D1IVjV6gmPxL8miz2AiH4MBgZyFgcxZcHHMPjU8+e5V9SOYzJI33/KghcZkha62VTsIjPCxr2OKvuAIL2hIbARRE0K1hyBDsJLW/zKFgXU0BnnSw58xo8Co/tpT5dZnZJAQ5FWDOX3oTh8/5IJV9xn2Yic/+hkHjW3VfoNzpe36EDiXkgU84AkP5iW68bz0rROdoCUh8GPujXCrurA+sHmt6I+OSx2HfuBci7gxfjSPKm1/RhYYKFsfWdc528fyXv3KmErv3BHXUVzAHmQvnfJGctN2c5MEdwezwAXvghjt4qoaPcHVvq3PBBUDEnUNplt5XAsdCb6fcLWNwMrvoRaCLfNj7QqO+QPdWhIHR4J2nXtrY7yF2DXxPDSfkeRSdkMgCW43vPaidsETJAwoJL1ZITjMAvJsTG8naJNIkcMunYRH4STzkAq7b/Csu2bt53RB32jXLM1DqDntMF+4epwVT7pbNkrMNeYchXm8tIFAJrScclLGCLBmwT8bgDE+S61JcEvoHKmPBe/1ggHX7zCj3e2hCc4To4GbKyiCCMEEXZaS65FMvjdswJuTM8meUwrBAhy8skpg/fcbGVfBJCKu0EcF/Exao/7aBh+vPcGZ+cwpjed6Vckz12/xQUXv22fmL1fRFbdvKfK2hEAS3BI6R+pjh1wXMc8EXb+hdDH0EXD2Lez04MHioLDr9h7/IS4U8PBbHrjwwyduAPbF+drGgQe4jOYI12N8d1srbJ7AV2x5hyfzr27m6J9901yTkf4gsBWBJLitSIUuCASBIBAELgqBJLiLcleUDQJBIAgEga0IJMFtRSp0QSAIBIEgcFEIJMFdlLuibBAIAkEgCGxFIAluK1KhCwJBIAgEgYtCIAnuotwVZYNAEAgCQWArAklwW5EKXRAIAkEgCFwUAklwF+WuKBsEgkAQCAJbEUiC24pU6IJAEAgCQeCiEEiCuyh3RdkgEASCQBDYikAS3FakQhcEgkAQCAIXhUAS3EW5K8oGgSAQBILAVgSS4LYiFbogEASCQBC4KASS4C7KXVE2CASBIBAEtiKQBLcVqdAFgSAQBILARSGQBHdR7oqyQSAIBIEgsBWBJLitSIUuCASBIBAELgqBJLiLcleUDQJBIAgEga0IJMFtRSp0QSAIBIEgcFEIJMFdlLuibBAIAkEgCGxFIAluK1KhCwJBIAgEgYtCIAnuotwVZYNAEAgCQWArAklwW5EKXRAIAkEgCFwUAklwF+WuKBsEgkAQCAJbEUiC24pU6IJAEAgCQeCiEEiCuyh3RdkgEASCQBDYikAS3FakQhcEgkAQCAIXhUAS3EW5K8oGgSAQBILAVgSS4LYiFbogEASCQBC4KASS4C7KXVE2CASBIBAEtiKQBLcVqdAFgSAQBILARSGQBHdR7oqyQSAIBIEgsBWBJLitSIUuCASBIBAELgqBJLiLcleUDQK7IfD169fdBlwZ9UO3/8rcubM5SXA7Q5YBQeDmhsDpbwse0t5VwEXOu3fvbn766aebjx8/blHx6miw/8WLFzd//vnnra+uzsAYtIpAEtwqRCE4JQI18PfnNbmHJgvlrckZ9T9//vzm0aNHt78PHz6MSP5rI9D+8MMPt7Q///zzf+2nevjnn39ukENy+/z586IYMaj14oB71ln19rmqSHJ/+vTpzcuXL5PkKjAP5DkJ7oE4+j6a+dtvv/2XJEgA/p49e3YbkP7++++p2iQNEswff/wxpVnrePLkyc3jx4/XyKb96IkOX758mdLYQRKE9tdff7Vpr5ogvlbQi+S2pBdJEPyhFXewIDFewolvF/3BwSS3hl36rwuBJLjr8ufFWcMVEoGfhEUhGPH8448/3iafT58+DW169erVbWD+66+/hv1rjZxsCOy//PLLGum0n2RAgthSTHBrp70lXib1pcRF0iJRLZ3cwBS9wbjiri/ue4LbR39sYp69f/9+CeL0XRkCSXBX5tBLM8eg2gOyAYlT1n0tu5zIPK12O3exDayW8CDxkdy4jpsVTj7wILn1zUAHMwAAABoFSURBVAGYk/TvczlEf5L6En732e7oth8CSXD74XZxo/rVlt8revtdG0bQmV0T0kcS6YFY3buuo3bbqp2jthmvGS06oZsnshmdfDntzeyEpo7nuRbef//991t58Km0lY7rWnRaOoGR/KDx5FbHX8LzIfqLT59Pl2B3dNwPgSS4/XC7mFEEw7dv397+NRlK885pgm8S7NbZ0c6uAU9tJCcOgi0nk1EZBTNPQozj2UKiwabaThvfomgjMViwmzZ+o8K3P2RzypG2jmeMwZITRZXd6eSPrJmdXR4+4S//LOigHj5Tc01bi6fh2lafOT2ix6WeYg7V301JnTcVnzxfHwLjFX59dj5Ii0hmBmp3rfyRAycJdvkkGILd0sliC3CzE8XaWK8hZwHHZNb/kAQbCNSeVKgJ9tjDqY8kwwlFviQ5n9HJxDpKRiR78Kh/gQhPfrWAK23U8u56Sb8UWNl8II9vgQRwdOO9JyETqjbLu9aM63rWfvXreFaa+/x8DP2ZNyO/32e7o9v+CCTB7Y/dvR/Zk9so0JpETIC7GkUi8TSIPE40o4Icgnctyp4FbQJRTWSONdB1fvQzhkDfTzeOpTbBwacWkwv2VN7o2fm5MXj9+vV/LGb2zJITyRT7+CcHtWh3bfN0VvWq/TyvBW90huaQ74Bd5l2+H0N/7O+bh7u0IbLuFoEkuLvF+86k8ddiLOa6WycB0VaDpMmiJpmtJzJk+OfoBE0CNcmlX3ly3QZdL6NAXmk4jYwCMu2zkwr06FBtrDx55koROupaTFBr36fgzXgSYS1i2TcLJqdKy7P214TDBgH9+3UmQXlms3zRCZ6jggz6LzW4H0t/MOCX8jAQiKev1M8Ewx7oR0HSIFsTgt+ySIhLBboanKHlRMP3IZIFwZpkR3LrQR9aAs0saM8CmsllpBtJGp79ZNZtMBF13cEHzNaKCbInQjHv40e4ax9j2FDwHQ5+0IJrxUvakc1VFrbPEpzYzPorn/v4fCz9wYhfysNAIJ6+Uj/3YGeQrMHfZNFPIrbX01+HCZoe4KUhAZDUSHSz5EYAR8eqj+OpTUIkylpmyQWa2VVgHc8zSYVEUos2b0kA6lYTpOP7yUvce3LSDv9oBB9wAsYGeNUi7QxvabvPbac2QXQ8K819fj6W/mDUfX+f7Y5uhyGQBHcYfvd2dA92Bklqy+xKzqDsKWLrlaV8t9Qmo1HQ9tsUQb8He5OLulVZXgX2MZVmlogMoLOEW3mQBPtJT3z7psD2bucMe+SAdy3ajI6WTkM7Os1OxNo3S3Ajfsq6D/Wx9O/r4j7YFh1Oh0AS3OmwPStnTk4sZv/owyBp8CeJEBCh64VgzFgLz8cODCajegpCHvoSpEf/EJl+kkvVTR2p2Zmv7c57oDSwm9T7SQu+0ihrhIUJyyTkmFm7ia9uOOCP/f0U2W3mOnOUiDud+lJrX+dNn3Oh+4I+7NCWyq8+j2i2tq3xsX9f/R1Pre87diNd67g8Xy4C/4til2tDNB8gwAmHJEGy4I88qPkRQDll0EdyM+FVFgR5AyGLn2TE9ZnBu9Lu8wxPkis/gwtBm2SAXv0bVJVhEGcc9J6MDICjBFXHG+TevHlz+52rXin6RzLoAn91qv85L8cjuxYTGViTtAyis3Zwx1bk4xPkMY42bZK/mwF4o9PMF56Ke9KUD3aQnPlH4+LOHwrhh9nGAPqKkbysxaNuOkZt0ENT6eRhvSZrH/3lTS0+9fQPVshlXo3WQh2f58tDIAnu8ny2WWMWMn/eTtBkEfO9h2cCRQ+ilSmJkADNzp56ibaO2/JMEEG+355qPfsGVfmiC2NIgjXJEKho71eEdSzPyCexQ4seNdjxjA7qBB38auDjnf6eRBgLtvz4QxvHzNrRhT51Gemj7tBhLzTox/uoIJNkNUvy9KMbOmojz+I+4mnwH/XRRjKTlzSjNvo6nfTWa7L20V/e1MxrcKxFf5LgjznPq4w8nw+BJLjzYX9nkgmIBI+14I9CBBGTIQFpafd+ZwZE0GYEPDFy0ju0eCreMm/uuyxOv8zrvjFRb3BLghON66mT4K7Hl1NLCFAs7tnOvw70ygZagyVJL+VyEOBU6HX0IVrjf3jdhf9PKYsrYE62/R/riw32cZpPuT4EkuCuz6ffWcQpjAW+pXAl6fcYd/B+r7qLnfwWHUOzjAABm2tHkhzXzPsUedxFcjulLK5LuYadJTdObSTXu7BzHz9kzGEIJMEdht9FjOa7w9YdKjv2Ggx45qqyfle6CKOj5O2VG/4kyD/EQvJi3j9U+x+iz7vNSXAdkbwHgSAQBILAVSCQBHcVbowRQSAIBIEg0BFIguuI5D0IBIEgEASuAoEkuKtwY4wIAkEgCASBjkASXEck70EgCASBIHAVCCTBXYUbY0QQCAJBIAh0BJLgOiJ5DwJBIAgEgatAIAnuKtwYI4JAEAgCQaAjkATXEcl7EAgCQSAIXAUCSXBX4cYYEQSCQBAIAh2BJLiOSN6DQBAIAkHgKhBIgrsKN8aIIBAEgkAQ6AgkwXVE8h4EgkAQCAJXgUAS3FW4MUYEgSAQBIJARyAJriOS9yAQBIJAELgKBJLgrsKNMSIIBIEgEAQ6AklwHZG8B4EgEASCwFUgkAR3FW6MEUEgCASBINARSILriOQ9CASBIBAErgKBJLircGOMCAJBIAgEgY5AElxHJO9BIAgEgSBwFQgkwV2FG2NEEAgCQSAIdASS4DoieQ8CQSAIBIGrQCAJ7ircGCOCQBAIAkGgI5AE1xHJexAIAkEgCFwFAklwV+HGGBEEgkAQCAIdgSS4jsgDev/69evNly9fNlv8999/b6YNYRAIAkHg3AgkwZ3bA2eS/+nTp5unT5/evHv3bpMGf/311y3927dvN9GHKAgEgSBwbgSS4M7tgTPIJ7k9fvx4mtw42fHrhST3448/3iTJdWTyHgSCwH1EIAnuPnrlhDpxzUhy++2334ZS6H/58uXtb0RAkvvhhx+S5EbgpC0IBIF7hUAS3L1yx+mVefbs2c3PP//8nSC+xZH0nj9/fvPo0aNpAmQgdCTJfJP7DsY0BIEgcI8QSIK7R844tSpcLZK8Pn78+J2oz58/33z48OG2nQQ4O+FBQDJ88uTJzYsXL77jk4YgEASCwH1BIAnuvnjiDvQgKY1Ob130WoKDngRIsuTKMiUIBIEgcB8RSIK7j145gU78tSQJactfTW5JcJzi4Mf3upQgEASCwH1EYK8E51/ZWd9Hwx6CTuI/+ovHbj/XiSSkLf/ubUuCg/8uPLs+W993sXErz9AFgSCwjMC5192WmLZswb+9Oye49+/f3/CHCvwlHT+CIX92ftdlzQFr/Xet77Hl8c0M/Ela/JaKpy18taVsTXB//PHHrewtp8Itcrnu7N/+/KMXbPQb4RZe56Rx7s10sH9Wz8b1dsf39q3vjl+rt/K7djpxuo92qlutD9GTfyNrbBl9s1fOITJmY4ltP/300+1tE3IOKcuRsXH2jxQIbBQCJ8mOv6jbcjJo7A56rcFdfWSILrW/B03pLr02ceGDpUJiYLKu4cBk4i8j/Vb3zz//DP89nLJISPA9xjUlSfL169ey/qbGPuQw8S+hgB/6zkpN2sxTf4wDS3DfUtbkrPGY6aE+2PDrr7+usXkw/eBNrLtPhTXLoYNNqX6j5r3HxV31dt2NYvupsWCtMz/RYet6GNk3X4WNmiw+mvD+scEoyzcWR301uKLTKMAb1K/5L/1McGtBiH5wWvNRXSA+M8mWCnyZ7IcUbgBGPpQni/VQGfI6de16AJdZwW/YA41JmzYSPG3Yu1a2yFnjMdLDMa73Y53O5XupNX5iTfzyyy/3ygTXNvMBf1KMfWtxYc0Q1uRoXd4lFuCNDtq2pnPvn6/CRjn73uJCWwuejd3BryY49TJQyNgFCt21Fm1cC0IEzBpMj4kHk28pmK/JMsgu2QD/S9ioMAfZ4TsnlxYldKOkbeJbwm0XOUt86JvpQR8Bva+rNX7pvzsE8A1rY3SD4hXfIdrMeB/Cc9exrCHm6L7JelOCQ8gsyBA8WQhLi3lXo7bQExCRTY1u/ThOO8AsFe+RrTvtqN02aottvlPbNqKrbXWMz3WstNbSULu5WEvi4ABGpygmz303ONiw5Cc3MtW/FZ9T2LQvTwINP/0yw8TANEra4rmkw1Y5Szzoc12PAuRsrNj3/t7e36G3jbqWUfuorY5Z4lfpKh+fa3/nM6KxjXpWKs2MTprKw7bZmErbn1kTrGvm21pRTqWzbSTbzXNdd4xdGiPvSuOzfbW2z7r21WdPqftstjZFPY+83Vi+yZHctgBcFT7GM4sSw12k/ShN/yiAIJs7XXTmQyr6899XhL7e9WKzH1q1jzZ2RkwqAhGFPsYTpOnHWeDkWIMH37Y4biMP2lnw63/Ew/jff//9P3kVO3RYSg7Sou+pEpyTb2aPOsxq9BfLEY2LuCZxT41gCeZrxQW0S73Gs/f/+eeft75lETInwHuGCTrT39cTPGlf8ukucrqO/X2khxh1Wt7rXK+2EQfQ2TVY6UZrou7GWXPMcXzJ6RX8bGNdjXDs6xfZrJFe3rx5898al7/rUVqux4kT9PtjPlb7aEcPfqPif94OfZWD7Fp2waSOW3p2c4/++G1W9pHNGOyt6w7+u2IBj9H6HmHG3B6V0SZ3RDdqG3usUXZjUY42jKXet7iYlpwz481iwMEUHAyQNcPTP9KNxcFCJFE5iakZ7+Tn/dWrV7fJE1ochCz5MZZnAhQ/kyxt8DBwqRfvfMtCP501crp/xKNs+MpD3Soea8kBWnUbyau89n3GZrDT5l34iLu4jsZiN/xrIUASTPR/7Rs9M093+a19dxzJYJ5oh5g4vzq9m4IePAjS2FoTQB+7i5w+tr+rJ0HetQjeo7lCooKe+cS8Y15STLjgyzdE5gE/5x1j6prQdn3Hpg+cTLZs8KAHG9rgWwu6srbr+kVfcKvrX9v0AX1gx1hLtYVnCrrDS/1o05YRLugDHtihfGXLYw0T13ifD+o5q9HLTQCYEdt6WZPd/eF4dfLdegsW+EYswJxfLW6IxAyeYFh9U+l53hLr+hjev40cI4qbm9sJb5DBaUw6lCbR7VtwBiDCC8OY2KNE5+SvcgTZCYFO6IczKfY7uetYF4Nj7bPdd2vaAZeks1SQjx1Oamid6Jzm0MkCbV8s6Eo7Tq9FHpUv/egPPf1LRb5dXh1jcFuqK319Vr81PeoYn7eMBVN1Rz9wYAF1/8nzHDWJiTmij7Wr+0zdWDvQi3fdMGKbfKS33lWO42Y1uDKHWIP+tswp1i36s4apZ76HV18TzseexG2HfuZbcKGfoF4xEm94WLDNeWMb/qjrS5l1HLSs16oDsrCl60w79vf1DQ/ou3zbscFYRdtIf9q3FPQ0ycGXhDsru/gDu0b674oFttXYiX7o0TeRyKJ9Vtb6Z+PmHMuI7ixAZSGyUPuuwUVbhg8fGWsAwNk4iV0I4y08Yxg7uVp4xwEWQYcnxd2g/dZOaOT0MgMQ25GFjFmRLzxqQQ7jtZM+dja0dR2U785HPk7+3g5m8OmL03HWM93spza4LdWVvj6rH/WuxbHdv/LRr1z3kAQIrOC25AvH3lWNLsyPar++qW3qo034ruLNeqpBT3rrXeU4bqlGhzpnnZtrc0r7WG8jG5E5m3e248tanAszftBKU9eTMYL4UQv8sW8p4BPHoIG2xp3Kh2fjSZ+nI30c27GlXdsr5rSr677zmnEkbmQyF2tyVp+ZbNurP8Rl5It9sFAH6lGcc6NUdahjZuM6zeh9NcEJQDfWxVADNRMPkPn1gFyFQ1fH0Qeo7IRYNFx7MOHYMdXsLw/aOxhOEuSia59EjGUcuvWJSh+7H/pq0fa+c6s0PDvRO19sgW8tYlQXqViOdKatJnN5gR/6ri0KbRjxltchtbb3+bGFJzphAzqOiovJRDDz3WhsbWMu7fqr45eemRv4h0WqDOYvuo4w0aa1OdVl7iqnj+/vzouuB1ivFcf2NVjHOS/6mtD+jo1zYWk+u0bBmh+8GEdyq+sJPYgnnmzYIEE/KsQXbDbujGjACH/2mMYY2rvOs/U8wwS74HVowRb0GfllJnvkD2NU9x36zbDAhlGcqjaJC7ayVti0Ioux/dRcx/Hs/JjFik7v+7cR3dZSu1sbGQuYCK7FhVjb+jO8+mSBhomCk5iY/Jb+4W9fIDoFfdGpL1z4C9LWCanta6CaXKudyACfnshNTNV+A0a3SR6jCcuk4LdW5N39tDZua78Lp+u+Zbxje3ByrIuJYKUdI79KP6tNkFvrfn0y4+t1y4gvvh9hok2j9XRMOTNetju3R3oQfJaK34qXfDFaE/Cc2U9wXJqjrgVwBW/iA34iRtS1VPWmHRroCaqzJAcWJkNs64Wxo7WGLqPEZCzq/h9hgo7wWcKy67P0Dq9RohnJhs/IH6MYpcwRFvpmyX+MZ67pP3xCUsM/4A+PpWLsXovFncdqghOYPol0TDcKcBhjcVfr+6G1YHZDbUcfHDxauILUdXCxU9ei7WvgI69PdANy5zmaIAb6viBs7zy0tSfPqnt9ZlJ1P9X+Q55nOm7h6dhut2O7v7Ymdcefuka/0fcXfV/Xgbpo09qckp56Hzl1/OjZud3XNbTMq5lPSADetPQ5X+XMEhY+JJnU4pXYUpAX09Gc7wm5v5tw6thOgz+wh7VSMXGtjXw5W1dii121jOKEuo3iVR3bn7v+9qPTKBnv4o9RjIL/DAt9s+Q/xrveRxvamT3a5bqpvrFvqV5NcLOgwrEfMHvwhd7FoUFka0A4RnEXMOLlxOqTVFpB8p2aXZ2Lrgedme11vIm+Lh76tb3bjW4uFp3qJBc3xnM6cMHJQ3onlPS2V73qMzL5naK4A1RHZaDTml7dDsdaozM+s4AxbU7yNf6OO0XNH0WhS7cbWdpVdafdALGUGLqu+8iBxxo2BDzmdy9cG43sgh/JjbVM4Nbvrpkqz4TV7fevRHvMcP6PAp/6iV0PoshlPTknqEd2jeZSlzdas/pytNY6T3T1VN9vAWZxos9p7aWumNZ29B7NIX3XP+vs6o9qV9VhhsXMtm6Dsbsnc2JwnyvVXp7pR69a0K3qV/t8/naErf9fqziMmZQwQxmemeh8GK9Fehc9wLJbY3K4ECr9rs/IN4nhzF5cKKMJDi16oDc1vODBREFHdK1FW3riqjQ8K7MvlpFDKk+dShs/9GIMeoEfCwQ9wZ5vOrQ5cZ1obDLQu/5lWNePd/Dok2NERxvy/c1oart26nP6nMj0rfkd3aDrRRvBwCLWBNotdjvu2DW+I0GAKVj1YiDHttpf2+GxVvaVI3Y9kSjPIMxa0tfUroc+V8SdOSpPfYwv4FcTj/TYD094GzOcw+pCvRTkKx1rovJELjGo8lQv5GmbsusaZd3zUz9t53Rai1i61tycQoM+rk9ksU6ZF+jU572YVB3g4frB1+jpfFcubb2IF/NJG5GNPUuyK3ZiUrFTDjahFzrVzYM6zbDA9oonttTYBCbEWnjCG93xF20dF3WxRifsq4V3ZDIHZ2UxwTlZmLwoygTnxyRA+e5E6TESoYwjeB+r4Ax1oO4FfWgfOU1a9AZQ6KihHemILdCMJpi8qPkGAF1PkPDuGwDoxZHFwYSx4GDG8IMntqCXbegp3tTwRi78RvrLl9pFtEYHLT6DL4tBeZVXf2aS9YDo4oHH2sTFH308MuSBHyzoA27oV/Gw/y5qMGT+owO//p0Ye+xTT/QCh9o+mhtV/33lwMNAhC69ME+dO1Wf+tx1A2v6Oz/54JM6VwzA0DN/GQtms7kAH+jWCjJcP/BkXJ0fjMd2aKqPoOuy8ZtzCV7IH8UCZGonfOs6p09stBGbKxbaNIsT6O8ar/gu+ZC1ge5iuyZ7V39oE3bXGLWEBbhUPBmLnh2LPv86puJVa8YQI+BXC/bjZ+LMrCwmOIIdjLcERgRID+CMq5NhpkDaT48ACwd/1Mk6koq/WHBMSibNGj084Nt3VvJGbg8s9lkjayn4SZd6OwJgfs71x9yZzYntVjxsymP68NL9YV6Z5SHix6wsJjh2/kvZsTNlUnP8JGhxdEQxCkfRlPMhQNIi4NUd4po27PrW6N1lQtsLc4Dd2Zbiro4xKYcjwLplF36Ogg+Za679c+hwDTKP5cNL94fX9KP5RF7hVNdP8dX/iwmOiVrvnOvA/iyQ7tgJeiQ5FES5BK+O2N296xsWzdayJcExuZgjfYIxB0iOu/icJDeaxFv1Dd2/CLDpOFdyQ4N9NlPx3bcIHNOHl+wPEhgxi+vOHks4zXHtu3bLNE1wDOSqsd97fuuK/70BJPQeIwlY3pGuKfE/Lnk6FQJsVEhGfaKM5C39Q+VKz+RjE7OFZx03e2burJ0aZ2PTfj8Q8FsT10ZZ9+f3yaX6wz/6GSW3XVCdJrhdmIT2/iPgTs4T9kxjrwS4ml5KNiQ1EuboenLGO+1BIAgEgTUE2OiS2NZi1Rof+pPgtqB0JTTsqteuKTnp8ddQJLelBEcfCY6TekoQCAJB4D4ikAR3H71yIp3YEZGUZldH/NMOrpXZQS19g+P0xtXk1j8iOZE5YRsEgkAQWEQgCW4Rnuvr5HQ2OsXxj2ZJWn5zhWZ2guMfepoIrw+hWBQEgsC1IJAEdy2e3GgHpzOSU01enMj4Jx7ce1u4quSvGrmC5PudhUTIHxMd435cnqmDQBAIAqdAIAnuFKjec57+hatJyr+0qleX/lUs15D+lSTJjeTouHtuZtQLAkHggSOQBPdAJwBJbpe/VII+/8WRBzpZYnYQuFAEkuAu1HHHUtvT2RZ+u9Bu4ReaIBAEgsApEUiCOyW64R0EgkAQCAJnQyAJ7mzQR3AQCAJBIAicEoH/A8/Yhh4hSN2sAAAAAElFTkSuQmCC\" width=\"440\" height=\"206\"\u003e\u003c/p\u003e\n\u003cp\u003eWhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(K\\)\u003c/span\u003e\u003c/span\u003e represents the three outlier signals (APA, Expression, and Splice), \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\beta }_{k}\\)\u003c/span\u003e\u003c/span\u003e are parameters defining the contribution of the 40 genomic features to the CRF of the three outlier signals, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\alpha\\)\u003c/span\u003e\u003c/span\u003e defines the intercept of the CRF for each outlier signal, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\theta\\)\u003c/span\u003e\u003c/span\u003e represent parameters defining the edge weights between pairs of the three outlier signals, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\phi }_{k}\\)\u003c/span\u003e\u003c/span\u003e are the paramters denoting the categorical distributions of each of the three outlier signal, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(C\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\lambda\\)\u003c/span\u003e\u003c/span\u003e are hyper-parameters.\u003c/p\u003e \u003cp\u003eTo train and evaluate aWatershed, we utilized all gene-individual pairs that have at least one of the three multi-tissue outlier signals, which are defined as the absolute value of Z-score greater than 3 or \u003cem\u003eP-value\u003c/em\u003e less than 0.0027 for splicing outliers, measured in GTEx v8 data. We also used a set of 38 binary and continuous genomic annotation features aggregated across all rare variants within the 10-kb region, flanking each gene. We then trained aWatershed to learn edge weights connecting the three transcriptomic outlier signals, weights representing the contribution of each genomic annotation for each type of outlier signal, and other parameters, as described previously\u003csup\u003e7,8\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eTo evaluate aWatershed, we selected pairs of individuals with the same set of rare variants associated with the same gene (known as \"N2pair\") from the training dataset. We estimated the posterior probability of a functional rare variant in the first individual of the pair and used the outlier status of the second individual as a label for evaluation. We also trained and evaluated the genomic annotation model on each layer of the three transcriptomic signals to determine whether the integration of transcriptomic outlier signals contributes to the prediction of rare functional variants. We compared the results to those obtained from the aWatershed model. After evaluation, we utilized the aWatershed prediction model to calculate posterior probabilities.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003e3\u0026prime;aQTL mapping across 49 GTEx tissues\u003c/h2\u003e \u003cp\u003eGenetic associations between GTEx common variants within 1 Mb of each gene and PEER-corrected APA usage were mapped by Matrix eQTL\u003csup\u003e67\u003c/sup\u003e, as described in our previous study\u003csup\u003e25\u003c/sup\u003e. Known covariates, including sex, RNA integrity number, platform, top five genotype principal components, and unobserved covariates inferred from PEER, were used during 3\u0026prime;aQTL mapping with Matrix eQTL. The number of PEER covariates for each tissue was used as suggested by the GTEx Consortium. We performed 1,000 rounds of permutation to obtain empirical \u003cem\u003eP\u003c/em\u003e-values for each gene, which were then adjusted using the R package qvalue.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003eColocalization analysis between GWAS summary statistics and 3\u0026prime;aQTL\u003c/h2\u003e \u003cp\u003eWe conducted colocalization analysis comparing GWAS summary statistics from the UKBB and literature and 3\u0026prime;aQTL summary data from 49 human tissues using the coloc v.5.1.0.1 package in R\u003csup\u003e68\u003c/sup\u003e. Only GWAS summary data for traits with at least 10,000 cases (binary traits) or 50,000 participates (continuous trait) and with at least 10 SNPs overlapped with aOutlier-associated RVs were kept, which resulted in 1,186 well-powered traits. We extracted the sentinel SNPs for each GWAS trait, defined as GWAS SNPs with \u003cem\u003eP\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;5 \u0026times; 10\u003csup\u003e‒8\u003c/sup\u003e, located at least 1 Mb away from more significant variants. We then searched for colocalizing signals within the 1-Mb region surrounding each sentinel SNP. The coordinates from 3\u0026prime;aQTL summary data were converted from human genome build 38 (hg38) to build 37 (hg19) by CrossMap software\u003csup\u003e69\u003c/sup\u003e to match the version used in all GWAS summary statistics. As defined by the coloc method, five posterior probabilities under five different null hypotheses were calculated. In detail, PP0 represents the null model of no association. PP1 and PP2 represent the probability that causal genetic variants are associated with disease signals or 3\u0026prime;aQTL only. PP3 represents the probability that the genetic effects of trait signals and 3\u0026prime;aQTL are independent, and PP4 represents the probability that trait signals and 3\u0026prime;aQTL data share causal SNPs. The current study classified colocalized events as those with PP4\u0026thinsp;\u0026gt;\u0026thinsp;0.75.\u003c/p\u003e \u003cp\u003e \u003cb\u003e3\u003c/b\u003e\u0026prime; \u003cb\u003eUTR APA transcriptome-wide association study (3\u003c/b\u003e\u0026prime;\u003cb\u003eaTWAS) analysis\u003c/b\u003e\u003c/p\u003e \u003cp\u003eWe used APA quantitative data that was previously used for 3\u0026prime;aQTL mapping\u003csup\u003e25,52,61\u003c/sup\u003e and genotype data of individual genomes from whole genome sequencing (WGS) of GTEx consortia to construct 3\u0026prime;aTWAS model using FUSION\u003csup\u003e70\u003c/sup\u003e for each of the 49 human tissues. To avoid the effects of confounders, well-established factors used in 3\u0026prime;aQTL mapping, including gender, sequencing platform, and other covariates, were incorporated to adjust APA usages. To build the TWAS model, four different models embedded in FUSION were used for weight calculation, including best linear unbiased predictor (blup), elastic-net regression (enet), lasso regression (lasso), and single best eQTL (top1). Subsequently, the cross-validation approach was employed to choose the optimal 3\u0026prime;aTWAS model for each gene. Of note, only genes exhibiting significant heritability estimates (\u003cem\u003ecis-h\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e) (Bonferroni-corrected P\u0026thinsp;\u0026lt;\u0026thinsp;0.05) were retained for subsequent analysis. The built models were then applied to GWAS summary statistics for gene-based association analysis, and a significant association was defined by the FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05. The disease risk genes identified by 3\u0026prime;aTWAS in two or more tissues were used for further analysis.\u003c/p\u003e \u003cdiv id=\"Sec23\" class=\"Section3\"\u003e \u003ch2\u003ePrioritization of trait-associated RVs\u003c/h2\u003e \u003cp\u003eTo determine the frequency with which randomly selected aWatershed-prioritized RVs exhibit larger GWAS effect sizes than matched non-prioritized RVs, we conducted a random sampling test (n\u0026thinsp;=\u0026thinsp;1000) on all RVs using posterior probabilities obtained from the aWatershed prediction model and effect sizes from UKBB GWAS summary statistics. We used aWatershed-prioritized RVs based on aOutlier signals, matched non-prioritized RVs, as well as GWAS effect sizes, gene IDs, and prioritized scores as input data. We defined matched non-prioritized RVs as those with a posterior probability of \u0026lt;\u0026thinsp;0.1 and MAF within \u0026plusmn;\u0026thinsp;0.001 of the selected prioritized RVs in the UKBB cohort.\u003c/p\u003e \u003cp\u003eFor each gene in each trait, we randomly selected one prioritized RV and one matched non-prioritized RV and then identified the one with the largest absolute GWAS effect size in the pair. By summarizing all genes in the trait, we computed the odds of observing a prioritized RV with a larger absolute effect size than a non-prioritized RV across all genes. To generate a null distribution of odds, we repeated this process for matched non-prioritized variants only and randomly selected and compared two non-prioritized RVs for each gene.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec24\" class=\"Section2\"\u003e \u003ch2\u003eCell culture\u003c/h2\u003e \u003cp\u003eHEK293T and MCF7 cells were purchased from the Cell Bank of the Type Culture Collection at the Shanghai Institute of Biochemistry \u0026amp; Cell Biology, Chinese Academy of Science. Cells were maintained in Dulbecco's modified Eagle medium (DMEM; Invitrogen, #11960044) supplemented with 10% fetal bovine serum (Gibco), 100-\u0026micro;g/ml streptomycin, and 100-units/ml penicillin at 37\u0026deg;C in a humidified incubator with 5% CO\u003csub\u003e2\u003c/sub\u003e.\u003c/p\u003e \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e \u003ch2\u003ePlasmid construction\u003c/h2\u003e \u003cp\u003eAll primers used in this study are listed in Supplementary Table\u0026nbsp;8. For intronic APA (IPA) minigenes, the candidate intron and its flanking exons were amplified from genomic DNA as wild-type fragments. For 3\u0026prime; UTR APA minigenes, the 3\u0026prime; UTR of each gene was amplified from genomic DNA as wild-type fragments, and mutations were introduced by PCR-based site-directed mutagenesis. In short, genomic DNA from HEK293T and MCF7 cells was amplified by PCR using primers to generate two 20\u0026ndash;25 bp overlapping fragments containing a mutant site. The IPA wild-type and mutant fragments were subcloned into the EcoRI and BamHI sites of the pcDNA3.1 vector, while 3\u0026prime; UTR APA wild-type and mutant fragments were subcloned into the XhoI and PmeI sites of the mpCHECK2 vector by the One Step Cloning Kit (Vazyme). Two sets of predesigned shRNAs from Sigma against \u003cem\u003eDDX18\u003c/em\u003e were used to clone into pLKO.1-puro vector.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e \u003ch2\u003eTransient transfection\u003c/h2\u003e \u003cp\u003eFor transient transfection, HEK293T and MCF7 cells were plated in a 2-ml culture medium at 6 \u0026times; 10\u003csup\u003e5\u003c/sup\u003e cells/well in six-well plates. After 24 h of culture, cells were transfected with 2 \u0026micro;g of wild-type or mutant minigene plasmid using Lipofectamine 2000 (Invitrogen), according to the manufacturer's instructions. The culture medium was replaced at 6 h post-transfection, and cells were harvested for RNA extraction at 48 h post-transfection. Total RNA was extracted using TRIzol reagent (Invitrogen), according to the manufacturer's instructions, and cDNA was synthesized using the FastKing RT Kit (Tiangen, KR116) with the S-CDS primer. All cDNA was diluted 4-fold in nuclease-free double-distilled H\u003csub\u003e2\u003c/sub\u003eO for further use.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec27\" class=\"Section3\"\u003e \u003ch2\u003e3\u0026prime; RACE\u003c/h2\u003e \u003cp\u003eThe total length of 3\u0026prime; UTR was identified and amplified from the total RNA of NCI-H1299 cells by 3\u0026prime; RACE using the HiScript-TS 5\u0026prime;/3\u0026prime; RACE Kit (Vazyme, RA101) following the manufacturer\u0026rsquo;s protocol. 3\u0026prime; RACE was performed using the S-PCR primer and pcDNA3.1-F or mpCHECK2-F primer to distinguish minigene RNA from endogenous RNA, respectively. The 3\u0026prime; RACE PCR products were separated by gel electrophoresis, and excised bands were purified for Sanger sequencing using the Zymoclean Gel DNA Extraction kit. Cleaned DNA fragments were cloned into the PCE2 vector using the 5 min TA/Blunt-Zero Cloning Kit (Vazyme, C601) and bidirectionally sequenced with M13 forward and reverse primers. At least five colonies were sequenced for every gel product that was purified. Primer sequences are listed in Supplementary Table\u0026nbsp;8.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec28\" class=\"Section2\"\u003e \u003ch2\u003eDual-luciferase reporter assay\u003c/h2\u003e \u003cp\u003eMCF-7 cells were seeded 1 day prior to transfection. The Renilla luciferase in the mpCHECK-2 vector was transfected into cells using Lipofectamine 3000 Transfection Reagent (Invitrogen, cat#: L3000015) according to the manufacturer\u0026prime;s instructions. Forty-eight hours post-transfection, firefly, and renilla luciferase activities were measured by Dual-Luciferase Assay System (Promega, #E1980) on a BioTek Synergy H1 plate reader with full waveband. Each assay was measured in three independent replicates.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec29\" class=\"Section2\"\u003e \u003ch2\u003eCell viability and proliferation assays for shRNA-mediated knockdown\u003c/h2\u003e \u003cp\u003eshRNA-expressing lentivirus was produced with the third-generation packaging system in human embryonic kidney (HEK) 293T cells. For lentivirus infection, target cells (MCF7) were seeded in a 6-well plate 16\u0026ndash;18 h before infection and were grown to 70\u0026ndash;80% confluency upon transduction. The culture medium was removed, and cells were incubated with virus supernatant along with 8 \u0026micro;g/ml polybrene. Puromycin was applied to kill non-infected cells 2days after infection. After two days of selection, when non-infected control cells were all dead, surviving cells were split and maintained with the same concentration of puromycin. Cells were trypsinized, resuspended at 1 \u0026times; 10\u003csup\u003e4\u003c/sup\u003e cells/ml, and seeded in 96-well plates, with each well containing 100ul medium of 1 \u0026times; 10\u003csup\u003e3\u003c/sup\u003e cells. Cell viability and proliferation were determined using CCK8 assays (Yeasen, cat#: 40203ES76) at designated time points (day 1, day 3, day 5, and day 7) by measuring the absorbance at 450 nm, following the manufacturer\u0026prime;s instructions. Values were obtained from three replicate wells for each treatment and time point. Results were representative of three independent experiments.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eThe comprehensive data portal for aOutliers\u003c/h3\u003e\n\u003cp\u003eWe have established a database along with a web interface called rareAPA (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://bioinfo.szbl.ac.cn/rareAPA/index.php\u003c/span\u003e\u003c/span\u003e) on a standard LAMP (Linux\u0026thinsp;+\u0026thinsp;Apache\u0026thinsp;+\u0026thinsp;MySQL\u0026thinsp;+\u0026thinsp;PHP) system, which serves as a comprehensive resource presenting detailed and comprehensive information on rare APA events and their associated RVs. All these data in the rareAPA were stored in MySQL (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ewww.mysql.com\u003c/span\u003e\u003c/span\u003e). The interactive web pages were implemented using HTML, CSS, JavaScript, and PHP languages (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ewww.php.net\u003c/span\u003e\u003c/span\u003e), with several JavaScript libraries (JQuery.js, DataTable.js, and IGV.js) and Bootstrap framework (a popular framework for developing interactive websites) on Red Hat Linux powered by an Apache server (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ewww.apache.org\u003c/span\u003e\u003c/span\u003e). This data portal is valuable for exploring aOutliers and their associated functional rare variants. With rareAPA, users can search, browse, and visualize important information on aOutliers in 49 human tissues. Users can search by gene or tissue name and scrutinize rare APA events among individuals in each tissue. Additionally, users can also visualize aOutliers using a scatter plot or explore them through a genome browser. Furthermore, rareAPA provides a curated list of prioritized RVs using the aWatershed algorithm, allowing users to examine rare variants and their aWatershed posterior scores. Additionally, rareAPA offers batch downloading of all single-tissue aOutliers and multi-tissue aOutliers. The rareAPA is freely available online without registration or login requirements.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eCode availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDaPars2 is available at https://github.com/3UTR/DaPars2, and IPAFinder can be accessed through https://github.com/ZhaozzReal/IPAFinder. The codes for mapping 3\u0026prime;aQTL are available at https://github.com/3UTR/3aQTL-pipe. The custom scripts and source codes for data analysis relevant to this study are available, under the MIT license, at Github repository: https://github.com/Xu-Dong/rareAPA and Zenodo: https://doi.org/10.5281/zenodo.10576656.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe raw data of whole transcriptome and genome sequencing data from the GTEx project V8 are available at the database of Genotypes and Phenotypes (dbGaP) under the accession number:\u0026nbsp;phs000424.v7.p2 [https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000424.v7.p2]\u003csup\u003e71\u003c/sup\u003e. All processed GTEx data, including gene expression outlier (eOutlier) and splicing outlier (sOutlier), are available via the GTEx portal (http://gtexportal.org). GWAS summary statistics used in this study were obtained from UK Biobank GWAS (https://www.nealelab.is/uk-biobank), Finn Gen (https://www.finngen.fi/en), and JENGER (http://jenger.riken.jp). The details about the GWAS summary statistics are listed in Supplementary Table 5. Genomic and functional annotations of rare variants are available via the Combined Annotation Dependent Depletion (CADD v1.5, https://cadd.gs.washington.edu/), and gnomAD v3.1(https://gnomad.broadinstitute.org/). The crosslinking and immunoprecipitation (CLIP) assay data for RNA binding proteins used in this study are available at The Encyclopedia of DNA Elements (ENCODE, https://www.encodeproject.org/). The data described in this study are freely available for querying, visualizing, and downloading at http://bioinfo.szbl.ac.cn/rareAPA/index.php, a website portal dedicated to rare APA.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor Contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eL.L. T.N., and W.L. conceived and supervised the project. X.Z., and Z.Z. performed the bioinformatics analysis with the help from K.X., and H.C. X.Z. constructed the website. Y.C. performed the experiments with the help from Z.W. and S.C., X.Z., T.N., W.L., and L.L. interpreted the data and wrote the manuscript. G.W. and S.X. reviewed and revised the manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting Interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe thank Dr. Jian Yang from Westlake University for providing feedback on the manuscript. We also thank members of the Li laboratory for helpful discussions. This work was supported by the National Natural Science Foundation of China (no. 32100533, 32370721, 32288101, 32030020) and startup funds from Shenzhen Bay Laboratory to L.L. We also thank Qin Wang at the Shenzhen Bay Laboratory Supercomputing Center for high-level computing support and the Medical Science Data Center of Fudan University.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eTaliun, D.\u003cem\u003e et al.\u003c/em\u003e Sequencing of 53,831 diverse genomes from the NHLBI TOPMed Program. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e590\u003c/strong\u003e, 290-299 (2021).\u003c/li\u003e\n\u003cli\u003eKeinan, A. \u0026amp; Clark, A.G. Recent explosive human population growth has resulted in an excess of rare genetic variants. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e336\u003c/strong\u003e, 740-3 (2012).\u003c/li\u003e\n\u003cli\u003eConsortium, U.K.\u003cem\u003e et al.\u003c/em\u003e The UK10K project identifies rare variants in health and disease. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e526\u003c/strong\u003e, 82-90 (2015).\u003c/li\u003e\n\u003cli\u003eNelson, M.R.\u003cem\u003e et al.\u003c/em\u003e An abundance of rare functional variants in 202 drug target genes sequenced in 14,002 people. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e337\u003c/strong\u003e, 100-4 (2012).\u003c/li\u003e\n\u003cli\u003eTennessen, J.A.\u003cem\u003e et al.\u003c/em\u003e Evolution and functional impact of rare coding variation from deep sequencing of human exomes. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e337\u003c/strong\u003e, 64-9 (2012).\u003c/li\u003e\n\u003cli\u003eWang, Q.\u003cem\u003e et al.\u003c/em\u003e Rare variant contribution to human disease in 281,104 UK Biobank exomes. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e597\u003c/strong\u003e, 527-532 (2021).\u003c/li\u003e\n\u003cli\u003eFerraro, N.M.\u003cem\u003e et al.\u003c/em\u003e Transcriptomic signatures across human tissues identify functional rare genetic variation. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e369\u003c/strong\u003e(2020).\u003c/li\u003e\n\u003cli\u003eLi, X.\u003cem\u003e et al.\u003c/em\u003e The impact of rare variation on gene expression across tissues. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e550\u003c/strong\u003e, 239-243 (2017).\u003c/li\u003e\n\u003cli\u003eHernandez, R.D.\u003cem\u003e et al.\u003c/em\u003e Ultrarare variants drive substantial cis heritability of human gene expression. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 1349-1355 (2019).\u003c/li\u003e\n\u003cli\u003eFresard, L.\u003cem\u003e et al.\u003c/em\u003e Identification of rare-disease genes using blood transcriptome sequencing and large control cohorts. \u003cem\u003eNat Med\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 911-919 (2019).\u003c/li\u003e\n\u003cli\u003eMayr, C. What Are 3\u0026apos; UTRs Doing? \u003cem\u003eCold Spring Harb Perspect Biol\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e(2019).\u003c/li\u003e\n\u003cli\u003eTian, B. \u0026amp; Manley, J.L. Alternative polyadenylation of mRNA precursors. \u003cem\u003eNat Rev Mol Cell Biol\u003c/em\u003e \u003cstrong\u003e18\u003c/strong\u003e, 18-30 (2017).\u003c/li\u003e\n\u003cli\u003eMayr, C. Regulation by 3\u0026apos;-Untranslated Regions. \u003cem\u003eAnnu Rev Genet\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 171-194 (2017).\u003c/li\u003e\n\u003cli\u003eBerkovits, B.D. \u0026amp; Mayr, C. Alternative 3\u0026apos; UTRs act as scaffolds to regulate membrane protein localization. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e522\u003c/strong\u003e, 363-7 (2015).\u003c/li\u003e\n\u003cli\u003eDi Giammartino, D.C., Nishida, K. \u0026amp; Manley, J.L. Mechanisms and consequences of alternative polyadenylation. \u003cem\u003eMol Cell\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, 853-66 (2011).\u003c/li\u003e\n\u003cli\u003eMitschka, S. \u0026amp; Mayr, C. Context-specific regulation and function of mRNA alternative polyadenylation. \u003cem\u003eNat Rev Mol Cell Biol\u003c/em\u003e \u003cstrong\u003e23\u003c/strong\u003e, 779-796 (2022).\u003c/li\u003e\n\u003cli\u003eSingh, I.\u003cem\u003e et al.\u003c/em\u003e Widespread intronic polyadenylation diversifies immune cell transcriptomes. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 1716 (2018).\u003c/li\u003e\n\u003cli\u003eZhao, Z.\u003cem\u003e et al.\u003c/em\u003e Cancer-associated dynamics and potential regulators of intronic polyadenylation revealed by IPAFinder using standard RNA-seq data. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e31\u003c/strong\u003e, 2095-2106 (2021).\u003c/li\u003e\n\u003cli\u003eMasamha, C.P.\u003cem\u003e et al.\u003c/em\u003e CFIm25 links alternative polyadenylation to glioblastoma tumour suppression. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e510\u003c/strong\u003e, 412-6 (2014).\u003c/li\u003e\n\u003cli\u003ePark, H.J.\u003cem\u003e et al.\u003c/em\u003e 3\u0026apos; UTR shortening represses tumor-suppressor genes in trans by disrupting ceRNA crosstalk. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e50\u003c/strong\u003e, 783-789 (2018).\u003c/li\u003e\n\u003cli\u003eMittleman, B.E.\u003cem\u003e et al.\u003c/em\u003e Alternative polyadenylation mediates genetic regulation of gene expression. \u003cem\u003eElife\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e(2020).\u003c/li\u003e\n\u003cli\u003eMariella, E., Marotta, F., Grassi, E., Gilotto, S. \u0026amp; Provero, P. The Length of the Expressed 3\u0026apos; UTR Is an Intermediate Molecular Phenotype Linking Genetic Variants to Complex Diseases. \u003cem\u003eFront Genet\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 714 (2019).\u003c/li\u003e\n\u003cli\u003eLi, L., Li, Y., Zou, X., Peng, F., Cui, Y., Wagner, E.J., Li, W. Population-scale genetic control of alternative polyadenylation and its association with human diseases. \u003cem\u003eQuantitative Biology\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 44-54 (2022).\u003c/li\u003e\n\u003cli\u003eGraham, R.R.\u003cem\u003e et al.\u003c/em\u003e Three functional variants of IFN regulatory factor 5 (IRF5) define risk and protective haplotypes for human lupus. \u003cem\u003eProc Natl Acad Sci U S A\u003c/em\u003e \u003cstrong\u003e104\u003c/strong\u003e, 6758-63 (2007).\u003c/li\u003e\n\u003cli\u003eLi, L.\u003cem\u003e et al.\u003c/em\u003e An atlas of alternative polyadenylation quantitative trait loci contributing to complex trait and disease heritability. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e53\u003c/strong\u003e, 994-1005 (2021).\u003c/li\u003e\n\u003cli\u003eFeng, X., Li, L., Wagner, E.J. \u0026amp; Li, W. TC3A: The Cancer 3\u0026apos; UTR Atlas. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e46\u003c/strong\u003e, D1027-D1030 (2018).\u003c/li\u003e\n\u003cli\u003eLiu, Z.\u003cem\u003e et al.\u003c/em\u003e Pan-cancer analysis identifies mutations in SUGP1 that recapitulate mutant SF3B1 splicing dysregulation. \u003cem\u003eProc Natl Acad Sci U S A\u003c/em\u003e \u003cstrong\u003e117\u003c/strong\u003e, 10305-10312 (2020).\u003c/li\u003e\n\u003cli\u003eAlsafadi, S.\u003cem\u003e et al.\u003c/em\u003e Genetic alterations of SUGP1 mimic mutant-SF3B1 splice pattern in lung adenocarcinoma and other cancers. \u003cem\u003eOncogene\u003c/em\u003e \u003cstrong\u003e40\u003c/strong\u003e, 85-96 (2021).\u003c/li\u003e\n\u003cli\u003eChen, E.Y.\u003cem\u003e et al.\u003c/em\u003e Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 128 (2013).\u003c/li\u003e\n\u003cli\u003eMcLaren, W.\u003cem\u003e et al.\u003c/em\u003e The Ensembl Variant Effect Predictor. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 122 (2016).\u003c/li\u003e\n\u003cli\u003eRentzsch, P., Witten, D., Cooper, G.M., Shendure, J. \u0026amp; Kircher, M. CADD: predicting the deleteriousness of variants throughout the human genome. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e47\u003c/strong\u003e, D886-D894 (2019).\u003c/li\u003e\n\u003cli\u003eBogard, N., Linder, J., Rosenberg, A.B. \u0026amp; Seelig, G. A Deep Neural Network for Predicting and Engineering Alternative Polyadenylation. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e178\u003c/strong\u003e, 91-106 e23 (2019).\u003c/li\u003e\n\u003cli\u003eZhao, Z.\u003cem\u003e et al.\u003c/em\u003e Comprehensive characterization of somatic variants associated with intronic polyadenylation in human cancers. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e49\u003c/strong\u003e, 10369-10381 (2021).\u003c/li\u003e\n\u003cli\u003eYeo, G. \u0026amp; Burge, C.B. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. \u003cem\u003eJ Comput Biol\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 377-94 (2004).\u003c/li\u003e\n\u003cli\u003eAlipanahi, B., Delong, A., Weirauch, M.T. \u0026amp; Frey, B.J. Predicting the sequence specificities of DNA- and RNA-binding proteins by deep learning. \u003cem\u003eNat Biotechnol\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 831-8 (2015).\u003c/li\u003e\n\u003cli\u003eJenal, M.\u003cem\u003e et al.\u003c/em\u003e The poly(A)-binding protein nuclear 1 suppresses alternative cleavage and polyadenylation sites. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e149\u003c/strong\u003e, 538-53 (2012).\u003c/li\u003e\n\u003cli\u003eDominguez, D.\u003cem\u003e et al.\u003c/em\u003e Sequence, Structure, and Context Preferences of Human RNA Binding Proteins. \u003cem\u003eMol Cell\u003c/em\u003e \u003cstrong\u003e70\u003c/strong\u003e, 854-867 e9 (2018).\u003c/li\u003e\n\u003cli\u003eLinder, J., Koplik, S.E., Kundaje, A. \u0026amp; Seelig, G. Deciphering the impact of genetic variation on human polyadenylation using APARENT2. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e23\u003c/strong\u003e, 232 (2022).\u003c/li\u003e\n\u003cli\u003eHamosh, A., Scott, A.F., Amberger, J.S., Bocchini, C.A. \u0026amp; McKusick, V.A. Online Mendelian Inheritance in Man (OMIM), a knowledgebase of human genes and genetic disorders. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, D514-7 (2005).\u003c/li\u003e\n\u003cli\u003eSlavotinek, A.M.\u003cem\u003e et al.\u003c/em\u003e Mutation analysis of the MKKS gene in McKusick-Kaufman syndrome and selected Bardet-Biedl syndrome patients. \u003cem\u003eHum Genet\u003c/em\u003e \u003cstrong\u003e110\u003c/strong\u003e, 561-7 (2002).\u003c/li\u003e\n\u003cli\u003eStone, D.L.\u003cem\u003e et al.\u003c/em\u003e Mutation of a gene encoding a putative chaperonin causes McKusick-Kaufman syndrome. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 79-82 (2000).\u003c/li\u003e\n\u003cli\u003eSlavotinek, A.M.\u003cem\u003e et al.\u003c/em\u003e Mutations in MKKS cause Bardet-Biedl syndrome. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 15-6 (2000).\u003c/li\u003e\n\u003cli\u003eKatsanis, N.\u003cem\u003e et al.\u003c/em\u003e Mutations in MKKS cause obesity, retinal dystrophy and renal malformations associated with Bardet-Biedl syndrome. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 67-70 (2000).\u003c/li\u003e\n\u003cli\u003eWuyts, W.\u003cem\u003e et al.\u003c/em\u003e Mutations in the EXT1 and EXT2 genes in hereditary multiple exostoses. \u003cem\u003eAm J Hum Genet\u003c/em\u003e \u003cstrong\u003e62\u003c/strong\u003e, 346-54 (1998).\u003c/li\u003e\n\u003cli\u003eStickens, D.\u003cem\u003e et al.\u003c/em\u003e The EXT2 multiple exostoses gene defines a family of putative tumour suppressor genes. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 25-32 (1996).\u003c/li\u003e\n\u003cli\u003eQuintas-Cardama, A. \u0026amp; Cortes, J. Molecular biology of bcr-abl1-positive chronic myeloid leukemia. \u003cem\u003eBlood\u003c/em\u003e \u003cstrong\u003e113\u003c/strong\u003e, 1619-30 (2009).\u003c/li\u003e\n\u003cli\u003eSalesse, S. \u0026amp; Verfaillie, C.M. BCR/ABL: from molecular mechanisms of leukemia induction to treatment of chronic myelogenous leukemia. \u003cem\u003eOncogene\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, 8547-59 (2002).\u003c/li\u003e\n\u003cli\u003eWeiner, D.J.\u003cem\u003e et al.\u003c/em\u003e Statistical and functional convergence of common and rare genetic influences on autism at chromosome 16p. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e54\u003c/strong\u003e, 1630-1639 (2022).\u003c/li\u003e\n\u003cli\u003eSchrode, N.\u003cem\u003e et al.\u003c/em\u003e Synergistic effects of common schizophrenia risk variants. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, 1475-1485 (2019).\u003c/li\u003e\n\u003cli\u003eSingh, T.\u003cem\u003e et al.\u003c/em\u003e Rare coding variants in ten genes confer substantial risk for schizophrenia. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e604\u003c/strong\u003e, 509-516 (2022).\u003c/li\u003e\n\u003cli\u003eCui, Y.\u003cem\u003e et al.\u003c/em\u003e Alternative polyadenylation transcriptome-wide association study identifies APA-linked susceptibility genes in brain disorders. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 583 (2023).\u003c/li\u003e\n\u003cli\u003eChen, H.\u003cem\u003e et al.\u003c/em\u003e A distinct class of pan-cancer susceptibility genes revealed by alternative polyadenylation transcriptome-wide association study. \u003cem\u003emedRxiv\u003c/em\u003e, 2023.02.28.23286554 (2023).\u003c/li\u003e\n\u003cli\u003eDong, G.\u003cem\u003e et al.\u003c/em\u003e DDX18 drives tumor immune escape through transcription-activated STAT1 expression in pancreatic cancer. \u003cem\u003eOncogene\u003c/em\u003e \u003cstrong\u003e42\u003c/strong\u003e, 3000-3014 (2023).\u003c/li\u003e\n\u003cli\u003eRedmond, A.M.\u003cem\u003e et al.\u003c/em\u003e Genomic interaction between ER and HMGB2 identifies DDX18 as a novel driver of endocrine resistance in breast cancer cells. \u003cem\u003eOncogene\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, 3871-80 (2015).\u003c/li\u003e\n\u003cli\u003eMcFarland, J.M.\u003cem\u003e et al.\u003c/em\u003e Improved estimation of cancer dependencies from large-scale RNAi screens using model-based normalization and data integration. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 4610 (2018).\u003c/li\u003e\n\u003cli\u003eTsherniak, A.\u003cem\u003e et al.\u003c/em\u003e Defining a Cancer Dependency Map. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e170\u003c/strong\u003e, 564-576 e16 (2017).\u003c/li\u003e\n\u003cli\u003eDemontis, D.\u003cem\u003e et al.\u003c/em\u003e Genome-wide analyses of ADHD identify 27 risk loci, refine the genetic architecture and implicate several cognitive domains. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e55\u003c/strong\u003e, 198-208 (2023).\u003c/li\u003e\n\u003cli\u003eWu, N.\u003cem\u003e et al.\u003c/em\u003e TBX6 null variants and a common hypomorphic allele in congenital scoliosis. \u003cem\u003eN Engl J Med\u003c/em\u003e \u003cstrong\u003e372\u003c/strong\u003e, 341-50 (2015).\u003c/li\u003e\n\u003cli\u003eDobin, A.\u003cem\u003e et al.\u003c/em\u003e STAR: ultrafast universal RNA-seq aligner. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 15-21 (2013).\u003c/li\u003e\n\u003cli\u003eQuinlan, A.R. \u0026amp; Hall, I.M. BEDTools: a flexible suite of utilities for comparing genomic features. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 841-2 (2010).\u003c/li\u003e\n\u003cli\u003eZou, X.\u003cem\u003e et al.\u003c/em\u003e Using population-scale transcriptomic and genomic data to map 3\u0026apos; UTR alternative polyadenylation quantitative trait loci. \u003cem\u003eSTAR Protoc\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 101566 (2022).\u003c/li\u003e\n\u003cli\u003eMa, X.\u003cem\u003e et al.\u003c/em\u003e ipaQTL-atlas: an atlas of intronic polyadenylation quantitative trait loci across human tissues. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e51\u003c/strong\u003e, D1046-D1052 (2023).\u003c/li\u003e\n\u003cli\u003eStegle, O., Parts, L., Piipari, M., Winn, J. \u0026amp; Durbin, R. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses. \u003cem\u003eNat Protoc\u003c/em\u003e \u003cstrong\u003e7\u003c/strong\u003e, 500-7 (2012).\u003c/li\u003e\n\u003cli\u003eGudmundsson, S.\u003cem\u003e et al.\u003c/em\u003e Variant interpretation using population databases: Lessons from gnomAD. \u003cem\u003eHum Mutat\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, 1012-1030 (2022).\u003c/li\u003e\n\u003cli\u003eWang, R., Zheng, D., Yehia, G. \u0026amp; Tian, B. A compendium of conserved cleavage and polyadenylation events in mammalian genes. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 1427-1441 (2018).\u003c/li\u003e\n\u003cli\u003eWang, R., Nambiar, R., Zheng, D. \u0026amp; Tian, B. PolyA_DB 3 catalogs cleavage and polyadenylation sites identified by deep sequencing in multiple genomes. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e46\u003c/strong\u003e, D315-D319 (2018).\u003c/li\u003e\n\u003cli\u003eShabalin, A.A. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 1353-8 (2012).\u003c/li\u003e\n\u003cli\u003eGiambartolomei, C.\u003cem\u003e et al.\u003c/em\u003e Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. \u003cem\u003ePLoS Genet\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, e1004383 (2014).\u003c/li\u003e\n\u003cli\u003eZhao, H.\u003cem\u003e et al.\u003c/em\u003e CrossMap: a versatile tool for coordinate conversion between genome assemblies. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 1006-7 (2014).\u003c/li\u003e\n\u003cli\u003eGrishin, D. \u0026amp; Gusev, A. Allelic imbalance of chromatin accessibility in cancer identifies candidate causal risk variants and their mechanisms. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e54\u003c/strong\u003e, 837-849 (2022).\u003c/li\u003e\n\u003cli\u003eConsortium, G.T. The GTEx Consortium atlas of genetic regulatory effects across human tissues. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e369\u003c/strong\u003e, 1318-1330 (2020).\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"nature-portfolio","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"","title":"Nature Portfolio","twitterHandle":"","acdcEnabled":false,"dfaEnabled":false,"editorialSystem":"ejp","reportingPortfolio":"","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-3907149/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3907149/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eAlthough rare non-coding variants (RVs) play crucial roles in human complex traits and diseases, understanding their functional mechanisms and identifying those most closely associated with diseases continue to be major challenges. Here, we constructed the first comprehensive atlas of alternative polyadenylation (APA) outliers (aOutliers) from 15,201 samples across 49 human tissues. Strikingly, these aOutliers exhibit unique characteristics markedly distinct from those of outliers based on transcriptional abundance or splicing. This is evidenced by a pronounced enrichment of RVs specifically within aOutliers. Mechanistically, aOutlier RVs frequently alter poly(A) signals and splicing sites, and experimental perturbation of these RVs indeed triggers APA events. Furthermore, we developed a Bayesian-based APA RV prediction model, which successfully pinpointed a specific set of RVs with significantly large effect sizes on complex traits or diseases. A particularly intriguing discovery was the observed convergence effect on APA between rare and common cancer variants, exemplified by the combinatorial regulation of APA in the \u003cem\u003eDDX18\u003c/em\u003e gene. Together, this study introduces a novel APA-enhanced framework for individual genome annotation and underscores the importance of APA in uncovering previously unrecognized functional non-coding RVs linked to human complex traits and diseases.\u003c/p\u003e","manuscriptTitle":"Impact of Rare Non-coding Variants on Human Diseases through Alternative Polyadenylation Outliers","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-03-07 18:54:12","doi":"10.21203/rs.3.rs-3907149/v1","editorialEvents":[],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"nature-communications","isNatureJournal":true,"hasQc":false,"allowDirectSubmit":false,"externalIdentity":"NCOMMS","sideBox":"Learn more about [Nature Communications](http://www.nature.com/ncomms/)","snPcode":"","submissionUrl":"https://mts-ncomms.nature.com/","title":"Nature Communications","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"ejp","reportingPortfolio":"Nature Communications","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"da086e35-e28a-441c-b114-8679abc7ed12","owner":[],"postedDate":"March 7th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":29155987,"name":"Biological sciences/Genetics/Gene regulation"},{"id":29155988,"name":"Biological sciences/Computational biology and bioinformatics/Data mining"}],"tags":[],"updatedAt":"2025-01-17T08:05:27+00:00","versionOfRecord":{"articleIdentity":"rs-3907149","link":"https://doi.org/10.1038/s41467-024-55407-3","journal":{"identity":"nature-communications","isVorOnly":false,"title":"Nature Communications"},"publishedOn":"2025-01-16 05:00:00","publishedOnDateReadable":"January 16th, 2025"},"versionCreatedAt":"2024-03-07 18:54:12","video":"","vorDoi":"10.1038/s41467-024-55407-3","vorDoiUrl":"https://doi.org/10.1038/s41467-024-55407-3","workflowStages":[]},"version":"v1","identity":"rs-3907149","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3907149","identity":"rs-3907149","version":["v1"]},"buildId":"7rjqhiLT3MXkJMwkYKINL","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}