Main

Autism spectrum disorder (ASD) is a common neurodevelopmental disorder (NDD), with a childhood prevalence of close to 2%18. The last decade of genetic studies has yielded hundreds of risk genes, consistent with extraordinary aetiological heterogeneity1,2,3,4. More than 100 high-confidence mutations have been associated with ASD in genetic studies2,5,6,7,8,9,10,11,12,13. These rare, usually de novo mutations with large effect sizes are expected to account for 10–15% of ASD cases19,20, whereas common genetic variation is predicted to explain at least 50% of genetic risk21,22,23. Overall, ASD shows a complex genetic architecture, with a substantial component derived from a collection of distinct, rare disorders with overlapping clinical features.

Despite its genetic heterogeneity, postmortem transcriptome analysis has revealed consistent changes in most individuals with idiopathic ASD, as well as individuals with a specific syndromic form of ASD, (dup)15q11–13 (refs. 14,15,16,17). However, the mechanisms by which distinct mutations can lead to convergent molecular pathology, and whether convergence occurs across rare forms of ASD remains unknown. Understanding how these processes develop is complicated by the fact that the expression of most ASD risk genes peaks during fetal development, yet gene expression studies in human brain are conducted after this critical developmental window has ended16,24,25,26,27. Several lines of evidence, including genetic28, genomic29,30, neuroimaging31 and neuropathology32, indicate that early neurodevelopment has an essential role in the development of ASD.

The advent of stem cell-based in vitro systems enables high-fidelity modelling of human brain development in NDDs33,34,35,36,37,38,39,40,41. Most studies have investigated relatively small numbers of lines from individuals with idiopathic ASD42, or focused on individual mutations, which have demonstrated the utility of human induced pluripotent stem (hiPS) cell-based systems to study the impact of ASD genetic risk on neurodevelopment43,44,45,46. In addition, recent work has highlighted the power of studying many genes in parallel, identifying evidence for molecular convergence using CRISPR-based perturbations of ASD risk genes on control genetic backgrounds47,48,49. However, studies of lines derived from affected individuals are a major gap in the field.

Here we profile a large cohort of hiPS cell lines, starting from 96 hiPS cell lines ascertained from individuals with 8 different mutations associated with ASD and 11 individuals with idiopathic ASD and 30 lines derived from 25 matched control participants. From each of these lines, we derived neural organoids using a guided differentiation approach to create human cortical organoids (hCOs)50,51. Using orthogonal analytic approaches, we find evidence for convergence during early neuronal differentiation in those with genetically defined forms, but no significant shared signal across the idiopathic cases. We identify and characterize a downregulated, chromatin and transcriptional network that contains several ASD risk genes, including members of the SWI–SNF complex. This network is predicted to drive the observed changes in downstream gene expression associated with these ASD susceptibility mutations. We use CRISPRi to validate the effects of many putative network drivers on downstream gene expression, which includes downregulation of important neurodevelopmental pathways including many ASD susceptibility genes.

Generation of hCOs from hiPS cell lines

We reprogrammed somatic cells (Methods) to generate hiPS cells from a cohort of individuals with eight different mutations associated with ASD: (1) 22q11.2 deletion, a 1.5–3-megabase (Mb) deletion that leads to a constellation of variably present symptoms including heart defects, craniofacial features, intellectual disability, ASD (roughly 20%) and psychosis (roughly 25%)52,53 (n = 18)44; (2) 22q13.3 deletion, known as Phelan–McDermid syndrome, a deletion spanning 130 kilobases (kb) to 9 Mb that includes SHANK3, among other genes and presents with developmental delay, hypotonia and impaired social interactions54 (*n *= 4); (3) 15q13.3 deletion, a 1.5-Mb deletion that presents variably with intellectual disability and epilepsy55 (n = 3); (4) 16p11.2 deletion, a roughly 600-kb deletion presenting variably with intellectual disability, motor impairments, communication deficits and ASD56 (n = 4); (5) 16p11.2 duplication, which can share intellectual disability and ASD phenotypes with the reciprocal deletion57 (n = 4); (6) Timothy syndrome, which is characterized by variants within the CACNA1C gene and presents with syndactyly, prolonged QT interval, ASD and intellectual disability58,59 (n = 2)45,60,61; (7) PCDH19-related disorder, which is associated with epilepsy and can also include intellectual disability and ASD62,63 (n = 2) and (8) the SHANK3 R522W mutation (n = 1), a point mutation associated with neurodevelopmental risk64; mutations in SHANK3 share behavioural phenotypes with the larger 22q13.3 deletion61,65,66. We also profiled individuals with idiopathic ASD with no known pathogenic variants (n = 11) and unaffected individuals (n = 25) (Fig. 1a and Supplementary Table 1) for a total of 74 individuals. For several individuals we used several hiPS cell lines to assess reproducibility, amounting to a total of 96 hiPS cell lines (Fig. 1a,b and Extended Data Fig. 1).

Fig. 1: Experimental workflow and validation.

a, Schematic workflow going from hiPS cells to cortical organoids to sequencing data. The number of hiPS cell lines and individuals (in parentheses) for each form of ASD and controls is indicated. b, Schematic representation of hCO differentiations derived from each hiPS cell line for two forms of ASD (16p11.2 deletion and 22q13.3 deletion). The other forms of ASD can be found in Extended Data Fig. 1. c, Spearman’s correlation of gene expression between samples from the same time point and form of ASD that were derived either from different lines (red) or from the same line (blue). The sample sizes for each group (number of differentiations) are as follows: day 25: control lines, n = 46; 15q13.3del, n = 7; 16p11.2del, n = 5; 16p11.2dup, n = 4; 22q11.2del, n = 13; 22q13.3del, n = 11; idiopathic, n = 12; PCDH19, n = 2; SHANK3, n = 2; and Timothy syndrome, n = 3. Day 50: control lines, *n *= 53; 15q13.3del, n = 7; 16p11.2del, n = 6; 16p11.2dup, n = 4; 22q11.2del, n = 15; 22q13.3del, n = 12; idiopathic, n = 15; PCDH19, n = 2; SHANK3, n = 2; and Timothy syndrome, n = 3. Day 75: control lines, n = 54; 15q13.3del, n = 7; 16p11.2del, n = 6; 16p11.2dup, n = 4; 22q11.2del, n = 23; 22q13.3del, n = 12; idiopathic, n = 15; PCDH19, n = 2; SHANK3, n = 2; and Timothy syndrome, n = 2. Day 100: control lines, n = 50; 15q13.3del, n = 6; 16p11.2del, n = 6; 16p11.2dup, n = 4; 22q11.2del, n = 15; 22q13.3del, n = 12; idiopathic, n = 15; PCDH19, n = 2; SHANK3, n = 2; and Timothy syndrome, n = 2. Boxplots show: centre, median; lower hinge, 25% quantile; upper hinge, 75% quantile; lower whisker, smallest observation greater than or equal to lower hinge −1.5× interquartile range; upper whisker, largest observation less than or equal to upper hinge +1.5× interquartile range. d, Genes within the CNVs are downregulated in deletions and upregulated in duplications, as exemplified by 16p11.2 deletion (del) and duplication (dup) and 15q13.3 deletion. 16p11.2del, 61.5% of genes; 16p11.2dup, 50% of genes and 15q13.3del, 52.9% of genes. *dreamlet P values of less than 0.005. Top illustration adapted from ref. 51, Springer Nature America; illustrations in the IP–MS and CRISPRi panels created in BioRender. Geschwind, D. (2025) (https://biorender.com/m2jkj03).

We differentiated hiPS cell lines into hCOs for 100 days51, which yields hCOs that contain developing radial glia and ultimately give rise to functional glutamatergic neurons in a process that closely resembles cerebral cortical development36,51. We performed 172 unique hCO differentiations of the 96 lines (replicates for 49% of lines; control n = 75, idiopathic ASD n = 21, 22q11.2 deletion n = 31, 22q13.3 deletion n = 15, 15q13.3 deletion n = 7, 16p11.2 deletion* n* = 9, 16p11.2 duplication n = 5, Timothy syndrome n = 3, PCDH19 mutation n = 4 and SHANK3 mutation n = 2). We performed whole-genome sequencing and removed any lines with large copy number variations (CNVs), viral integration or in which the mutation of interest was absent or further mutations were present, which yielded 70 verified lines from 55 individuals passing quality control for downstream functional genomic analyses (Methods, Extended Data Figs. 1–3 and Supplementary Table 2).

We next analysed the transcriptomes of pooled hCOs at days 25, 50, 75 and 100 of hCO differentiation (Fig. 1a,b and Extended Data Fig. 1), which provides a quantitative, unbiased metric67,68. Day 25 corresponds to an early period in cortical development36,50, consisting primarily of cycling progenitors and radial glia50,69. As differentiation progresses, the cell composition evolves to consist of intermediate progenitors, deep layer postmitotic neurons and astroglial lineage cells44,50. After stringent RNA sequencing (RNA-seq) quality control and outlier analysis to remove low-quality samples (Fig. 1a, Extended Data Figs. 1–3 and Methods), 464 samples from 55 individuals (70 lines) were retained (Supplementary Table 1).

We observed high reproducibility both between (overall mean Spearman correlation 0.96, range 0.88–0.99) and within (overall mean Spearman correlation 0.97, range 0.92–0.98) individuals, similar to published metrics36,44,51 (Fig. 1c and Methods). The largest driver of variation was stage of differentiation (Extended Data Fig. 4). Genes within the CNV boundaries were downregulated in the deletion carriers and upregulated in the duplication carriers as expected (Fig. 1d and Extended Data Fig. 2). Although several of the genes affected by point mutations showed a trend towards downregulation (PCDH19, log fold change (logFC) of −0.77, false discovery rate (FDR) of 0.45; CACNA1C: logFC = −0.22, FDR = 0.78), these changes were not significant, as expected70,71,72 (Extended Data Figs. 2 and 3).

Transcriptomic relationships across forms and time

We observed several reliable, distinct, gene-expression clusters across forms (Methods and Supplementary Table 3), with participants harbouring the 16p11.2 deletion clustering with those with PCDH19-related disorder, 16p11.2 duplication clustering with Timothy syndrome and 22q13.3 and SHANK3 mutation clustering together (Fig. 2a). We performed bootstrapping to show that these clusters were robust (Extended Data Fig. 5 and Methods). Despite the strong effect of differentiation day on transcriptional profiles, most samples from individual genetic mutations clustered together across time points, showing a consistent effect of mutation over developmental progression (Fig. 2a).

Fig. 2: Clustering and overlap of ASD forms points to increasing convergence across development.

a, Heatmap of hierarchical clustering based on correlation of logFC between different conditions (ASD form and differentiation day). Annotation bars (shades of blue) represent differentiation day and colourful boxes represent ASD form. b, The number of DEGs in each ASD form during hCO differentiation. Colours and sizes of the circles represent the number of DEGs. c, Overlap in DEGs between forms of ASD at day 25. Number and percentage of total genes are indicated. Statistics were derived from Fisher’s exact test (two-sided). ***FDR < 0.005. d, Gene ontology terms significantly enriched in unique and intersecting DEGs. Statistics were derived from clusterProfiler enrichGO. e, Overlaps between DEGs at day 25 and ASD2,8, NDD and intellectual disability risk genes91. Colour represents OR and size represents −log10FDR. Only positive significant overlaps (OR > 1 and FDR < 0.05) are shown. Statistics were derived from Fisher’s exact test (two-sided). f, Network representation of the correlation (Spearman’s rho) between logFC of different ASD forms at days 25 and 100. Colour corresponds to rho value and line thickness corresponds to the rho absolute value. g, All correlations (Spearman’s rho) between logFC of the different forms. ANOVA followed by posthoc Tukey’s HSD; ANOVA F3,140 = 10.97, P = 1.6 × 10−6, Tukey HSD: day 100 to day 25 P = 2.67 × 10−5, day 100 to day 50 P = 8 × 10−6, day 100 to day 75 P = 7 × 10−4. Boxplots show: centre, median; lower hinge, 25% quantile; upper hinge, 75% quantile; lower whisker, smallest observation greater than or equal to lower hinge −1.5× interquartile range; upper whisker, largest observation less than or equal to upper hinge +1.5× interquartile range. Statistics were derived from n = 36 correlations per time point. ***Tukey HSD for all contrasts less than 0.001. h, The number of differentially expressed meta-significant genes comparing all affected forms with control participants. Colour and size of the circle represent the number of DEGs. BMP, bone morphogenetic protein; ID, intellectual disability.

Next, we found that the largest changes in gene expression for individual forms were observed at early stages of hCO differentiation, with the largest number of differentially expressed genes (DEGs) found at day 25 in four forms (16p11.2 deletion, 979 genes; 16p11.2 duplication, 1,556 genes; PCDH19, 514 genes and Timothy syndrome, 832 genes), or at day 50 in 2 of the remaining forms (15q13.3 deletion, 616 genes and 22q13.3 deletion, 373 genes; Fig. 2b and Supplementary Table 3). Significant similarity was seen between genes differentially expressed at day 25 in samples harbouring deletion and duplication of 16p11.2 (375 out of 2,160 genes, P < 10−16, hypergeometric test) (Fig. 2c). The intersecting genes were significantly enriched for chromatin remodelling related terms (Fig. 2d and Supplementary Table 4), which have previously been associated with ASD risk genes8,24,73. Overall, genes downregulated at 25 days in these four forms of ASD were enriched for known high-confidence ASD and NDD risk genes (16p11.2del: SFARI genes odds ratio (OR) 3.0, FDR = 3.8 × 10−9; 16p11.2dup: SFARI genes OR = 3.3, FDR = 1.4 × 10−13; PCDH19: SFARI genes OR = 3.1, FDR = 1.5 × 10−5; Timothy syndrome: SFARI genes OR = 2.0, FDR = 4.3 × 10−2; 16p11.2del: brain disorder genes OR = 2.4, FDR = 6.6 × 10−3; 16p11.2dup: brain disorder genes OR = 2.2, FDR = 1.2 × 10−2; PCDH19: brain disorder genes OR = 3.2, FDR = 2.2 × 10−3; Fisher’s exact test, Fig. 2e, Extended Data Fig. 6 and Methods). Gene set enrichment analysis (GSEA) found downregulation of terms related to neural precursor proliferation (for example, ‘Beta catenin binding’, ‘Neural precursor cell proliferation’ and ‘Canonical WNT signalling pathway’) in all 4 of these forms (Extended Data Fig. 6 and Supplementary Table 5). We detected only two significantly differentially expressed genes in the individuals with idiopathic ASD (PRRC2C at day 50 and the long non-coding RNA RP11-114H21.2 at day 100; FDR < 0.05).

Although the largest mutation-specific changes were observed early, the shared differential expression signature of mutational effects increased with maturation. Examination of the correlation of log2FCs over time revealed that the correlations between hCO harbouring different genetic events were significantly stronger at day 100 (analysis of variance (ANOVA) F3,140 = 11, P = 1.6 × 10−6; Tukey’s honestly significant difference (HSD): day 100 to day 25 P = 2.7 × 10−5, day 100 to day 50 P = 8.0 × 10−6, day 100 to day 75 P = 7.0 × 10−4; Fig. 2f,g) suggesting that the changes in gene expression converge as cortical differentiation progressed in vitro, with distinct mutational forms showing more consistent differences from control participants and less mutational specificity. Meta-analysis across forms found the highest number of significant genes at day 100, more than 50% greater than at day 25 (day 25, 4,034; day 50, 2,606 genes; day 75, 3,957 genes and day 100, 6,680 genes; Fig. 2h). This increase in convergence of gene-expression changes could not be explained by alterations in variability in control participants or the distinct genetic forms across development (Extended Data Fig. 4, η2 = 0.04).

Dysregulated biological processes and cell types

We observed that six different forms of ASD profiled were enriched for upregulated genes related to neuronal cell fate and protein translation-related terms by day 75 (ref. 74) (Fig. 3a and Supplementary Table 5). We also observed broad downregulation of synapse and ion channel related terms such as ‘Synaptic membrane’ and ‘Gated channel activity’, across many time points and all the main mutational forms and in meta-analyses, particularly at day 75 and 100 (Fig. 3a, Extended Data Fig. 6 and Supplementary Table 5). The early disruption in pathways related to neurogenesis, neural progenitors and WNT signalling described above (for example, Fig. 2), coupled to the later emergence (day 75 and 100) of enrichment in synaptic pathways, provides a potential link between early mutational effects on neurogenesis that converge on signalling and synaptic biology during cortical neurogenesis across different genetic forms of ASD.

Fig. 3: Pathways and cell types affected across forms of ASD.

a, Select gene ontology terms enriched in upregulated (red) or downregulated (blue) genes in each ASD form. Colour shows normalized enrichment score (NES); point size shows −log10(FDR). Statistics were derived from GSEA. b, Cell proportion changes of ASD forms. Top annotation, differentiation day (shades of blue); lower annotation, ASD form (colourful). Statistics were derived from linear mixed-effects models and contrasts comparing each form with control participants with two-sided Wald tests. *FDR < 0.05. c, Adjusted R of WGCNA modules associated with either a form of ASD (Form) or with the presence of ASD (Dx) at day 25 and day 100 (left). Changes (linear model beta) in module eigengene in each ASD form (right). +FDR < 0.1, *FDR < 0.05, **FDR < 0.01, ***FDR < 0.005. d, Cell-type enrichment of modules using EWCE (red), genome-wide association study (GWAS) enrichment using a stratified linkage disequilibrium score (blue) and ASD-related gene set enrichment using a two-sided Fisher’s test (green). Only FDR < 0.05 are shown. e, Select gene ontology terms enriched in modules. Black lines represent FDR = 0.05. Statistics were derived from clusterProfiler enrichGO. f, The number of modules associated with either ASD forms (specific mutations) or with Dx. Significance was tested using two-sided Fisher’s exact test. g, Magnitude of module association with ASD forms or ASD broadly (Dx). Significance was tested using Welch’s two sample t-test (sample size: n = 22 modules per time point). Boxplots show: centre, median; lower hinge, 25% quantile; upper hinge, 75% quantile; lower whisker, smallest observation greater than or equal to lower hinge −1.5× interquartile range; upper whisker, largest observation less than or equal to upper hinge +1.5× interquartile range. ****P *< 0.005. ADHD, attention-deficit/hyperactivity disorder; BD, bipolar disorder; CycPro, cycling progenitor; Dx, diagnosis; ExNeu-1, excitatory neurons 1 (early); ExNeu-2, excitatory neurons 2 (late); InNeu, interneurons; IPC, intermediate progenitor cells; MDD, major depressive disorder; ME, module eigengene; ORG, outer radial glia; RG-1, radial glia 1 (early); RG-2, radial glia 2 (late); SCZ, schizophrenia.

To further explore the relationships of gene-expression changes to specific cell types, we deconvolved bulk RNA-seq into cell-type transcriptomes using single-cell RNA-seq (scRNA-seq) from hCOs44 (Extended Data Fig. 7). Predicted changes in cell-type proportions were subtle and largely seen at one or two time points per form (Fig. 3b and Extended Data Fig. 8). The least mature progenitors (CycPros) and more mature neuronal cell types (ExNeu-2) showed predicted reductions in 16p11.2del and 22q13.3del, suggesting changes in the timing of neuronal maturation (ExNeu-2 decreased at day 25 and CycPro decreased at day 75), consistent with recent observations in a 16p11.2 organoid model75. Intermediate progenitors showed changes in 15q13.3 deletion, 16p11.2 duplication and Timothy syndrome, suggesting convergent effects on this cell type during development (Fig. 3b and Extended Data Fig. 8). Correlation of deconvoluted cell-type compositional shifts increased at day 100 compared with day 25, demonstrating convergence in cell-type changes, similar to what we observed with differential expression and changes in biological process (Extended Data Fig. 8). By day 100, although the correlation across forms increased, only InNeu showed any significant shifts from control participants in any given form (22q13.3del, 15q13.3del and PCDH19). This indicates that observed transcriptional convergence is probably related to both altered timing of differentiation of cell types and altered gene expression signatures within cell types, the latter appearing more prominent here.

We used an orthogonal method, independent components analysis to determine whether pathway changes were robustly identified; we observed high concordance between independent components analysis and expression-based clustering (Rand index 0.58, Extended Data Fig. 9). Upregulated independent components across 16p11.2 deletion and duplication were downregulated for processes involved in neuronal precursor renewal and differentiation, such as ‘Beta catenin-TCF complex assembly’ and ‘Neuron fate specification’76 (Extended Data Fig. 9 and Supplementary Table 6). The broad role of these pathways in ASD has previously been suggested42,43,77, highlighting their potential role as a point of convergence. Independent component 31 (IC31), enriched for synaptic terms, was downregulated in 16p11.2 duplication (IC31 beta = −0.13, FDR = 1.11 × 10−3) and Timothy syndrome (IC31 beta = −0.14, FDR = 2.9 × 10−3; Extended Data Fig. 9 and Supplementary Table 6), similar to what was observed in the differential expression analysis.

Network analysis supports convergence over time

As transcripts do not act independently, but as part of highly regulated transcriptional networks78, we constructed co-expression networks using weighted gene co-expression network analysis (WGCNA)78,79,80. We tested whether modules derived from the day 25 samples changed in their association with individual forms and with all ASD-associated forms (diagnosis) compared with control participants over time (Fig. 3c, Extended Data Fig. 10 and Supplementary Table 7). Notably, we observed that the 22 identified co-expression modules from day 25 data were preserved at all other time points, suggesting they represent robust, continuing biological processes over this 100-day period. Moreover, their conservation in vivo80 and in vitro81,82,[83](#re