A microexon regulatory code reveals autism-linked genetic variation

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:

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:

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:

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:

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

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:

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:

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.