RESEARCH ARTICLE | GENETICS OPEN ACCESS Genome-wide detection of human variants that disrupt intronic branchpoints Peng Zhanga,1 , Quentin Philippotb,c, Weicheng Rend, Wei-Te Leia, Juan Lia, Peter D. Stensone, Pere Soler Palacınf,g,h, Roger Colobranh,i,j, Bertrand Boissona,b,c , Shen-Ying Zhanga,b,c, Anne Puela,b,c , Qiang Pan-Hammarstro € md,2, Qian Zhanga,b,c,2, David N. Coopere,2 , Laurent Abela,b,c,3, and Jean-Laurent Casanovaa,b,c,k,1,3 Contributed by Jean-Laurent Casanova; received June 30, 2022; accepted September 14, 2022; reviewed by William Fairbrother and Stylianos Antonarakis Pre-messenger RNA splicing is initiated with the recognition of a single-nucleotide intronic branchpoint (BP) within a BP motif by spliceosome elements. Forty-eight rare Significance variants in 43 human genes have been reported to alter splicing and cause disease by disrupting BP. However, until now, no computational approach was available to effi- The search for candidate variants ciently detect such variants in massively parallel sequencing data. We established a com- underlying human disease in prehensive human genome-wide BP database by integrating existing BP data and massive parallel sequencing data generating new BP data from RNA sequencing of lariat debranching enzyme DBR1- typically focuses on coding regions mutated patients and from machine-learning predictions. We characterized multiple and essential splice sites, mostly features of BP in major and minor introns and found that BP and BP-2 (two nucleoti- ignoring noncoding variants. The des upstream of BP) positions exhibit a lower rate of variation in human populations RNA spliceosome recognizes and higher evolutionary conservation than the intronic background, while being com- intronic branchpoint (BP) motifs at parable to the exonic background. We developed BPHunter as a genome-wide compu- the beginning of splicing and tational approach to systematically and efficiently detect intronic variants that may operates mostly within introns to disrupt BP recognition. BPHunter retrospectively identified 40 of the 48 known patho- genic BP variants, in which we summarized a strategy for prioritizing BP variant candi- define the exon–intron dates. The remaining eight variants all create AG-dinucleotides between the BP and boundaries; however, BP variants acceptor site, which is the likely reason for missplicing. We demonstrated the practical have been paid little attention. We utility of BPHunter prospectively by using it to identify a novel germline heterozygous established a comprehensive BP variant of STAT2 in a patient with critical COVID-19 pneumonia and a novel genome-wide database and somatic intronic 59-nucleotide deletion of ITPKB in a lymphoma patient, both of knowledgebase of BP and which were validated experimentally. BPHunter is publicly available from https:// developed BPHunter for hgidsoft.rockefeller.edu/BPHunter and https://github.com/casanova-lab/BPHunter. systematic and informative genome-wide detection of intronic branchpoint j splicing j intronic variant j disease genetics j software variants that may disrupt BP and Pre-mRNA (messenger RNA) splicing is a necessary step for protein-coding gene splicing, together with an effective expression in eukaryotic cells. The biochemical process is initiated by the recognition strategy for prioritizing BP variant of the branchpoint (BP) in an intron, followed by the identification and ligation of the candidates. BPHunter not only 50 splice site (50 ss = donor site) and 30 splice site (30 ss = acceptor site) to join two constitutes an important resource exons and the removal of the intervening intron as a circular lariat, to yield a mature for understanding BP, but should mRNA of coding sequence for protein translation (Fig. 1A). Splicing takes place in the also drive discovery of BP variants nucleus and is regulated by an array of cis-acting elements (the RNA sequences with in human genetic diseases and their splicing codes) and trans-acting elements (proteins and small nuclear RNAs traits. (snRNAs) that bind to the cis-acting elements), which together constitute a garden around the splice sites promoting the formation of the complex and dynamic spliceo- some (1–3). When splicing is initiated, the BP motif is recognized by a spliceosomal snRNA, with the cooperation of splicing factors (4, 5). BP recognition is achieved dif- ferently between the two types of spliceosome (Fig. 1B): the major (U2-type, in >99% of human introns) and the minor (U12-type). The main difference between these spli- Reviewers: W.F., Brown University; and S.A., University ceosomes is the use of different snRNAs. The major spliceosome recruits U1, U2, U4, of Geneva Medical School. U5, and U6 snRNAs, in which the U1 and U2 snRNAs recognize 50 ss and BP inde- The authors declare no conflict of interests. pendently (1). In the interaction between BP and U2 snRNA, the BP nucleotide bulges Copyright © 2022 the Author(s). Published by PNAS. out to bind into a pocket formed by SF3B1 and PHF5A (6), and its flanking sequences This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY). base-pair with U2 snRNA (5, 7), while U2 snRNA is stabilized by SF3B6 (6). The 1 To whom correspondence may be addressed. Email: minor spliceosome recruits U11, U12, U4atac, U5, and U6atac snRNAs, in which the pzhang@rockefeller.edu or casanova@rockefeller.edu. U11 and U12 snRNAs first form a di-snRNA complex before recognizing the 50 ss and 2 Q.P.-H., Q.Z., and D.N.C. contributed equally to this work. BP simultaneously (7, 8). 3 L.A. and J.-L.C. contributed equally to this work. BP constitutes a single-nucleotide cis-acting element within the intron. It is typically an adenine and generally located [20, 40] (from 20 to 40) nucleotides (nt) away This article contains supporting information online at http://www.pnas.org/lookup/suppl/doi:10.1073/pnas. from the 30 ss. The BP motif is thought to follow the consensus sequence YTNAY (4) (A 2211194119/-/DCSupplemental. denotes the BP site, Y = C/T, n = any nucleotide). There are also distal BP deeper inside Published October 28, 2022. PNAS 2022 Vol. 119 No. 44 e2211194119 https://doi.org/10.1073/pnas.2211194119 1 of 12 Spliceosome iCLIP experimentation identifies cross-links of spli- ceosomal factors on pre-mRNA sequences at nucleotide-level res- olution to pinpoint BP positions (18). Some machine-learning models have also been developed to predict BP positions based on one or two experimentally identified BP datasets (19–21). However, these individual BP studies have yet to be integrated into a comprehensive database and knowledgebase of human BP. The functional role of splicing makes it a critical process for the correct assembly of gene products, while its complexity ren- ders it vulnerable to deleterious variants that can perturb any part of the splicing machinery. A variant may affect a cis- or trans-element, leading to altered splicing with alternative 50 ss/30 ss, exon skipping, or intron retention, any of which could result in the synthesis of a defective gene product (22, 23). As previously reported (24–26), at least 10% of disease-causing variants are associated with missplicing [of which ∼65% are in essential exonic/intronic splice sites, ∼25% are in the remainder of the intron, and ∼10% are in the remainder of the exon (25)]; this may, however, be a serious underestimate (24) owing to our rather limited ability to detect variants that impair different aspects of the spliceosome. Currently, the search for candidate pathogenic variants in massively parallel sequence (MPS) data generally focuses on nonsynonymous variants located within cod- ing sequences or variants in essential splice sites, while mostly ignoring synonymous coding exonic variants and noncoding intronic variants. Data from RNA-seq or other studies of mRNA occasionally lead to the detection of intronic or exonic variants that disrupt splicing (27). However, the spliceosome operates mostly within intronic regions to define the exon–intron bound- aries and hence the coding sequences. It follows that introns probably harbor a substantially larger number of pathogenic var- iants than has so far been appreciated (28). The BP, involved in the initiation of splicing through its interactions with spliceosome elements, constitutes a key vulnerability of splicing by virtue of its potential variants. Employing the Human Gene Mutation Database (HGMD) Fig. 1. Schematics of a splicing event, the two types of spliceosome, and Professional version (25) in concert with an in-depth literature typical consequences of a BP variant. (A) Schematic of the biological pro- review, we identified 48 BP variants reported in 43 genes cesses of transcription, splicing, and translation from DNA to pre-mRNA to underlying 33 human diseases whose pathogenicity has been mature mRNA and protein. In pre-mRNA splicing, the BP is first recognized, supported by accompanying experimental evidence. These BP two exons are then joined, and the intervening intron is finally released as a circular lariat. (B) Schematic of the major (U2-type, on the left) and minor variants resulted in complete/partial exon skipping or intron (U12-type, on the right) spliceosomes, with an illustration of the interaction retention, or a combination of these consequences (Fig. 1C). between pre-mRNA sequence and U2/U12 snRNA. (C) Schematic of the The misspliced mRNAs may be rapidly degraded by nonsense- potential molecular consequences of a BP variant, including complete/ partial exon skipping and complete/partial intron retention. The percentage mediated decay (NMD), resulting in loss of expression, and the in parentheses refer to the observed fraction of each category from the translated truncated proteins could exhibit loss of function, published pathogenic BP variants (note that one variant could result in both being potentially disease-causing. However, since the first more than one missplicing consequence). report of a disease-causing BP variant in the L1CAM gene in a family with X-linked hydrocephalus in 1992 (29), the discovery introns: recursive splicing within the intron for multistep intron rate of pathogenic BP variants appears to have been rather low removal (9), stem-loop RNA structures bring distal BP closer to (only 48 variants discovered in the last 30 y). This is probably 30 ss (10), stochastic splice site selection leads to kinetic variation because 1) genome-wide knowledge of BP has not been system- in intron removal (11) (SI Appendix, Fig. S1). So far, consider- atically integrated and characterized, 2) genome-wide detection able efforts have been made to identify BP positions, including a software has not been available for the efficient identification of lariat sequencing method to study BP in Saccharomyces cerevisiae BP variants from MPS data, and 3) in-depth understanding of (12), a spliceosome profiling method to sequence RNA copurify- all reported pathogenic BP variants has not been analyzed, ing with spliceosomal complex in Schizosaccharomyces pombe (13), which could potentiate the strategic prioritization of candidate and analyzing large-scale experimental data in humans [RNA variants that could functionally disrupt BP. We have therefore sequencing (RNA-seq) (5, 14–17) and iCLIP (18)]. To locate aimed to fill these gaps both for biological and medical reasons. BP positions from RNA-seq data, the sequencing reads from the circular lariats that traverse the junctions between 50 ss and Results BP are identified by aligning to the intronic sequences of 50 ss and 30 ss, respectively. Since lariats are transient and only exist Collection of Experimentally Identified BP Data. We first col- at low abundance, their detection is technically difficult, lated five datasets of experimentally identified BP (eBP), based on hence large-scale data are required for identifying BP positions. RNA-seq and iCLIP data: eBP_Mercer (59,359 BP from total 2 of 12 https://doi.org/10.1073/pnas.2211194119 pnas.org Table 1. A comprehensive collection of eBP and cBP, including the additional BP data generated in this study (eBP_BPHunter and cBP_BPHunter), followed by consensus-guided positional adjustment No. of No. of Δ no. BP dataset raw BP A, % C, % G, % T, % adj. BP A, % C, % G, % T, % of BP eBP_Mercer 59,359 78.5 8.5 4.7 8.4 54,531 84.9 8.0 3.5 3.6 4,828 eBP_Taggart 27,795 67.8 15.3 5.5 11.4 27,635 73.7 13.1 4.1 9.1 160 eBP_Pineda 138,313 75.0 7.6 9.5 7.9 133,652 77.0 7.3 8.9 6.9 4,661 eBP_Talhouarne 240 4.6 84.2 4.6 6.7 240 8.3 81.7 3.3 6.7 0 eBP_Briese 43,630 100.0 0.0 0.0 0.0 43,630 100.0 0.0 0.0 0.0 0 eBP_BPHunter 8,682 45.9 22.4 18.2 13.5 7,990 56.6 21.8 15.7 5.8 692 eBP 210,986 73.4 9.1 8.5 9.1 198,256 76.9 8.6 7.7 6.9 12,730 cBP_BPP 232,828 98.8 1.1 0.0 0.2 232,820 98.8 1.1 0.0 0.2 8 cBP_Branchpointer 352,322 98.8 0.5 0.0 0.7 325,304 99.5 0.6 0.0 0.0 27,018 cBP_LaBranchoR 203,914 98.3 1.6 0.0 0.0 203,907 98.3 1.6 0.0 0.0 7 cBP_BPHunter 20,718 100.0 0.0 0.0 0.0 20,718 100.0 0.0 0.0 0.0 0 cBP 479,366 97.9 1.5 0.0 0.6 451,215 98.3 1.6 0.0 0.1 28,151 Total BP 546,559 90.40 4.25 2.78 2.57 mBP 102,912 99.15 0.83 0.01 0.01 RNA-seq) (5), eBP_Taggart (27,795 BP from total RNA-seq) machine-learning model: a set of three algorithms (gradient (15), eBP_Pineda (138,313 BP from total RNA-seq) (16), boost machine, random forest, and logistic regression) that are eBP_Talhouarne (240 BP from cytoplasmic RNA-seq) (17), tuned for high precision and integrated with majority voting and eBP_Briese (43,630 BP from iCLIP) (18) (Table 1). In for the final prediction (Methods). These three models were addition to the published BP datasets, we revisited our RNA-seq independently trained on the 198,256 consensus-guided posi- data derived from DBR1-mutated patients with brainstem viral tion-adjusted eBP positions (discussed in the next section; SI infection (30). DBR1 encodes the only known lariat debranching Appendix, Fig. S2 B and C) versus 1 million random intronic/ enzyme, and the patients with DBR1 variants, who were charac- exonic positions and then optimized by tuning parameters and terized by a residual DBR1 activity of 3 to 10%, have been thresholds through cross-validation (SI Appendix, Fig. S2D). We reported to exhibit an elevated cellular level of lariats (30). We used the model to predict BP in [3, 40] nt of 30 ss, and it designed a computational program to detect BP positions from yielded 20,718 BP sites (cBP_BPHunter). In these four cBP RNA-seq reads (Methods), which identified 8,682 BP sites datasets, the most abundant BP nucleotides were cBP_BPP (eBP_BPHunter) from 15 RNA-seq samples from the fibroblasts (99% A), cBP_Branchpointer (99% A), cBP_LaBranchoR (98% of three DBR1-mutated patients (SI Appendix, Fig. S2A and A), and cBP_BPHunter (100% A) (Table 1), showing that the Table S1). In these six eBP datasets, the most abundant BP prediction models favored adenine-BP. In total, we obtained nucleotides were eBP_Mercer (79% A), eBP_Taggart (68% A), 479,366 unique cBP. Please note that cBP does not denote eBP_Pineda (75% A), eBP_Talhouarne (84% C), eBP_Briese cytosine-BP. When we mention the nucleotides of BP, we fully (100% A), and eBP_BPHunter (46% A) (Table 1). This showed spell out the nucleotide (e.g., adenine-BP, cytosine-BP). that adenine-BP were the most frequent in lariats from total RNA- seq, whereas cytosine-BP were the most frequent in the stable Consensus-Guided Positional Adjustment of BP and Integration lariats in cytoplasm albeit in very low amounts (17). DBR1- of BP Datasets. We tested the established BP consensus sequence mutated patients had 22% stable cytosine-BP. The iCLIP- YTNAY (4) in eBP but found only 23.5% eBP matched this pat- identified BP were invariably adenine. In total, we obtained tern precisely. Considering the noncanonical 20 -to-50 linkage 210,986 unique eBP. between 50 ss and BP, and the possibility of mutations being intro- duced by reverse transcriptase when traversing the 50 ss-BP junc- Collection of Computationally Predicted BP Data. We also col- tion in RNA-seq (5, 15), we anticipated that a number of eBP lated three datasets of computationally predicted BP (cBP): may have been mislocated in the raw dataset (SI Appendix, Fig. cBP_BPP (232,828 BP) (19), cBP_Branchpointer (352,322 S2B), which could also have led to some mislocated cBP predic- BP) (20), and cBP_LaBranchoR (203,914 BP) (21) (Table 1). tions. We therefore screened a window of [2, +2] nt from each These studies trained machine-learning models based on one/ BP position for consensus sequence matching and adjusted the two eBP datasets and performed predictions in the region raw BP position to its closest neighbor that perfectly matched [21, 34], [18, 44] and [1, 70] nt of 30 ss respec- YTNAY (SI Appendix, Fig. S2C). As a result, about 6% of raw tively. These cBP datasets supplemented the BP positions that BP positions were adjusted (SI Appendix, Table S2), and multiple were not detected from experimental data and hence increased neighboring raw BP positions were merged into a single adjusted the BP coverage across the human genome. These computa- position, yielding 198,256 eBP (Δ = 12,730) and 451,215 cBP tional methods overlapped the predictions in [21, 34] nt of (Δ = 28,151) after adjustment (Table 1). By integrating eBP 30 ss, with less attention being paid to in the region closer to 30 ss, and cBP, we assembled a comprehensive collection of 546,559 BP which is generally considered to comprise polypyrimidine tracts in the human genome, with 102,912 mutually shared BP (mBP) (PPT). However, the use of a hard distance cutoff for PPT or (Fig. 2A, Table 1, SI Appendix, Fig. S3 and Table S3, and the use of an algorithm to estimate PPT to skip BP prediction Dataset S1). This BP database contained 90.4% adenine-BP, 4. could miss potentially genuine BP sites. For this reason, we 25% cytosine-BP, 2.78% guanine-BP, and 2.57% thymine- postulated that there could be additional BP candidates in the BP. We further decomposed adenine-BP by adopting increas- region closer to 30 ss. Therefore, we developed a high-precision ingly relaxed consensus sequences: YTNAY (35.8%), YTNA PNAS 2022 Vol. 119 No. 44 e2211194119 https://doi.org/10.1073/pnas.2211194119 3 of 12 found that these BP were largely located in proximity. By taking the first BP in each intron as the reference (i.e., the one closest to the 30 ss), 31% and 59% of nonfirst BP were found to be located within 5 nt and 10 nt upstream of the first BP (SI Appendix, Fig. S7). With multiple BP, the first BP could be the most important, as one study showed that BP closer to 30 ss were much more likely to influence splicing efficiency (31). The BP usage is also expected to be context-specific in alternatively spliced transcripts, but we do not have the requisite data to dif- ferentiate the usage. Thus, in this study, we considered the BP–intron associations in a general setting. Nucleotide Composition, BP-30 ss Distance, and BP-U2/U12 snRNA Binding Energy in Major and Minor Introns. In major introns, BP motifs exhibited high conservation of A (92%) at pos:0 and T (72%) at pos:2 and high pyrimidine content (C/T) at pos:2 (84%), pos: 3 (69%), and pos:+1 (66%). In minor introns, BP showed high conservation of A (95%) at pos:0 and T (78%) at pos: 2 and high pyrimidine content at pos:[2, 5] and pos:+1 with frequencies 89%, 75%, 67%, 66%, and 67%, respectively (Fig. 3A and SI Appendix, Table S4). We found 2% and 67% of BP in major introns were distributed in regions [3, 15] and [16, 30] nt of 30 ss, whereas 27% and 43% of BP in minor introns were found within these two regions (Fig. 3B and SI Appendix, Table S5). BP in minor introns were generally located closer to 30 ss than BP in major introns, which is consistent with the major spliceosome requiring the binding of U2AF65 with PPT that needs a larger spacing, which Fig. 2. Integration of all BP datasets and mapping of BP onto introns. (A) The nucleotide composition displays the motif [9, +3] of BP, where the locations of BP are marked with a red background. (B) BP were mapped to 30 -proximal intron. (C) Distance from the mapped BP to their corresponding 30 ss. (D) The number of BP mapped to each 30 -proximal intron. (53.65%), TNA (74.32%), and YNA (77.24%) (SI Appendix, Fig. S4). We then examined the region [50, +20] nt of BP for motif enrichment and reported 14 enriched motifs (P < 0.05, occurrence >5%) in the 50-nt upstream region of BP, which could represent cis-acting elements with a role in facilitating splicing (SI Appendix, Fig. S5). Mapping BP to Introns. We identified 200,059 unique introns from 41,975 transcripts of 17,372 protein-coding genes and classified introns into 199,393 (99.7%) major and 666 (0.3%) minor types (SI Appendix, Fig. S6 and Methods). Because the major and minor spliceosomes regulate splicing differently (7) (Fig. 1B), it is important to distinguish intron types when study- ing BP and splicing. Major introns were dominated by GT-AG (98.9%) as the 50 ss–30 ss terminal dinucleotides, whereas 71% and 24% of minor introns were GT-AG and AT-AC, respec- tively. We mapped 546,559 BP onto 200,059 introns, retained the BP within [3, 100] nt of 30 ss (30 -proximal introns), and finally obtained 386,209 BP associated with 200,059 introns (Fig. 2B and Datasets S2 and S3). We excluded the distal BP in this study, as more complicated splicing mechanisms may occur (SI Appendix, Fig. S1). Among the intron-associated BP, 57% of Fig. 3. Characterization of BP. (A) Nucleotide composition of BP, 50 ss, and 30 ss motifs in major and minor introns. (B) The distance from BP to 30 ss in them were located within a 10-nt window [21, 30] nt of major and minor introns. (C) The binding energy between BP motif and U2/ 30 ss, and 88% were located within a 25-nt window [15, 40] U12 snRNA in major and minor introns (the lower the energy, the higher the nt of 30 ss (Fig. 2C); 32%, 41%, and 18% of 30 -proximal introns binding affinity). (D) The proportion of each genomic position harboring pop- had one, two, and three BP mapped, respectively, while 9% of ulation variants (Upper), and the MAF distribution of population variants (Lower). (E) The distribution of the conservation scores GERP (Left) and PhyloP 30 -proximal introns harbored more than three BP (Fig. 2D). (Right) in each genomic position. (F) The cross-compared conservation scores With 68% of 30 -proximal introns containing multiple BP, we between BP and BP2 positions. 4 of 12 https://doi.org/10.1073/pnas.2211194119 pnas.org is absent in the minor spliceosome (7). We further estimated the positions (fourth quadrant). Thus, if we consider at least one binding energy between BP and snRNA and found that BP in positive score for BP or BP2 position, 86% of BP regions major introns displayed weaker binding with U2-snRNA than were evolutionarily conserved. Taken together, these findings that between BP and U12-snRNA in minor introns (Fig. 3C and suggested that the BP region is more highly conserved than SI Appendix, Table S5). This may imply that the BP–U2 interac- the intronic background but less conserved than the exonic tion in the major spliceosome can be complemented by other background and splice sites. It may also hint at a possible col- spliceosomal elements (e.g., SF3B1, SF3B6, PHF5A, SFA2, laborative role between BP and BP2 to ensure a certain SFA3) (6), and hence a weaker wobble BP–U2 base pairing could level of conservation. An extended window around splice sites is be acceptable. Since the minor spliceosome forms a U11–U12 shown in SI Appendix, Fig. S9. A recent study of the bovine di-snRNA complex to bind to 50 ss and BP simultaneously, it genome showed that bovine BP are under evolutionary con- requires more perfect BP–U12 base pairing, which is also straint (35). It is worth noting that splicing evolves rapidly reflected in more conserved 50 ss and BP motifs in minor introns between species (36), with one study showing significant splicing (Fig. 3A). differences between primates, including human and chimpanzee (37). This suggests that the conservation scores of splicing regu- Population Variants in BP. We obtained human population var- latory elements may not be invariably correlated with their iants from the gnomAD database r3.1.1 (32), with which to functional importance. In other words, BP with low/negative compare the occurrence of population variants in the BP region, conservation scores could be functionally important. essential/flanking splice sites, and intronic/exonic background. We focused on the [3, +1] positions of 386,209 BP, took a Comparative Characterization of BP. We made seven compari- 4-nt window for each splice site ([1, +3] of 50 ss and [3, +1] sons of BP in the major introns, at intron, transcript, gene, and of 30 ss) from 200,059 introns and randomly sampled 1 million chromosome levels: 1) 353,031 adenine-BP vs. 14,476 cytosine- positions from the intron and exon backgrounds. In this way, we BP vs. 9,626 guanine-BP vs. 7,381 thymine-BP (SI Appendix, Fig. generated 15 sets of genomic locations for comparison with their S10). We found that adenine-/thymine-BP had much lower rates population variants (Fig. 3D and SI Appendix, Table S6). We of population variation than cytosine-/guanine-BP, and they had found that 17.4% of intronic and 16% of exonic background two spikes of high PhyloP conservation scores. 2) 94,581 mBP vs. positions harbored population variants; 13.9 to 15.1% of the 60,450 exclusively eBP vs. 229,483 exclusively cBP (SI Appendix, flanking splice sites had population variants, whereas the essential Fig. S11). mBP displayed the most conserved nucleotide composi- splice sites exhibited a much lower proportion (6.4 to 10.7%). tion, the lowest rate of population variation, and the largest pro- Among the BP regions, 14.5% of BP positions, 15.4% of BP1 portion of positive conservation scores, followed by exclusively cBP positions, and 12.7% of BP2 positions harbored population and then exclusively eBP. 3) 180,766 first BP vs. 203,748 nonfirst variants, which were lower than those for the intron and exon BP (SI Appendix, Fig. S12) had no distinct feature noted. 4) backgrounds and similar to the flanking positions of splice sites. 57,254 BP in 30 -proximal introns with single BP vs. 327,260 BP The other two positions in the BP region (BP3 and BP+1) had in 30 -proximal introns with multiple BP (SI Appendix, Fig. S13). 17.1% and 18.7% of population variants respectively, which BP in introns of single BP showed an absolute dominance (96%) were close to or even higher than the intronic background. The of T in the BP2 position, a lower rate of population variation, minor allele frequencies (MAF) of these population variants at and a larger proportion of positive conservation scores. 5) 335,617 each position also suggested a similar pattern (Fig. 3D). An BP in introns without alternative splicing vs. 48,897 BP in introns extended window around splice sites ([2, +5] of 50 ss and [5, with alternative splicing (SI Appendix, Fig. S14) had no notable +2] of 30 ss) is shown in SI Appendix, Fig. S8. These findings sug- difference. 6) 148,571 BP in genes encoding a single isoform vs. gested that the BP and BP2 positions exhibit a lower rate and a 235,943 BP in genes encoding multiple isoforms (SI Appendix, lower frequency of population variations than the intronic back- Fig. S15) had no notable difference. 7) 13,876 BP in genes located ground and are similar to the exonic background and flanking on the sex (X and Y) chromosomes vs. 370,638 BP in genes positions of essential splice sites. located on the autosomes (SI Appendix, Fig. S16). The sex chro- mosomes had a slightly higher proportion of 30 -proximal Cross-Species Conservation of BP. We evaluated genome-wide introns with single BP (36% vs. 31%) than the autosomes. conservation scores, GERP (33) and PhyloP (34), in the above- We found that 8.7% of BP on sex chromosomes and 15% of mentioned 15 sets of genomic locations. A positive score BP on autosomes harbored population variants, revealing a denotes that the position is likely to be evolutionarily con- lower rate of population variation of BP on sex chromosomes served, whereas a negative score indicates that the position is than autosomes. This may be due to the three population probably evolving neutrally (Fig. 3E and SI Appendix, Table copies of X chromosomes (two copies in females, but one S7). Both scores were found to be comparable for each position copy in males) compared to four population copies of auto- studied. The essential splice sites and their exonic flanking somes (two copies in both females and males) and some other regions displayed an absolute positive distribution of conserva- factors (e.g., effective population size, demography) (38). tion scores, whereas their intronic flanking sites exhibited a reduced level of conservation. Intronic positions were mostly Pathogenic BP Variants. Variants in BP may lead to reduced, around zero, whereas exonic positions were mostly positive. altered, or abolished binding to spliceosome elements during the The BP region displayed a larger proportion of positive scores initiation of splicing, leading to the disruption of splicing and giv- (evolutionarily conserved) than negative scores (evolving neu- ing rise to misspliced gene products with the consequent loss of trally), with two peaks of high PhyloP scores at BP and BP2 expression/function (Fig. 1C). Through the HGMD database and positions. We cross-compared the scores between BP and primary literature, we identified 48 pathogenic BP variants (42 BP2 (Fig. 3F) and found that ∼23% of negatively scored BP single-nucleotide variants [SNVs] and 6 microdeletions) reported positions had positive scores for their corresponding BP2 in 43 human genes underlying 33 different disorders that had positions (second quadrant), while ∼22% of negatively scored been experimentally confirmed (Table 2 and SI Appendix, Table BP2 positions had positive scores for their corresponding BP S8), from the first report of a BP variant in the L1CAM gene PNAS 2022 Vol. 119 No. 44 e2211194119 https://doi.org/10.1073/pnas.2211194119 5 of 12 Table 2. The 48 reported pathogenic BP variants underlying human inherited disorders, with experimentally confirmed molecular consequences Pathogenic BP variants BPHunter detection Distance Gene Variant to 30 ss Disease Consequence Rank Hit position ABCC8 g.17452526T > C 20 Hyperinsulinemic hypoglycemia Partial retention of intron-11 (73nt) #1/2 0 ALPL g.21896765_21896784del 33 Hypophosphatasia Complete skipping of exon-8 and exon-7/8 #1/1 2j1j0 BBS1 g.66287067A > T 21 Retinitis pigmentosa Complete skipping of exon-8 and exon-7/8, #1/1 0 partial skipping of exon-8 (30 nt) BTK g.100609705T > C 23 Agammaglobulinemia Complete skipping of exon-16 #1/1 0 C21orf2 g.45750232T > A 23 Axial spondylometaphyseal dysplasia Complete retention of intron-6 #2/3 0 CAPN3 g.42684808del 29 Calpainopathy Partial retention of intron-6 (389 nt) #1/1 2 CD40LG g.135736500_135736507del 32 X-linked hyper-IgM syndrome Complete skipping of exon-3 #1/1 2j1j0 CDT1 g.88873665A > G 24 Meier-Gorlin syndrome Complete retention of intron-8, complete #1/1 0 skipping of exon-9 COL4A5 g.107849932A > G 40 Alport syndrome Complete skipping of exon-29, partial skipping #1/1 0 of exon-29 (43 nt) COL5A1 g.137686903T > G 25 Ehlers-Danlos syndrome type II Partial skipping of exon-33 (45 nt) #1/2 2 COL7A1 g.48616971T > C 23 Dystrophic epidermolysis bullosa Complete retention of intron-58 and #1/1 0 intron-58/59, complete skipping of exon-59 CPS1 g.211452758A > G 24 Hyperammonemia Complete skipping of exon-7 #1/2 0 DYSF g.71817308A > G 33 Limb-girdle muscular dystrophy Complete retention of intron-31 #2/2 0 ENG g.130578354A > G 22 Pulmonary arterial hypertension Complete skipping of exon-13 #1/2 2 F8 g.154130469T > C 27 Hemophilia A Complete skipping of exon-19 #1/1 0 F9 g.138619496A > G 25 Hemophilia B Partial retention of intron-2 (25 nt) #1/1 0 FAS g.90770494A > G 16 Autoimmune lymphoproliferative Complete skipping of exon-6 N.D. syndrome FBN2 g.127670562A > C 26 Congenital contractural arachnodactyly Complete skipping of exon-31 #1/1 2 FBN2 g.127671284T > C 15 Congenital contractural arachnodactyly Complete skipping of exon-29 N.D. FGD1 g.54476769del 35 Aarskog–Scott syndrome Complete skipping of exon-13 #1/1 0 HEXB g.74014605A > G 17 Sandhoff disease Partial retention of intron-10 (37nt) N.D. IKBKG g.153788599A > T 23 Ectodermal dysplasia with primary Complete skipping of exon-5 and exon- #1/1 0 immunodeficiencies 3/4/5/6, complete retention of intron-4 ITGB2 g.46321660A > C 12 Leukocyte adhesion deficiency Partial skipping of exon-6 (149 nt) N.D. ITGB4 g.73732344T > A 25 Epidermolysis bullosa with pyloric Complete retention of intron-14 and #1/1 2 atresia intron-14/15 ITGB4 g.73748508T > A 19 Epidermolysis bullosa with pyloric Complete retention of intron-31, partial #1/1 2 atresia skipping of exon-32 (38 nt) KCNH2 g.150646165T > C 28 Long QT syndrome Partial retention of intron-9 (147 nt) #1/1 0 L1CAM g.153131293T > G 19 X-linked hydrocephalus Complete skipping of exon-19, partial #2/3 0 retention of intron-18 (69 nt) LCAT g.67976512A > G 22 Fish-eye disease Complete retention of intron-4 #1/2 2 LIPC g.58830518A > G 14 Hypertriglyceridemia and Complete retention of intron-1, partial N.D. cardiovascular disease retention of intron-1 (13 nt, 78 nt) LMX1B g.129377625_129377641del 37 Nail patella syndrome Complete skipping of exon-2 #1/1 2j1j0 MLH1 g.37090369T > G 26 Inherited cancer Complete skipping of exon-17 #2/3 2 MLH1 g.37090371A > G 24 Inherited cancer Complete skipping of exon-17 #2/3 0 MSH2 g.47709894A > G 24 Inherited cancer Complete skipping of exon-16 and 30 #1/1 0 untranslated region, partial retention of intron-15 (85 nt, 141 nt) MSH6 g.48032731T > G 26 Inherited cancer Complete skipping of exon-5 #2/3 2 NPC1 g.21137182T > C 28 Niemann-Pick type C disease Complete skipping of exon-7 #1/1 0 NTRK1 g.156843392T > A 33 Congenital insensitivity to pain with Partial retention of intron-7 (137nt) #1/1 2 anhidrosis RB1 g.49039315A > T 26 Retinoblastoma Complete skipping of exon-24 #1/2 1 SLC25A20 g.48921567A > C 10 Carnitineacylcarnitine translocase Complete skipping of exon-3 and exon-3/4 N.D. deficiency SLC5A2 g.31499327_31499349del 31 Familial renal glycosuria Complete skipping of exon-8 #1/2 2j1j0 TH g.2187017A > T 24 Extrapyramidal movement disorder Complete skipping of exon-12, partial #3/3 2 retention of intron-11 (36 nt) TSC2 g.2138031A > G 18 Tuberous sclerosis Complete retention of intron-38, partial #1/2 0 skipping of exon-39 (74 nt) UROS g.127477605A > C 31 Congenital erythropoietic porphyria Partial retention of intron-9 #1/1 2 (81 nt, 246 nt, 358 nt, 523 nt) USH2A g.216040529T > C 17 Usher syndrome Partial skipping of exon-44 (39 nt) N.D. VMA21 g.150572076A > C 27 Autophagic vacuolar myopathy Showed significant reduction of expression #1/1 0 and activity VMA21 g.150572076A > T 27 Autophagic vacuolar myopathy Showed significant reduction of expression #1/1 0 and activity VWF g.6101204T > A 20 von Willebrand disease Complete skipping of exon-38 #1/2 0 XPC g.14209889A > T 9 Xeroderma pigmentosum Complete skipping of exon-4 N.D. XPC g.14209904T > C 24 Xeroderma pigmentosum Complete skipping of exon-4 #1/1 0 All references are available in SI Appendix, Table S8. Full information on these variants with BPHunter’s annotations are available in Dataset S4. (The genomic positions are on the hg19/ GRCh37 human genome assembly. Hit position: BP positions disrupted by the variant. N.D.: not detected.) causing X-linked hydrocephalus in 1992 (29) to a recent study of oncogenic, metabolic, neurological, immunological, and circula- BP variants in hereditary cancer genes in 2021 (39). These disor- tory. These BP variants comprised 16 homozygotes, 23 heterozy- ders were attributed to different classes of disease: developmental, gotes, 6 hemizygotes, and 3 without zygosity information and 6 of 12 https://doi.org/10.1073/pnas.2211194119 pnas.org were located within [9, 40] nt of 30 ss. Taking the reported Box 1: Input and Output of BPHunter missplicing consequences together (one variant may have multiple BPHunter supports both hg19/GRCh37 and hg38/GRCh38 human missplicing consequences), we identified 29 (60%) complete and reference genomes and provides an option to focus on all or canonical 7 (15%) partial exon skippings and 10 (21%) complete and 11 transcripts. Given a VCF file of variants with the first five columns as (23%) partial intron retentions, suggesting that complete exon mandatory fields (CHROM, POS, ID, REF, ALT), BPHunter outputs a skipping may be the predominant molecular consequence from text file of the variants that disrupt [2, 0] nt region of BP, with the fol- the disruption of BP recognition. The discovery narrative of these lowing informative annotations: pathogenic BP variants generally occurred in a similar context (SI • variant type (snv, x nt-deletion, x nt-insertion) Appendix, Box S1). Their discovery timeline showed that 11, 15, • gene symbol and 22 pathogenic BP variants were published in the 1990s, the • BP name (m/e/cBP_chrom_pos_strand_nucleotide) 2000s, and since 2010, respectively. Thus, neither the advent of • BP ranking in this 30 -proximal intron (#rank/total) MPS technology (since 2010) nor the publication of large-scale • BP position hit (2, 1, 0) human BP datasets (since 2015) has resulted in a major increase • distance from BP to 30 ss in the identification of pathogenic BP variants. In fact, paradoxi- • level of consensus (1:YTNAY, 2:YTNA, 3:TNA, 4:YNA, 0:none) cally, the proportion of pathogenic BP variants among all patho- • number of sources, and list of sources genic variants discovered since the 1990s has steadily decreased. • population variant MAF at BP and BP-2 positions (noted as “.” if The detection efficiency of pathogenic BP variants is rather low no variant found) when compared with other types of variant. Therefore, there is an • conservation scores GERP and PhyloP at BP and BP2 positions urgent need to bridge BP data with MPS data, thereby extending • intron type (major, minor), and intron length the utility of the MPS data, making more practical use of the BP • transcript ID and intron number (e.g., ENST123456789_ data and contributing to pathogenetic discovery. IVS10) • BPHunter score (SI Appendix, Fig. S20) BPHunter: Detection of Intronic Variants That Disrupt the BP. Variants that affect BP recognition may be individually rare, In addition, the BPHunter webserver provides a negative BP reporting platform, which requires the submitters to provide variant information, but their impact on splicing, gene expression, and protein func- cell type, assay type, and experimental evidence following the template tion can be highly deleterious. Therefore, in searching for (Dataset S6). disease-causing variants, we should routinely consider BP var- iants as a discrete category of potentially damaging variants, especially when the search within coding regions and essential respectively, 22 SNVs that disrupted BP sites, 11 SNVs that dis- splice sites has turned out to be unproductive. We therefore rupted BP2 sites, and 1 SNV that disrupted a BP1 site (Fig. developed BPHunter for the genome-wide detection of human 4B). In 24 cases, the variants disrupted the only BP in their 30 - intronic variants that disrupt the BP in MPS data, thereby proximal introns, in 9 cases, the variants disrupted the first BP, and enabling the identification of BP variant candidates that could in 7 cases, the variants disrupted the nonfirst BP (Fig. 4C). The potentially result in aberrant splicing at a biochemical level and affected BP comprised 19 mBP, 1 exclusively eBP and 20 exclu- hence disease at a physiological level. Given a VCF file of var- sively cBP (Fig. 4D). There were 24 BP that matched the YTNAY iants, BPHunter detects those variants that disrupt BP with consensus sequence, whereas the other 8, 7, and 1 BP matched informative outputs (Box 1). We developed the BPHunter as a the increasingly relaxed YTNA, TNA, and YNA consensus motifs one-line command that can be easily implemented into MPS (Fig. 4E). There were 39 BP positions supported by more than analysis and also as a user-friendly webserver for users with less one evidence (Fig. 4F). No population variants were found for computational expertise. We believe that there may have been 38 BP positions, whereas very rare variants were noted in the many negative findings involving BP variants that were tested other two BP positions (Fig. 4G). The conservation scores GERP but never reported. To improve our understanding of functional and PhyloP were both positive for 33 BP positions and both neg- BP, their mutational consequences, their utilization of alternative ative for 5 BP positions (Fig. 4H). We also used two missplicing BP, and their cell-type specificity, we ought to respect the poten- prediction tools [SpliceAI (40) and MMSplice (41)] and a muta- tial value to be derived from these negative data. Therefore, we created a reporting platform in the BPHunter webserver for tion deleteriousness prediction score [CADD (42)] to evaluate researchers to submit their tested BP variants that do not display these pathogenic BP variants (Fig. 4I). SpliceAI provides a proba- aberrant splicing. These data will be periodically collated and bility score for altering the acceptor sites (the higher the score, analyzed, and the report will become available on our website the higher the likelihood of missplicing), and it predicted no vari- and as a preprint, with all contributing researchers acknowledged ant above its high-precision cutoff of 0.8 (the highest was 0.74), as members of the BPHunter Effort Group. As we cannot vali- and five variants above its medium cutoff of 0.5. MMSplice pro- date the submitted negative BP variants centrally, we shall vides a score for acceptor disruption (the lower the score, the require all investigators to take responsibility for their submitted higher the likelihood of missplicing), and it predicted no variant data. Therefore, BPHunter constitutes a computational method beyond its suggested cutoff of 2 (the lowest was 1.944). capable of detecting candidate BP-disrupting variants and a plat- CADD provides a mutation damage score which is heavily reliant form for collating the negative results from BP genetics studies on conservation-related metrics. It identified 3 variants having that should not be wasted. scores >20 (genome-wide top 1%), 26 variants having scores in 10 to 20 (top 1 to 10%), and 11 variants with scores <10. With Validation of BPHunter on the Reported Pathogenic BP Variants. the use of BPHunter score (introduced in a later section) we We deployed BPHunter retrospectively on the 48 reported patho- found 38/40 variants had score ≥3 (Fig. 4J). genic BP variants and successfully captured 40 of them (Fig. 4A, Table 2, SI Appendix, Fig. S17, and Dataset S4). There were 4 dele- Pathogenic BP Variants That BPHunter Did Not Detect. tions (8/17/20/22 nt) that removed the entire BP motif, 2 single- BPHunter detected 40/48 reported pathogenic BP variants, nucleotide deletions that removed one BP site and one BP2 site, leaving 8 variants undetected (Table 2), all of which were PNAS 2022 Vol. 119 No. 44 e2211194119 https://doi.org/10.1073/pnas.2211194119 7 of 12 Comparison of the Pathogenic BP Variants with Common BP Variants in the General Population. Focusing on SNVs at BP and BP2 positions, we compared 33 pathogenic variants (22 variants on BP and 11 variants on BP2) against 659 com- mon population variants (386 variants on BP and 273 variants on BP2, with MAF ≥1%) (SI Appendix, Fig. S19A). We found that common variants occurred significantly less frequently in 30 -proximal introns with single BP (6%) versus pathogenic var- iants (58%, P < 0.00001 by Fisher’s exact test) and were signifi- cantly more frequent in BP of rank ≥3 (28% vs. 2%, P value < 0.00001). Thirty-five percent of common variants were located within BP motifs that did not match any level of consensus, which was 0% in pathogenic variants (P = 0); 63% of common variants were in BP sites supported by only one source, which was 3% in pathogenic variants (P < 0.00001). By contrasting the positively GERP/PhyloP-scored BP positions with the neg- atively GERP/PhyloP-scored BP positions, common variants (36% vs. 54%) showed significant enrichment (P < 0.00001) in having negative conservation scores against pathogenic var- iants (79% vs. 14%). There were 66% of common variants obtained low BPHunter scores = 0∼2 (the next section), whereas only one pathogenic BP variant was scored ≤2 (P < 0.00001). Furthermore, we analyzed the nucleotide changes in BP and BP2 positions separately. We compared 22 patho- genic variants against 386 common variants at BP positions (SI Appendix, Fig. S19B). A > G was the most frequent (68%) among the pathogenic BP variants (47) but not statistically Fig. 4. Detection of pathogenic BP variants by BPHunter. (A) Timeline of significant when compared with common BP variants (51%, the 48 published pathogenic BP variants (Left) and the distance to their 30 ss P = 0.128). We then compared 11 pathogenic variants against (Right). (B) Disrupted positions of BP. (C) Ranking of the disrupted BP in its 30 -proximal intron. (D) Type of disrupted BP. (E) Level of consensus of the 273 common variants at BP2 positions (SI Appendix, Fig. disrupted BP motif (1: YTNAY, 2: YTNA, 3: TNA, and 4: YNA). (F) Number of S19C). The nucleotide changes T > A and T > G were both data sources supporting the disrupted BP. (G) Population variation of the significantly enriched in pathogenic BP2 variants against com- disrupted BP. (H) Conservation scores of the disrupted BP. (I) Missplicing scores and deleteriousness score of the pathogenic BP variants. (J) mon BP2 variants (37% vs. 10% with P = 0.025 for T > A BPHunter scores of the pathogenic BP variants. and 45% vs. 10% with P = 0.0038 for T > G). The nucleotide change T > C was the least frequent (18%) in pathogenic located very close ([9, 17] nt) to 30 ss and 6 of which were BP2 variants although the most frequent (38%) in common published before 2004 (Fig. 4A). The eight introns harboring BP2 variants, albeit not significant (P = 0.2198). These statis- these variants all have BP documented in the BPHunter, but tics were based on the small number of reported pathogenic var- none of these pathogenic variants disrupted the BP. Interest- iants, and with additional BP variants reported in the future this ingly, by extracting and aligning their wild-type (wt) and analysis should be refreshed. mutated (mt) sequences using SeqTailor (43), we found that all of them had created an AG-dinucleotide in the region between Strategy for Prioritizing BP Variant Candidates. The analysis BP and 30 ss (SI Appendix, Fig. S18). By comparing the MaxEnt of the published pathogenic BP variants, in concert with the 30 ss strength (44) between the constitutional AGs and the cre- analysis of common BP variants, helped us to formulate a strat- ated AGs, we found that these eight variants created one stron- egy to prioritize BP variant candidates: 1) deletion of the entire ger 30 ss, six slightly weaker 30 ss, and one unlikely 30 ss. These BP motif, or disruption of the BP or BP2 positions; 2) variant eight original studies all reported complete/partial exon skip- of the only BP in the 30 -promixal intron or the first/second BP ping or intron retention, rather than de novo 30 ss at the created in the 30 -promixal intron with no more than three BP; 3) variant AGs. We therefore hypothesize that the AG-dinucleotides cre- in a BP motif matching the consensus YTNAY, while also allow- ated in the BP-30 ss region may have promoted binding with ing a slightly relaxed consensus; 4) variant at a BP position with U2AF35 that recognizes the acceptor site, thereby sterically more than one supporting source, regardless of its type (mBP/ interfering with the recognition of PPT by U2AF65 and the eBP/cBP); 5) variant at a BP position with no or rare population formation of the spliceosome complex (45, 46) (Fig. 1B). variants; 6) variant at a BP position with positive-scored GERP Thus, these eight variants should not be considered as BP var- or PhyloP may be preferred; and 7) variant C > T or T > C at iants. Rather, they should be regarded as causing aberrant BP2 could be deprioritized. We implemented a BPHunter splicing through the disruption of spliceosome organization scoring scheme to take account of these important features (Fig. on the 30 side of the intron. BPHunter therefore successfully 4J and SI Appendix, Fig. S20) to help to prioritize the candidate detected all 40 published pathogenic variants that truly dis- variants. The score ranges from 0 to 10, with a higher score indi- rupted BP. When these 40 variants were tested with existing cating a higher likelihood of deleteriousness. We recommend a cBP datasets, we found that 36, 30, and 37 can be detected threshold score ≥3 and suggest focusing on SNVs and deletions, using BPP, Branchpointer, and LaBranchoR, respectively. as the impact of insertions on splicing is difficult to predict using This confirms that the detection output from BPHunter is an the current tool. The SpliceAI, MMSplice, and CADD scores inclusive outcome that outperforms the use of any single could provide additional reference but should not be fully relied BP dataset. upon, particularly when evaluating newly identified BP variants 8 of 12 https://doi.org/10.1073/pnas.2211194119 pnas.org (see our case study of a STAT2 variant in a critical COVID-19 This BP matched the YTNAY consensus sequence, had a MAF patient). Additionally, homozygous and hemizygous variants of 6.98e-6, exhibited a GERP of 0.257 and a PhyloP of 0.152, should be prioritized, and candidate variants in genes associated and was supported by five data sources. These features satisfied with the disease under study should be prioritized. The sequenc- our prioritization strategy requirements for a promising candidate ing quality and the complexity/repeats of neighboring nucleoti- BP variant, and it had a BPHunter score = 6, even though the des need to be carefully checked. This prioritization strategy will missplicing/deleteriousness predictions were not suggestive of its be updated as the number of discovered BP variants increases. pathogenicity (SpliceAI = 0, MMSplice = 0.4, CADD = 10.1). We therefore hypothesized that this intronic variant might dis- A Germline Heterozygous SNV of STAT2 Disrupts BP and Splicing rupt BP recognition, leading to missplicing of STAT2 transcripts. in a Patient with Critical COVID-19. We performed BPHunter To assess its impact on mRNA splicing, we performed an exon on whole-exome sequencing (WES) data from a cohort of 1,035 trapping assay. Compared to the wild-type control, mRNA patients with life-threatening COVID-19 (48), focusing on extracted from cells transfected with this variant displayed a 13 genes that have been reported to carry deleterious variants majority of abnormally spliced STAT2 transcripts: 6% normal impairing type-I interferon (IFN) production and underlying transcripts in the mutant vs. 70% in the control, 82% of com- critical illness (49). BPHunter detected one heterozygous variant plete exon-6 skipping in the mutant vs. 25% in the control, and g.56749159T > A (IVS5-24A > T) in STAT2 that disrupted the 4% of complete intron-5 retention in the mutant vs. none in the only BP (mBP) located 24 nt from 30 ss of intron-5 (Fig. 5A). control (Fig. 5B). We extracted mRNA from whole blood of the patient and a healthy control and amplified the STAT2 comple- mentary DNA (cDNA) from exon-3 to exon-7/8 boundary. Sanger sequencing of the PCR products following TOPO-TA cloning showed that the patient’s cells had more misspliced tran- scripts (exon skipping and intron retention) than the healthy control (Fig. 5C). The cDNA from the patient contained 35% unannotated transcripts in GTEx, whereas the control contained 9% of unannotated transcripts (Fig. 5D). We also assessed STAT2 mRNA by RT-qPCR to estimate the amount of canoni- cal mRNA transcripts extracted from whole blood and observed about half the level in the patient as compared with the control (Fig. 5E). If we also consider the misspliced gene products that could have been degraded by NMD (50), the patient’s cells should carry more aberrantly spliced products. This case study effectively demonstrated the identification of a novel BP variant in STAT2 from a life-threatening COVID-19 patient and the experimental validation of its biochemical consequences. It also showed that a BP position with a very low level of evolu- tionary conservation can nevertheless be functional in splicing regulation, and it is for this reason that existing missplicing/ deleteriousness predictions are inadequate to identify promis- ing candidate BP variants. A Somatic Intronic 59-nt Deletion of ITPKB Removes All Three Potential BP Sites and Leads to the Retention of the Intron in a Lymphoma Patient. We performed BPHunter on somatic tumor-specific variant data derived from the whole-genome sequencing (WGS) of 53 diffuse large B-cell lymphoma (DLBCL) tumor samples (51–53). Focusing on a set of 212 lymphoma-associated genes (51), BPHunter detected one intronic deletion in ITPKB that removed 59 nucleotides from intron-3 (g.226835095_226835153del) (Fig. 6A). This intron contains three potential BP sites, and the deletion eliminated all of them, suggesting that the ITPKB gene product is highly likely to be misspliced in this tumor sample. ITPKB has been shown to be a negative regulator of the BCR/NF-κB signaling path- way which is very important for the growth and survival of lymphoma cells (54). It has also been shown that truncated ITPKB protein boosts PI3K/Akt signaling, which is a key growth-promoting pathway in lymphoma cells (55). We there- Fig. 5. Detection and validation of a germline heterozygous STAT2 variant fore hypothesized that this intronic deletion would lead to which disrupts BP and splicing in a life-threatening COVID-19 patient. (A) Detection of a heterozygous STAT2 variant disrupting BP. (B) STAT2 transcripts aberrant splicing, resulting in a truncated/decayed ITPKB gene and proportions from exon trapping assay in COS-7 cells. (C) STAT2 tran- product with loss of function/expression, leading to promotion scripts and proportions from RNA extracted from whole blood. (D) Ratio of of the growth and proliferation of lymphoma cells in the the GTEx-annotated transcripts among the total transcripts. (E) Estimation affected patient. We studied RNA-seq data of this tumor sam- of STAT2 canonical transcripts based on RT-qPCR, measuring STAT2 mRNA levels in whole blood using two probes spanning intron 5 (probe 5-6) and ple derived from the patient, as well as another three tumor intron 6 (probe 6-7). biopsies from random DLBCL patients who do not carry this PNAS 2022 Vol. 119 No. 44 e2211194119 https://doi.org/10.1073/pnas.2211194119 9 of 12 inhibitor was used to retain misspliced transcripts. This analysis demonstrated additional BP variants that lead to missplicing consequences at the biochemical level. Discussion Intronic sequences comprise ∼38% of the human genome, whereas exonic sequences comprise ∼3%. In the search for can- didate variants underlying human disease, current approaches mostly focus on nonsynonymous variants in coding regions and variants residing in essential splice sites. If an investigator wishes to screen for candidate variants among the vast number of intronic variants, a combination of approaches may be applied, e.g., testing for enrichment of intronic variants in a case-control study, computing and ranking the deleteriousness/missplicing scores (40–42) of intronic variants, referencing the MAF of popu- lation variations (32, 60), and clustering genetic heterogeneity (presumably deleterious intronic variants together with coding/ splice-site variants) underlying physiological homogeneity (61). However, enrichment of intronic variants relies on having a group of carriers in cases versus controls, so individual variants (albeit functionally impactful) in sporadic cases are rarely cap- tured. Moreover, many variant deleteriousness scores (42) are heavily reliant upon evolutionary conservation, and splicing has been shown to evolve rapidly with significant differences Fig. 6. Detection and validation of a private somatic intronic deletion of between species (36, 37). This suggests that the application of ITPKB in a lymphoma patient. (A) Detection of a somatic intronic 59-nt dele- tion of ITPKB that removed all three BP sites in intron-3. (B) RNA-seq reads conservation-derived deleteriousness scores to intronic splicing alignment on ITPKB that showed intronic retention in the lymphoma patient regulatory elements may not be particularly effective, a conclusion carrying this variant but not in the patients without this variant. (C) Sashimi supported by our STAT2 case study. We also showed that mis- plot from the RNA-seq data. splicing scores (40, 41) are underpowered to prioritize promis- deletion, and observed complete retention of intron-3 specifi- ing BP variants. Further, even though some intronic variant cally in this patient (Fig. 6 B and C). Hence, we concluded candidates may be identified, they do not come with the bio- that this intronic deletion, which would have been filtered out chemical interpretations that other types of variants do (e.g., in the absence of BPHunter, resulted in aberrant splicing of missense, stop-gain, essential splice site) to guide investigators ITPKB, which consequently contributed to the lymphomagene- in assessing their pathogenicity. sis in this patient. This case study successfully demonstrated With BPHunter, we are now able to detect intronic variants the use of BPHunter for the detection and exploration of that may disrupt BP recognition, and hence splicing and gene somatic variants. In addition, we screened the noncoding expression/function, in a systematic, efficient, and informative somatic variants in COSMIC (56), focused on 1,148 genes manner. It provides an opportunity for medical genetics studies to known to be associated with cancer formation and progression explore patients’ MPS data more deeply, by identifying candidate (Methods), and identified 129 candidate variants that may dis- BP variants in known disease-associated genes, or alternatively in rupt BP (106 SNVs and 23 deletions, in 101 genes, from 121 new and functionally promising candidate genes. It also provides patients involving 19 cancer types) (Dataset S5). an opportunity to systematically investigate the biomolecular effect of BP variants, at least in some genes, particularly in those 30 - Genome-Wide Detection of BP Variants and Validation of proximal introns containing multiple BP. Such studies will delin- Their Consequences from Paired WES and RNA-Seq Data. We eate the intron-specific molecular associations between BP and identified 14 individuals’ RNA-seq data that we previously pub- splicing. We expect an increasing number of characterized BP lished in different studies (30, 48, 57–59), whose paired WES variants, which alongside the characterization of their biochemi- data were also available in-house (SI Appendix, Table S9). We cal consequences and pathophysiological impact which will help used BPHunter for genome-wide detection, retained SNVs and us to comprehend BP function and dysfunction, from the posi- deletions that had BPHunter scores ≥3, removed variants that tive results standpoint. BPHunter also provides a platform to failed quality checking, and obtained 69 variants to be validated. report BP variants that do not give rise to missplicing, from the We checked their corresponding RNA-seq data (in different cell standpoint of negative results, which cannot be published but types, under nonstimulation) and found that 39 variants in genes should nevertheless not be neglected. Lessons learned from neg- that have no or very low expression, which we discarded for fur- ative BP variants will feed back valuable input that can be used ther analysis. Of the remaining 30 variants that had sufficient to determine functional BP in different contexts, which will in RNA-seq read coverage, we observed 11 (37%) variants that dem- turn improve our understanding of genuine BP positions and onstrated evidence of missplicing (retained intronic reads, exon alternative BP usages, while improving our ability to discover skipping, and reduced splicing efficiency) (SI Appendix, Fig. S21). pathogenic BP variants. We have also noted a study of 120 Cell type-specific missplicing was observed with variant #5, human BP variants (some were their own data, and some were where intronic retention was evident in neurons but not in fibro- collected from the published literature), in which the authors blasts. It is worth noting that the RNA-seq data used for this val- showed 82 BP variants (from 79 genomic positions) having no idation were generated by different configurations (protocols, missplicing effect (62). We found that 72/79 of these positions cell types, read lengths/depths, and machines), and no NMD were not documented as BP by BPHunter. Moreover, it should 10 of 12 https://doi.org/10.1073/pnas.2211194119 pnas.org be noted that many of the variants in this paper were misquoted datasets of computationally predicted BP (19–21). We created eBP_BP- from their cited origins and had no clear experimental validation. Hunter, by analyzing RNA-seq data of three DBR1-mutated patients (30). In this study, we built a comprehensive database of human BP We created cBP_BPHunter, by developing a high-precision machine-learn- which represents a qualitative and quantitative knowledgebase of ing model to scan the region [3, 40] nt of 30 ss. We screened a window BP in the human genome. This BP database is, however, not yet of [2, +2] nt from each raw BP position for consensus sequence matching optimal. The eBP are likely to include false-positive BP, caused and adjusted the position to its closest neighbor that perfectly matched by mislocation from RNA-seq reads. The cBP would have missed the consensus. false-negative BP (some real non-adenine-BP for example), as Human Population Variants. We obtained human genetic variants from machine-learning methods always learn the major distinguishing the gnomAD database (32) v3.1, which contained 76,156 WGS samples. We features from training data to differentiate between positives and extracted the total allele count (AC) and MAF of population variants, which negatives. Regarding the BP consensus sequence, we noted that were categorized into singleton (AC = 1), rare (MAF < 1%), and common while only ∼40% of BP motifs matched the long-established (MAF ≥ 1%). YTNAY, most of the remaining BP motifs fitted with relaxed Variant Deleteriousness and Missplicing Prediction Scores. We used consensus patterns. Our hypothesis I is that BP motifs that fit the CADD v1.6 (42), which predicts variant deleteriousness. We used two missplicing relaxed consensus (mostly displaying low binding affinity with scores, SpliceAI (40) and MMSplice (41), to evaluate the predicted probability of snRNA) may have other cis-elements (e.g., intronic splicing BP variants in disrupting splicing. enhancers) in the vicinity that interact with trans-elements with stronger affinity, to guarantee the collaborative organization of A Cohort of Patients with Critical COVID-19 and a Cohort of Lymphoma the spliceosome necessary for correct splicing (6). Our hypothesis Patients. We used the WES data from 1,035 patients with critical COVID- II is that the definition of the BP consensus sequence could be 19, which were recruited through The COVID Human Genetic Effort (48). context-specific (e.g., tissue-, cell-type-, developmental-stage-, or We also used WGS data from 53 diffuse large B-cell lymphoma patients, whose paired RNA-seq data from the tumor tissues were avail- intron-specific or RNA-structure-dependent) to some extent, to able (51–53). meet the highly context-specific physiological requirements of Written informed consent was obtained in the country of residence of each different tissues and cells. A variant that creates a new BP may patient, in accordance with local regulations and with institutional review board disrupt splicing as well, but BPHunter cannot detect this in its (IRB) approval. Data analysis was conducted under the approval of the IRB of the current form, which suggests one opportunity for improvement. Rockefeller University Institutional Review Board in New York. In addition, we are aware of some splicing mechanism studies The full methods section is available in SI Appendix. involving BP: recursive splicing within the intron allowing multi- step intron removal (9), stem-loop structures bringing distal BP Data, Materials, and Software Availability. BPHunter datasets and soft- closer to 30 ss (10), stochastic splice site selection leads to kinetic ware are publicly accessible from https://hgidsoft.rockefeller.edu/BPHunter (66) variation in intron removal (11), closer-to-50 ss BP leading to and https://github.com/casanova-lab/BPHunter (67), under CC BY-NC-ND 4.0 mutually exclusive splicing (63), and BP haplotypes associated License. All other study data are included in the article and/or supporting with alternative splicing (64). Taken together, these studies pre- information. sent a fascinating glimpse of the complexity of the spliceosome ACKNOWLEDGMENTS. We thank Y. Nemirovskaya, M. Woollett, L. Lorenzo, that calls for greater effort in understanding splicing and in inves- D. Liu, and E. Williams for administrative support; A. Cobat, E. Jouanguy, and tigating splicing-related pathogenesis. S. Boisson-Dupuis for assistance; and Z. Yang and Y. Nemirovskaya for artwork BPHunter enables the efficient genome-wide detection of design. The Laboratory of Human Genetics of Infectious Diseases is supported by candidate BP variants in MPS data. We hope that this tool will the Howard Hughes Medical Institute, the Rockefeller University, the St. Giles help to improve our ability to diagnose human genetic disease. Foundation, the NIH (R01AI088364), the National Center for Advancing Transla- We also anticipate the application of BPHunter to the study of tional Sciences (NCATS), NIH Clinical and Translational Science Award (CTSA) pro- human somatic variants, in cancers for example (65). Our timely gram (UL1 TR001866), the Yale Center for Mendelian Genomics and the GSP recruitment of BPHunter to screen our COVID-19 cohort Coordinating Center funded by the National Human Genome Research Institute directly led to the identification of a novel and deleterious BP (NHGRI) (UM1HG006504 and U24HG008956), the French National Research variant in the STAT2 gene from a patient with life-threatening Agency (ANR) under the “Investments for the Future” program (ANR-10- COVID-19, which was immediately tested experimentally to IAHU-01), ANR grants (ANR-14-CE14-0008-01, ANR-18-CE15-0020-02, and demonstrate its biochemical consequences. Our subsequent appli- ANR-20-CO11-0001), the Integrative Biology of Emerging Infectious Diseases cation of BPHunter to lymphoma patients also directly led to the Laboratory of Excellence (ANR-10-LABX-62-IBEID), the French Foundation for detection of an intronic deletion of ITPKB that removed all Medical Research (FRM) (EQU201903007798), a Fast Grant from Emergent three potential BP sites in that intron leading to its complete Ventures, Mercatus Center at George Mason University, the Fisher Center for Alzheimer’s Research Foundation, the Meyer Foundation, the JPB Foundation, retention as evidenced by RNA-seq data. Without BPHunter, the FRM and ANR GENCOVID project (ANR-20-COVI-0003), ANRS Nord-Sud neither of these deleterious variants that impaired the function (ANRS-COV05), the Square Foundation, Grandir-Fonds de solidarite pour of these key genes in the disease studied would have been picked l’enfance, the SCOR Corporate Foundation for Science, The Foundation du Souffle, up by preexisting conventional methods (e.g., missplicing predic- Institut National de la Sante et de la Recherche Medicale (INSERM) and the Paris tions, deleteriousness scores, enrichment tests). BPHunter therefore Cite University. We also thank the Swedish Cancer Society (Cancerfonden), the constitutes not only an important resource that provides us with Swedish Research Council, Radiumhemmets, the Center for Innovative Medicine an improved genome-wide understanding of BP but also a (CIMED), and the Knut and Alice Wallenberg Foundation. powerful tool that should drive us to discover new BP variants underlying human genetics diseases and traits. Methods Author affiliations: aSt. Giles Laboratory of Human Genetics of Infectious Diseases, Rockefeller Branch, The Rockefeller University, New York, NY 10065; bLaboratory of Collection of Eight Published BP Datasets, Creation of Two New BP Human Genetics of Infectious Diseases, Necker Branch, INSERM UMR1163, 75015 Paris, France; cParis Cite  University, Imagine Institute, 75015 Paris, France; Datasets, and Consensus-Guided Positional Adjustment of BP. We col- d Department of Biosciences and Nutrition, Karolinska Institutet, 14183 Huddinge, lected five datasets of experimentally identified BP (5, 15–18) and three e Sweden; Institute of Medical Genetics, School of Medicine, Cardiff University, Cardiff PNAS 2022 Vol. 119 No. 44 e2211194119 https://doi.org/10.1073/pnas.2211194119 11 of 12 CF14 4XN, United Kingdom; fInfection in Immunocompromised Pediatric Patients Barcelona Hospital Campus, Vall d’Hebron University Hospital, 08035 Barcelona, Spain; j Research Group, Vall d’Hebron Research Institute, Vall d’Hebron Barcelona Hospital Immunology Division, Genetics Department, Vall d’Hebron University Hospital, Vall Campus, Vall d’Hebron University Hospital, 08035 Barcelona, Spain; gPediatric d’Hebron Barcelona Hospital Campus, Autonomous University of Barcelona, 08035 Infectious Diseases and Immunodeficiencies Unit, Vall d’Hebron University Hospital, Barcelona, Spain; and kHHMI, New York, NY 10065 Vall d’Hebron Research Institute, Vall d’Hebron Barcelona Hospital Campus, Author contributions: P.Z., L.A., and J.L.-C. designed research; P.Z., Q.P., and W.R. Autonomous University of Barcelona, 08035 Barcelona, Spain; hJeffrey Modell performed research; P.Z., W.-T.L., J.L., P.D.S., P.S.P., R.C., B.B., S.-Y.Z., A.P., Q.P.-H., Q.Z., Diagnostic and Research Center for Primary Immunodeficiencies, 08035 Barcelona, and D.N.C. contributed new reagents/analytic tools; P.Z., Q.P., and W.R. analyzed data; Spain; iDiagnostic Immunology Group, Vall d’Hebron Research Institute, Vall d’Hebron and P.Z., D.N.C., L.A., and J.L.-C. wrote the paper. 1. M. M. Scotti, M. S. Swanson, RNA mis-splicing in disease. Nat. Rev. Genet. 17, 19–32 (2016). 37. J. Xiong et al., Predominant patterns of splicing evolution on human, chimpanzee and macaque 2. G. S. Wang, T. A. Cooper, Splicing in disease: Disruption of the splicing code and the decoding evolutionary lineages. Hum. Mol. Genet. 27, 1474–1485 (2018). machinery. Nat. Rev. Genet. 8, 749–761 (2007). 38. H. Ellegren, The different levels of genetic diversity in sex chromosomes and autosomes. 3. P. A. Sharp, Splicing of messenger RNA precursors. Science 235, 766–771 (1987). Trends Genet. 25, 278–284 (2009). 4. K. Gao, A. Masuda, T. Matsuura, K. Ohno, Human branch point consensus sequence is yUnAy. 39. D. M. Canson et al., The splicing effect of variants at branchpoint elements in cancer genes. Nucleic Acids Res. 36, 2257–2267 (2008). Genet. Med. 10.1016/j.gim.2021.09.020 (2021). 5. T. R. Mercer et al., Genome-wide discovery of human splicing branchpoints. Genome Res. 25, 40. K. Jaganathan et al., Predicting splicing from primary sequence with deep learning. Cell 176, 290–303 (2015). 535–548.e24 (2019). 6. J. Tholen, M. Razew, F. Weis, W. P. Galej, Structural basis of branch site recognition by the human 41. J. Cheng et al., MMSplice: Modular modeling improves the predictions of genetic variant effects spliceosome. Science 375, 50–57 (2022). on splicing. Genome Biol. 20, 48 (2019). 7. J. J. Turunen, E. H. Niemel€a, B. Verma, M. J. Frilander, The significant other: Splicing by the minor 42. P. Rentzsch, D. Witten, G. M. Cooper, J. Shendure, M. Kircher, CADD: Predicting the spliceosome. Wiley Interdiscip. Rev. RNA 4, 61–76 (2013). deleteriousness of variants throughout the human genome. Nucleic Acids Res. 47 (D1), 8. R. Bai et al., Structure of the activated human minor spliceosome. Science 371, eabg0879 (2021). D886–D894 (2019). 9. C. R. Sibley et al., Recursive splicing in long vertebrate genes. Nature 521, 371–375 (2015). 43. P. Zhang et al., SeqTailor: A user-friendly webserver for the extraction of DNA or protein sequences 10. C. L. Lin, A. J. Taggart, W. G. Fairbrother, RNA structure in splicing: An evolutionary perspective. from next-generation sequencing data. Nucleic Acids Res. 47 (W1), W623–W631 (2019). RNA Biol. 13, 766–771 (2016). 44. G. Yeo, C. B. Burge, Maximum entropy modeling of short sequence motifs with applications to 11. Y. Wan et al., Dynamic imaging of nascent RNA reveals general principles of transcription RNA splicing signals. J. Comput. Biol. 11, 377–394 (2004). dynamics and stochastic splice site selection. Cell 184, 2878–2895.e20 (2021). 45. K. Wimmer et al., AG-exclusion zone revisited: Lessons to learn from 91 intronic NF1 30 splice site 12. A. R. Awan, A. Manfredo, J. A. Pleiss, Lariat sequencing in a unicellular yeast identifies mutations outside the canonical AG-dinucleotides. Hum. Mutat. 41, 1145–1156 (2020). regulated alternative splicing of exons that are evolutionarily conserved with humans. Proc. Natl. 46. S. J. Bryen et al., Prevalence, parameters, and pathogenic mechanisms for splice-altering acceptor Acad. Sci. U.S.A. 110, 12762–12767 (2013). variants that disrupt the AG exclusion zone. HGG Adv 3, 100125 (2022). 13. W. Chen et al., Transcriptome-wide interrogation of the functional intronome by spliceosome 47. J. Kralovicova, H. Lei, I. Vorechovsky, Phenotypic consequences of branch point substitutions. profiling. Cell 173, 1031–1044.e13 (2018). Hum. Mutat. 27, 803–813 (2006). 14. A. J. Taggart, A. M. DeSimone, J. S. Shih, M. E. Filloux, W. G. Fairbrother, Large-scale mapping of 48. J. L. Casanova, H. C. Su, COVID Human Genetic Effort, A global effort to define the human genetics branchpoints in human pre-mRNA transcripts in vivo. Nat. Struct. Mol. Biol. 19, 719–721 (2012). of protective immunity to SARS-CoV-2 infection. Cell 181, 1194–1199 (2020). 15. A. J. Taggart et al., Large-scale analysis of branchpoint usage across species and cell lines. Genome 49. Q. Zhang et al.; COVID-STORM Clinicians; COVID Clinicians; Imagine COVID Group; French COVID Res. 27, 639–649 (2017). Cohort Study Group; CoV-Contact Cohort; Amsterdam UMC Covid-19 Biobank; COVID Human 16. J. M. B. Pineda, R. K. Bradley, Most human introns are recognized via multiple and tissue-specific Genetic Effort; NIAID-USUHS/TAGC COVID Immunity Group, Inborn errors of type I IFN immunity in branchpoints. Genes Dev. 32, 577–591 (2018). patients with life-threatening COVID-19. Science 370, eabd4570 (2020). 17. G. J. S. Talhouarne, J. G. Gall, Lariat intronic RNAs in the cytoplasm of vertebrate cells. 50. F. H€auser et al., A non-invasive diagnostic assay for rapid detection and characterization of aberrant Proc. Natl. Acad. Sci. U.S.A. 115, E7970–E7977 (2018). mRNA-splicing by nonsense mediated decay inhibition. Mol. Genet. Metab. 130, 27–35 (2020). 18. M. Briese et al., A systems view of spliceosomal assembly and branchpoints with iCLIP. Nat. Struct. 51. W. Ren et al., Genetic landscape of hepatitis B virus-associated diffuse large B-cell lymphoma. Mol. Biol. 26, 930–940 (2019). Blood 131, 2670–2681 (2018). 19. Q. Zhang et al., BPP: A sequence-based algorithm for branch point prediction. Bioinformatics 33, 52. X. Ye et al., Genome-wide mutational signatures revealed distinct developmental paths for human 3166–3172 (2017). B cell lymphomas. J. Exp. Med. 218, e20200573 (2021). 20. B. Signal, B. S. Gloss, M. E. Dinger, T. R. Mercer, Machine learning annotation of human 53. K. Georgiou et al., Genetic basis of PD-L1 overexpression in diffuse large B-cell lymphomas. branchpoints. Bioinformatics 34, 920–927 (2018). Blood 127, 3026–3034 (2016). 21. J. M. Paggi, G. Bejerano, A sequence-based, deep learning model accurately predicts RNA splicing 54. R. Schmitz et al., Genetics and pathogenesis of diffuse large B-cell lymphoma. N. Engl. J. Med. branchpoints. RNA 24, 1647–1658 (2018). 378, 1396–1407 (2018). 22. T. A. Cooper, L. Wan, G. Dreyfuss, RNA and disease. Cell 136, 777–793 (2009). 55. E. Tiacci et al., Pervasive mutations of JAK-STAT pathway genes in classical Hodgkin lymphoma. 23. R. K. Singh, T. A. Cooper, Pre-mRNA splicing in disease and therapeutics. Trends Mol. Med. 18, Blood 131, 2454–2465 (2018). 472–482 (2012). 56. J. G. Tate et al., COSMIC: The catalogue of somatic mutations in cancer. Nucleic Acids Res. 47 (D1), 24. S. Bao, D. F. Moakley, C. Zhang, The splicing code goes deep. Cell 176, 414–416 (2019). D941–D947 (2019). R V 25. P. D. Stenson et al., The Human Gene Mutation Database (HGMD ): Optimizing its use in a clinical 57. F. G. Lafaille et al., Human SNORA31 variations impair cortical neuron-intrinsic immunity to HSV-1 diagnostic or research setting. Hum. Genet. 139, 1197–1207 (2020). and underlie herpes simplex encephalitis. Nat. Med. 25, 1873–1884 (2019). 26. T. Sterne-Weiler, J. Howard, M. Mort, D. N. Cooper, J. R. Sanford, Loss of exon identity is a 58. P. Bastard et al., Herpes simplex encephalitis in a patient with a distinctive form of inherited common mechanism of human inherited disease. Genome Res. 21, 1563–1571 (2011). IFNAR1 deficiency. J. Clin. Invest. 131, e139980 (2021). 27. L. S. Kremer et al., Genetic diagnosis of Mendelian disorders via RNA sequencing. Nat. Commun. 59. M. Ogishi et al., Impaired IL-23-dependent induction of IFN-γ underlies mycobacterial disease in 8, 15824 (2017). patients with inherited TYK2 deficiency. J. Exp. Med. 219, e20220094 (2022). 28. D. N. Cooper, Functional intronic polymorphisms: Buried treasure awaiting discovery within our 60. P. Zhang et al., PopViz: A webserver for visualizing minor allele frequencies and damage genes. Hum. Genomics 4, 284–288 (2010). prediction scores of human genetic variations. Bioinformatics 34, 4307–4309 (2018). 29. A. Rosenthal, M. Jouet, S. Kenwrick, Aberrant splicing of neural cell adhesion molecule L1 mRNA 61. P. Zhang et al., A computational approach for detecting physiological homogeneity in the midst of in a family with X-linked hydrocephalus. Nat. Genet. 2, 107–112 (1992). genetic heterogeneity. Am. J. Hum. Genet. 108, 1012–1025 (2021). 30. S. Y. Zhang et al., Inborn errors of RNA lariat metabolism in humans with brainstem viral infection. 62. R. Leman et al., Assessment of branch point prediction tools to predict physiological branch points Cell 172, 952–965.e18 (2018). and their alteration by variants. BMC Genomics 21, 86 (2020). 31. A. B. Rosenberg, R. P. Patwardhan, J. Shendure, G. Seelig, Learning the sequence determinants of 63. C. W. Smith, B. Nadal-Ginard, Mutually exclusive splicing of alpha-tropomyosin exons enforced by alternative splicing from millions of random sequences. Cell 163, 698–711 (2015). an unusual lariat branch point location: Implications for constitutive splicing. Cell 56, 749–758 32. K. J. Karczewski et al.; Genome Aggregation Database Consortium, The mutational constraint (1989). spectrum quantified from variation in 141,456 humans. Nature 581, 434–443 (2020). 64. J. Kralovicova, S. Houngninou-Molango, A. Kr€amer, I. Vorechovsky, Branch site haplotypes that 33. G. M. Cooper et al.; NISC Comparative Sequencing Program, Distribution and intensity of control alternative splicing. Hum. Mol. Genet. 13, 3189–3202 (2004). constraint in mammalian genomic sequence. Genome Res. 15, 901–913 (2005). 65. C. Calabrese et al.; PCAWG Transcriptome Core Group; PCAWG Transcriptome Working 34. K. S. Pollard, M. J. Hubisz, K. R. Rosenbloom, A. Siepel, Detection of nonneutral substitution rates Group; PCAWG Consortium, Genomic basis for RNA alterations in cancer. Nature 578, 129–136 on mammalian phylogenies. Genome Res. 20, 110–121 (2010). (2020). 35. N. K. Kadri, X. M. Mapel, H. Pausch, The intronic branch point sequence is under strong 66. P. Zhang et al., BPHunter. https://hgidsoft.rockefeller.edu/BPHunter/. Accessed 17 October evolutionary constraint in the bovine and human genome. Commun. Biol. 4, 1206 (2021). 2022. 36. H. Keren, G. Lev-Maor, G. Ast, Alternative splicing and evolution: Diversification, exon definition 67. P. Zhang et al., BPHunter, Genome-wide detection of human variants that disrupt intronic and function. Nat. Rev. Genet. 11, 345–355 (2010). branchpoints. GitHub. https://github.com/casanova-lab/BPHunter. Deposited 8 January 2022. 12 of 12 https://doi.org/10.1073/pnas.2211194119 pnas.org