Nature 论文：微外显子调控密码揭示与自闭症相关的遗传变异

Main

Autism spectrum disorder (ASD) is a highly heterogeneous condition. Common characteristics of affected individuals include altered social and communication abilities, and repetitive behaviours. Hundreds of genetic alterations have been linked to ASD^5,6, and approximately 130 genes have been confidently associated with ASD susceptibility, with functions that converge on transcriptional and chromatin regulation, synaptic signalling and translational control^7,8,9. However, high-confidence ASD-associated variants are detected in only around 15% of affected individuals. Among the approximately 85% of idiopathic individuals, missing genetic aetiology may be attributable, in part, to variants overlapping non-coding genomic regions with important yet difficult-to-predict functions¹⁰.

Alternative splicing, a critical process by which exon and intron sequences are differentially included in transcripts to generate functional and regulatory diversity, is widespread in the mammalian nervous system^11,12. Neuronally enriched, 3–27-nucleotide-long microexons represent a highly conserved class of alternative splicing^1,2. Most annotated microexons preserve reading frames, and characterized examples function in modulating protein–protein interactions required for nervous system development and function^3,4. Consistently, mice haploinsufficient of the major activator of neuronal microexons, Srrm4, display autistic-like neuroanatomical, electrophysiological and behavioural features¹³. Moreover, perturbations of individual microexons also display ASD-related transcriptomic, neurodevelopmental and behavioural phenotypes^{3,4,14,15,16,17}.

Transcriptome-wide profiling of post-mortem brain samples from ASD-affected individuals has revealed recurring, distinct patterns of altered gene and splice isoform changes^18,19,20, with the most prominent splicing dysregulation affecting microexons¹. This dysregulation may arise in part through factors that affect neuronal activity, as the altered transcriptomic patterns correlate with those detected in activity-stimulated neural cells^13,19. However, whether and to what extent microexon splicing is directly affected by human genetic variation linked to neurodevelopmental disorders is not understood. Moreover, the extent to which microexons remain undetected is unclear, as previous studies focused primarily on analysing RNA-sequencing (RNA-seq) data from intact brain tissues, potentially missing the identification of microexons spliced in specialized neuronal subtypes or those that trigger transcript instability.

Here we address these questions through an integrative analysis encompassing the identification of an expanded repertoire of neuronal microexons, as well as the development of a machine learning model that delineates sequence elements critical for microexon splicing, which can be used to predict the impact of genetic variants on microexon regulation. Through these approaches, we reveal predominantly rare human genetic variation predicted to convergently affect a subset of microexons in genes with known and unknown genetic and functional links to nervous system biology and neurodevelopmental disorders.

An expanded landscape of neuronal microexons

To define a more comprehensive set of microexons, we applied a de novo detection method²¹ using a precomputed set of cDNA sequences representing the inclusion and skipping of all putative microexons, of which there are approximately 2 million across human introns (Fig. 1a). Putative microexons are defined here as intronic regions of 3–27 nucleotides flanked by sequences resembling splice sites, but which were not previously detected in spliced transcripts. From analysing 4,790 whole tissue RNA-seq samples from the human brain and other sources, and data from 47,380 single nuclei from human cortical regions²², we detect 1,514 microexons not previously annotated in GENCODE or VastDB (Methods), of which 292 display differential splicing between neural and non-neural cells and tissues. Approximately 43% of these are detected in both bulk and single-nucleus RNA-seq (snRNA-seq) datasets (Methods, Extended Data Fig. 1a–d and Supplementary Tables 1 and 2). Similarly, for mouse, by analysing RNA-seq data from 297 whole tissue samples, 133 targeted translating ribosome affinity purification sequencing (TRAP-seq) datasets and two single-cell RNA-seq (scRNA-seq) profiling datasets, the latter comprising 25,479 cells derived from primary visual and anterior lateral motor cortex^23,24, 1,115 new microexons were identified, of which 95 are detected as neurally spliced, often in multiple datasets (Extended Data Fig. 1e and Supplementary Tables 1 and 2).

Compared with putative microexons, the newly detected microexons are associated with higher 5′ and 3′ splice site strength scores and higher levels of surrounding sequence conservation (Fig. 1b–d; P < 3 × 10⁻³, two-sided Mann–Whitney U-tests), and are in genes expressed, on average, at lower steady-state levels compared with those containing previously annotated microexons, but more highly expressed than genes with putative microexons (Fig. 1e; P < 2 × 10⁻⁴, two-sided Mann–Whitney U-test with Bonferroni correction). Consistent with these observations, compared with annotated microexons, the new microexons more often introduce premature termination codons (PTCs) predicted to induce nonsense-mediated mRNA decay (NMD), and those that preserve reading frame less often overlap protein domains and intrinsically disordered regions (Fig. 1f and Extended Data Fig. 1f; P = 1.6 × 10⁻¹⁴, two-sided χ² test with Yates’ correction).

Hierarchical clustering of the new microexons across the analysed samples on the basis of percentage spliced-in (PSI) measurements reveals that they generally have higher inclusion levels in snRNA-seq compared to bulk tissue RNA-seq data (Fig. 1g, Methods and Extended Data Fig. 1b,c,g–i), including those with NMD-inducing potential, consistent with NMD requiring cytoplasmic translation (Fig. 1g). Supporting a regulatory role, genes containing predicted NMD-inducing known and new microexons, relative to frame-preserving microexons, have significantly reduced mRNA expression in brain RNA-seq data (Fig. 1h; P = 0.0068, two-sided Mann–Whitney U-test). Moreover, cultured mouse glutamatergic neurons treated with cycloheximide (which inhibits NMD by blocking translation) display increased splicing of candidate NMD-inducing microexons, on the basis of RNA-seq analysis, an effect that is strongly enriched among newly identified microexons, and assays using PCR with reverse transcription (RT–PCR) validated increased levels for 8 out of 11 tested examples (Fig. 1i,j, and Extended Data Fig. 2a). Similarly, treatment of cultured human neurons with the NMD inhibitor SMG1i resulted in increased detection of 4 out of 7 analysed new, predicted NMD-inducing microexons (Extended Data Fig. 2b). Consistent with the larger proportions of new human and mouse microexons (relative to known microexons) predicted to elicit NMD, and evidence that NMD-inducing alternative exons often are less conserved than regulated coding exons²⁵, a relatively small percentage (12.0%) of the new neural-differential human microexons have clear 1:1 orthologues in mouse. Collectively, these results reveal numerous examples of new microexons with neural-differential regulation and potential to control gene expression through activation of NMD, thus broadening the functional scope of this class of alternative splicing events.

Functional characteristics of new microexons

A subset of the new microexons analysed in mouse display marked neuronal subtype-specific splicing (Fig. 2a and Extended Data Fig. 3a). These include examples in genes with functions in cell migration and embryonic brain development (Actn4), vesicle fusion and synaptic transmission (Scfd2) and transcriptional control of neuronal development and survival (Klf7) (Fig. 2a and Extended Data Fig. 3a–d). A subset of new microexons displays dynamic splicing changes during embryonic and postnatal mouse brain development, including examples in genes such as Ralgapa1, Gtf2i and Psd2 that are crucial for neurodevelopment and implicated in brain disorders (Fig. 2b–d).

To assess whether developmentally regulated microexons predicted to trigger NMD affect gene expression, we systematically compared correlations between changes in microexon PSI and levels of the corresponding transcripts, across mouse embryonic stages and brain tissues (Methods). Consistent with NMD limiting detection of changes in transcript levels containing PTC-introducing microexons, increased detection of NMD-inducing microexons is, on average, modestly yet significantly correlated with decreased transcript levels, as compared to frame-preserving microexons with comparable developmental changes in PSI (mean ρ = −0.10 versus +0.11; P = 0.0046, two-sided Mann–Whitney U-test; Extended Data Fig. 4a). Overall, we identified 34 mouse microexons predicted to elicit NMD that are regulated during brain development, including 7 newly identified in this study (examples are shown in Extended Data Fig. 4b,c and Supplementary Table 3).

The expanded landscapes of microexons detected in human and mouse are significantly enriched in genes associated with Gene Ontology (GO) terms representing critical nervous-system-specific processes (Fig. 2e (data for mouse) and Extended Data Fig. 5 (data for human)). The newly identified microexons are additionally enriched in terms related to the Rho family of small GTPases, focal adhesion and lamellipodia (Fig. 2e, Extended Data Fig. 5 and Supplementary Table 4). Collectively, these results reveal more extensive repertoires of microexons with neuronal subtype and developmental splicing patterns, as well as NMD-activating potential, than previously appreciated. These repertoires serve as a basis for investigating the impact of genetic variation on microexon regulation and function.

Code governing neuronal microexon splicing

To define a sequence code that controls neuronal microexon splicing, we developed a machine learning approach that discriminates sequences surrounding true-positive microexons (detected by RNA-seq) from those surrounding a background, negative set of microexon-resembling sequences (not detected by RNA-seq) (Fig. 3a). For training, we compiled a positive dataset comprising around 1,800 neural-differential microexons from human, mouse and additional mammalian species, and a negative dataset consisting of 91,920 intronic sequences matching known microexons in terms of length, flanking splice site strength score distributions, and additional contextual features, such as their locations within genes and surrounding intronic GC content (Methods and Extended Data Fig. 6a). As the number of true-positive microexons was insufficient to train effective deep neural networks, we adopted a transfer-learning approach. Specifically, we trained gradient boosted decision trees that incorporate, as input features, the outputs of previously described deep-learning models of general splicing regulation and RNA-binding protein (RBP) binding^26,27,28. Additional features used by the code comprise characteristics of 5′ and 3′ splice sites, branch sites, nucleotide composition biases, and experimentally derived binding profiles and motifs of diverse RBPs (Extended Data Fig. 6b,c and Supplementary Tables 5 and 6). Cis-regulatory elements that control tissue-dependent alternative exons are typically concentrated within sequences surrounding the exons, but also overlapping the immediate upstream and downstream exons and their flanking intron sequences^29,30. Accordingly, the model learns features within these regions to discriminate the positive from the background set of microexon-associated sequences (Fig. 3a). The output of the resulting ‘microexon code-model’ is a score that estimates the likelihood that an input sequence represents a neuronally spliced microexon, on the basis of the composition of cis-RNA features overlapping these regions.

A train–test routine distinguished the positive set with a classification performance of 0.72 area under the precision-recall curve (AUPRC) (Extended Data Fig. 6d), which is a 45-fold increase in AUPRC over the baseline expected from random classification (AUPRC = 0.016). This performance was not significantly affected by excluding splicing-related features or orthologous sequences from test data that are present in training datasets (Methods and Extended Data Fig. 6e). Further supporting its performance, the model assigns significantly higher scores to known microexons analysed in Fig. 1 that have neuronal splicing patterns, compared to microexons that are constitutive, or spliced in subsets of non-neural samples (Extended Data Fig. 6f; P ≤ 0.001, Mann–Whitney U-test with Bonferroni correction). Moreover, although the model was not trained to predict PSI levels, its scores correlate more highly with PSI measurements for microexons in neural versus non-neural tissues (Extended Data Fig. 6g; Mann–Whitney U-test, P < 10⁻⁹).

To assess the accuracy of our model in pinpointing critical regulatory elements, we quantified the impact on prediction scores of in silico saturation mutagenesis (ISM) of each microexon and its neighbouring sequences described above (Fig. 3a,b and Extended Data Fig. 7a). Consistent with previous results indicating that microexons are on average associated with relatively strong 5′ splice sites and suboptimal 3′ splice sites^1,31,32, in silico mutation of 5′ splice site sequences has the largest impact on prediction scores, followed by mutations at 3′ splice sites (Fig. 3b). Moreover, the average predicted impact of mutations across the surrounding intronic regions strongly correlates with vertebrate genomic conservation of the corresponding sequences (Fig. 3b and Extended Data Fig. 7b; Pearson’s r = 0.92, P < 1 × 10⁻¹⁰⁰). For example, intronic sequences between 0 and −90 are predicted to have greater regulatory importance, consistent with this region encompassing higher evolutionary conservation, as well as concentrating cis-motifs critical for microexon splicing^1,32 (Fig. 3b,c; see below). Applying the model to the new neuronal microexons reveals that 59 of these have ISM profiles resembling those of known microexons (Extended Data Fig. 7c–f). This subset may therefore be spliced using similar rules as known neuronal microexons, whereas other new microexons may be regulated by rules not captured by our model.

To further assess performance of the microexon code-model, we compared PSI changes observed after experimental mutagenesis of splicing reporters with changes in code scores resulting from corresponding in silico mutations. PSI changes resulting from mutagenesis of motifs corresponding to putative intronic SRRM4- and PTBP1-binding sites proximal to microexons correlate with code score changes (Spearman’s ρ = 0.62, P = 1.36 × 10⁻⁶, 10 microexons and 49 mutations analysed; Extended Data Fig. 8a and Supplementary Table 7). Furthermore, using data from massively parallel human microexon reporter splicing assays performed in HEK293 cells engineered to inducibly express SRRM4 (to mimic a neuronal-like context)³³, we observed significant correlations between changes in PSI measurements and in silico code scores at three levels of SRRM4 expression (Spearman’s ρ = 0.049–0.115, all P < 4 × 10^–10, 128 microexons and 16,258 analysed mutations; Supplementary Table 8).

To assess how changes in code scores compare to changes in PSI in response to variants in vivo, we analysed variants within 300-nucleotide windows surrounding microexons detected in whole-genome sequencing (WGS) data, using patient-matched RNA-seq data from PsychENCODE brain samples^20,34,35. The top and bottom 5th percentiles of variant-dependent code score changes correlate significantly with changes in PSI measurements (Extended Data Fig. 8b; Pearson’s r = 0.25, P = 2.6 × 10⁻⁴). Moreover, our code-model compares favourably in predicting variant-dependent PSI changes from the reporter and in vivo brain datasets to Pangolin-Brain³⁶ (Pearson’s r = 0.24, P = 4.8 × 10⁻⁴) and DeltaSplice³⁷ (Pearson’s r = 0.23, P = 1.3 × 10⁻³), two previously described convolutional neural network (CNN) models trained on brain RNA-seq data and designed to predict PSI measurements (Methods, Extended Data Fig. 8b,c and Supplementary Table 9). Illustrating specificity, predictions from Pangolin models trained on heart, liver and testis RNA-seq data do not significantly correlate with variant-dependent ΔPSI measurements in patient brain RNA-seq data (Extended Data Fig. 8b). However, the Pangolin-Brain and DeltaSplice correlations are largely driven by a single (+1) 5′-splice-site-disrupting variant associated with the NRG1 microexon; removing this datapoint abolishes significance for all but the microexon code-model (Pearson’s r = 0.17, P = 0.034). This suggests that the microexon code-model is effective at capturing regulatory features relatively distal to core splice site dinucleotides, pinpointing cis-regulatory elements that are critical for microexon splicing and predicting effects of variants surrounding these exons.

Features of the microexon sequence code

To characterize the regulatory code learned by our model, we measured the impact of in silico mutagenesis of all hexamer sequences overlapping all known and the 59 new human microexons with similar ISM profiles. Mutations with a relatively large impact are enriched in subsets of related hexamer sequences resembling binding sites of known RBPs important for microexon regulation, including SRRM4, SRSF11, PTBP1 and RBFOX^1,2,32 (Fig. 3c). The code further correctly maps the location of these motifs. For example, UGC and C/U repeat-containing motifs upstream of the 3′ splice site represent experimentally defined bipartite microexon splicing enhancer sequences bound by SRRM4 and SRSF11, respectively³² (Fig. 3d).

Binding sites for undetermined or less-well-established regulators are also predicted by the code. For example, a prominent role for CUAAY motifs, primarily residing between −20 to −40, is inferred (Fig. 3c). While this motif corresponds to one of multiple possible branch site sequences³⁸, it also matches recognition sites of Quaking (QKI), a multifunctional RBP that controls alternative splicing, and that is linked to brain disorders³⁹. To investigate a role for QKI in microexon regulation, we depleted it in mouse neuroblastoma (N2A) cells using small interfering RNAs (siRNAs; Extended Data Fig. 9a,b) and performed RNA-seq. Microexons proximal to consensus QKI-recognition motifs^40,41,42 (that is, ACUAAY) within upstream introns show increased inclusion after Qki knockdown compared with those without QKI motifs^40,41,42 (P < 0.01, Wilcoxon test; Fig. 3e and Extended Data Fig. 9c). Moreover, these motifs peak within −20 to −40 nucleotides upstream of QKI-repressed microexons, consistent with the in silico predictions, including those from a targeted ISM analysis of QKI-binding motifs (Extended Data Fig. 9d,e). To link the ACUAAY motif to QKI regulation, we generated minigene splicing reporters containing the QKI-repressed mouse microexons Ppfia1 and Abi2, which contain upstream intronic QKI consensus binding motifs (Fig. 3f). Knockdown of Qki reveals that microexons from the reporters with the wild-type (WT) motif (ACUAAC) are upregulated after Qki knockdown, whereas mutation of the motif to ACUGAC (bold indicates the mutation), which preserves predicted overlapping branchpoints but disrupts the QKI-binding site, abrogates this effect (Fig. 3f). These results further support the efficacy of the code in inferring known and undefined regulatory elements surrounding microexons.

Effects of autism-linked variation on splicing

We next used the code-model to assess whether common and ASD-linked human genetic variation directly affects neuronal microexons. WGS data were initially analysed from 2,666 ASD-affected and unaffected sibling pairs from the Simons Simplex Collection (SSC) and Simons Foundation Powering Autism Research for Knowledge (SPARK) datasets⁴³. For the 369 code-selected microexons (310 known microexons and the subset of 59 new microexons described above), we calculated the differences in code score values (Δcode scores) between discordant sibling pairs as a consequence of individual or combinations of variants within the sequence intervals analysed above (Fig. 4a), focusing on those predicted to have the largest impact (Fig. 4b).

Variants with the largest sibling-pair Δcode scores were classified by their allele frequency (AF) across reference populations⁴⁴ as common (AF ≥ 1 × 10⁻²), rare (AF < 1 × 10⁻², ≥1 × 10⁻⁵) or ultrarare (AF < 1 × 10⁻⁵) and individually analysed. Within the analysed genomic regions, 43% contain a single dominant variant, and 57% are associated with two or more predicted impactful variants (Extended Data Fig. 10a–c). Notably, there is a marked inverse relationship between the frequency of these variants, and the magnitude of their predicted impact on sibling-pair Δcode scores (Extended Data Fig. 10d,e). Specifically, rare and ultrarare variants have a significantly greater predicted impact than common variants (Fig. 4c, P = 4.4 × 10⁻⁸, two-sided Mann–Whitney U-test). Consistently, DeltaSplice also yields a significant difference when comparing the predicted impact of (ultra)rare versus common variants on microexon PSI (Extended Data Fig. 10f). Variants overlapping highly conserved intronic regions exhibit more negative Δcode scores than those overlapping less-conserved intronic sequences (Fig. 4d; P = 1.6 × 10⁻⁸, two-sided Mann–Whitney U-test). When analysing combined annotation dependent depletion (CADD) scores—a genome-wide metric of variant deleteriousness that integrates various genomic features⁴⁵—we observed that 22.7% of the microexon-associated variants in ASD-affected individuals have predicted pathogenic potential (Methods), with significantly higher sibling-pair Δcode scores relative to 77.3% of variants scored as lacking pathogenic potential (P = 0.0013, two-sided Mann–Whitney U-test). These results indicate that the microexon code-model predicts impactful variants, with concordant support from the independent variant effect predictors, DeltaSplice and CADD.

To further assess reproducibility and test whether a subset of microexons is recurrently affected by rare variants, we analysed WGS data from an independent cohort of 5,112 ASD-affected and 263 unaffected siblings from the MSSNG resource⁸. Within each cohort (SSC, SPARK and MSSNG), we focused on microexons with the largest (that is, top decile) Δcode scores between affected and unaffected individuals, that is, those linked to (ultra)rare variants. This analysis revealed a significant overlap between microexons predicted to be impacted by variants across the three cohorts (Fig. 4e; 216 out of 369, 58.5% in the three-way overlap; Bonferroni-adjusted empirical P = 4 × 10⁻⁷). When restricting to Δcode scores inferred from ultrarare variants only, we still observe a significant three-cohort overlap (Extended Data Fig. 10g; Bonferroni-adjusted empirical P = 4 × 10⁻⁷). These results reveal that a common subset of microexons is enriched in predicted code-disrupting (ultra)rare variants derived from independent cohorts of ASD-affected individuals.

We next evaluated differences in the predicted cumulative impact of microexon-disrupting variants detected in ASD-affected versus unaffected siblings, by computing for each sibling a cumulative score that represents the summation of code score differentials (with respect to the hg38 reference genome assembly), for all predicted high-impact variants associated with all (n = 369) analysed microexons (Methods). Using data from the aforementioned massively parallel microexon splicing reporter assay data, we defined a lower-bound variant-effect threshold corresponding to a mean change in PSI ≥ 5, thereby identifying a set of variants referred to below as predicted impactful (Extended Data Fig. 8d). ASD-affected siblings were observed to have significantly higher cumulative differential code scores, relative to their unaffected siblings, above this threshold (Fig. 4f; P < 0.001, one-sided Wilcoxon signed-rank test). Moreover, cumulative sibling-pair Δcode scores predict ASD status (P = 0.0083, sibling-pair-matched conditional logistic regression controlling for sex; Supplementary Table 10). Consistently, when performing the same cumulative-impact analysis using CNN-based models, DeltaSplice yields comparable results (Fig. 4f).

To further assess the association between predicted impactful microexon variants and ASD status, we performed a rare-variant generalized disequilibrium test (RV-GDT)⁴⁶, which evaluates whether affected relatives carry a higher aggregate burden of specific variants compared to unaffected relatives, using phenotype permutations within pedigrees to obtain empirical P values. Using predictions from our code-model and DeltaSplice, predicted impactful microexon-disrupting variants were collected and analysed across 25,132 individuals, representing 7,391 families from SSC (n = 2,419), SPARK (n = 2,033) and MSSNG (n = 2,939), while stratifying variants by AF. Consistent with the cumulative-impact analysis, ultrarare variants prioritized by the code-model show a significant within-family association with ASD status (Fig. 4g,h; mean empirical P < 0.05 after Benjamini–Hochberg correction). Notably, while DeltaSplice predictions alone did not reach significance, an intersection set prioritized by both models did (Fig. 4g,h; mean empirical P < 0.05, after Benjamini–Hochberg correction). From these analyses, we defined a superset of 72 predicted impactful variants significantly associated with ASD, comprising 28 ultrarare variants prioritized by the code-model and 44 additional variants highlighted by both models (Fig. 4g,h). These variants are detected in 2.5% of ASD-affected individuals and absent from their unaffected siblings (Extended Data Fig. 10h). Finally, minigene reporter assays in human induced neurons show that representative rare and ultrarare variants included among the RV-GDT superset alter microexon splicing (Extended Data Fig. 10i and Supplementary Table 11). Notably, these examples include a variant that directly targets a microexon (but not its splice site consensus sequences) within the MAP4K4 gene. These results provide evidence that variants that are predicted to impact microexon splicing are significantly over-represented in ASD-affected individuals and associated with ASD status.

We next examined whether microexons associated with predicted impactful variants detected in ASD-affected individuals converge on specific processes and functions. Focusing on microexons predicted to be recurrently affected by rare variants (SSC/SPARK/MSSNG triple overlap) and associated with ASD status (RV-GDT superset), we observe convergence on genes that function in the formation and function of synapses (Fig. 5). These encode proteins that are localized to the presynaptic active zone, such as members of the PPFIA (liprin-α), ERC and RIMS families of proteins, and that form trans-synaptic connections, such as the LAR receptor tyrosine phosphatase family members PTPRS and PTPRD^47,48. Members of the Shank and Neurexin families of synaptic genes, as well as multiple genes involved in calcium signalling pathways (such as CADPS, CACNA1A, CACNA2D1, CACNB1 and RYR1), are also targets of the microexon-associated variants. Other predicted microexon-disrupting variants are in genes involved in signal transduction (such as SRC, PTK2, VAV2 and GIT2), vesicle trafficking (such as DCTN2, SYNJ1 and DNM1) and ubiquitination (such as NEDD4L, UBE2D3 and OS9). Overall, genes containing predicted impactful and ASD-associated variants, at the protein level, form an extended interaction network (Fig. 5), consistent with growing evidence that microexons often control protein–protein interactions that are critical for nervous system development and cognitive functioning^3,4,49,50,51.

Discussion

Neuronal microexons represent a critical hub of alternative splicing regulation, with characterized examples having defined roles in nervous system development and higher-order cognitive-functioning, as well as multiple links to brain disorders^3,4. Here we reveal an expanded landscape of this important class of alternative splicing events that are directly affected by human variation. This variation is predicted to contribute to genetic burden underlying ASD, through disruption of microexon-dependent functions that converge on synaptic biology and additional neuronal pathways.

Previous studies using machine learning have provided evidence for an impact of genomic variants on splicing patterns in neurological disorders^37,52,53,54. More-recent studies have further highlighted exons that are differentially spliced between individuals with ASD and unaffected individuals that are predicted to be affected by ASD-associated variants^55,56,57. An important feature of our model is that it enables the inference of the collective impact of multiple variants on splicing and ASD association. Through analysing ASD-affected and unaffected individuals from the SSC, SPARK and MSSNG cohorts, our study highlights an extensively overlapping subset of microexons associated with variants that have a relatively large inferred disruption of splicing, with rare and ultrarare variants contributing the most to this signal. Moreover, predominantly the latter variants are associated with ASD status and enriched in a network of genes with previously established, as well as less well-defined or unknown links, to neurological disorders.

We observe that more than a third of individuals, independent of ASD status, have one or more variants predicted to affect microexon splicing, whereas a smaller proportion (2.5%) of ASD-affected individuals have variants predicted to disproportionately contribute to genetic burden associated with this disorder (Extended Data Fig. 10h). Analogous to other variant-associated mechanisms in ASD⁵⁸, rarer variants are predicted to be more disruptive to microexon splicing than common variants, and are more often linked to ASD status. The majority (82%) of individuals carrying the ASD-associated variants (that is, from the RV-GDT superset) are idiopathic, whereas a smaller subset (18%) harbour previously detected⁸ candidate disruptive variants (Methods). Notably, similar to previously defined LGD variants⁵⁹, a subset of the ultrarare variants predicted to disrupt microexon regulation concentrates in genes that function in ubiquitin-dependent proteasome degradation pathways (Fig. 5).

Our study offers a resource for a previously underappreciated source of rare genetic variation with potential to contribute to convergent mechanisms that underlie ASD, as well as common genetic variation that may, more broadly, affect individual-specific neurobiological differences. To facilitate further exploration of regulatory mechanisms and genotype–phenotype relationships involving microexons, we provide a repository of code predictions for all single-nucleotide variants and indels evaluated in this work (Supplementary Table 12). In the short term, this resource may inform the redesign of targeted genotyping approaches, such that they can capture microexon-disrupting variants (Extended Data Fig. 10j), and also facilitate the interpretation of future WGS data collected from individuals affected not only by autism, but also by other brain disorders in which dysregulation of microexons may have a role²⁰. We expect this resource and insights from the present study to enable a deeper understanding of the multifaceted and critical functions of microexons in nervous system development, function and disorders.

Methods

Ethics statement

This study analysed previously collected, de-identified human data. No new participants were recruited and no new human samples were collected.

The study protocols, which were approved by the University of Toronto Research Ethics Board, encompassed the analysis of WGS data (RIS Human Protocol, 00039747) and PsychENCODE/NIMH Data Archive datasets (RIS human protocol 00037785).

Controlled access to these datasets was granted by the relevant data access committees under the terms of the respective data use agreements. Informed consent was obtained from all of the participants in the original studies that generated each dataset.

The use of human genetic and phenotypic data complied with all applicable ethical regulations.

Microexon discovery and quantitative analyses using RNA-seq data

We used a modified version of MicroExonator²¹ (https://github.com/hemberg-lab/MicroExonator/tree/blencowe) to detect known and new microexons across a large collection of human and mouse bulk RNA-seq, TRAP-seq, snRNA-seq and scRNA-seq (Extended Data Fig. 1d,e and Supplementary Table 1). Human snRNA-seq data generated by the BRAIN Initiative were retrieved from the NeMO Archive using gsutil v.4.59; the remaining human and mouse RNA-seq data were obtained from the Short Read Archive (SRA) using fastq-dump v.2.11.0. MicroExonator was used to perform a directed search of new microexons from a pre-computed list of putative microexon sites that were not previously annotated in GENCODE⁶⁰ or VastDB⁶¹. For detection and quantification purposes, nuclei and cells derived from snRNA-seq and scRNA-seq data were randomly allocated into three pseudobulks per cell cluster. Microexon inclusion was calculated as PSI across bulk RNA-seq and pseudobulk pools, for which only PSI measurements based on at least five spanning reads were considered for downstream quantitative analyses. Bulk RNA-seq experiments and pseudobulk pools were grouped into sample groups according to their corresponding biological condition or cell type of origin. We only considered putative microexons as detected if at least half of the bulk RNA-seq experiments or pseudobulk pools associated with the same sample group had a PSI equal to or greater than 10. Newly detected microexons were classified as neural if in at least one neural sample group exhibited an SSMD equal or greater than 1 and mean PSI difference equal or greater than 20 with respect to a set of non-neural bulk RNA-seq samples.

In our RNA-seq analyses, we aimed to estimate the proportion of both known and newly discovered microexons. To achieve this, we computed the count of microexons identified across incrementally larger fractions of randomly selected input data. This procedure was replicated 1,000 times to ensure robustness of our findings. Subsequently, we fit the microexon detection curves obtained from each replication to a nonlinear regression model, using the following logarithmic function:

\[{N}_{\mathrm{MIC}}(D)=\frac{\log (m\times D\times a+1)}{a}\]

where N_MIC(D) is the number of unique microexons, D is the percentage of data used, and m and a are free parameters that were fitted within the nonlinear regression model. The percentage of data (D) contributed by each sample was estimated by considering only the amount of reads from each sample that span across known and new neuronal microexon sites.

We assessed the potential of each known and newly detected microexon to trigger NMD by directly examining the implications of their inclusion across annotated isoforms. A microexon was classified as NMD-inducing if its inclusion resulted in a PTC situated more than 55 nucleotides upstream of the final splice junction in a minimum of 50% of the annotated isoforms. For comparisons of host-gene expression (Fig. 1e,h), each gene was represented by its median log₁₀[TPM] across 13 brain tissues from GTEx v8, and a gene was classified as NMD-inducing if any of its neuronal microexons was predicted to induce NMD. We further annotated the genomic context of each known and newly detected microexon by intersecting microexon coordinates with GENCODE untranslated region (UTR) and coding sequence (CDS) annotations to classify events as impacting UTRs or CDS. To assess predicted protein-level impacts, we intersected microexon coordinates (mapped to coding transcripts where applicable) with UniProt⁶² feature annotations to identify events overlapping annotated protein domains or intrinsically disordered regions. For clustering purposes, we imputed missing PSI values using a probabilistic principal component analysis function from pcaMethod R package (v.1.94.0)⁶³. We clustered PSI values applying Ward’s minimum variance criterion⁶⁴ using hclust R function (ward.D2 method).

To identify microexons that are differentially included during mouse development, we analysed embryonic brain RNA-seq data from ENCODE^65,66 and whole-cortex data across various embryonic and postnatal developmental stages⁶⁷. PSI measurements were grouped by neural tissue type (neural tube, hindbrain, midbrain, forebrain and whole cortex). For each microexon detected in at least five different developmental stages of a single tissue, we calculated the Spearman correlation of PSI across developmental time and associated P values. These P values were adjusted for multiple comparisons using the Bonferroni correction. A microexon was considered to be regulated through development if the adjusted P value was equal to or lower than 0.05 and the Spearman’s correlation was greater than 0, so that this definition selects microexons with increasing inclusion during development. Furthermore, we sought microexons with robust neuron-type-specific inclusion through downstream quantitative analyses of microexon PSI across pseudobulk pools derived from neuronal types extracted from the primary visual cortex and anterior lateral motor cortex²⁴. Pseudobulk pools were grouped by their associated neuronal cell type. A microexon was defined as having an inclusion specific to a given cell type if at least two pseudobulk pools associated with that cell type exhibited a PSI with a lower bound binomial 95% confidence interval 20 units greater than the median PSI across all pseudobulk pools.

We analysed RNA-seq data from mouse embryonic and postnatal brain tissues (forebrain, midbrain, hindbrain, neural tube and whole cortex) spanning multiple developmental stages^66,67 to identify microexons that are dynamically regulated during development. For each microexon and tissue, PSI measurements were modelled as a function of developmental time (days post coitum) using three approaches: (1) linear regression (ordinary least squares); (2) binomial generalized linear models (GLMs); and (3) piecewise binomial regression to capture non-linear relationships between PSI and time (implemented using the statsmodels module v0.14.2 in Python). A microexon was classified as developmentally regulated if at least one model yielded a significant time effect (Benjamini–Hochberg adjusted P < 0.05) in at least one regression in three or more tissues, a criterion that captures changes in either direction. Microexons classified in this way are those reported as developmentally regulated in Supplementary Table 3 and analysed in Extended Data Fig. 4. Gene expression TPM values were obtained with Whippet (v.1.6.1). For each developmentally regulated microexon, we calculated a single Spearman correlation between microexon PSI and host-gene expression (TPM), pooling measurements across developmental stages and brain tissues and requiring at least three measurements for which both values were available (Supplementary Table 3). Microexons predicted to induce NMD and frame-preserving microexons were then compared using a two-sided Mann–Whitney U-test (Extended Data Fig. 4a).

Functional enrichment analysis

To assess functional enrichment of genes harbouring microexons, we first obtained enriched terms using the G:GOSt module of g:Profiler⁶⁸ comparing all genes with known or new microexons with a custom gene space comprising all genes expressed in brain at greater than 0.1 transcripts per million, with the three GO domains and Reactome as annotation sources. Only terms comprising at least 5 and at most 1,000 genes, and with an intersection with the input of at least 3 genes, were considered. When there were multiple enriched categories with identical intersections, only the most enriched category was retained. The results were further processed with Cytoscape (v.3.10.0)⁶⁹, using EnrichmentMap (v.3.3.6)⁷⁰ together with the gene space GMT file from g:Profiler. A network was constructed using the default parameters. The node size was mapped to intersection size and edge width was mapped to overlap represented as the combined Jaccard-overlap coefficient (50%/50% weights) between pairs of nodes. Nodes with an FDR below 0.05 and edges with a combined coefficient above 0.375 were displayed. Nodes that did not reach this FDR cut-off in a parallel search using only genes with known microexons were highlighted. The AutoAnnotate (v.1.4)⁷¹ app in Cytoscape was then used to cluster nodes with default settings, and cluster names were manually curated from automatically generated three-word cluster labels.

Identification of newly detected neuronal microexons conserved in mouse

To assess conservation of newly detected neuronal microexons between human and mouse, we extracted genomic intervals spanning each human microexon either alone or including ±100 nucleotides of flanking intronic sequence. Human intervals were converted to mouse coordinates using UCSC liftOver, and the corresponding mouse sequences were retrieved from the reference genome. Human and mouse sequences were then aligned pairwise using the pairwise2 module in Biopython v.1.83 with the following scoring parameters: match = 1, mismatch = −0.25, gap open = −2 and gap extend = −0.25. For each interval, we computed a length-normalized alignment score (S_norm) as:

\[{S}_{\mathrm{norm}}=\frac{{S}_{\mathrm{align}}}{L}\]

where S_align is the raw pairwise alignment score and L is the length (in nucleotides) of the evaluated interval. Human microexons were classified as conserved in mouse if S_norm ≥ 0.50 across the evaluated intervals.

Cell culture

Mouse neuroblastoma (N2A) cells (ATCC CCL-131; sex not reported by the repository) were grown in DMEM (high glucose with l-glutamine, Wisent) supplemented with 10% FBS, sodium pyruvate, non-essential amino acids and penicillin–streptomycin. Human induced pluripotent stem (iPS) cell line 37 (WT), derived from a female patient, was previously described and characterized⁷². An independent human iPS cell line (WTC11; Coriell Institute for Medical Research, GM25256; derived from a male donor) was used for microexon validation experiments shown in Extended Data Fig. 2b. Neural precursor cells (NPCs) derived from iPS cell lines 37 (WT) were grown as a monolayer and split every 5–7 days in NPC medium containing DMEM/F12, N2, B27, non-essential amino acids, penicillin–streptomycin, 2 μg ml⁻¹ heparin, 1 μg ml⁻¹ laminin and 20 ng ml⁻¹ FGF2. Cells were maintained under subconfluent conditions. All cell lines were maintained at 37 °C with 5% CO₂. Cells were tested for absence of mycoplasma contamination by PCR using the e-Myco VALiD commercial kit (iNtRON Biotechnology). N2A and WTC11 cells were characterized by their source repositories, and iPS cell line 37 was previously characterized as described above.

In vitro neuron differentiation and NMD inhibition treatments

In vitro differentiation of CGR8 mouse embryonic stem cells (ECACC 07032901; derived from a male embryo) into glutamatergic neurons was performed using a protocol as described previously³². At 7 days in vitro, mature neurons were treated with 100 μg ml⁻¹ cycloheximide (Sigma-Aldrich) for 12 h to block translation and indirectly inhibit NMD, after which total RNA was collected. NPCs were plated on poly-l-ornithine + laminin-coated plates at a density of 750,000 cells per well in a six-well plate and cultured for 4 weeks in neural differentiation medium (Neurobasal, N2, B27, 1 μg ml⁻¹ laminin, penicillin–streptomycin, 10 ng ml⁻¹ BDNF, 10 ng ml⁻¹ GDNF, 200 ng ml⁻¹ ascorbic acid and 1 μM cAMP). After 4 weeks of neural differentiation, cells were treated with 1 μM SMG1i (Sigma-Aldrich), or with DMSO as a control for 12 h and RNA was collected using the QIAGEN RNeasy kit.

Transfection of siRNA and minigene reporters

N2A cells were transfected with 40 nM of siGENOME siRNA pools (Dharmacon) targeting mouse Qki or a non-targeting (siNT) siRNA pool (D-001206-13), using Lipofectamine RNAiMax (Thermo Fisher Scientific) as recommended by the manufacturer. Cells were seeded at 100,000 cells per well in 12-well plates, reverse transfected and total RNA was collected at 72 h after transfection. Plasmid constructs were transfected into N2A cells using Lipofectamine 3000 (Thermo Fisher Scientific) as recommended by the manufacturer. The day before transfection, cells were seeded in 12-well plates at a density of 100,000 cells per well (N2A). The next day, 250 or 500 ng of plasmid was transfected per well, and total RNA was collected at 48 h after transfection. For co-transfection of siRNA and plasmids, 100,000 N2A cells were seeded and reverse-transfected with siRNA in 12-well plates. Then, 24 h after seeding, 250 ng of minigene plasmid was transfected using Lipofectamine 3000 (Thermo Fisher Scientific). Total RNA was collected 48 h after plasmid transfection.

For transfection of minigene reporters in human cultured neurons, NPCs were plated at a density of 375,000 cells per well in a 12-well plate and cultured for 4 weeks in neural differentiation medium (see details above). Differentiated neurons were transfected with 1 μg of minigene plasmid using Lipofectamine 2000 (Thermo Fisher Scientific) and, at 48 h after transfection, total RNA was collected.

RNA extraction and RT–PCR assays

Total RNA was extracted from cells using the RNeasy Mini Kit (QIAGEN) as recommended by the manufacturer, with on-column DNase treatment. RNA was eluted in nuclease-free water and stored at −80 °C. RT–PCR assays were performed using the OneStep RT-PCR kit (QIAGEN) according to the manufacturer’s recommendations, in a 10 μl final volume with 40 ng of total RNA and 6 pmol primers per reaction. The reaction products were separated on 3–4% agarose gels and imaged using Bio-Rad Image Lab v6.0.1. PSI values were calculated from band intensities quantified using ImageJ v.1.54p or Fiji v.2.9.0 as the relative intensity of the exon-included (upper) band divided by the sum of the exon-included and exon-skipped (lower) bands. Uncropped gel images are provided as Supplementary Figures.

To validate new microexon sequences, RT–PCR bands of interest were either directly excised from the gel and purified using the GeneJET Gel Extraction Kit (Thermo Fisher Scientific), or the entire reaction was cloned as a pool using the Topo TA subcloning kit (Invitrogen). Purified products were directly Sanger sequenced using the RT–PCR primers, or multiple colonies from the pooled cloning reaction were sequenced from the M13Rev primer on the Topo vector backbone.

RT–qPCR

Two-step quantitative PCR with reverse transcription (RT–qPCR) was performed by generating cDNA using the Maxima H Minus First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) according to manufacturer’s recommendations. Then, 250 ng of total RNA was added per 20 μl reaction with 100 pmol random hexamer and 60 pmol Oligo(dT). The cDNA product was diluted tenfold. qPCR was performed using the SensiFAST SYBR No-ROX Kit (BIOLINE) with the CFX96 Real-Time PCR Detection System (Bio-Rad) according to the manufacturer’s recommendations. In total, 4 μl of the diluted cDNA was used per qPCR reaction in a final volume of 10 μl, with 500 nM of gene-specific forward and reverse primers. The ΔΔC_t method was used to calculate the fold change, with mouse Gapdh used as the housekeeping gene.

Immunoblotting

Cells were trypsinized, washed with PBS, lysed in Tris lysis buffer (10 mM Tris-HCl pH 8.0, 150 mM NaCl, 1% NP-40, and 10% glycerol) and sonicated. The lysates were further centrifuged at 4 °C, 10,000g for 10 min. The supernatant was mixed with Laemmli buffer and boiled for 5 min at 95 °C. The samples were separated by SDS–PAGE on 10% Tris-glycine gels, transferred to a PVDF membrane and blocked with 5% milk/PBST. The membranes were incubated with primary antibodies overnight: anti-pan-QKI (N147/6, Sigma-Aldrich, MABN624; 1:500) and anti-alpha-tubulin (B-5-1-2, Sigma-Aldrich, T6074, 1:10,000). The membranes were then incubated with HRP-conjugated secondary antibodies and imaged using ECL chemiluminescence (Perkin-Elmer). The specificity of the anti-pan-QKI antibody was supported by loss of the expected signal after Qki knockdown (Extended Data Fig. 9b); the alpha-tubulin antibody was used as a loading control on the basis of manufacturer validation for immunoblotting. Uncropped blot images are provided as Supplementary Figures.

Generation of minigene splicing reporters

Minigene reporters were based on the ExonTrap pET01 mammalian expression vector (MoBiTec). Each alternative exon of interest and around 300 nucleotides of its flanking intronic sequences were amplified from human or mouse gDNA, then restriction cloned in between two constitutive exons on the vector backbone. The Q5 site-directed mutagenesis kit (NEB) was used to generate the indicated mutations. After transfection and total RNA extraction, microexon inclusion levels for each minigene reporter was assayed using RT–PCR with primers annealing to the constitutive exons of the vector backbone.

RNA-seq

RNA was extracted from N2A cells after Qki knockdown, or from mouse embryonic-stem-cell-derived glutamatergic neurons after cycloheximide (CHX) treatment, using the RNeasy Mini Kit (QIAGEN) as described above. Extracted total RNA was submitted to the Donnelly Sequencing Centre for further processing. In total, 500 ng of total RNA was used for generation of polyA+ libraries using the NEB Next Ultra II Directional kit. Sequencing was performed by the Donnelly Sequencing Centre using the Nova-Seq 6000 instrument and S1 flow cell to an average depth of 50 million paired-end 150-nucleotide reads.

RNA-seq analysis after Qki knockdown and CHX treatment

Alternative splicing analysis of RNA-seq data from Qki knockdown experiments was performed using vast-tools (v.2.5.1)⁶¹. From the primary output, the set of mouse neural-high microexons were filtered for coverage (quality score 1 of VLOW or higher) and junction balance (quality score 4 of OK/B1/B2/Bn). Differential splicing was scored using vast-tool’s diff module requiring MV.dPsi._at_0.95 > 0 and a point estimate of |dPSI| > 10. Matt (v.1.3.0)⁷³ was used to identify occurrences of ACUAA[C|U] motifs within the 150 nucleotides surrounding intronic sequences of mouse neural-high microexons. Quantification of microexons in RNA-seq data obtained for CHX experiments was performed using MicroExonator, providing mouse genome assembly GRCm38 (mm10), GENCODE transcript annotation vM16 and splice site position frequency matrices from SpliceRack⁷⁴. To enable microexon detection, we provided the list of putative microexons previously detected through our mouse microexon discovery analyses described above. As part of the MicroExonator computational workflow, we assessed differential microexon inclusion using Whippet⁷⁵, defining microexons with |DeltaPsi| >0.1 and probability > 0.9 as differentially spliced.

Datasets for microexon code-model development

We constructed the foreground dataset from microexons annotated in the VastDB database⁶¹ downloaded on 27 August 2019 for human (hg38) and mouse (mm10), as well as pre-release data for cow (bosTau6), opossum (monDom5), chimpanzee (panTro4), rhesus macaque (rheMac2) and rat (rn6). To define neural-included microexons, we grouped RNA-seq samples into neural (adult brain and differentiated neural tissues and cell types) and non-neural, and excluded incompletely differentiated neural samples, such as neural precursors and embryonic tissues, resulting in 4–34 neural samples and 26–123 non-neural samples per species. For every microexon with a length between 3 and 27 nucleotides, we required that the mean PSI in neural samples be at least 10% greater than in non-neural samples and, if more than one sample per group was present (which was the case in 74–96% of microexons), a strictly standardized mean difference (SSMD) of at least 1:

\[\mathrm{SSMD}=\frac{({\mu }_{{\rm{n}}}-{\mu }_{{\rm{c}}})}{\sqrt{{\mathrm{var}}_{{\rm{n}}}+{\mathrm{var}}_{{\rm{c}}}}}\]

where μ_n and μ_c are the PSI means, and var_n and var_c the variances of neural and non-neural samples, respectively. This procedure yielded 344 human microexons that were classified with high confidence as neurally spliced, as well as 329 from mouse, 188 from cow, 178 from opossum, 253 from chimpanzee, 261 from rhesus macaque and 236 from rat, for a total of 1,789 neurally spliced microexons.

Next, to train a classifier that detects neurally spliced microexons, we constructed a background set of intronic sequences that closely resemble annotated neuronal microexons in both length (3–27 nucleotides) and splice site strength.

Splice site strength can be used to train an extremely highly accurate discriminator of microexons versus length-matched, randomly selected intronic sequences.

Therefore, to incentivize the model to learn informative microexon splicing signals beyond splice site strength alone, such as RBP-binding sites and branch-point strength, we developed a sampling procedure to select a set of intronic sequences that match the bivariate 3′/5′ splice site strength distribution of known microexons.

To build this background set, we first downloaded genome annotations for all seven species in our data and extracted all intronic regions. To avoid bias from different levels of gene expression we only retained intronic sequences from genes that contain at least one neurally spliced microexon. We then identified a set of candidate microexon-like sequences in these introns, defined as all sequences between 3 and 27 nt that are flanked by 3′ and 5′ splice sites of which the strengths, as measured by MaxEntScan⁷⁶, fall within the ranges observed for known microexons (−10 to +14 and −5 to +12, respectively), and that are not in the known microexon set. We then sampled from these candidate sequences to create the background set. We selected an initial set of 2 million candidate sequences for human, 1.6 million for mouse and slightly fewer for each of the other species.

To ensure that the background set matched the known microexons not only in the range of splice site strengths represented but also in their bivariate (that is, two dimensional) distributions of 3′ versus 5′ splice site strength, we used the known microexons to construct a two-dimensional target distribution for 3′/5′ splice site strength and then sampled the sequences for the background set from the candidate set according to this target distribution. To create the target distribution, we first separately estimated the one-dimensional marginal 3′ and 5′ splice site strength distributions of the known microexons using a kernel density estimator using the KernelDensity function from scikit-learn v.1.5.0 with bandwidth = 1. From these one-dimensional marginal distributions, we then created a discretized smoothed bivariate splice site strength distribution on a 20 by 20 grid that we used as the target distribution. These grid points were equally spaced within the splice site strength ranges indicated above (Extended Data Fig. 6a). We then binned the 3′ and 5′ splice site strengths of the candidate sequences on their respective grids to obtain normalized one-dimensional histograms of each marginal distribution. The outer product of these two histograms defined the normalized two-dimensional histogram used as the target distribution.

Given the two-dimensional target distribution represented by the two-dimensional normalization histogram, we next constructed a background set with a similar distribution of splice site strengths as the target distribution. Owing to the large number of candidate sequences, for efficiency reasons, we used a two-stage procedure to select the background set. First, we sampled a fixed number of candidate set sequences (n = 500) in each of the 400 (20 × 20) grid cells (that is, bins) in the two-dimensional histograms. If there were fewer than n candidate sequences in a bin, we sampled with replacement. Otherwise, we sampled without replacement. We selected the background sequences from this initial set of 200,000 candidates. Next, for each bin, we sampled, without replacement, m(x,y) sequences from the (x,y) bin:

\[m \sim \mathrm{Binomial}(n,a(x,y))\]

where, given the target distribution’s probability p(x,y), the acceptance probability for a candidate sequence in cell (x,y) is

\[a(x,y)=\frac{p(x,y)}{{\max }_{x,y}p(x,y)}.\]

These sampled sequences constitute the background set. Note that, even though the sequences were sampled without replacement, it is possible for the background set to contain multiple versions of the same sequence in bins where there were fewer than n candidate sequences.

We used the above procedure to separately generate background sequence sets for all seven mammalian species included in the training set, but always using the same target distribution originally made using the known microexons from all species combined.

Finally, we set apart a test set consisting of all known human neural microexons and human background sequences from chromosomes 1, 3, 5, 7 and 9. To prevent any sequence overlap between the training and test sets, for example, due to orthologous genes, we also used BLAT to remove sequences from the training set that have 95% or more sequence identity with a microexon and ±40 nt of its surrounding intronic sequence in the test set.

Features

For each foreground and background sequence, we computed a collection of features on the basis of a curated selection of pretrained machine learning models and tools covering splice site strength and features, RBP binding, branch point strength and nucleotide composition biases (Supplementary Tables 5 and 6).

Seqweaver

Seqweaver²⁸ is a neural network that is trained to predict 231 RNA CLIP-Seq RNA-binding profiles, as well as alternative exons and splice sites from sequence alone. We applied Seqweaver by obtaining prediction scores from the model for 217 of these RNA profiles on seven sequence intervals for each microexon:

1. (1)

   100 nucleotides upstream of the constitutive donor site upstream of the microexon.

2. (2)

   300 nucleotides downstream of the constitutive donor site upstream of the microexon.

3. (3)

   300 nucleotides upstream of the microexon’s acceptor site.

4. (4)

   The microexon itself.

5. (5)

   300 nucleotides downstream of the microexon’s donor site.

6. (6)

   300 nucleotides upstream of the constitutive acceptor site downstream of the microexon.

7. (7)

   100 nucleotides downstream of the constitutive acceptor site downstream of the microexon.

For each predicted feature we computed the maximum within each of the above intervals and concatenated them to obtain a total of 1,519 (=217 × 7) features.

COSSMO

COSSMO²⁷ is a neural network trained to predict alternative 3′/5′ splice site selection. There are two versions of the model—one to predict alternative acceptor selection and one to predict alternative donor selection. We used the former model to predict the choice between the microexon’s acceptor and the downstream constitutive exon’s acceptor. We used the latter model to predict the choice between the upstream exon’s donor site or the microexon’s donor site.

After applying each model, we obtained the predicted probabilities to choose either the microexon or the flanking exon’s splice site, the logits (unnormalized splice site strengths) of these probabilities along with values of the model’s final hidden layer (200 values for each splice site). These final-layer values were used to capture a higher-order representation of the splice site learned by the model. This yielded a total of 808 ((2 × 200 + 2 + 2) × 2) features. COSSMO feature generation was performed using TensorFlow (v.2.2.0).

MaxEntScan

To obtain splice site strength features, we used the MaxEntScan algorithm⁷⁶. Specifically, we used MaxEntScan to predict the strengths of the acceptor sites of both the microexon and the downstream flanking exon, as well as the donor sites of the microexon and the upstream flanking exon. Moreover, we used the difference in the strengths between the microexon and downstream acceptor scores and between the microexons and upstream donor scores.

Branch point scores

We used svm-bp⁷⁷ to predict the location and strength of both the microexon’s branch point and the downstream flanking exon’s branch point. We extracted 100 nucleotides of sequence upstream of both the microexon’s and the downstream flanking exon’s acceptor site and used svm-bp to identify the branch point candidate with the maximum strength in each of these sequences, and used their strengths, as well as their distance from their respective acceptor sites, as features (four total features).

CisBP PWMs

We used the CisBP-RNA compendium of binding motifs for 117 human RBPs. We extracted the sequence of each microexon, as well as 300 nucleotides of flanking intronic sequence on both sides, and used the FIMO tool to predict scores for each motif at every position in the three intervals (5′ intron flank, microexon, 3′ intron flank). We next computed the maximum match within each of the three intervals and concatenated the values to obtain 351 total features.

Known microexon splicing factor features

In addition to PWM features from CisBP we also computed some motif-matching features for RBPs that have known roles in alternative microexon splicing, namely NOVA (motif: YCAY), CELF (motif: YGCUYK), STAR (motif: YUAAY), SRRM4 (motif: UGC) and SRSF11 (motifs: UCUCU or CUCUCU).

To obtain our features, we separately counted the number of matches of each of these motifs in the microexon, as well as the ±300 nucleotides of flanking intronic sequence. Note that Y indicates C or U, and K indicates G or U.

SRRM4 is a known microexon splicing factor with a well-defined binding region.

Thus, we computed additional SRRM4 features as binary indicators of the presence of any of the motifs UGCUGC, CUGCUG, GUGCC, GUGCU or UGGACG in each 3 nucleotide bin between 60 nucleotides and 6 nucleotides upstream of the start of the microexon.

Enriched k-mers

To obtain additional sequence features that are predictive of known microexons, we identified RNA hexamers enriched or depleted in the 100 nucleotide intronic flanking sequence upstream of known microexons versus the background sequences.

We focused only on this 100 nucleotide interval because it is known to be the most critical for microexon splicing and extending the interval would have added too many features.

We also limited this feature set to only the 100 most enriched and 100 most depleted hexamers.

Specifically, for each known microexon or background sequence, we counted the number of occurrences of each hexamer in the 100 nucleotides upstream of the microexon in the training set.

We next identified the 100 hexamers that are the most enriched in known microexons versus background sequences and vice versa.

To compute features, we divided the 100 nucleotide upstream intronic sequence into five equal bins of 20 nucleotides and computed a binary value for each hexamer that indicates whether it is present in each bin, giving 1,000 (200 × 5) additional features.

GC content

We computed the fraction of G and C nucleotides within the microexon, within the 100 nucleotides of upstream intron immediately flanking the microexon, and within the 100 nucleotides of downstream intron immediately flanking the microexon. This gave a total of three features. The total number of features used by the model was 3,724.

Information leakage

All machine-learning models used in the input features had their trainable parameters frozen during training of the microexon code, so no information leakage between our train and test sets should be expected even in cases in which our test overlaps the training set of the component models. However, to rigorously assess potential information leakage from pretrained models, we evaluated overlap between our test sequences and training data of the pretrained models used in the features. COSSMO did not have any sequence overlap with our training set. While Seqweaver training regions overlapped 99.5% of our sequences, this does not constitute information leakage because Seqweaver was trained to predict general RBP profiles, not to classify microexons specifically. Our model learned to combine these general features into microexon-specific patterns. SVM-BP overlapped 11.0% (93) of foreground and 14.4% (2,583) of background. Excluding sequences overlapping SVM-BP training regions did not significantly affect classification performance (P = 0.42, Wilcoxon signed-rank test). Moreover, we used VastDB exon orthology annotations to identify test-set microexons whose orthologues were still present (after BLAT filtering) in the training set; removal of these microexons did not significantly affect performance (P = 0.11, Wilcoxon signed-rank test).

Model training

We trained the neural microexon classifier using gradient-boosted decision trees with a logistic objective function implemented in XGBoost v.1.5.0. The model outputs the probability that a candidate exon is a neurally expressed microexon.

In addition to the held-out test set described above, which was used only to assess the model, we set aside 10% of the training data as a validation set to fit model hyperparameters and determine when to stop training.

The model was trained for a maximum of 1,000 boosting rounds.

Four metrics were computed on the validation set after every boosting round (classification error, area under the curve, AUPRC and log-loss); following the XGBoost convention that the last metric supplied in eval_metric governs early stopping, validation log-loss served as the early-stopping criterion.

Training stopped when the validation log-loss failed to improve for 100 consecutive boosting rounds, and predictions were generated from the best-performing iteration.

Remaining hyperparameters were a learning rate (eta) of 0.1, a maximum tree depth (max_depth) of 4, maximum Δstep (max_delta_step) of 7, minimum loss reduction to make a further partition on a leaf node (gamma) of 0.5, L2 regularization (lambda) of 0.1 and L1 regularization (alpha) of 0.8.

To account for the imbalance between the positive and background sets, scale_pos_weight was set to the ratio of negative to positive training examples. After training, we evaluated performance on the held-out human test set using the AUPRC metric.

SHAP analysis

To interpret the contributions of individual features and feature groups to the microexon classification model, we performed SHAP (Shapley additive explanations) analysis⁷⁸ using TreeExplainer with a background dataset of 2,000 randomly sampled instances from the background set. SHAP values were computed for all positively labelled foreground microexons against this background distribution. Feature-level importance was visualized using a SHAP summary plot displaying the top 30 features ranked by their impact on model output, with individual datapoints coloured by feature value to reveal the direction and magnitude of each feature’s effect. To provide a higher-level view, SHAP values were aggregated into interpretable feature groups (splice site strength, GC content, branch point features, RBP motifs, CisBP, COSSMO, Seqweaver, nSR100 and hexamers) by summing the per-instance SHAP contributions of constituent features within each group, and the resulting local group SHAP values were visualized as a grouped summary plot to illustrate the distribution and directionality of each feature group’s contribution across microexons.

ISM

To simulate the effect of genetic variation on the recognition of microexons, we performed ISM, evaluating every possible single nucleotide substitution (relative to the human or mouse reference) in each of the 344 human and 329 mouse neuronal microexons across all regions evaluated by the model. For each substitution, we recomputed the input features from the altered sequence and reapplied the microexon classification model. We then computed the predicted score differences as:

\[\Delta \mathrm{code}\,\mathrm{score}=\frac{{\mathrm{logit}(p}_{\mathrm{alt}})-\mathrm{logit}({p}_{\mathrm{ref}})}{S}\]

where p_alt is the model prediction using the alternative allele, p_ref is using the reference allele and S is a scaling factor computed for each microexon. Scaling factors were defined as the maximum absolute difference between logit(p_alt) and logit(p_ref) across all computed variants for a given microexon. All code-model Δcode scores reported throughout the study were computed using this definition.

To obtain an overall mutational profile, for each microexon, we summarized the ISM results within a ±150 nucleotide flanking intronic window by computing, at each position, the maximum and minimum Δcode score observed across all single-nucleotide substitutions. These position-wise extrema profiles were used as input for dimensionality reduction by principal component analysis and clustering using Ward’s minimum variance criterion (hclust R function with method = “ward.D2”). As the majority of known microexons used to train the code-model clustered within cluster 1 (Extended Data Fig. 7c), we defined code-selected microexons as known and newly detected microexons assigned to cluster 1 with a silhouette score ≥ 0.1. Microexons located within uncharacterized non-coding RNA genes (clone-based RP11 identifiers) were also excluded at this step. This yielded 369 code-selected microexons, comprising 310 of the 344 known human neuronal microexons and 59 of the 292 newly detected human neuronal microexons.

Hexamer motif analyses

We extracted all hexamers that overlap each variant across microexons and 300 nucleotides of their flanking 5′/3′ intronic sequences, grouping them into coordinate interval sets. Microexon sequences were divided into quartiles and flanking intronic regions into 20 nucleotide bins. For each of the 4,096 possible hexamers across each interval, we estimated a difference of their associated Δcode scores with background values as z scores. We selected the top 100 hexamers with higher absolute z-score values across all analysed intervals for further analysis. We used the Biopython⁷⁹ implementation of the Smith–Waterman algorithm (pairwise2.align.localms function) to perform pairwise alignments across the selected hexamers, providing a matching score of 1 and penalties of −1, −2 and −1 for mismatches, gap openings and gap extensions, respectively. We used the obtained alignment scores to compute pairwise sequence distances (D_ij) as:

\[{D}_{{ij}}=1-\frac{{S}_{{ij}}}{{S}_{\max }}\]

where S_ij is the alignment score and S_max is the maximum alignment score among hexamer pairs.

Moreover, we calculated Euclidean distances between the mean Δcode score values of hexamers. To cluster hexamers, we applied Ward’s minimum variance criterion over a compound distance metric—a weighted sum of both hexamer sequence and Δcode score distances—with respective weights of 0.4 and 0.6. To derive motifs from selected hexamers within clusters, we used clustalw⁸⁰ to perform multiple-sequence alignments, and we represented the frequency of each nucleotide within the alignments as sequence logos, which were generated using the seqLogo R library.

Mouse minigene analyses

We compiled a set of minigene assay results for ten different mouse microexons from previously published data and additional experiments^29,81,82. The tested variants were introduced within or near the microexons and their flanking intronic regulatory regions, including UGC and C/U-rich motifs. We used our model to evaluate the effect of these variants in silico, and we conducted correlation analyses to evaluate the prediction performance.

Complementary variant effect predictors

To complement microexon code–model predictions, we evaluated variants within ±300 nucleotides of neuronal microexons present across analysed cohorts and represented on massively parallel splicing assay (MaPSy) data using the deep-learning splicing (CNN-based) predictors DeltaSplice (v.1.0.1)³⁷ and Pangolin (v.1.0.1)³⁶. DeltaSplice predictions were obtained using the brain-trained model and Pangolin predictions were obtained using tissue-specific models trained on brain, heart, liver and testis. For each variant, we extracted predicted changes in splice-site usage (ΔSSU) at both the 5′ and 3′ splice sites flanking the microexon. For downstream analyses, including correlation analyses with microexon MaPSy measurements, we considered the ΔSSU prediction with the largest absolute magnitude across the two splice sites as a single per-variant summary. Using these models, variants were classified as impactful if this maximum absolute ΔSSU was ≥0.05. Moreover, we used the Combined Annotation–Dependent Depletion (CADD; v.1.6)⁸³ framework to obtain PHRED-scaled CADD scores, a genome-wide measure of predicted variant deleteriousness that integrates diverse annotations (including sequence conservation and regulatory features), and classified variants as having predicted pathogenic potential if their PHRED-scaled CADD score was ≥10.

Analysis of variant-dependent PSI changes in human brain data

We processed matched WGS and bulk brain RNA-seq data from 276 individuals in the PsychENCODE BrainGVEX cohort (Synapse ID: syn4590909). RNA-seq samples were processed using MicroExonator, using analogous parameters used for processing the rest of human and mouse samples. Single-nucleotide variants were extracted from a ±300 nucleotide window around microexons. For each variant–microexon pair, we calculated an average ΔPSI by subtracting the mean PSI of samples without the variant from the mean PSI of samples carrying the variant. We assessed concordance between predicted variant effects and observed splicing changes by computing Pearson’s correlation between per-variant prediction scores and average ΔPSI, restricting to variants above the 95th percentile or below the 5th percentile according to predicted effects on splicing for each model (Extended Data Fig. 8b). We repeated the correlation analyses after removing a single outlier corresponding to a canonical splice-site (+1) mutation (Extended Data Fig. 8b).

Processing of microexon massively parallel splicing assays

We used publicly available microexon massively parallel splicing assay (MaPSy) data³³ to calculate correlations with Δcode scores and to derive a data-driven threshold for classifying predicted impactful variants. The MaPSy dataset comprises reporter-derived ΔPSI measurements for variants introduced within a fixed window surrounding neuronal microexons and was generated across multiple SRRM4 induction levels; we used the highest SRRM4 induction condition (most neuronal like) for calibration. For each assayed variant, we computed the corresponding absolute Δcode score (|Δcode score|) by scoring the reference and alternate sequences using the microexon code-model, and computing predicted score differences as described above for the ISM analyses.

To map code-model Δ scores to experimentally observed splicing effects, we scanned thresholds over |Δcode score| between 0 and 0.5 in increments of 0.01. For each threshold t, variants were partitioned into those with |Δcode score| ≥ t versus <t, and we computed the mean absolute MaPSy effect size (mean |ΔPSI|) for each group. We summarized the separation between groups as the mean difference in |ΔPSI| (Extended Data Fig. 8d). Statistical significance of the difference in |ΔPSI| between groups at each threshold was assessed using a two-sided Mann–Whitney U-test, with P values corrected for multiple testing across thresholds using a Bonferroni correction. We defined the impactful-variant Δ threshold as the |Δcode score| cut-off at which the MaPSy measurements indicate a mean difference in splicing effect of at least 5 percentage points (difference in mean |ΔPSI| ≥ 5.0; |Δcode score| = 0.29).

ASD variant analyses

To assess the impact of microexon-linked genetic variation on ASD risk, we processed whole-genome sequence data provided by SSC, SPARK and MSSNG. We filtered the original VCF according to the following BCF expression to obtain all high-quality variant calls: FILTER=“PASS” & DP>=10 & ((GT=“het” & (TYPE=“snp” | TYPE=“mnp”) & GQ>=99 & AD[:1]/DP>=0.3 & AD[:1]/ DP<0.8) | (GT=“het” & TYPE=“indel” & GQ>=90 & AD[:1]/ DP>=0.3 & AD[:1]/ DP<0.8) | (GT=“AA” & GQ>=25 & AD[:1]/ DP>=0.8)).

We used BCFtools (v.1.20)⁸⁴ to compute individualized consensus genomes from the hg38 reference genome and the filtered VCF files, which were then used as input to the microexon model, resulting in microexon code scores across all analysed genomes. Moreover, we extracted variants overlapping with microexon regions to compute single variant Δcode score predictions, which were evaluated with the same procedure described for ISM analyses. Only code-selected microexons (defined above) were considered for these analyses.

Discordant sibling-pair processing

We computed the differences in microexon code scores between 2,666 discordant sibling pairs from SSC and SPARK (ASD affected − unaffected sibling). To reduce the contributions of variants that were evenly distributed between ASD-affected and unaffected individuals, we implemented an imbalance filter: for each microexon, we grouped sibling-pair Δcode scores by their absolute value and required that at least two-thirds of sibling pairs sharing a given absolute Δcode score had the same sign (that is, ≥2/3 positive or ≥2/3 negative). We ranked the obtained microexon code score differences and extracted genomic variants associated with microexon code score differences within the top decile of the most positive and negative scores. For each variant, we calculated the maximum observed AF across general population cohorts sourced from gnomADv3 (ref. ⁴⁴) and used this value to classify variants as ultrarare (max AF < 1 × 10⁻⁵) or rare (1 × 10⁻⁵ ≤ max AF < 1 × 10⁻²). Variants without frequency information were considered as common for downstream analyses. Furthermore, we assessed the evolutionary conservation of sites impacted by selected variants by calculating their mean phyloP score from the 100-vertebrate track within a ±5 nucleotide range and classified variants as overlapping highly conserved sequence if the mean phyloP score was ≥1.5. This comparison was restricted to intronic variants. Of the 3,028 variant entries analysed in Fig. 4c, 76 exonic entries were removed and 13 variants appearing more than once were collapsed to the entry with the largest absolute Δcode score, leaving 2,939 unique intronic variants. For visualization purposes, absolute Δcode scores of 1 were omitted from Fig. 4b.

Assessing recurrent impact of rare and ultrarare microexon variants across cohorts

To determine whether a subset of microexons is recurrently predicted to be impacted by rare or ultrarare variants across independent ASD cohorts, we applied a uniform variant-processing and scoring workflow to SSC, SPARK and MSSNG. Within each cohort, we computed microexon-level Δcode score values for each ASD-affected individual relative to a cohort-specific baseline, defined as the mean code score across unaffected siblings for the same microexon. We then decomposed each individual’s microexon Δcode score into the relative contributions of rare and ultrarare variants by aggregating single-variant |Δcode scores|. For each cohort, microexons were ranked by absolute Δcode score magnitude and we selected those in the upper decile (top 10%) of the most extreme differences, additionally requiring that the observed microexon Δ signal was predominantly driven by either rare or ultrarare variants (≥50% of the aggregated |Δcode scores| attributable to variants in the corresponding frequency class). The overlap across cohorts was quantified by direct set intersection of the selected microexon lists, including the three-way (SSC ∩ SPARK ∩ MSSNG) overlap. Statistical significance of the observed triple overlap was assessed by permutation testing (10⁷ permutations), in which set sizes were preserved, and empirical P values were computed as (r + 1)/(n + 1), where r is the number of permutations yielding an overlap equal to or greater than that observed and n is the total number of permutations. The reported empirical P values were adjusted for multiple testing using the Bonferroni method.

Cumulative impact of predicted microexon-disrupting variants

For each individual in discordant sibling pairs from SSC and SPARK, we computed a cumulative impact score as the sum of |Δcode score| values for variants classified as impactful, restricting to microexons that passed the ‘imbalance filter’ described above. Statistical significance of paired differences in cumulative impact scores between affected and unaffected siblings was assessed using a one-sided Wilcoxon signed-rank test. To assess robustness to the choice of impactful-variant threshold, we repeated the analysis across a range of |Δcode score| cut-offs expressed as multiples of the s.d. of the |Δcode score| distribution for common variants (Fig. 4f). We performed analogous cumulative-impact analyses using the CNN-based splicing predictors described above. For each model, variants were classified as impactful using the threshold corresponding to a ≥5% predicted splicing change: |Δcode score| = 0.29 for the microexon code-model, calibrated against MaPSy measurements (Extended Data Fig. 8d), and maximum |ΔSSU| ≥ 0.05 for DeltaSplice. These thresholds are indicated by vertical dashed lines in Fig. 4f. Moreover, we tested whether cumulative code score impact predicts ASD status using a sibling-pair-matched conditional logistic regression model, with ASD status as the outcome and sex as a covariate.

Generalized disequilibrium test

We performed generalized disequilibrium tests using RV-GDT (v.1.0)⁴⁶ on sets of impactful variants derived from microexons that passed the ‘imbalance filter’ across analysed discordant sibling pairs. We evaluated three variant sets: (1) impactful variants prioritized by the microexon code-model; (2) impactful variants prioritized by DeltaSplice; and (3) the intersection of variants prioritized by both models. RV-GDT analyses were performed across a larger dataset comprising 25,132 individuals from 7,391 families drawn from SSC (n = 2,419), SPARK (n = 2,033) and MSSNG (n = 2,939). For code-model-derived variant sets, we incorporated variant weights to prioritize microexons of which the regulatory profiles most closely matched the neuronal microexon code, using the cluster 1 silhouette score derived from the ISM profile clustering described above. Empirical P values were obtained through within-pedigree phenotype permutations, as implemented in RV-GDT. For each tested variant set and allele-frequency stratum, we performed 10,000 phenotype permutations per run and repeated the procedure across 10 independent runs (total of 100,000 permutations per condition), summarizing results as the mean empirical P value across runs. Multiple-testing correction across RV-GDT tests was performed using the Benjamini–Hochberg procedure, with significance assessed at FDR = 0.05.

Classification of superset variant carriers by previous genetic findings

Individuals carrying at least 1 variant from the 72-variant RV-GDT superset were identified across the SSC, SPARK and MSSNG cohorts (n = 629 carriers). Carriers were classified as idiopathic or as harbouring one or more previously detected candidate disruptive variants using participant-level annotations available for the SSC and MSSNG cohorts⁸; equivalent annotations were not available for SPARK participants. A participant was classified as idiopathic only if no previous finding was recorded in any of nine categories: SNV/indel findings with dominant inheritance; SNV/indel findings with recessive inheritance; chromosomal abnormalities; genomic disorders; large or gene-rich CNVs; structural variants disrupting ASD-associated genes; uniparental isodisomy; tandem repeat expansions; and mitochondrial variants. Of the 629 carriers, 162 had annotations available, of which 133 (82%) were classified as idiopathic and 29 (18%) had at least one previous genetic finding. Applying the same classification to all annotated individuals (n = 7,519) yields 85.7% idiopathic and 14.3% with previous genetic findings.

Protein–protein interaction network analyses

We used the STRING database (v.11.5)⁸⁵ to identify potential physical interactions between human proteins encoded by genes containing neuronal microexons. We configured our query to consider only experiments and databases as knowledge sources. The obtained network represents all found interactions with at least medium confidence, allowing a maximum of five interactions within the first shell. We clustered the obtained network using the multilevel clustering algorithm from clusterMaker2 (v.2.3.4)⁸⁶. To label clusters, we used AutoAnnotate (v.1.4) to process GO terms associated to each protein. Both clusterMaker2 and AutoAnnotate were operated as apps within Cytoscape.

Computational workflows and figure preparation

Computational workflows were implemented with Snakemake (v.7.32.4) and custom Python and R notebooks run in JupyterLab (v.4.0.8). Some figures were prepared using GraphPad Prism (v.10).

Statistics and reproducibility

No statistical methods were used to predetermine sample size. Experiments were not randomized, and investigators were not blinded to allocation during experiments and outcome assessment.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

RNA-seq data generated in this study are publicly available in the NCBI Gene Expression Omnibus (GEO) under accession GSE268269. A complete inventory of previously deposited and publicly available RNA-seq datasets reanalysed in this study is provided in Supplementary Table 1. In summary, human snRNA-seq data generated by the BRAIN Initiative were obtained from the NeMO Archive under accession dat-ptqn5sg. Mouse cortical scRNA-seq data from refs. ^23,24 were obtained from the NCBI Sequence Read Archive under accessions SRP061902 (primary visual cortex) and SRP150473 (primary visual and anterior lateral motor cortex). Microexon annotations were obtained from VastDB (https://vastdb.crg.eu) and GENCODE (https://www.gencodegenes.org). Gene expression data were obtained from the GTEx Portal (v8; https://gtexportal.org). Mouse embryonic and postnatal brain RNA-seq data were obtained from the ENCODE portal (https://www.encodeproject.org), and mouse whole-cortex developmental RNA-seq data were obtained from the NCBI Sequence Read Archive under accession SRP055008. Matched WGS and bulk brain RNA-seq data from the PsychENCODE BrainGVEX cohort were obtained from Synapse under accession syn4590909. Microexon massively parallel splicing assay (MaPSy) measurements were obtained from the NCBI Gene Expression Omnibus under accession GSE276143. RBP motifs were obtained from CisBP-RNA (http://cisbp-rna.ccbr.utoronto.ca). Population allele frequencies were obtained from gnomAD v.3 (https://gnomad.broadinstitute.org). The human (hg38) and mouse (mm10/GRCm38) reference genome assemblies were obtained from the UCSC Genome Browser (https://hgdownload.soe.ucsc.edu). WGS and associated phenotype data from the SSC and SPARK are available through SFARI Base (https://base.sfari.org) after registration and approval by the Simons Foundation. Data from the MSSNG resource are available through the MSSNG portal (https://research.mss.ng) after registration and approval. These datasets are subject to controlled access under the terms of their respective data-use agreements and cannot be redistributed by the authors.

Code availability

All custom code used for the analyses and to generate the figures in this study is available at GitHub (https://github.com/geparada/microexon-regulatory-code). The repository release corresponding to this publication is archived at Zenodo⁸⁷ (https://doi.org/10.5281/zenodo.21940006).

References

1. Irimia, M. et al. A highly conserved program of neuronal microexons is misregulated in autistic brains. Cell 159, 1511–1523 https://doi.org/10.1016/j.cell.2014.11.035 (2014).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

2. Li, Y. I., Sanchez-Pulido, L., Haerty, W. & Ponting, C. P. RBFOX and PTBP1 proteins regulate the alternative splicing of micro-exons in human brain transcripts. Genome Res. 25, 1–13 https://doi.org/10.1101/gr.181990.114 (2015).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

3. Gonatopoulos-Pournatzis, T. & Blencowe, B. J. Microexons: at the nexus of nervous system development, behaviour and autism spectrum disorder. Curr. Opin. Genet. Dev. 65, 22–33 https://doi.org/10.1016/j.gde.2020.03.007 (2020).

   Article  CAS  PubMed  Google Scholar 

4. Mackensen, T. & Irimia, M. From tiny exons to big insights: the expanding field of microexons. Annu. Rev. Genom. Hum. Genet. 26, 77–102 https://doi.org/10.1146/annurev-genom-121323-103648 (2025).

   Article  CAS  Google Scholar 

5. Eyring, K. W. & Geschwind, D. H. Three decades of ASD genetics: building a foundation for neurobiological understanding and treatment. Hum. Mol. Genet. 30, R236–R244 https://doi.org/10.1093/hmg/ddab176 (2021).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

6. Vorstman, J. A. S. & Scherer, S. W. Contemplating syndromic autism. Genet. Med. 25, 100919 https://doi.org/10.1016/j.gim.2023.100919 (2023).

   Article  CAS  PubMed  Google Scholar 

7. Satterstrom, F. K. et al. Large-scale exome sequencing study implicates both developmental and functional changes in the neurobiology of autism. Cell 180, 568–584 https://doi.org/10.1016/j.cell.2019.12.036 (2020).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

8. Trost, B. et al. Genomic architecture of autism from comprehensive whole-genome sequence annotation. Cell 185, 4409–4427 https://doi.org/10.1016/j.cell.2022.10.009 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

9. Fu, J. M. et al. Rare coding variation provides insight into the genetic architecture and phenotypic context of autism. Nat. Genet. 54, 1320–1331 https://doi.org/10.1038/s41588-022-01104-0 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

10. Dominguez-Alonso, S., Carracedo, A. & Rodriguez-Fontenla, C. The non-coding genome in autism spectrum disorders. Eur. J. Med. Genet. 66, 104752 https://doi.org/10.1016/j.ejmg.2023.104752 (2023).

    Article  CAS  PubMed  Google Scholar 

11. Raj, B. & Blencowe, B. J. Alternative splicing in the mammalian nervous system: recent insights into mechanisms and functional roles. Neuron 87, 14–27 https://doi.org/10.1016/j.neuron.2015.05.004 (2015).

    Article  CAS  PubMed  Google Scholar 

12. Vuong, C. K., Black, D. L. & Zheng, S. The neurogenetics of alternative splicing. Nat. Rev. Neurosci. 17, 265–281 https://doi.org/10.1038/nrn.2016.27 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

13. Quesnel-Vallieres, M. et al. Misregulation of an activity-dependent splicing network as a common mechanism underlying autism spectrum disorders. Mol. Cell 64, 1023–1034 https://doi.org/10.1016/j.molcel.2016.11.033 (2016).

    Article  CAS  PubMed  Google Scholar 

14. Gonatopoulos-Pournatzis, T. et al. Autism-misregulated eIF4G microexons control synaptic translation and higher order cognitive functions. Mol. Cell 77, 1176–1192 https://doi.org/10.1016/j.molcel.2020.01.006 (2020).

    Article  CAS  PubMed  Google Scholar 

15. Parras, A. et al. Autism-like phenotype and risk gene mRNA deadenylation by CPEB4 mis-splicing. Nature 560, 441–446 https://doi.org/10.1038/s41586-018-0423-5 (2018).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

16. Polinski, P. et al. A highly conserved neuronal microexon in DAAM1 controls actin dynamics, RHOA/ROCK signaling, and memory formation. Nat. Commun. 16, 4210 https://doi.org/10.1038/s41467-025-59430-w (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

17. Dupas, S. J. et al. Splice isoform-perturbation coupled to single cell transcriptome profiling reveals functions of microexons in neurogenesis and autism-linked pathways. Nat. Commun. 17, 1217 https://doi.org/10.1038/s41467-025-67931-x (2026).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

18. Voineagu, I. et al. Transcriptomic analysis of autistic brain reveals convergent molecular pathology. Nature 474, 380–384 (2011).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

19. Parikshak, N. N. et al. Genome-wide changes in lncRNA, splicing, and regional gene expression patterns in autism. Nature 540, 423–427 https://doi.org/10.1038/nature20612 (2016).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

20. Gandal, M. J. et al. Transcriptome-wide isoform-level dysregulation in ASD, schizophrenia, and bipolar disorder. Science 362, eaat8127 https://doi.org/10.1126/science.aat8127 (2018).

21. Parada, G. E. et al. MicroExonator enables systematic discovery and quantification of microexons across mouse embryonic development. Genome Biol. 22, 43 https://doi.org/10.1186/s13059-020-02246-2 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

22. Jorstad, N. L. et al. Transcriptomic cytoarchitecture reveals principles of human neocortex organization. Science 382, eadf6812 https://doi.org/10.1126/science.adf6812 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

23. Tasic, B. et al. Adult mouse cortical cell taxonomy revealed by single cell transcriptomics. Nat. Neurosci. 19, 335–346 https://doi.org/10.1038/nn.4216 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

24. Tasic, B. et al. Shared and distinct transcriptomic cell types across neocortical areas. Nature 563, 72–78 https://doi.org/10.1038/s41586-018-0654-5 (2018).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

25. Pan, Q. et al. Quantitative microarray profiling provides evidence against widespread coupling of alternative splicing with nonsense-mediated mRNA decay to control gene expression. Genes Dev. 20, 153–158 https://doi.org/10.1101/gad.1382806 (2006).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

26. Olier, I. et al. Transformational machine learning: Learning how to learn from many related scientific problems. Proc. Natl Acad. Sci. USA 118, e2108013118 https://doi.org/10.1073/pnas.2108013118 (2021).

27. Bretschneider, H., Gandhi, S., Deshwar, A. G., Zuberi, K. & Frey, B. J. COSSMO: predicting competitive alternative splice site selection using deep learning. Bioinformatics 34, i429–i437 https://doi.org/10.1093/bioinformatics/bty244 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

28. Zhou, J. et al. Whole-genome deep-learning analysis identifies contribution of noncoding mutations to autism risk. Nat. Genet. 51, 973–980 https://doi.org/10.1038/s41588-019-0420-0 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Barash, Y. et al. Deciphering the splicing code. Nature 465, 53–59 https://doi.org/10.1038/nature09000 (2010).

    Article  ADS  CAS  PubMed  Google Scholar 

30. Ule, J. & Blencowe, B. J. Alternative splicing regulatory networks: functions, mechanisms, and evolution. Mol. Cell 76, 329–345 https://doi.org/10.1016/j.molcel.2019.09.017 (2019).

    Article  CAS  PubMed  Google Scholar 

31. Raj, B. et al. A global regulatory mechanism for activating an exon network required for neurogenesis. Mol. Cell 56, 90–103 https://doi.org/10.1016/j.molcel.2014.08.011 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

32. Gonatopoulos-Pournatzis, T. et al. Genome-wide CRISPR-Cas9 interrogation of splicing networks reveals a mechanism for recognition of autism-misregulated neuronal microexons. Mol. Cell 72, 510–524 https://doi.org/10.1016/j.molcel.2018.10.008 (2018).

    Article  CAS  PubMed  Google Scholar 

33. Bonnal, S., Bajew, S., Martinez-Corral, R. & Irimia, M. Core splicing architecture and early spliceosomal recognition determine microexon sensitivity to SRRM3/4. Nat. Struct. Mol. Biol. 32, 2022–2034 https://doi.org/10.1038/s41594-025-01634-1 (2025).

    Article  CAS  PubMed  Google Scholar 

34. Wang, D. et al. Comprehensive functional genomic resource and integrative model for the human brain. Science 362, aat8464 https://doi.org/10.1126/science.aat8464 (2018).

35. Roussos, P. et al. PsychENCODE at 10: from genomic maps to mechanistic insights in mental illness. Neuron 113, 4099–4102 https://doi.org/10.1016/j.neuron.2025.11.027 (2025).

    Article  CAS  PubMed  Google Scholar 

36. Zeng, T. & Li, Y. I. Predicting RNA splicing from DNA sequence using Pangolin. Genome Biol. 23, 103 https://doi.org/10.1186/s13059-022-02664-4 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

37. Xu, C. et al. Reference-informed prediction of alternative splicing and splicing-altering mutations from sequences. Genome Res. 34, 1052–1065 https://doi.org/10.1101/gr.279044.124 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

38. Corioni, M., Antih, N., Tanackovic, G., Zavolan, M. & Kramer, A. Analysis of in situ pre-mRNA targets of human splicing factor SF1 reveals a function in alternative splicing. Nucleic Acids Res. 39, 1868–1879 https://doi.org/10.1093/nar/gkq1042 (2011).

    Article  CAS  PubMed  Google Scholar 

39. Darbelli, L. & Richard, S. Emerging functions of the Quaking RNA-binding proteins and link to human diseases. Wiley Interdiscip. Rev. RNA 7, 399–412 https://doi.org/10.1002/wrna.1344 (2016).

    Article  CAS  PubMed  Google Scholar 

40. Galarneau, A. & Richard, S. Target RNA motif and target mRNAs of the Quaking STAR protein. Nat. Struct. Mol. Biol. 12, 691–698 https://doi.org/10.1038/nsmb963 (2005).

    Article  CAS  PubMed  Google Scholar 

41. Ray, D. et al. A compendium of RNA-binding motifs for decoding gene regulation. Nature 499, 172–177 https://doi.org/10.1038/nature12311 (2013).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

42. Hayakawa-Yano, Y. et al. An RNA-binding protein, Qki5, regulates embryonic neural stem cells through pre-mRNA processing in cell adhesion signaling. Genes Dev. 31, 1910–1925 https://doi.org/10.1101/gad.300822.117 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

43. Fischbach, G. D. & Lord, C. The Simons Simplex Collection: a resource for identification of autism genetic risk factors. Neuron 68, 192–195 https://doi.org/10.1016/j.neuron.2010.10.006 (2010).

    Article  CAS  PubMed  Google Scholar 

44. Karczewski, K. J. et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature 581, 434–443 https://doi.org/10.1038/s41586-020-2308-7 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

45. 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 https://doi.org/10.1093/nar/gky1016 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

46. He, Z. et al. The rare-variant generalized disequilibrium test for association analysis of nuclear and extended pedigrees with application to Alzheimer disease WGS data. Am. J. Hum. Genet. 100, 193–204 https://doi.org/10.1016/j.ajhg.2016.12.001 (2017).

    Article  CAS  PubMed  Google Scholar 

47. Kim, S. et al. Alternatively spliced mini-exon B in PTPδ regulates excitatory synapses through cell-type-specific trans-synaptic PTPδ-IL1RAP interaction. Nat. Commun. 16, 4415 https://doi.org/10.1038/s41467-025-59685-3 (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

48. Han, K. A. et al. Specification of neural circuit architecture shaped by context-dependent patterned LAR-RPTP microexons. Nat. Commun. 15, 1624 https://doi.org/10.1038/s41467-024-45695-0 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

49. Roth, J. F. et al. Systematic analysis of alternative exon-dependent interactome remodeling reveals multitasking functions of gene regulatory factors. Mol. Cell 83, 4222–4238 https://doi.org/10.1016/j.molcel.2023.10.034 (2023).

    Article  CAS  PubMed  Google Scholar 

50. Choi, Y. et al. SALM5 trans-synaptically interacts with LAR-RPTPs in a splicing-dependent manner to regulate synapse development. Sci. Rep. 6, 26676 https://doi.org/10.1038/srep26676 (2016).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

51. Cornejo, F., Cortes, B. I., Findlay, G. M. & Cancino, G. I. LAR receptor tyrosine phosphatase family in healthy and diseased brain. Front. Cell Dev. Biol. 9, 659951 https://doi.org/10.3389/fcell.2021.659951 (2021).

    Article  PubMed  PubMed Central  Google Scholar 

52. Xiong, H. Y. et al. The human splicing code reveals new insights into the genetic determinants of disease. Science 347, 1254806 https://doi.org/10.1126/science.1254806 (2015).

    Article  ADS  MathSciNet  CAS  PubMed  Google Scholar 

53. Jaganathan, K. et al. Predicting splicing from primary sequence with deep learning. Cell 176, 535–548 https://doi.org/10.1016/j.cell.2018.12.015 (2019).

    Article  ADS  CAS  PubMed  Google Scholar 

54. Cheng, J., Celik, M. H., Kundaje, A. & Gagneur, J. MTSplice predicts effects of genetic variants on tissue-specific splicing. Genome Biol. 22, 94 https://doi.org/10.1186/s13059-021-02273-7 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

55. Yang, Y. et al. Single-cell long-read sequencing in human cerebral organoids uncovers cell-type-specific and autism-associated exons. Cell Rep. 42, 113335 https://doi.org/10.1016/j.celrep.2023.113335 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

56. Chau, K. K. et al. Full-length isoform transcriptome of the developing human brain provides further insights into autism. Cell Rep. 36, 109631 https://doi.org/10.1016/j.celrep.2021.109631 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

57. Litman, A. et al. Variant-resolved prediction of context-specific isoform variation with a graph-based attention model. Cell Genom. 6, 101126 https://doi.org/10.1016/j.xgen.2025.101126 (2026).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

58. Margolis, M. P. et al. From variants to mechanisms: neurogenomics in the post-GWAS era. Neuron 113, 3509–3529 https://doi.org/10.1016/j.neuron.2025.10.014 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

59. Wilfert, A. B. et al. Recent ultra-rare inherited variants implicate new autism candidate risk genes. Nat. Genet. 53, 1125–1134 https://doi.org/10.1038/s41588-021-00899-8 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

60. Frankish, A. et al. Gencode 2021. Nucleic Acids Res. 49, D916–D923 https://doi.org/10.1093/nar/gkaa1087 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

61. Tapial, J. et al. An atlas of alternative splicing profiles and functional associations reveals new regulatory programs and genes that simultaneously express multiple major isoforms. Genome Res. 27, 1759–1768 https://doi.org/10.1101/gr.220962.117 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

62. UniProt Consortium. UniProt: the Universal Protein Knowledgebase in 2025. Nucleic Acids Res. 53, D609–D617 https://doi.org/10.1093/nar/gkae1010 (2025).

    Article  CAS  Google Scholar 

63. Stacklies, W., Redestig, H., Scholz, M., Walther, D. & Selbig, J. pcaMethods—a Bioconductor package providing PCA methods for incomplete data. Bioinformatics 23, 1164–1167 https://doi.org/10.1093/bioinformatics/btm069 (2007).

    Article  CAS  PubMed  Google Scholar 

64. Murtagh, F. & Legendre, P. Ward’s hierarchical agglomerative clustering method: which algorithms implement Ward’s criterion? J. Classif. 31, 274–295 https://doi.org/10.1007/s00357-014-9161-z (2014).

65. ENCODE Project Consortium. The ENCODE (ENCyclopedia Of DNA Elements) Project. Science 306, 636–640 https://doi.org/10.1126/science.1105136 (2004).

    Article  ADS  CAS  Google Scholar 

66. He, P. et al. The changing mouse embryo transcriptome at whole tissue and single-cell resolution. Nature 583, 760–767 https://doi.org/10.1038/s41586-020-2536-x (2020).

    Article  ADS  CAS  PubMed  Google Scholar 

67. Weyn-Vanhentenryck, S. M. et al. Precise temporal regulation of alternative splicing during neural development. Nat. Commun. 9, 2189 https://doi.org/10.1038/s41467-018-04559-0 (2018).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

68. Raudvere, U. et al. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 47, W191–W198 https://doi.org/10.1093/nar/gkz369 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

69. Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504 (2003).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

70. Merico, D., Isserlin, R., Stueker, O., Emili, A. & Bader, G. D. Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS ONE 5, e13984 https://doi.org/10.1371/journal.pone.0013984 (2010).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

71. Kucera, M., Isserlin, R., Arkhangorodsky, A. & Bader, G. D. AutoAnnotate: a Cytoscape app for summarizing networks with semantic annotations. F1000Res. 5, 1717 https://doi.org/10.12688/f1000research.9090.1 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

72. Cheung, A. Y. et al. Isolation of MECP2-null Rett syndrome patient hiPS cells and isogenic controls through X-chromosome inactivation. Hum. Mol. Genet. 20, 2103–2115 https://doi.org/10.1093/hmg/ddr093 (2011).

    Article  MathSciNet  CAS  PubMed  PubMed Central  Google Scholar 

73. Gohr, A. & Irimia, M. Matt: Unix tools for alternative splicing analysis. Bioinformatics 35, 130–132 https://doi.org/10.1093/bioinformatics/bty606 (2019).

74. Sheth, N. et al. Comprehensive splice-site analysis using comparative genomics. Nucl. Acids Res. 34, 3955–3967 (2006).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

75. Sterne-Weiler, T., Weatheritt, R. J., Best, A. J., Ha, K. C. H. & Blencowe, B. J. Efficient and accurate quantitative profiling of alternative splicing patterns of any complexity on a laptop. Mol. Cell 72, 187–200 https://doi.org/10.1016/j.molcel.2018.08.018 (2018).

    Article  CAS  PubMed  Google Scholar 

76. Yeo, G. & Burge, C. B. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. J. Comput. Biol. 11, 377–394 https://doi.org/10.1089/1066527041410418 (2004).

    Article  CAS  PubMed  Google Scholar 

77. Corvelo, A., Hallegger, M., Smith, C. W. & Eyras, E. Genome-wide association between branch point properties and alternative splicing. PLoS Comput. Biol. 6, e1001016 https://doi.org/10.1371/journal.pcbi.1001016 (2010).

    Article  ADS  MathSciNet  CAS  PubMed  PubMed Central  Google Scholar 

78. Lundberg, S. M. & Lee, S.-I. A unified approach to interpreting model predictions. In Adv. Neural Information Processing Systems Vol. 30 (eds Guyon, I. et al.) 4765–4774 (Curran Associates, 2017).

79. Cock, P. J. et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics 25, 1422–1423 https://doi.org/10.1093/bioinformatics/btp163 (2009).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

80. Larkin, M. A. et al. Clustal W and Clustal X version 2.0. Bioinformatics 23, 2947–2948 https://doi.org/10.1093/bioinformatics/btm404 (2007).

    Article  CAS  PubMed  Google Scholar 

81. Nakano, Y. et al. A mutation in the Srrm4 gene causes alternative splicing defects and deafness in the Bronx Waltzer mouse. PLoS Genet. 8, e1002966 https://doi.org/10.1371/journal.pgen.1002966 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

82. Ohnishi, T., Shirane, M. & Nakayama, K. I. SRRM4-dependent neuron-specific alternative splicing of protrudin transcripts regulates neurite outgrowth. Sci. Rep. 7, 41130 https://doi.org/10.1038/srep41130 (2017).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

83. Kircher, M. et al. A general framework for estimating the relative pathogenicity of human genetic variants. Nat. Genet. 46, 310–315 https://doi.org/10.1038/ng.2892 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

84. Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10, giab008 https://doi.org/10.1093/gigascience/giab008 (2021).

85. Szklarczyk, D. et al. The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49, D605–D612 https://doi.org/10.1093/nar/gkaa1074 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

86. Utriainen, M. & Morris, J. H. clusterMaker2: a major update to clusterMaker, a multi-algorithm clustering app for Cytoscape. BMC Bioinform. 24, 134 https://doi.org/10.1186/s12859-023-05225-z (2023).

    Article  Google Scholar 

87. Parada, G. E. et al. Microexon code model and resources accompanying ‘A microexon regulatory code reveals autism-linked genetic variation’. Zenodo https://doi.org/10.5281/zenodo.21940006 (2026).

Download references

Acknowledgements

We thank J. Ule, S. Cordes and A. Packer, as well as members of the Blencowe laboratory, including W. Winick-Ng, E. Hoover, K. Rebelo and M. Wu, for discussions and assistance with experiments and analyses; and the Donnelly Sequencing Centre for generating sequencing data.

Funding

This work was supported by grants from the Simons Foundation (SFARI AN-AR-Gene Therapies-00893981), Canadian Institutes of Health Research, Wellcome Trust (227465/Z/23/Z) and the University of Toronto McLaughlin Centre to B.J.B.; and a grant from Autism Speaks to S.W.S.

Q.D.M. was supported by NIH/NCI Cancer Center Support (P30 CA008748) and NIH/NHGRI R01 (HG013328) grants, and a Canada CIFAR AI Chair through the Vector Institute; and G.E.P. by a Charles H. Best Postdoctoral Fellowship. H.B.

was partially supported by a Vector Institute Research Grant. J.D.L. and S.J.D. were supported by NSERC Canada Doctoral Graduate Scholarships; and J.L.F. by an NSERC Vanier Canada Graduate Scholarship. S.W.S.

holds the Northbridge Chair in Paediatric Research at the Hospital for Sick Children and University of Toronto. B.J.B. holds the University of Toronto Banbury Chair in Medical Research and the Canada Research Chair in RNA Biology and Genomics.

Author information

Author notes

1. These authors contributed equally: Guillermo E. Parada, Hannes Bretschneider

Authors and Affiliations

1. Donnelly Centre, University of Toronto, Toronto, Ontario, Canada

   Guillermo E. Parada, Hannes Bretschneider, Jack Daiyang Li, Jacob L. Fine, Steven J. Dupas, Jonathan D. Ellis, Ulrich Braunschweig, Quaid D. Morris & Benjamin J. Blencowe

2. Vector Institute, Toronto, Ontario, Canada

   Guillermo E. Parada, Hannes Bretschneider, Jacob L. Fine, Quaid D. Morris & Benjamin J. Blencowe

3. Centre for Developmental Neurobiology, King’s College London, Guy’s Campus, London, UK

   Guillermo E. Parada & Benjamin J. Blencowe

4. The Francis Crick Institute, London, UK

   Guillermo E. Parada & Benjamin J. Blencowe

5. Department of Molecular Genetics, University of Toronto, Toronto, Ontario, Canada

   Jack Daiyang Li, Jacob L. Fine, Steven J. Dupas, James Ellis, Brett Trost, Stephen W. Scherer, Quaid D. Morris & Benjamin J. Blencowe

6. Centre for Genomic Regulation, Barcelona Institute of Science and Technology, Barcelona, Spain

   Sophie Bonnal & Manuel Irimia

7. Computational Biology and Medicine Program, Weill-Cornell Graduate School, New York, NY, USA

   Taykhoom Dalal & Quaid D. Morris

8. Computational and Systems Biology Program, Memorial Sloan Kettering Cancer Center, New York, NY, USA

   Taykhoom Dalal & Quaid D. Morris

9. The Centre for Applied Genomics and Genetics and Genome Biology Program, The Hospital for Sick Children, Toronto, Ontario, Canada

   Worrawat Engchuan, Mehdi Zarrei, Brett Trost & Stephen W. Scherer

10. Developmental, Stem Cell & Cancer Biology Program, The Hospital for Sick Children, Toronto, Ontario, Canada

    Wei Wei & James Ellis

11. Department of Biochemistry, Donnelly Centre, University of Toronto, Toronto, Ontario, Canada

    Nadia Zafar & Liliana Attisano

12. The Gene Lay Institute of Immunology and Inflammation, Brigham and Women’s Hospital, Massachusetts General Hospital, and Harvard Medical School, Boston, MA, USA

    Martin Hemberg

13. Universitat Pompeu Fabra, Barcelona, Spain

    Manuel Irimia

14. ICREA, Barcelona, Spain

    Manuel Irimia

15. Prosserman Centre for Population Health Research, Lunenfeld-Tanenbaum Research Institute, Sinai Health, Toronto, Ontario, Canada

    Michael Wainberg

16. Department of Psychiatry, University of Toronto, Toronto, Ontario, Canada

    Michael Wainberg

17. Department of Computer Science, University of Toronto, Toronto, Ontario, Canada

    Michael Wainberg & Quaid D. Morris

18. Biostatistics Division, Dalla Lana School of Public Health, University of Toronto, Toronto, Ontario, Canada

    Michael Wainberg

19. Molecular Medicine Program, The Hospital for Sick Children, Toronto, Ontario, Canada

    Brett Trost

Authors

1. Guillermo E. Parada
2. Hannes Bretschneider
3. Jack Daiyang Li
4. Jacob L. Fine
5. Steven J. Dupas
6. Jonathan D. Ellis
7. Ulrich Braunschweig
8. Sophie Bonnal
9. Taykhoom Dalal
10. Worrawat Engchuan
11. Mehdi Zarrei
12. Wei Wei
13. Nadia Zafar
14. Martin Hemberg
15. Liliana Attisano
16. James Ellis
17. Manuel Irimia
18. Michael Wainberg
19. Brett Trost
20. Stephen W. Scherer
21. Quaid D. Morris
22. Benjamin J. Blencowe

Contributions

Project conceptualization: B.J.B., H.B. and Q.D.M. Study design: H.B., G.E.P., U.B., Q.D.M. and B.J.B. Machine learning model development: H.B. and Q.D.M. Data analysis: G.E.P., J.L.F., S.J.D., U.B., H.B., J.D.L., B.T. and T.D.

with input from W.E., M.Z., M.I., M.H., M.W., S.W.S., Q.D.M. and B.J.B. Experiments: J.D.L., J.D.E., S.B., W.W. and N.Z. Manuscript figure preparation: G.E.P., J.D.L., J.L.F., U.B., H.B. and B.J.B. Manuscript writing: B.J.B., G.E.P., H.B.

and J.D.L., with input from all of the authors. Supervision: B.J.B., Q.D.M., M.I., M.H., J.E. and L.A. Funding acquisition: B.J.B., Q.D.M. and S.W.S.

Corresponding authors

Correspondence to Quaid D. Morris or Benjamin J. Blencowe.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature thanks Neelroop Parikshak and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data figures and tables

Extended Data Fig. 1 Discovery and characterization of new microexons.

A) Detection of neuronal microexons in human RNA-seq data as a function of the percentage of total data used. Lines represent a non-linear model fitted to the average number of microexons detected across 1,000 iterations (Methods). Shaded areas show distribution ranges between the 5^th and 95^th percentiles of the detection rate, representing 90% of the data obtained across all iterations. Dotted lines represent an extrapolation of the model, estimating the number of microexons that would be detected if 10 times more data were analysed. The detection rate for known microexons plateaus after utilizing ~10% of the total data, whereas the detection rate for newly discovered microexons, as extrapolated from the curve, suggests the potential for the discovery of additional microexons with increased data analysis. B) Enrichment analysis of newly detected neuronal microexons in bulk and single-nucleus RNA-seq data. For each microexon, mean PSI was compared across four contexts: bulk non-neural datasets, bulk neural datasets, sn (single-nucleus) non-neural cell types and sn neural cell types. Boxplots summarize per-event mean PSI values (red, neural contexts; grey, non-neural), showing significant enrichment in bulk neural samples but no significant difference between sn-neural and sn-non-neural cell types (n = 292 new neuronal microexons, quantified across 3,524 neural and 228 non-neural bulk RNA-seq samples, and across 330 snRNA-seq neural and 27 non-neural pseudobulks). Brackets indicate two-sided Mann-Whitney U tests with Bonferroni-correction; ***, P ≤ 0.001; ns, not significant. C) Enrichment analysis of newly detected and code-selected microexons (i.e., the subset of new microexons whose ISM profiles resemble those of known neuronal microexons) across the sample contexts described in panel b (n = 59 code-selected microexons, quantified across the same samples as in b). This subset of new microexons shows higher mean PSI in both bulk neural samples and sn-neural cell types (Bonferroni-adjusted; *, P ≤ 0.05; ***, P ≤ 0.001). In b and c, box plots show the median (centre line), the 25th and 75th percentiles (bounds of the box) and whiskers extending to the most extreme value within 1.5 × the interquartile range; minima and maxima are the most extreme points plotted. D, E) Venn diagrams showing overlaps between neuronal microexons detected in the different RNA-seq datasets indicated from human (d) and mouse (e). F) Genomic and protein-context features of known neuronal microexons (blue) and new microexons identified as neuronal-differential (orange). Bars indicate the percentage of events located in UTRs, predicted to elicit nonsense-mediated decay (NMD), overlapping UniProt-annotated protein domains, overlapping UniProt-annotated intrinsically disordered regions, or located in other regions of coding sequence (CDS). G) Heatmap showing hierarchical clustering of PSI measurements obtained from total detected microexons across various human RNA-seq samples. Rows and columns represent microexons and samples, respectively, each classified as indicated by the colour keys. Sample clusters are labelled A to L. PSI measurements derived from snRNA-seq data, after pooling nuclei to generate pseudobulk-based PSI measurements, are indicated. The colour scale from blue to red represents PSI values ranging from 0 to 100, with blue indicating low PSI values and red indicating high PSI values. Grey tiles represent missing PSI values. Only PSI measurements based on a minimum of five reads are displayed. H) Unrooted dendrogram representing the distances between various sample clusters in panel g, labelled ‘A’ to ‘L’. Each cluster is associated with a word cloud highlighting the most common terms related to the sample’s metadata. The branches connecting the clusters are colour-coded to indicate their membership with each sample cluster. I) Heatmap showing hierarchical clustering based on PSI values of newly detected neuronal microexons representing three data source categories: “Non-neural RNA-seq”, “Neural RNA-seq”, and “Neural snRNA-seq”. Each row in the heatmap corresponds to a specific microexon, and these are divided into two categories: NMD-inducing and Not NMD-inducing. Blue, indicates low PSI values; white, intermediate PSI values; red, high PSI values; grey, missing PSI values.

Extended Data Fig. 2 Validation of candidate NMD-inducing microexons detected in mouse and human neurons.

A) Cultured mouse glutamatergic neurons were treated with the translation inhibitor cycloheximide (CHX) and candidate NMD-inducing microexons were quantified by RT-PCR. Lanes show untreated control and +CHX conditions. B) Cultured human neurons were treated with SMG1i (an inhibitor of the NMD factor SMG1) or with DMSO as a control, and predicted NMD-inducing microexons were quantified by RT-PCR. Lanes show control and +SMG1i conditions. In both panels, schematics to the left of each gel indicate the expected RT-PCR products corresponding to isoforms with and without microexon inclusion, and PSI values from band quantification are shown below each gel. In a, two independent experiments were performed with similar results. The experiment in b was performed once; two of the events shown, BCO1 and AFG1L, were validated in an independent hiPSC line (WTC11).

Extended Data Fig. 3 Human and mouse microexons with cell type specific splicing patterns.

A) Heatmap showing new human microexons clustered based on differential splicing patterns across cortical neuronal sub-types. Rows indicate microexon-containing gene and columns indicate cell types. A microexon was considered specific to a given cell type if at least two pseudobulk pools associated with that cell type exhibited a PSI with a lower bound binomial 95% confidence interval 20 units greater or lower than the median PSI across all pseudobulk pools. *, microexons with significant cell type-specific inclusion. B-D) Schematics illustrating open reading frames (ORFs) affected by inclusion of new neuronal microexons in mouse with detected cell type-specific inclusion patterns. Annotated UniProt domains are displayed for each ORF. Triangles indicate microexon location within ORFs, triangle colours indicate cell type-specific inclusion patterns. Red asterisks denote microexons predicted to trigger NMD (when included individually). B, Actn4, in which new microexons overlap a Vinculin (VCL)-interacting domain. C, Scfd2, in which 11 new microexons are detected, 3 of which exhibit neuronal subtype-specific splicing and one of which has the potential to alter the reading frame and induce NMD; a subset of these microexons is predicted to modify the Sec1 domain, a common feature of a family of proteins that regulate vesicle fusion and synaptic transmission. D, Klf7, in which a 3 nt microexon is predicted to directly affect one of the three C2H2-type zinc finger domains of this Kruppel-like transcription factor.

Extended Data Fig. 4 Predicted NMD-inducing microexons show developmental anti-correlation with gene expression levels and temporally-dynamic inclusion patterns.

A) Spearman correlations between microexon PSI and host-gene expression (TPM, transcripts per million) across stages of mouse brain development (Methods). Predicted NMD-inducing microexons exhibit a significantly lower mean correlation between PSI and TPM than frame-preserving/non-NMD microexons with comparable developmental PSI changes (mean ρ = −0.104 and mean ρ = 0.113, respectively; n = 34 NMD-inducing and 360 non-NMD microexons; **, P = 0.0046, two-sided Mann-Whitney U test). Box plots as in Extended Data Fig. 1b. B) Examples of developmental trajectories of predicted NMD-inducing microexons that were identified to be dynamically regulated through mouse embryonic development. These include known (top row) and newly detected microexons (bottom row), across mouse brain tissues (forebrain, midbrain and hindbrain). PSI is plotted against developmental time (days post coitum); points indicate mean PSI at each time point with binomial 95% confidence intervals. Lines show fitted models: linear regression (blue), binomial GLM (red) and piecewise binomial GLM (green). C) Representative predicted NMD-inducing microexons showing strong anti-correlations between microexon inclusion and host-gene expression during mouse forebrain development. Blue traces indicate microexon PSI values and red traces indicate host-gene expression as TPM across developmental time, measured in days post coitum (dpc). Points show the mean and error bars the 95% confidence interval across biologically independent samples. Spearman correlation coefficients (ρ) and associated P values are shown for each microexon-host gene pair.

Extended Data Fig. 5 Functional enrichment analysis of genes harbouring known and new human neuronal microexons.

Each node within the network represents an enriched gene function category (FDR < 0.05) from Gene Ontology and Reactome; edge width indicates overlap between categories (combined Jaccard–overlap coefficient, edges > 0.375 are shown).

Extended Data Fig. 6 Training sets, performance, and characterization of microexon code model.

A) Bivariate distribution of 3′ and 5′ microexon splice strengths based on MaxEnt Scores (MES) for known microexons (foreground) and candidate sequences resembling microexons (background). The background sequences were sampled from the larger set of intronic microexon-resembling sequences to match the splice site strength distribution of the foreground sequences as closely as possible to the foreground set. The scatter plot with overlaid contour lines represents the density of occurrences, for both foreground and background microexons. The marginal distributions of MES3 and MES5 are displayed at the top and right sides of the plot, respectively. B) SHapley Additive exPlanation (SHAP) analysis of feature importance. Summary of how different feature categories contribute to model predictions. Each dot represents one microexon from the test set. The x-axis shows the SHAP value, which quantifies each feature’s contribution to the model output: positive values shift the prediction toward classifying the sequence as a true neuronal microexon, while negative values shift toward classifying it as background. Dot colour indicates the feature value for that example (red: high/above average; blue: low/below average). For example, high Seqweaver feature values (red dots) tend to have positive SHAP values, indicating that strong predicted RBP binding increases the likelihood of microexon classification. C) SHAP feature importance plot for the top 30 most-important features. An explanation of feature names is provided in Supplementary Table 6. The plot reveals which features contributed positively or negatively to the model predictions. For example, for positively contributing features, high feature values correlate with positive SHAP values (e.g. alt_ss_bp_score, the microexon’s branch point score), and for negatively contributing features, high feature values correlate with negative SHAP values (e.g. alt_ss_bp_dist, the microexon’s branch point’s distance). D) Precision-recall curve for code model predictions while training on multi-species or human-only datasets. The multi-species model (blue line) exhibits an area under the precision-recall curve (AUPRC) of 0.72, indicating increased performance in estimating the likelihood of neuronally-spliced microexons compared to the human-only model (orange line), which has an AUPRC of 0.50. E) Controls for potential information leakage impacting code model classifier performance. Boxplot showing the distribution of area under the precision-recall curve (AUPRC) values obtained for the microexon code-model when evaluated on held-out chromosomes (n = 5 held-out chromosomes). AUPRC values were compared for all microexons (“All events”), after excluding flanking sequences in the test set that overlap SVM-BP training regions (“Without SVM-BP”), or after excluding microexon sequences in the test set that are annotated in VastDB (using a liftOver-based approach) as orthologs in the training-set (“Without orthologs”). Two-sided Wilcoxon signed-rank tests assess whether differences in AUPRC values between the different test-set comparisons are significant; ns, non-significant (All events vs Without SVM-BP, P = 0.42; All events vs Without Orthologs, P = 0.11). F) Code-model scores relate to neuronal splicing patterns. Comparison of code-model prediction scores between microexons defined by RNA-seq analysis as not neuronal (i.e. mostly constitutively spliced) and neuronally spliced (n = 491 constitutively non-neuronal and 493 neuronal microexons, restricted to those falling on the held-out test chromosomes used in e). Brackets indicate two-sided Mann-Whitney U tests with Bonferroni correction (***, P ≤ 0.001). G) Correlation between code-model scores and microexon splicing levels across human tissues. Boxplots show per-sample Pearson correlation (r) between code-model prediction scores and microexon percent spliced-in (PSI) values, grouped by tissue or cell-type categories from VastDB. Each point represents one RNA-seq sample (n = 19 neuronal and 77 non-neuronal samples). Brackets indicate two-sided Mann-Whitney U tests comparing correlations between Embryonic Brain and Neuronal groups, and all other analysed cell and tissue types (***, P < 10⁻⁹). Box plots in e-g as in Extended Data Fig. 1b.

Extended Data Fig. 7 Analysis of in silico mutagenesis profiles across human known and newly identified microexons.

A) In silico saturation mutagenesis profiles (ISMs) across regions used to train and test the microexon code model. The model was trained and tested using ±300 nt regions (Fig. 3a), however only splice site flanking windows up to ±150 nt are shown as these capture most of the variance in mean Δ code scores. Red line, mean ∆ code score changes in response to all three possible nucleotide substitutions at each analysed position associated with all analysed known and new microexons (n = 369); light grey trace indicates Δ code score distributions representing the 5^th−95^th percentile range; dark grey trace indicates the 25^th−75^th percentile interquartile range. C1 and C2, flanking constitutive exons. B) Magnitude of ISM ∆ code scores correlates with vertebrate sequence conservation. Scatterplot showing a comparison of ISM ∆ code scores (representing the mean of the maximum absolute score changes upon mutation of each nucleotide to every other nucleotide, for each position and for all analysed microexons) versus the mean vertebrate conservation score (phyloP scores derived from 100 vertebrate genomes for the same nucleotide positions used to generate ISM measurements). Each point represents a given position within 300 nt windows directly upstream and downstream of microexon splice sites. Red line, least-squares linear regression fit. Pearson’s r and P values are indicated. C) Principal Component Analysis of microexon ISM profiles reveals two distinct clusters: Cluster 1 (red) primarily represents annotated microexons, whereas Cluster 2 (blue) mostly represents new microexons. D) Principal Component Analysis of density map highlighting new microexons that belong to Cluster 1. The colour scale indicates the Silhouette scores for these microexons, which provide a measure of their clustering assignment validity. E) Microexon code score distributions for annotated and new microexons that form Cluster 1 and Cluster 2, where the code scores are based on the hg38 reference genome. The colour scale indicates Silhouette scores (two-sided Mann-Whitney U test, P < 1 × 10⁻¹⁵). F) ISM profiles of annotated and new microexons that are represented by Cluster 1 and Cluster 2 defined in panel c. Traces as in a.

Extended Data Fig. 8 Validation and benchmarking of microexon code model scores.

A) Correlation analysis directly comparing microexon Δ code scores from our model with experimentally determined changes in Percent Spliced In (PSI) for variants introduced within or near mouse microexons and their flanking intronic regulatory regions. Blue line, least-squares linear regression fit; shaded area, 95% confidence interval of the fit (n = 49 variants). Spearman and Pearson correlation coefficients are shown. Details for minigene experiments are provided in Supplementary Table 7. B) Correlation between predicted impact of variants on splicing using different models and measurements of variant impact on splicing (∆PSI) from analysing patient brain RNA-seq data. ∆PSI values were determined as the average PSI for a microexon in patient RNA-seq data harbouring a specific variant, compared to the average PSI for the same microexon in individuals without that variant. For each model, correlations were calculated only using variant effect scores from predicted high-impact variants (top and bottom 5th percentiles of each model’s score distribution, red dots). Light blue circles indicate a variant that disrupts the 5′-splice site of a microexon in the NRG1 gene. Red lines indicate least-squares linear regression fits, and shaded bands indicate 95% confidence intervals. Two-sided Pearson’s r and Bonferroni-corrected P values are indicated within each panel (n = 128 variants per model). C) Change in correlation between variant-dependent effects on microexon splicing (PSI) measured using a massively parallel splicing assay (MaPSy)³³, following induction of SRRM4 expression, and predicted effects of the variants on splicing using the different models indicated. Measurements of variant-dependent changes in percent spliced in (ΔPSI) were generated using MaPSy data for microexons in response to extensive mutagenesis of surrounding sequences, and changes in model predictions for splicing effects were calculated for the same variants and microexons using the microexon code model (this study), DeltaSplice and Pangolin models. MaPSy experiments were performed in HEK293 cells expressing GFP (baseline control) or expressing SRRM4 at increasing levels (Low, Mid, High) following induction with doxycycline, to establish a neuronal-like splicing context. For the code-model, effects are ∆ code scores; for DeltaSplice and Pangolin, effects are predicted ΔPSI levels. Log₂-fold changes in Pearson correlation for each condition are shown relative to GFP, which is normalized to zero. Raw correlations and corresponding P values are reported in Supplementary Table 8. D) Calibration of microexon Δ code scores with MaPSy-based measurements of ∆PSI. Using measurements from the MaPSy, we computed—across a sliding threshold representing absolute Δ code scores (x-axis)—the mean absolute ΔPSI difference associated with variants above or below the threshold (y-axis), compared using two-sided Mann-Whitney U tests. Point size indicates −log10(adjusted P) and colour indicates significance after Bonferroni correction (red, significant; grey, not significant). Vertical dashed line indicates the threshold at which the mean |ΔPSI| difference exceeds 5.0.

Extended Data Fig. 9 Analysis of the role of QKI in microexon regulation.

A) Validation of QKI knockdown efficiency in N2A cells by qRT-PCR 72 h following transfection of siRNAs targeting QKI transcripts. The ddCT method was used to calculate fold change with Gapdh as a normalization control. Bars show the mean of three biologically independent experiments (n = 3), and error bars indicate the standard deviation. B) Validation of QKI knockdown efficiency in N2A cells by Western blot 72 h following siRNA transfection. Detection of alpha-tubulin is used as a loading and recovery control. The blot shown is from a single experiment. Knockdown efficiency was additionally quantified across three biologically independent experiments by qRT-PCR (a). C) RT-PCR assays monitoring inclusion levels of five microexons harbouring upstream intronic QKI motifs (ACUAAY), following knockdown of QKI. Bars show the mean of three biologically independent experiments (n = 3), with the individual replicate values overlaid as points. Error bars indicate the standard deviation. Statistical significance was assessed using two-sided Student’s t tests, with P values adjusted for multiple comparisons across the five microexons using the Bonferroni method; **, P ≤ 0.01; ***, P ≤ 0.001. D) RNA motif map showing the occurrence of QKI motifs (ACUAAY) within surrounding introns of alternative microexons that are repressed by QKI (n = 41, blue line), or alternative microexons that do not change upon QKI knockdown (n = 89, red line). E) Positional impact of in silico mutagenesis of predicted QKI binding sites on microexon code scores. Hexamers matching QKI motifs surrounding all code-analysed microexons were identified. For each QKI-motif hexamer, the largest absolute Δ code score from in silico mutagenesis was determined. Boxplots show the distribution of these scores across positional bins surrounding microexons (x-axis). Box colour indicates the median z-score from mutating QKI-related hexamer binding sites relative to mutating all other hexamers at the same position bin. The largest effects cluster 20–50 nt upstream of microexons, in agreement with our experimental data supporting QKI-dependent repression of microexons mediated through binding to these regions. Box plots show the 25th and 75th percentiles (bounds of the box) and whiskers extending to 1.5× the interquartile range; minima and maxima are the most extreme points plotted.

Extended Data Fig. 10 Analysis of variants in ASD affected individuals predicted to impact microexon alternative splicing, and microexon coverage by exome capture panels.

A) Distributions of variants contributing the sibling-pair absolute ‘genome-level’ Δ code scores for microexons. ‘Relative contribution’ (x-axis) represents the fraction of the genome-level Δ score attributable to different numbers of variant, computed from single-variant predictions and restricted to variants present in only one sibling (y-axis). Analyses shown in A (and panels B, C) focus on the top decile genome-level ∆ code scores between siblings. B) Histogram of the numbers of variants that contribute to the Δ code scores for each sibling pair. C) Distribution of the maximum relative contribution of any single variant associated with any microexon within each sibling pair and microexon. Values near 1 indicate cases where one variant explains all or nearly all of the genome-level ∆ scores. D) Significant differences in Δ code score values calculated for individual variants. Groups of variants with different allele frequency intervals. *, P < 0.05; **, P < 0.01; ***, P < 0.0001; ns, not significant, two-sided Mann-Whitney U test corrected by multiple comparisons using the Benjamini-Hochberg method. E) Significant correlation between Δ code scores and variant rarity (−log₁₀[allele frequency]) (Spearman’s ρ = −0.30, P < 1.0 × 10⁻¹⁵). F) Predictions for effects of common vs rare variants across different tested splicing models. Violin plots compare distributions of model-predicted effects for variants near microexons, stratified by allele frequency (‘rare’ set includes ultrarare variants). n (rare, common): code model 2,067 and 961; DeltaSplice 3,727 and 2,519; Pangolin brain 3,380 and 1,830; Pangolin liver 1,843 and 932; Pangolin heart 2,354 and 1,168; Pangolin testis 2,251 and 1,116 microexon–variant predictions. The microexon code model and DeltaSplice show significantly larger effects for (ultra-)rare variants; Pangolin models do not. Significance determined using two-sided Mann–Whitney U tests with Bonferroni correction across the six models (**, P ≤ 0.01; ****, P ≤ 0.0001; ns, P > 0.05). Violins show the full data distribution; box plots in d and f as in Extended Data Fig. 1b. G) Overlaps of microexons in SSC (red), SPARK (blue) and MSSNG (green) WGS data that have top decile Δcode scores (affected vs. unaffected individuals) driven by ultrarare variants. Numbers indicate counts of microexons in each set overlap, with areas approximately proportional to set sizes. Three-way cohort overlap is significantly greater than expected by chance (1.85 times the expected overlap; Bonferroni-adjusted empirical P = 4 × 10⁻⁷, one-sided permutation test). H) Numbers of ASD-affected individuals carrying at least one variant with a predicted ∆PSI > 5.0, and the subset of these individuals with ASD-associated variants as scored using Rare-Variant Generalized Disequilibrium Test (RV-GDT) analyses. Indicated variant groups (x-axis) were stratified by the highest allele frequency observed across general population cohorts⁴⁴: All variants; Rare + Ultrarare (AF < 1 × 10⁻²); Ultrarare (AF < 1 × 10⁻⁵). Y-axis indicates the numbers and percentages of ASD-affected individuals carrying at least one predicted impactful variant that is absent from their unaffected sibling(s). I) Minigene assays validating predicted effects of rare and ultrarare variants on microexon splicing in cultured human neurons. Top, schematic summarizing three different minigene reporter constructs containing microexons and native surrounding intronic sequences from the human LYSMD4, MAP4K4, and UGGT2 genes. C1 and C2, flanking constitutive exons; μ, microexon; blue arrows, positions of variants relative to microexons. Bottom, RT-PCR analysis of the splicing of minigene transcripts: WT, wild type reporter; Mut, reporter harbouring mutation as indicated above. Exon diagrams beside each gel indicate expected (based on size) spliced isoforms. Asterisks indicate unexpected RT-PCR products, which were not confirmed by sequencing and were excluded from PSI quantification. PSI, percent-spliced-in. PSI values are averages from three biologically independent replicates. PSI measurements from individual replicate experiments are provided in Supplementary Table 11. J) Coverage of known (green) and newly identified (orange) human neuronal microexons by three commercial exome capture panels: xGen Exome Research Panel v2 (IDT), Twist Human Core Exome Panel (Twist Bioscience) and TruSeq DNA Exome Panel (Illumina). Bars show the percentage of microexons whose genomic coordinates overlap panel target regions; absolute numbers of microexons covered in each category are indicated above the bars.

Supplementary information

Reporting Summary (download PDF )

Supplementary Table 1 (download CSV )

List of RNA-seq datasets analysed for the discovery and quantification of microexons.

sample_ID, ENCODE, SRA and NEMO accession codes; species, human or mouse; experiment type, the class of sequencing data analysed; source, the repository holding each accession; condition, the tissue, cell type or experimental condition described by the original study; sample_type, whether the sample was analysed as neural or as part of the non-neural background.

Bulk RNA-seq and TRAP-seq rows correspond to individual sequencing runs, whereas single-cell and single-nucleus RNA-seq rows correspond to complete datasets.

Supplementary Table 2 (download CSV )

Newly identified human and mouse microexons. ME-ID, MicroExonator event ID (chromosome, strand, start and end); MIC_length, length of each microexon; MIC_type, new microexon with or without neuronal inclusion pattern. Coordinates for human and mouse microexons are based on the hg38 and mm10 genome assemblies, respectively.

Supplementary Table 3 (download CSV )

List of microexons dynamically regulated during mouse brain development. Rows correspond to quantified microexons identified by ME-ID, the MicroExonator event ID, which comprises the genomic coordinates (chromosome, strand, start and end). Other columns indicate number of observations pooled across all analysed tissues (n_points), types and number of tissues in which the microexon was detected (tissues, n_tissues), Spearman correlation between PSI and TPM computed across all tissues and timepoints for each given microexon with nominal and FDR-corrected P values (rho_psi_tpm, pval_psi_tpm, pval_psi_tpm_fdr) and associated significance labels (psi_tpm_significant, is_neg_sig, is_pos_sig), whether the microexon is predicted to trigger nonsense-mediated decay (is_NMD); microexon length in nucleotides (length), associated gene identifiers (gene_id, common_name) and whether the microexon is classified as developmentally-regulated according to regression model criteria (is_dev_regulated). Microexon coordinates are based on the mm10 genome assembly.

Supplementary Table 4 (download CSV )

GO categories significantly enriched among genes with human and mouse microexons. GO.ID, GO IDs; log2Enr, log₂[enrichment] for each category; FDR, corresponding adjusted P value; delta_log2Enr, differential log₂[enrichment] before and after including the set of new microexons; gained_0.05, the latter category was found to be significant after including new microexons.

Supplementary Tables 5 and 6 (download XLSX )

Supplementary Table 5: summary of features used to train the microexon code model. Descriptions and the number of features associated with each category are indicated. Supplementary Table 6: comprehensive list and descriptions of features used to train the microexon code model. Columns indicate the following. feature_name, internal feature identifier; Publication Name, human-readable label used in figures; Source, tool or resource from which each feature was derived; dtype, feature data type; Description, information on how each feature was computed.

Supplementary Table 7 (download XLSX )

Mouse minigene microexon splicing assays: tested variants, splicing effect measurements (ΔPSI) and corresponding code model Δ scores. The following are indicated. The microexon identifiers (ME-ID and VastDB ID), minigene, the gene from which the tested microexon and its flanking sequences were derived, tested variants and their corresponding genomic positions (variant, chrom, position, ref and alt); assembly, genome assembly (mm10); the experimentally measured outcomes: ref_psi and alt_psi, PSI percentages for the reference and alternative alleles, respectively, on a 0–100 scale; dpsi, the alternative-allele PSI minus the reference-allele PSI, reported in percentage points; code model Δ scores, the corresponding microexon code model prediction for each variant, computed as the difference between model predictions for the alternative and reference alleles, scaled by the maximum absolute difference observed across all variants evaluated for that microexon (Methods); pmid, PubMed unique identifier number referencing the source of each minigene experiment; cell_line, the cell line used for each assay; change_direction_experimental_15_cutoff, the direction of the measured splicing change at a 15-percentage-point cutoff (increased inclusion, decreased inclusion or no change); number_UGC_mutated and number_UCUC_mutated, the numbers of UGC and UCUC motifs mutated in each construct, with the substitutions indicated.

Supplementary Table 8 (download CSV )

Correlations between MaPSy-measured variant effects on ΔPSI and corresponding microexon Δ code scores. Correlation statistics from comparing MaPSy-derived ΔPSI values and predicted variant-effect scores from the indicated models across control (GFP) and SRRM4 induction levels (LOW, MID, HIGH). The columns indicate SRRM4 induction level (SRRM4_induction), model, correlation method (Pearson or Spearman), correlation coefficient (cor) and corresponding P value (p.value).

Supplementary Table 9 (download CSV )

Correlation statistics from comparisons between variant-dependent Δ code scores and in vivo microexon ΔPSI measurements from PsychENCODE BrainGVEX brain RNA-seq data. Analyses were restricted to variants within ±300 nucleotides of neuronal microexons and the top/bottom 10% of variants with the largest absolute predicted effects for each model (Threshold). Rows indicate the model used (Model_name) and whether canonical splice-site mutations were excluded or included (SS_mutations). For each condition, the number of analysed variants (N_variants) is indicated along with the Pearson’s correlation coefficient and P value (Pearson_R, Pearson_p) and corresponding Bonferroni-adjusted P value and significance label (Pearson_p_bonf, Pearson_signif). Spearman correlation statistics are similarly indicated (Spearman_rho, Spearman_p, Spearman_p_bonf, Spearman_signif).

Supplementary Table 10 (download CSV )

Conditional logistic regression of ASD status on cumulative-impact scores in discordant sibling pairs. Conditional logistic regression results from testing whether cumulative-impact scores predict ASD status in discordant sibling pairs when controlling for sex. Columns indicate model terms, odds ratios, standard errors (on the log-odds scale), Z statistics, two-sided P values and 95% confidence intervals.

Supplementary Table 11 (download XLSX )

Human iNeuron microexon minigene splicing assays. Columns indicate the following: Reporter, gene name and whether construct is WT or mutant (mut) with information on which exon(s) are quantified for PSI levels; Variant, the substitution introduced in the mutant construct (chromosome:position:reference:alternative), blank for wild-type constructs; PSI R1–R3, PSI measurements for each of three replicates; Average, mean PSI across three replicates; SD, standard deviation for the three replicate experiments.

Supplementary Table 12 (download CSV )

Microexon Δ code scores and additional effect predictions for analysed variants from SSC, SPARK and MSSNG cohorts across analysed human neuronal microexons. ME-ID, MicroExonator event ID for each microexon (chromosome, strand, start and end).

Other columns indicate: event, VastDB event ID; variant, variant identifier (chromosome:position:reference:alternative); geneName, gene symbol; MAX_AF, maximum AF observed across general population cohorts; var_freq_type, variant frequency class (common, rare or ultra-rare); var_type, variant class (SNV or indel); impactful_code_model and impactful_DeltaSplice, whether the variant was classified as impactful by the microexon code model and by DeltaSplice, respectively; RV_GDT_super_set, whether the variant belongs to the set of variants found to be significantly associated with ASD status by RV-GDT; MIC_group, whether the microexon was previously annotated (known_mic) or newly detected in this study (new_mic); MIC_dist, distance in nucleotides between the variant and the microexon; var_MIC_region, the region in which the variant is located relative to the microexon: ex, microexon; up, upstream intronic region; dn, downstream intronic region; c1, upstream constitutive exon splice site; c2, downstream constitutive exon splice site; mt_code_score and wt_code_score, microexon code model scores for the mutant and wild-type sequences; logit_scaling_factor, per-microexon normalization factor used for the ∆ code scores reported in this study; delta_code_model, the corresponding ∆ code score; logit_scaling_factor_all_vars, factor calculated across all released variants, including MSSNG; delta_code_model_all_vars, corresponding all-variant ∆ code score; delta_DeltaSplice, corresponding variant effect predicted by DeltaSplice; delta_PangolinBrain, delta_PangolinLiver, delta_PangolinHeart and delta_PangolinTestis, variant effects predicted by Pangolin in the indicated tissues; CADD_phred, PHRED-scaled CADD score; assembly, genome assembly.

The 0.29 threshold used to define impactful_code_model applies only to delta_code_model; delta_code_model_all_vars has no calibrated impactful threshold. Both ∆ code score columns are bounded to [-1, 1].

Scores cover 632 microexons; missing values denote unavailable predictions. Coordinates are based on the hg38 genome assembly.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Parada, G.E., Bretschneider, H., Li, J.D. et al. A microexon regulatory code reveals autism-linked genetic variation. Nature (2026). https://doi.org/10.1038/s41586-026-11119-w

Download citation

• Received: 29 May 2024

• Accepted: 03 September 2026

• Published: 14 October 2026

• Version of record: 14 October 2026

• DOI: https://doi.org/10.1038/s41586-026-11119-w