Rare coding variation illuminates the allelic architecture, risk genes, cellular expression patterns, and phenotypic context of autism

Patterns of coding variants in a large ASD exome dataset

We aggregated whole exome sequencing (WES) data across 33 ASD research cohorts, totalling 63,237 individuals: 15,036 affected offspring, 5,492 unaffected offspring, 28,522 parents, as well as 5,591 affected and 8,597 unaffected individuals from case-control studies (one sample exists in our dataset as both mother and affected offspring, Fig. 1a, Supplementary Tables 1, 2, 3, 4.1). Data from 35% of these individuals had not previously been analyzed or published. All accessible DNA sequence reads were aligned to the GRCh38 human reference genome and coding SNVs, indels, and CNVs were identified using GATK haplotype caller⁴² while a new method, GATK-gCNV⁴³, was used for CNV delineation. Variant counts were consistent across cohorts, with an average of 1.64 (1.66/affected, 1.57/unaffected) de novo SNVs, 0.18 (0.18/affected, 0.16/unaffected) de novo indels, and 0.035 (0.042/affected, 0.014/unaffected) de novo CNVs per offspring.

• Download figure
• Open in new tab

Fig. 1 Overview of SNV/indel and CNV rates in ASD by mode of inheritance.

a, The ASD cohort for SNV/indel analysis consisted of 49,049 family-based samples (15,036 cases) and 14,188 case/control samples (5,591 cases). b, The relative difference in PTV frequency between cases and unaffected controls (top) and average per-sample variant count in unaffected controls (bottom) across inheritance classes (color) and LOEUF deciles (5,446 genes in top 3 deciles of LOEUF). Using a binomial test, we found significant enrichment in cases for PTVs among the most constrained genes (lower LOEUF deciles), which weakened as negative selection against PTVs is relaxed (higher LOEUF deciles). c, Equivalent information and statistical significance for missense and synonymous variants. Synonymous variants were not enriched in cases or controls, as evaluated via binomial tests, in keeping with well-calibrated variant calling. d, The GATK-gCNV exome CNV detection pipeline when compared against WGS on overlapping samples achieved a sensitivity of 82% and positive predictive value of 89% for rare CNVs (<1% site frequency) that spanned 3 or more callable exons (red line). e, The relative difference in variant frequency between cases and controls for deletions. Using binomial tests, we found that deletions overlapping genes in the lowest LOEUF decile were observed at an elevated rate in affected individuals. f, Equivalent analysis for duplications, demonstrating a similar pattern of enrichment compared to deletions but with a weaker effect.

Abbreviations: PTV: protein truncating variant; CNV: copy number variant; WGS: whole genome sequencing; WES: whole exome sequencing; misB: missense variants with MPC score >=2; misA: missense variants with MPC score >=1 and <2; Del: Deletion CNV; Dup: Duplication CNV.

Statistical tests: b,c,e,f: binomial test.

Consistent with prior studies, we observed more PTVs and damaging missense variants in individuals with ASD compared to unaffected individuals (Fig. 1b-c). Enrichment of de novo and inherited PTVs in ASD cases was greatest in genes with the strongest evidence of selective constraint, represented by low LOEUF scores³⁹ (Supplementary Tables 4.2-4.7). Both de novo and inherited PTVs were enriched in the top three deciles of LOEUF (binomial test; 5,446 total genes; Fig. 1b). We grouped missense variants into three bins based on MPC score, which we termed MisB (MPC ≥ 2) and MisA (2 > MPC ≥ 1), with MPC < 1 for all remaining missense variants. MisB variants and, to a lesser degree, MisA variants were significantly enriched in ASD individuals (Fig. 1c). Our aggregated cohort demonstrated the greatest risk from de novo variation, with modest but significant ASD risk observed in rare case/control (for which de novo status cannot be determined) and inherited PTVs in the top three deciles of constrained genes, as well as in MisB and MisA variants.

Detection and analysis of rare and de novo CNVs in ASD

Previous studies using chromosomal microarrays have established a clear etiological role for large, rare and de novo CNVs in ASD cases compared to their unaffected siblings and the general population^{21,22,25–27,44–46}. Despite their considerable impact on ASD risk, CNV discovery at the resolution of individual exons from exome sequencing has represented a significant technical challenge. To overcome these challenges, we applied and benchmarked a pipeline built around GATK-gCNV, a novel Bayesian read depth-based CNV discovery tool for short-read sequencing data⁴². Our pipeline includes a PCA-based method for inferring large systematic differences, such as different exome capture kits (Supplementary Fig. 1), and we performed extensive benchmarking of this method from 7,165 samples from ASD quartet families for which matching gold-standard CNV calls from microarray and WGS were available^47,48. We measured sensitivity and positive predictive value (PPV) as a function of the number of captured exons from canonical protein-coding transcripts. Our application of GATK-gCNV achieved a sensitivity of 82% and a PPV of 89% for rare CNVs (site frequency < 1%) that spanned 3 or more such exons (Fig. 1d, Supplementary Fig. 2). When we considered de novo CNVs detected by WGS, GATK-gCNV achieved 82% sensitivity and 97.5% PPV for de novo CNV discovery at three-exon or greater resolution (Supplementary Fig. 3).

We applied GATK-gCNV to the 57,294 samples from the ASD cohort for which raw sequencing data were available (Supplementary Methods, Supplementary Table 3), including 13,694 affected and 5,007 unaffected offspring with data from both parents and 608 cases and 11,312 controls from non-family based cohorts. We focused our analysis on high quality CNVs spanning at least three captured protein-coding exons (Supplementary Methods), resulting in 16,889 rare (frequency < 1%) and 650 de novo autosomal CNVs, corresponding to an average of 0.90 (0.91/affected, 0.89/unaffected) rare and 0.035 (0.042/ affected, 0.014/unaffected) de novo CNVs per offspring.

Consistent with de novo coding CNVs mediating substantial risk for ASD, we observed that 3.93% of ASD cases and 1.40% of unaffected siblings harbored at least one de novo coding CNV, respectively (odds ratio [OR]: 2.88, p=9.6 × 10⁻²¹, Fisher’s exact). Splitting the CNVs by dosage (deletions versus duplications) and assigning each event to a LOEUF decile based on the most constrained gene within its boundaries, we observed that de novo deletions with at least one PTV intolerant gene showed the greatest degree of enrichment in ASD cases across all variant classes evaluated (OR: 9.10, p=2.81×10⁻²⁰, binomial test), with a relative difference approximately three-fold higher than that observed for de novo PTVs in the same constraint decile (p=2.9×10⁻⁴, simulation test). Enrichment of deletions was also observed in case/ control but not in rare inherited analyses. A similar but less significant pattern was observed for duplications across constraint deciles and modes of inheritance (Fig. 1e-f). Additionally, duplications with breakpoints residing within a gene have been suggested to harbor risk differently than duplications of an entire gene, likely through a loss of function mechanism^28,36. Leveraging our de novo duplication data, we found that partially duplicated genes in cases harbored more de novo PTV evidence compared to partially duplicated genes in controls (1.3 fold, p=1.7×10⁻³, binomial test).

To characterize the impact of exon and gene-level resolution in CNV discovery, we first considered 75 large genomic segments (Supplementary Table 4.8) associated with NDDs from curated literature^28,49–54. Of the 650 de novo CNVs discovered from WES here, 223 (34.3%) matched one of these loci with at least 50% reciprocal overlap, including 109 de novo GD deletions in cases compared to 6 in controls (OR: 6.7, p=3.9×10⁻⁹, Fisher’s exact) and 100 de novo GD duplications in cases compared to 8 in controls (OR: 4.6, p=5.4×10⁻⁷, Fisher’s exact, Fig. 2b). Considering rare inherited CNVs, we observed a non-significant trend towards enrichment of GD CNVs in cases (OR:1.3, p=0.07, Fisher’s exact). Excluding the GD loci, the remaining 427 de novo CNVs were also enriched in cases but with more modest effect sizes (de novo non-GD deletion: OR = 2.8, p=1.5×10⁻⁸, Fisher’s exact; de novo non-GD duplications: OR = 1.8, p=5.1×10⁻³, Fisher’s exact). A greater degree of enrichment was observed for non-GD deletions that overlapped a constrained gene (OR: 5.10, p=9.7×10⁻¹¹, Fisher’s exact, Fig. 2c) than for those that did not (OR: 1.32, p=0.32, Fisher’s exact, Fig. 2c) or for de novo PTVs in constrained genes (OR: 2.74, p=8.5×10⁻³⁴, Fisher’s exact, Fig. 2c). Effect sizes for de novo duplications were consistently smaller, although similar patterns were observed in relation to GD loci and constrained genes (Fig. 2c). Finally, when we considered de novo, non-GD CNVs that altered one of 71 genes associated with ASD risk (FDR ≤ 0.001, see below), we observed 46 variants in cases (29 deletions, 17 duplications) compared to none in siblings (p=8.5×10⁻⁷, Fisher’s exact). Considering the list of 373 genes implicated in both ASD and DD risk (see below) this increased to 130 de novo CNVs in cases and one in siblings (OR: 48.0, p=1.13×10⁻¹⁶, Fisher’s exact).

• Download figure
• Open in new tab

Fig. 2 Contribution of CNVs to ASD by mechanism and genomic location

a, CNVs involve deletions (DEL) or duplications (DUP) of genomic segments, and include a subset of sites known as genomic disorder (GD) loci that reoccur in human disorders. GDs that are mediated on both ends by non-allelic homologous recombination (NAHR) have recurrent breakpoints in the population, whereas non-NAHR GDs are located at sites susceptible to rearrangements. CNVs outside of known GD loci can still contribute to ASD risk as more unique events. b, Using Fisher’s exact tests, we find that de novo CNVs are more frequent in affected than unaffected offspring, and show a higher elevation than that observed in de novo PTVs (p=1.3×10⁻⁵, logistic regression). c, High ORs are observed for GD CNVs, a subset of which have no events observed in controls in this cohort (e.g., 16p11.2 deletions, 15q11.2-q13 duplications). Excluding GD CNVs, having a CNV that overlaps a constrained gene (LOEUF < 0.4) confers greater risk than having CNVs that overlap only unconstrained genes (p=3.3×10⁻⁴, logistic regression). d, Using Fisher’s exact tests, we find that de novo risk-mediating variants are observed at a higher frequency in affected females than affected males. e, Parent-of-origin analysis of de novo CNVs using binomial tests shows maternal bias for NAHR-mediated CNVs at GD regions, which is especially pronounced for 16p11.2.

Abbreviations: CNV: copy number variant; DEL: deletion CNV; DUP: duplication CNV; NAHR: non-allelic homologous recombination; PTV: protein truncating variant; GD: genomic disorder. Statistical tests: b,c,d: Fisher’s exact test, logistic regression; e: binomial test.

Consistent with a female protective effect in ASD, which proposes that a higher burden of risk factors is required for an ASD diagnosis in females^27,55, the burden of de novo CNVs overall was higher in female cases than male cases (5.94% vs. 3.47%, OR:1.76, p=3.4×10⁻⁸, Fisher’s exact). This effect was more pronounced for CNVs than PTVs (p=2.6×10⁻², logistic regression) or missense variants (p=2.0×10⁻⁵, logistic regression, Fig. 2d). While de novo SNVs and indels frequently arise on the paternal allele^9,44,56, a maternal allele bias has been observed from de novo CNVs in ASD⁵⁷. Using parental SNP data, we estimated the allelic parent-of-origin for 330 de novo CNVs (Supplementary Methods) and observed no overall bias (Fig. 2e, 130 maternal vs. 144 paternal, p=0.43, binomial test). However, restricting analyses to de novo CNVs at NAHR-mediated GD loci, we recapitulated prior findings with 70% arising on the maternal allele (70 maternal vs. 30 paternal, p = 7.9×10⁻⁵, binomial test). Including 20 additional samples from the Simons Searchlight project⁵⁸, we find that 94.3% of 16p11.2 CNVs occur on the maternal allele (Fig. 2e, 33 maternal vs. 2 paternal, p=3.7×10⁻⁸, binomial test). In contrast, CNVs at non-NAHR-mediated GD loci showed a 65.5% paternal bias (Fig. 2e, 60 maternal vs. 114 paternal, p=5.2×10⁻⁵, binomial test), consistent with a mechanistic bias in CNV formation and previous findings of a paternal origin across all non-recurrent/non-pathogenic CNVs⁴⁸.

Integrated discovery of ASD-associated genes across variant types and inheritance classes

Our analyses identified enrichment of rare protein-coding variants in ASD cases with effect sizes varying by mode of inheritance, variant class (PTV, misB, misA, deletion, duplication), and evolutionary constraint. We sought to utilise these insights and the size of our cohort (63,237 samples) to refine gene discovery in ASD by extending the TADA analytic framework^6,10,22,41 to include: (1) rare and de novo CNVs, (2) variants present in unaffected offspring, and (3) a continuous measure of evolutionary constraint from gnomAD (LOEUF³⁹, Supplementary Methods, Supplementary Fig. 4). A Bayes Factor (BF) was calculated to represent evidence of association for each autosomal protein-coding gene across variant types and modes of inheritance, taking into account null mutation rates and prior relative risks, when appropriate (Fig. 3).

• Download figure
• Open in new tab

Fig. 3 Integrating variant types and inheritance classes significantly boosts association power and reveals mutational biases

a, Our implementation of the TADA model includes de novo, case/control, and rare inherited modules for each variant type: PTV, MisB, MisA, deletion, and duplication. We leverage information from ASD probands as well as unaffected siblings in evaluating the effect of de novo variants. b, The evidence of ASD association contributed by each variant type for each of the 71 ASD genes with FDR ≤ 0.001. Some genes are predominantly associated through missense variants and duplications (e.g., PTEN, SLC6A1), suggesting mechanisms such as gain-of-function in contrast to haploinsufficient loss-of-function for genes with large contributions from PTVs and deletions. c, Applying TADA to our aggregated ASD dataset yields 71 genes at FDR ≤ 0.001, compared to 32 and 19 genes at the same threshold in previous studies on a subset of the samples (Satterstrom et al. 2020 and Sanders et al. 2015, respectively). Our expanded TADA model makes more efficient use of the available evidence of association by integrating more information into our statistical modeling. d, Across 71 ASD genes, the majority of evidence for ASD association is derived from PTVs and MisB variants and, e, de novo variants.

Abbreviations: BF: Bayes factor; PTV: protein truncating variant; misB: missense variants with MPC score >=2; misA: missense variants with MPC score >=1 and <2; Del: deletion CNV; Dup: duplication CNV; Inh: inherited; CC: case/control; DN: de novo. Statistical tests: b Extended TADA model.

Applying this model to the aggregated ASD data (TADA-ASD) we identify 71 genes associated with ASD at a FDR ≤ 0.001 and 183 genes associated at FDR ≤ 0.05 (Supplementary Table 4.9). Notably, the FDR ≤ 0.001 threshold is approximately equivalent to Bonferroni significance when back-calculating a p-value and correcting for 18,128 autosomal genes (p < 2.8×10⁻⁶), making it comparable to other recent studies of schizophrenia¹³ and DD¹¹. We focus on these 71 genes here, and illustrate in Fig. 3 that 93% of these genes have contributions from multiple variant classes (Fig. 3b). They also represent a meaningful increase in ASD discovery -a prior study from Satterstrom et al. included a subset of 11,986 of the cases aggregated here, producing 32 ASD risk genes at this same FDR ≤ 0.001 threshold using the TADA parameters applied in that study, which increases to 51 genes when incorporating the updated model parameters and the CNVs discovered here (Fig. 3c). The BFs within TADA allow us to assess genomic architecture through the relative contributions across variant types and modes of inheritance (Fig. 3d-e, Fig. 4a-b). In keeping with haploinsufficiency as the primary pathogenic mechanism, PTVs and deletions account for over 90% of the evidence in 20 of the 71 ASD-associated genes (28.2%). However, for 10 genes (14.1%), over 90% of the evidence comes from missense variants and duplications (e.g., DEAF1, SLC6A1; Fig. 4a). Included within these 10 genes is PLXNA1, for which the evidence comes primarily from inherited missense variation localized within the Plexin domain of the encoded protein (Fig. 4b-c).

• Download figure
• Open in new tab

Fig. 4 Relative contribution of evidence types in ASD risk genes

a, For 71 ASD-associated genes (FDR ≤ 0.001) the relative evidence of ASD association in the extended TADA model (log10BF) is shown for variants acting via a likely loss-of-function (LOF) mechanism (PTVs and deletions) on the x-axis versus variants that may act via alternative mechanisms (missense variants and duplications) on the y-axis. b, Equivalent plot of the relative evidence from de novo (x-axis) versus inherited (y-axis) variation for the same 71 ASD-associated genes. c, Evidence for ASD-association for the gene PLXNA1 is derived from de novo and transmitted missense variants, especially in the Plexin domain at the C-terminus of the Plexin-A1 protein.

Abbreviations: BF: Bayes factor. Statistical tests: c: Transmission disequilibrium test.

We next sought to estimate the magnitude of the ASD risk imparted by the genes discovered here and described in Fig. 3. Within the 71 genes, de novo PTVs or damaging missense variants were detected in 4.0% of cases and 0.5% of controls, giving a combined OR of 8.19. Using cross validation (Supplementary methods), we refined this estimate by variant type and mode of inheritance with OR estimates ranging from 16.0 for de novo PTVs to 1.02 for inherited MisA (Supplementary Table 4.10). We also applied this cross validation approach to estimate the yield of future gene discovery efforts, predicting a greater than linear increase in gene discovery at this same statistical threshold from a two-fold increase in the current sample size (150 genes; 95% CI: 108-278) or a three-fold increase (344 genes; 95% CI: 218-500).

Comparing the genetic architectures of ASD and general NDDs

ASD is frequently comorbid with other syndromic and non-syndromic NDDs, including intellectual disability (ID), and neurological disorders, including seizures^59,60. Prior analyses have clearly demonstrated both overlapping and distinct genetic architectures across ASD and DD, with significant heterogeneity across these broadly defined NDDs. To explore commonalities in genes that impact such NDDs, we integrated our ASD data with an independent cohort of 91,605 family-based samples encompassing 31,058 offspring with DD, the vast majority of which (>95%) were diagnosed with at least one NDD¹¹. Exome sequencing from this DD cohort was recently analyzed by the Deciphering Developmental Disorders (DDD) project using DeNovoWest, a permutation-based frequentist method, and reported exome-wide significant associations between rare SNVs and indels in 252 autosomal genes. However, these analyses did not include case-control samples, rare inherited SNVs and indels, or CNVs. We sought to re-analyze these DDD data to enable direct comparisons between the ASD and DD cohorts using uniform statistical models and significance thresholds between studies.

We applied the extended TADA framework to de novo variants in the DD cohort (TADA-DD) and found 309 autosomal genes associated at an FDR threshold of ≤ 0.001, including 94% of the initial 252 autosomal genes in Kaplanis et al. 2020¹¹ (Supplementary Table 4.9, TADA-DD found 477 at FDR ≤ 0.05). We observed the FDR values generated by our implementation of TADA (TADA-DD) to be highly correlated to those derived from the DeNovoWest significance values in Kaplanis et al. 2020¹¹ (r=0.95, p<10×10-22, Supplementary Fig. 5). Given the enrichment of cases with severe and syndromic disorders in this DD cohort compared to the ASD cohort^10,11,61, the de novo PTV, MisB, and MisA counts in this DD cohort showed the expected stronger but overall similar trends in variant enrichment in probands across the top three deciles of LOEUF (Fig. 5a-b).

• Download figure
• Open in new tab

Fig. 5 Integration of ASD and DD dataset identifies discovery of general neurodevelopmental disorder genes and heterogeneity

We analyzed 31,058 DD trios from Kaplanis et al. 2020 alongside our assembled ASD cohort to compile a combined dataset of 46,094 NDD cases. a, De novo PTVs are enriched in DD cases, with the effect diminishing as constraint relaxes. The enrichment in DD cases is higher than in the ASD cases, consistent with expectation. b, De novo missense variants are also enriched in DD cases above and beyond the enrichment observed in our ASD cases. c, 477 TADA-DD genes with FDR ≤ 0.05 have non-uniform p-value distributions from TADA-ASD that suggest 70.1% are associated with ASD. d, 183 TADA-ASD genes with FDR ≤ 0.05 have non-uniform p-value distributions from TADA-DD that suggest 86.3% of these genes are associated with DD. e, Using PTV and MisB variant data, we devised a chi-squared statistic, denoted the C statistic, to measure if a gene has more observed variants in one cohort relative to the other. A mixture model was used to deconvolve the commingled distributions. f, We transformed the fitted mixture distribution into posterior probability for ASD enrichment. Using a cutoff of <0.01, we found 82 DD-predominant genes, while using a cutoff of >0.99 we found 36 ASD-predominant genes.

Abbreviations: PTV: protein truncating variants; BF: Bayes factor; DD: developmental delay; FDR: false discovery rate. Statistical tests: a-b: binomial test, c-d: limma:: propTrueNull:, e-f: mixture model.

Since the first unbiased whole-genome analyses in ASD⁶², a notable overlap has been observed between genes affecting ASD and those affecting development more broadly, including intellectual development. Thus, it is reasonable to conjecture that integrating evidence for association from ASD and DD cohorts could yield additional genes influencing risk for ASD and broader development. Still, a cardinal rule of such meta-analysis is that the data should not be too heterogeneous. To determine whether the genes identified in the ASD cohort were also associated in the DD cohort, and vice versa, we converted the distribution of TADA q-values to p-values for each study. If the genes identified as associated in the DD cohort were also associated with ASD, the distribution of their ASD association p-values would be skewed toward zero; if they were not associated, the p-values would follow a uniform distribution; and, if some were associated and others were not, the p-values would follow a mixture distribution. Selecting the genes associated in the DD cohort at FDR ≤ 0.05, we evaluated the distribution of their p-values in the ASD cohort. The estimated proportion of associated genes (i.e., mixing parameter) was 0.701 (limma::propTrueNull, Fig. 5c), indicating that 70.1% of these genes affect risk for ASD. The converse analysis conditioning on the ASD-associated genes suggests that 86.3% of ASD risk genes have some broad effects on development (limma::propTrueNull, Fig. 5d).

Having observed that the results from the ASD and DD cohorts are somewhat complementary, we used the Bayesian framework underpinning TADA to integrate the genetic evidence for each gene across the cohorts by combining the BFs, which is conceptually similar to a frequentist meta-analysis. This combined analysis (TADA-NDD) revealed 373 genes associated with general NDDs at FDR ≤ 0.001 and 662 genes at FDR ≤ 0.05 (Supplementary Table 4.9). We observed that 54/373 and 124/662 genes did not meet the matching threshold in either the TADA-ASD or TADA-DD analyses alone, demonstrating a 15-19% increase in yield by combining evidence from these cohorts.

Heterogeneity of mutation patterns between ASD and DD associated genes

Due to the frequent comorbidity of ASD and DD phenotypes, isolating genes that exert a greater effect on ASD than they do on other DDs has remained challenging. Still, as documented above, 13.7% of the TADA-ASD FDR ≤ 0.05 genes show little evidence for association in the DD cohort at the same threshold (Fig. 5d). Moreover, while the remaining 86.3% of the genes are likely pleiotropic⁶³, some could have a far greater impact on risk for ASD while others have greater impact on risk for other features of development. To evaluate the overall gene-set heterogeneity between the ASD and DD cohorts, we retained only de novo SNVs/indels for independent gene-level BF calculations. For the 373 genes at TADA-NDD FDR ≤ 0.001, we observed a Pearson’s correlation of 0.78 (p=9.1×10⁻⁷⁸) of the gene-level log BF between the two major subcomponents of the ASD cohort (ASC+SSC versus SPARK), compared to only 0.42 (p=2.4×10⁻¹⁷) between the BFs of the ASD and DD cohorts. These observations signal much more consistent gene-level evidence between the ASD subcohorts than between the ASD and DD cohorts (Supplementary Fig. 6) and reflect both differential and shared genetic architecture underlying these often-comorbid phenotypes.

Next, we determined if certain genes were observed to be more commonly mutated in one cohort or the other. For these analyses, we first identified 464 “signal genes,” defined as any gene with FDR ≤ 0.05 in either TADA-ASD or TADA-DD from de novo PTV and MisB variants (which showed the greatest association and therefore were less likely to be benign). Of these signal genes, 120 belonged to TADA-ASD, 428 to TADA-DD, and 84 to both. Notably, even for the 84 genes significant in both cohorts, we also observed significant variant count heterogeneity (X2=317.6, DF=83, p=3.8×10⁻²³) between the cohorts.

Given this heterogeneity in mutational patterns, we sought to capture which of the 464 genes have more variation in either the ASD or the DD cohorts. A common way to assess this is by a standardized chi-squared test statistic (C statistic), but its power to discriminate is abrogated by the much higher burden of risk variants in the DD than ASD samples (Fig. 5a-b). We therefore adjusted for the difference in mutational burden between the cohorts by downsampling the DD mutations to be equal to the count of ASD mutations, then computed the C statistic. Genes with fewer than expected mutations in the ASD cohort relative to the DD cohort achieved a negative C statistic; a positive C statistic indicated overrepresentation in the ASD cohort. Next, a mixture model was adopted to disentangle the two commingled distributions, which assigns posterior probabilities that a gene is from the ASD or DD component of the statistical distribution (Fig. 5e-f, Supplementary Table 4.11). Using a posterior probability cutoff of greater than 0.99, we find 36 genes to be a part of the ASD mixture component and 82 genes to be a part of the DD component (Fig. 5f).

Mutational enrichment identifies differential neuronal layers impacted by ASD and DD genes

To gain insight into differences in expression of genes identified from the ASD versus DD cohorts, we examined single-cell gene expression patterns from human fetal brains. Two studies provided more than 37,000 cortical cells ranging from 6-27 weeks post-conception^64,65. To integrate these data, we first removed strong batch effects due to different protocols using cFIT and then obtained integrated factor loadings and gene expression for all measured cells⁶⁶. UMAP plots visualizing all cells from two sources, after removal of batch effects, showed that similar cell types from the different batches group together, while cells unique to either batch were also preserved⁶⁶ (Fig. 6a, Supplementary Fig. 7).

• Download figure
• Open in new tab

Fig. 6 Single-cell data reveals differential neuronal layers impacted by ASD and DD genes

a, A UMAP plot visualization after integrating two studies (Nowakowski et al. 2017; Polioudakis et al. 2019) that provided single-cell gene expression of human fetal brains consisting of 37,000 cortical cells 6-27 weeks post-conception. Similar cell types from the two batches group together while preserving cells unique to either study. b, Both ASD- and DD-predominant genes (right and left, respectively) were found to be enriched in interneurons and excitatory neurons compared to glial cells. Compared to DD-predominant genes, ASD-predominant genes are relatively more neuron-enriched than progenitor-enriched.

Abbreviations: UMAP: Uniform Manifold Approximation and Projection Statistical tests: b: Fisher’s exact

With the integrated data, we applied unsupervised clustering to identify cell subtypes in the context of a hierarchical tree to illustrate the relationships of major and minor cell type clusters. Using the MRtree method⁶⁷, we observed that cells of each labeled type were merged across both datasets into common clusters. Visualizing the tree, the major branches corresponded to glial and progenitor cells, excitatory neurons, deep layer enriched excitatory neurons, and inhibitory neurons. Likewise, minor splits reflected the expected relationships between cell types (Supplementary Fig. 8a). Based on the trajectory analysis of Polioudakis et al.⁶⁵, the ExN clade is less mature than the ExM clade, which in turn is less mature than the ExMU clade.

We subsequently assessed any differentials in strength of enrichment of the DD- and ASD-predominant risk genes within cell clusters for genes meeting the posterior probability 0.99 threshold (Fig. 5f). Of the 36 genes classified as ASD-predominant, 18 were expressed in these cell types; and, of the 82 classified as DD-predominant genes, 51 were expressed. Using odds ratio to reflect the strength of signal, we observed that both DD- and ASD-predominant genes were enriched in both interneurons and excitatory neurons compared to glial and progenitor cells (Fig. 6b, Supplementary Fig. 9, Supplementary Tables 4.12-4.13). While the enrichment was broadly distributed, the signal for DD-predominant enrichment tended to be found for cell types appearing early in the lineages of neurons. DD-predominant genes highlight four cell types, namely intermediate progenitors (IP), precursors of interneurons from the caudal ganglionic eminence (InCGE), migrating excitatory neurons (ExN3), and maturing excitatory neurons (ExM2). By contrast, ASD-predominant genes are strongly enriched in only one cell type, maturing excitatory neurons (ExMU1) and its clade. In keeping with a more mature neuron type, the ASD-predominant genes play roles in cytoskeleton organization, synaptic signaling, and activity-dependent programs. Our results regarding ASD agree with the conclusion of Polioudakis et al.⁶⁵, which was based on a subset of these data. If we judge enrichment solely by significance after Bonferroni correction for 21 cell types, these conclusions still hold: ExMU1 remained significant for enrichment of ASD-predominant genes; likewise ExN3 remained significant for enrichment of DD-predominant genes. Our results are consistent with DD-predominant genes being expressed earlier in development and in less mature cells than ASD-predominant genes.

Emergence of shared ASD schizophrenia risk genes

The genetic risk for ASD has also been suggested to overlap with that of other neuropsychiatric disorders, especially schizophrenia⁶⁸, such that the joint study of ASD and schizophrenia might reveal additional insights into both disorders. The Schizophrenia Exome Meta-Analysis (SCHEMA) Consortium recently analyzed exome sequencing from 24,248 schizophrenia cases and 97,322 controls, identifying 10 genes where ultra-rare coding variants were associated with schizophrenia after Bonferroni correction, and 244 genes at p < 0.01¹³. We compared our ASD- and NDD-associated genes to the SCHEMA findings to determine if there was any overlap between ASD and schizophrenia at the level of individual risk genes and, if so, whether it was related to ASD-DD overlap. Among the 71 ASD genes we discovered at an FDR ≤ 0.001, we found that 63 show an association with DD (using FDR ≤ 0.05, based on TADA-DD), and 8 show an association with schizophrenia (using p ≤ 0.01, based on values reported by Singh et al.). If the two associations were independent, we would expect ∼7 of the 8 ASD-schizophrenia genes to also show an association with DD (based on 63/71 = ∼88% of the ASD genes overlapping DD). However, we find 4 of the ASD-schizophrenia genes (NRXN1, ANK2, BRSK2, and DSCAM) lack an association with DD, which is a significant overrepresentation compared to random chance alone (p = 0.008, binomial test). This suggests that one subset of ASD risk genes may overlap DD, while a different subset overlaps schizophrenia.

We also analyzed these findings with an alternative approach: the distribution of genes in Fig. 5e (the 464 genes with FDR ≤ 0.05 for ASD and/or DD in our heterogeneity analysis) gives 36 ASD-enriched genes and 82 DD-enriched genes. Of the 244 genes identified by SCHEMA as schizophrenia-associated at p ≤ 0.01, 6 genes (ANK2, ASH1L, BRSK2, CGREF1, DSCAM, and NRXN1) overlap the 36 ASD genes, while only 3 (ATP2B1, GRIN2A, and HIST1H1E) overlap the 82 DD genes. If we compare to the null hypothesis that each of the genes in our TADA model has an equal chance of being schizophrenia-associated, then the ASD-schizophrenia overlap is significantly enriched (p = 8.5×10⁻⁶, binomial test), while the DD-schizophrenia overlap is not (p = 0.10, binomial test). The two outcomes (6/36 vs. 3/82) are also significantly different when compared to each other (p = 0.023, Fisher’s exact test), reinforcing the idea that the shared genetic risk between ASD and schizophrenia may be distinct to that observed between ASD and DD.