Transcriptional patterns of coexpression across autism risk genes converge on established and novel signatures of neurodevelopment

CRISPR perturbation functional models

We collected data from 17 newly generated and previously published CRISPR experiments targeting 15 different ASD-associated genes. The transcriptomic data from all these CRISPR experiments were generated from isogenic iPSC-derived glutamatergic neurons induced by Neurogenin 2 (Ngn2) overexpression. Data of previously published cellular models generated by CRISPR perturbation were ascertained from the NCBI’s Gene Expression Omnibus (GEO) and Synapse. For the CRISPR-mediated LoF gene editing models described in Deneault et al., gene counts and covariate information were ascertained from GEO (Accession # GSE107878). These consisted of 10 ASD-relevant genes (ANOS1, ASTN2, ATRX, CACNA1C, CHD8, DLGAP2, AFF2/FMR2, KCNQ2, SCN2A, TENM1). Gene counts and covariate information for MBD5 CRISPR-mediated LoF gene editing model were ascertained from GEO (Accession #: GSE144279). For the CRISPR activation (CRISPRa) targets (TSNARE1, SNAP91, CLCN3, and FURIN), gene counts and covariate information were ascertained from Synapse (syn20502314)¹⁸. CRISPR interference (CRISPRi) for SCN2A was ascertained from Synapse (syn26970716). From those datasets, only neuronal lines were included, and experiments were removed if there was any evidence of ineffective perturbation such as discordant RNA and protein levels of the perturbed gene.

Additionally, we used unpublished transcriptomic data from SCN2A and KCTD13 CRISPR-mediated LoF gene editing models which leveraged a dual-gRNA strategy to promote gene deletions. SCN2A sgRNAs targeted intron 4 and intron 11 (NM_001040142) to generate a 7.5kb partial gene deletion, and the crRNA sequences used were: 5’-tatcgtagggggaccaacc-3’ and 5’-gcgtggatctagtgaactt-3’. KCTD13 sgRNAs targeted 5’ and 3’ regions from KCTD13 (NM_178863) to generate a 23.3kb full gene deletion, and the crRNA sequences used were: 5’-taaaaaggatggatgtaggc-3’ and 5’-tgcctgtgttaggaggtatc-3’.

The deletion lines were generated in the male control human iPSC line GM08330-8 using the Human Stem Cell Nucleofector Kit 1 (Lonza), transfecting 1μg CRISPR/Cas9 PX459 plasmid and 1μg of each gRNA using the Amaxa Nucleofection II device (Lonza), according to the manufacturer’s instructions. Cells were subsequently plated on Matrigel plates in mTeSR1 or Essential 8 medium supplemented with ROCK inhibitor for 24 hours. For clonal isolation of SCN2A models, puromycin selection was started 24 hours after transfection and resistant colonies were picked and expanded 48 hours after selection. For clonal isolation of KCTD13 models, cells were separated by fluorescence-activated cell sorting (FACS) 48h after transfection. Genotyping of the resultant colonies for SCN2A and KCTD13 was performed by Sanger sequencing of the deletion-specific region and ddPCR assays for copy number. A total of 4-6 successfully edited clones with heterozygous deletions plus 2-6 unedited (i.e. WT-Cas9 exposed) clonal colonies were expanded per target. Prior to neuronal differentiation, iPSC clones were split into multiple replicates and each was manipulated in parallel during subsequent experiments.

For differentiation of SCN2A and KCTD13 human iPSC models into glutamatergic neurons, Ngn2-neuronal induction was performed as previously described²⁵. Briefly, iPSCs were seeded at a density of 1e6 cells/mL and transduced with a lentivirus expressing TetO-Ngn2-GFP-Puro or TetO-Ngn2-Puro along with rtTA. Twenty-four hours after transduction, doxycycline was added to initiate Ngn2 expression, and then 24h later the cells were selected with puromycin. Ngn2-glutamatergic iPSC-derived neurons were cultured in neuronal maintenance medium supplemented with BDNF and GDNF growth factors for an additional 22 days. Subsequent experiments were performed with 24 day-old Ngn2-glutamatergic neurons, using 6-34 total replicates per genotype (i.e. WT, heterozygous deletion) per target gene.

RNA sequencing of CRISPR perturbations

SCN2A RNAseq libraries were prepared from 200 ng of total RNA using a TruSeq stranded mRNA Sample Prep Kit (Illumina cat# RS-122-2102). Libraries were multiplexed, pooled and sequenced on multiple lanes of the Illumina NovaSeq platform, generating an average of 30.7M paired-end 150 bp-cycle reads for 30 samples (20 SCN2A^+/-, 10 SCN2A^+/+). RNAseq data was processed using a standard workflow, which includes quality assessment of fastq reads using FastQc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Raw sequence reads were trimmed against Illumina adapters using Trimmomatic²⁶ (v. 0.36) with parameters ILLUMINACLIP:adapter.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75. Sequence reads were aligned to the human reference genome (GRCh37, Ensembl build 75) using STAR²⁷ (v. 2.5.2.a) with parameters ‘–outSAMunmapped Within –outFilterMultimapNmax 1 – outFilterMismatchNoverLmax 0.1 –alignIntronMin 21 – alignIntronMax 0 –alignEndsType Local –quantMode GeneCounts –twopassMode Basic’. STAR aligner also generated gene level counts for all libraries relying on the human genome annotation provided for Ensembl GRCh37, build 75. Quality checking of alignments was assessed by custom scripts utilizing PicardTools (https://broadinstitute.github.io/picard/), RNASeQC²⁸, RSeQC²⁹ and SamTools³⁰. The deletion introduced by CRISPR on SCN2A loci was validated generating exon-level coverage using DEXseq³¹ and visually investigating target loci using IGV³². KCTD13 was included in gene edited as a potential contributor to 16p11.2 genomic disorder³³, and these lines are described elsewhere (Tai et al., in preparation). In brief, RNAseq libraries were prepared using a TruSeq stranded mRNA library kit (Illumina) and were multiplexed, pooled and sequenced on multiple lanes of Illumina HiSeq 2500 platform, generating an average of 46.5M paired-end reads of 75bp for 11 samples (5 KCTD13^+/-, 6 KCTD13^+/+). The same RNAseq data processing pipeline without trimming step as above was applied to KCTD13 RNAseq libraries. Processing of MBD5 RNAseq libraries were previously described elsewhere³⁴.

For the CRISPR-mediated LoF gene editing models described in Deneault et al., gene counts and covariate information were ascertained from NCBI’s Gene Expression Omnibus (GEO accession: GSE107878). The processing of samples followed standard RNA sequencing pipelines as previous described in Deneault et al. (2018)¹⁰.

Differential expression analysis for CRISPR perturbations

Gene counts were input into DESeq2 for differential expression analysis. Genes with normalized counts < 10 in at least half of samples were excluded. Covariates including RNA integrity number (RIN), batch, and processing date were included when available. A Wald’s test was used for differential expression and Z-scores were calculated by dividing the fold change by the standard error. Each CRISPR perturbation was analyzed separately.

Post-mortem brain cohorts

Two separate cohorts of post-mortem brain samples of the dorsolateral prefrontal cortex (DLPFC) were used. The CommonMind Consortium (CMC) included tissue samples from Mount Sinai NIH Brain Bank and Tissue Repository, The University of Pittsburgh NIH NeuroBioBank Brain and Tissue Repository, and the University of Pennsylvania Brain Bank of Psychiatric Illnesses and Alzheimer’s Disease Core Center. The DLPFC was dissected at each bank and sent to a centralized centre, the Icahn School of Medicine at Mount Sinai (ISMMS) for RNA extraction. Tissues from bipolar disorder (BD) or schizophrenia cases were included if they met the DSM-IV diagnostic criteria for schizoaffective disorder or schizophrenia, or for BD, which were determined in consensus conferences after reviewing of medical records, direct clinical assessments, and care provider interviews. Samples were excluded if donors had a history of Alzheimer’s disease, Parkinson’s disease, were on ventilators near time of death, or had acute neurological insults (anoxia, strokes and/or traumatic brain injury) before death.

The Human Brain Core Collection (HBCC) cohort consisted of DLPFC samples from the NIMH HBCC. The samples were analyzed clinically, neuropathologically, and toxicologically. The DSM-IV clinical diagnosis was determined through review of medical records by two psychiatrists and family interviews. Non-psychiatric controls did not have a history of substance use disorder or psychiatric conditions. Across all 933 samples, there were 345 females and 588 males. The self-reported ethnicities were 637 Europeans, 249 African, 33 Hispanic, 13 Asian and 1 other. Nearly fifty percent of samples were nonpsychiatric controls (N=462) and the remainder had a psychiatric diagnosis (113 BD, 350 schizophrenia, and 8 affective disorder)³⁵. The research abided by ethical regulations and was approved by the Vanderbilt University Medical Center Review Board (IRB: 220287).

RNA sequencing of post-mortem samples

Approximately 50 mg of homogenized tissue from the DLPFC was used to isolate RNA. The two cohorts were processed separately. Samples with age <18 were excluded prior to analysis. The RNA-seq processing is identical to that described in Han et al. (2020)³⁶ except for not using surrogate variable analysis (SVA) here to avoid removing trans-regulatory effects. Briefly, STAR was used to align RNA sequencing reads to GRCh37. FeatureCounts (v1.5.2) was used to count uniquely mapped reads that overlapped genes using the Ensembl v75 annotations. Fixed/mixed effects modelling were used for library normalization and covariate adjustments. Genes that were expressed at levels > 1 counts per million (CPM) in at least half of the samples in each study were retained for analysis. Conditional quantile normalization was done to account for variation in GC-content and gene length. A weighted-linear model using voom-limma was used to assess the sampling abundance confidence. Normalized log₂(CPM) values were used for hierarchical clustering and principal component analysis to detect outlier samples. Samples were removed if deemed outliers using either method. For the CMC cohort, covariates were identified using a stepwise fixed/mixed effect regression model to identify covariates significantly associated with gene expression. The covariates were added if there was an association with principal components explaining greater than 1% of expression residual variance. For the HBCC cohort, model selection was determined through Bayesian information criteria (BIC). BIC were used to find fixed effect covariates that improved the model for most genes. Covariate adjustment was done using with a fixed/mixed effect linear regression variant, choosing mixed effect models when several samples were available per donor: gene expression ∼ covariates + sex + diagnosis + (1|Donor). Observation weights were calculated using voom-limma to adjust for the mean-variance relationship. The covariate-adjusted expression was generated after adding back the diagnostic component. This was done for both HBCC and CMC cohorts.

Generating pairwise coexpression of genes in the DLPFC

For each of the two cohorts, the covariate-adjusted expression was used to calculate coexpression values across all 16,992 genes. A pairwise Pearson’s correlation was calculated for each pair of genes and the correlation coefficient was subsequently transformed into a Z-score using a Fisher transformation. The coexpression Z-scores for each pair of genes were subsequently meta-analyzed across the two cohorts using Stouffer’s weighted Z-score method.

Assessing relationship between CRISPR perturbation and coexpression

For each experiment, we calculated Pearson correlation between the rank normalized differential expression and the perturbed gene’s coexpression profiles from the post-mortem brain tissue. ASD CRISPR perturbations were meta-analyzed using Stouffer’s weighted Z-score and compared to the meta-analyzed coexpression of the same genes. The relationship was subsequently assessed with a Pearson’s correlation.

ASD convergent coexpression meta-analysis

We included 71 genes implicated in risk of ASD from a cross-consortia exome sequencing study that combined de novo and inherited single-nucleotide variant (SNV), indel, and CNV analyses (FDR < 0.001, approximating exome-wide Bonferroni correction)⁵. After filtering out genes due to low expression level, we performed a meta-analysis of 70 (out of 71) genes. Coexpression Z-scores were meta-analyzed using Stouffer’s weighted Z-score method to generate convergent effect sizes. To assess whether convergence of ASD risk genes is tissue specific, we used Genotype Tissue Consortium (GTEx) RNA-seq counts of skeletal muscle (N=803) and liver (N=251) as negative controls. We selected these as they are not derived from the ectoderm or previously implicated in ASD. Convergence was calculated using the same methods as described previously. To assess the null distribution for ASD convergence, we conducted 100,000,000 permutations where for each permutation a meta-analysis was performed using 70 randomly-selected genes, excluding the 70 ASD genes. The empirical p-values were calculated, as shown here: (# of absolute convergent Z-scores greater than or equal to the absolute convergent ASD Z-score + 1) / (Total # of permutations + 1). A Bonferroni-correction was then applied to the empirical p-values and a threshold of p < 0.01 after correction was used to increase stringency.

To account for the potential effects of confounding from LoF observed/expected upper bound fraction (LOEUF), the permutation was repeated by matching on LOEUF scores within −/+ 0.05 of each gene. Genes that were Bonferroni-significant and large effect (absolute Z-score > 2) were input into ToppGene for pathway enrichment, with the background gene set being all unique genes with coexpression values. Afterwards, we sought to determine whether the convergent coexpression relates to differentially expressed ASD genes. Summary statistics from the ASD vs control post-mortem brain PsychENCODE dataset were ascertained and correlated against the convergent coexpression.

Dissecting relationship between convergence, intolerance and ASD association

Since ASD risk genes are strongly intolerant and intolerant genes are more coexpressed with each other, we wanted to assess whether among the significant convergent genes were also tolerant genes associated with ASD. To define intolerance, we divided the genes into two sets, tolerant (LOEUF >0.35) and intolerant (LOEUF <0.35). The correlation was assessed between transcriptional convergence and the −log₁₀(p) significance of the exome data. Next, we sought to identify whether transcriptional convergence can identify novel ASD genes that have not been previously implicated due to limitations of genetic studies. We defined ASD association as having a total Bayes factor (BF) > 2 for the exome data which is includes the presence of even a weak association. A BF of greater than 3 is typically used to represent meaningful significance. We assessed the correlation between LOEUF and transcriptional convergence for both associated and non-associated ASD genes.

Classification of 71 ASD associated genes

We categorized 71 ASD associated genes into three main functional groups including chromatin, transcription and synaptic by their association with chromatin function related GO terms (“GOMF_CHROMATIN_BINDING”, “GOMF_CHROMATIN_DNA_BINDING”, “GOBP_REGULATION_OF_CHROMATIN_ASSEMBLY_OR_DISASSEMBLY”,”GOBP_CHROMATIN_ ORGANIZATION”,”GOBP_CHROMATIN_REMODELING”,”GOBP_REGULATION_OF_CHROMATIN_ BINDING”,”GOBP_CHROMATIN_MEDIATED_MAINTENANCE_OF_TRANSCRIPTION”,”GOBP_CHR OMATIN_MAINTENANCE”,”GOBP_REGULATION_OF_CHROMATIN_ORGANIZATION”) and transcription related GO terms (“GOBP_MRNA_TRANSCRIPTION”, “GOBP_REGULATION_OF_TRANSCRIPTION_REGULATORY_REGION_DNA_BINDING”, “GOBP_MRNA_TRANSCRIPTION_BY_RNA_POLYMERASE_II”,”GOBP_CHROMATIN_ORGANIZATI ON_INVOLVED_IN_REGULATION_OF_TRANSCRIPTION”,”GOMF_RNA_POLYMERASE_II_TRANS CRIPTION_FACTOR_BINDING”,”GOMF_TRANSCRIPTION_COREGULATOR_BINDING”,”GOMF_DN A_BINDING_TRANSCRIPTION_FACTOR_ACTIVITY”,”GOMF_TRANSCRIPTION_FACTOR_BINDIN G”,”GOMF_TRANSCRIPTION_REGULATOR_ACTIVITY”) and synaptic function based on manual curation reported in SynGO v1.1 database³⁷. The above listed GO terms and their associated genes were retrieved from MSigDB (v7.4) database. Genes that were not identified in any of these categories were manually further classified based on https://www.genecards.org/.

Characterization of CRISPR-mediated perturbations

RNA-sequencing data from 17 CRISPR perturbation experiments in neurons comprising 15 unique genes were included from three independent sources (see Methods). In total, the CRISPR experiments included: twelve CRISPR-mediated heterozygous or homozygous LoF mutational models (AFF2, ANOS1, ASTN2, ATRX, CACNA1C, CHD8, KCTD13, KCNQ2, MBD5, SCN2A [2x], and TENM1)^10,34, four CRISPR activation models (FURIN, SNAP91, TSNARE1, and CLCN3)³⁸, and one CRISPR interference model (SCN2A). We tested for differential expression between the CRISPR edited cell lines and the unedited cell lines for each gene passing quality control including observing the expected effect in the perturbed gene (see Methods). Each experiment was analyzed separately, and the differential expression effect sizes were converted to Z-scores representing “CRISPR perturbation.” The experiments had variable impact on global expression patterns, with the total number of significantly differentially expressed genes identified ranging from 164 to 4,857 and lambda inflation factors ranging from 0.29 – 4.13 (Table 1).

One gene (SCN2A) was perturbed in three independent experiments allowing us to quantify the variability of global gene expression changes due to CRISPR perturbation and genetic background. We found significant but modest Pearson correlations of differential expression across experiments that varied from 0.202 between the two CRISPR-mediated gene editing experiments to 0.343 and 0.335 between the CRISPRi experiment and the two gene edit experiments (Supplementary Figure 1).

Gene coexpression correlates with downstream transcriptional consequences of CRISPR perturbation

We hypothesized that gene coexpression would proxy differential expression driven by CRISPR perturbations in the experiments described above. We leveraged 993 post-mortem brain samples from the dorsolateral prefrontal cortex (DLPFC) and calculated pairwise gene coexpression using Pearson correlation (see Methods). We identified consistent and significant negative correlation between differential expression due to CRISPR-mediated gene knockdown and the corresponding perturbed gene’s normalized coexpression profile ranging from −0.089 to −0.45 (Figure 1A). The exception was the ANOS1 LoF model which had a nominally significant positive correlation (R=0.018, P=0.046). Overall, genes that are more significantly downregulated by a specific CRISPR perturbation are more highly and positively coexpressed with the perturbed gene. This holds true when CRISPR targets are downregulated via LoF mutation or CRISPRi as well as when they are upregulated using CRISPRa (Supplementary Figure 2).

• Download figure
• Open in new tab

Figure 1. Coexpression proxies CRISPR-induced differential expression.

(A) Correlation between differential expression profiles for ten CRISPR knockdown models (y-axis) with respective coexpression profiles (x-axis) for each gene (points). (B) Meta-analysis across 10 CRISPR perturbations with respective meta-analyzed coexpression. Coexpression is represented as a Fisher transformed Pearson’s correlation Z-score. A Pearson correlation was done to assess the correlation between coexpression and differential expression. The curve for each panel was fitted with a locally weighted smoothing (LOESS) regression.

With a variably sized but consistently significant negative relationship seen across single gene CRISPR perturbation and coexpression, we next asked whether meta-analyzing differential expression from multiple CRISPR perturbations could also be proxied by a meta-analysis of the coexpression profiles of the perturbed genes. In other words, we wanted to assess whether genes that were consistently differentially expressed in the same direction across multiple gene perturbations (convergent genes) could be inferred from a similar convergence of coexpression of the respective genes from post-mortem brain tissue. When we separately meta-analyzed all CRISPR-mediated gene edits and corresponding coexpression profiles for the perturbed genes, we found that convergent CRISPR perturbation was significantly negatively correlated with convergent coexpression (Pearson’s R = −0.44, P < 1×10⁻³⁰⁰, Figure 1B).

Given the significant correlation between CRISPR perturbation and coexpression, we next assessed how that compared to the correlation of differential expression from different CRISPR perturbations of the same gene (SCN2A) across three different experiments. Significant Pearson correlation was observed between SCN2A coexpression and differential expression in each CRISPR experiment of SCN2A ranging from −0.25 – −0.43 (Figure 2). The correlation statistics across CRISPR SCN2A experiments (0.2 – 0.34) were similar to those seen when compared to coexpression suggesting that coexpression provided a similar proxy to transcriptional dysregulation as an independent CRISPR experiment. After meta-analyzing the differential expression profiles from the three SCN2A CRISPR experiments, the correlation with coexpression was stronger than each individual experiment alone (Pearson’s R = −0.45, P < 1×10⁻³⁰⁰).

• Download figure
• Open in new tab

Figure 2. Coexpression consistently proxies differential expression for SCN2A across three independent perturbations.

Coexpression represented as a Fisher transformed Pearson’s correlation Z-score (x-axis). A Pearson correlation was calculated to assess the correlation between coexpression and differential expression (y-axis) for each gene (points). The curve for each panel was fitted with a locally weighted smoothing (LOESS) regression.

In silico functional convergence of ASD risk genes

We next sought to leverage coexpression to test convergence of risk genes implicated from the most recent ASD exome sequencing study⁵. Among the 71 implicated genes, 70 had coexpression data and could be meta-analyzed to calculate the convergent coexpression. After calculating meta-analysis effect sizes of Fisher’s transformed correlations, we assessed significance by running 100,000,000 permutations where the coexpression of 70 randomly selected genes (among those with coexpression) were meta-analyzed (see Methods). Significant overall convergence was identified by comparing the variance of the distribution of convergent coexpression of the ASD genes to the variance of the permuted random genes (ASD variance = 1.70, mean permuted null variance = 0.131, p = 9.9×10⁻⁷). A total of 3,798 genes were found to be significantly convergently coexpressed at a Bonferroni-adjusted empirical p-value less than 0.01. Nearly 60% of those (n=2,271) were positively convergent suggesting that their expression was consistently in the same direction as the ASD genes.

The analysis that identified the 70 ASD genes directly incorporated constraint against LoF variation (LOEUF³⁹) into estimates of prior relative risk, which could inflate coexpression convergence among intolerant genes. To account for this potential confounding, we reran the permutation procedure randomly selecting 70 genes with matching LOEUF scores to the ASD genes. The p-values between intolerance-matched permutation and not were strongly correlated (Spearman’s Rank Correlation, rho = 0.93, P < 1×10⁻³⁰⁰), suggesting that the models of ASD risk gene discovery did not significantly inflate discovery of convergent genes (Supplementary Figure 3). Further, an analysis based on a Fisher’s combined test using only de novo and missense variants without incorporating intolerance identified 49 ASD genes (FDR < 0.001) that showed no significant difference in intolerance scores compared to the 70 ASD genes (Two-sided Wilcoxon Test, W = 1825.5, p = 0.648). Given that intolerance did not seem to confound our analyses, we meta-analyzed the initial 70 genes for downstream analyses.

Convergent genes were overrepresented in genesets related to neuronal and synaptic function with the most significant pathways across categories including the synapse (p = 6.00×10⁻⁴²), synaptic signaling (p = 2.77×10⁻²⁶), neuronal system (p = 7.66×10⁻¹⁶) and abnormal CNS synaptic transmission (p = 2.06×10⁻¹³) (Supplementary Table 1). Convergent genes also represented 126 of the 485 genes implicated in neurodevelopmental disorders from DisGeNet⁴⁰ (p = 1.72×10⁻⁹) and 72 of the 239 genes in the cation channel complex (p = 5.23×10⁻⁷). Next, we partitioned the significantly convergent genes by direction of effect and found that positively convergent genes were largely enriched in synaptic pathways (Supplementary Table 2) whereas negatively convergent genes had much less clear pattern of enrichment across pathways (Supplementary Tables 3).

Convergence captures ASD signal from exome and post-mortem brain studies

Given the strong enrichment of synaptic functions observed across convergent genes and the neurodevelopmental deficits observed in ASD, we reasoned that convergent genes that contribute to ASD etiology would display disease-relevant tissue specificity. We tested this hypothesis and found no meaningful convergence using either liver or muscle coexpression data from GTEx (Figure 3A). The mean absolute convergence Z-score for the liver and muscle were 0.254 and 0.234 respectively, while the mean absolute Z-score from the DLPFC was 1.069 and was significantly different than both the liver and the muscle (Two-sided Wilcoxon Test, p < 2.2×10⁻¹⁶), suggesting tissue specificity of ASD convergence. To determine if the enrichment of convergent genes associated with synaptic terms was driven by ASD risk genes associated with synaptic functions, we assigned 60 of the 70 ASD genes to relevant pathway terms associated with synaptic, transcription factor, and chromatin terms. The convergence permutations were rerun for each gene set and we found that convergence was strongly driven by synaptic genes, whereas the chromatin genes and transcription factors displayed progressively weaker convergence (Synaptic variance: 2.57, chromatin variance = 0.39, transcription variance = 0.18) (Supplementary Figure 4).

• Download figure
• Open in new tab

Figure 3. Convergence is tissue specific and associated with ASD.

(A) Distributions of transc iptional convergence of 70 high-confidence ASD genes in different tissues. (B) Smoothed relationship with co fidence interval between ASD exome significance and significantly convergent coexpression effect sizes (P_bonf < 0.01). (C) Correlation between convergently coexpressed genes (P_bonf < 0.01) and differential expression of ASD post-mortem DLPFC compared to controls. The curve was fitted with a locally weighted smoothing (LOESS) regression.

Finally, we asked whether the convergent genes could inform on the biology of ASD. Excluding the 70 ASD genes used to calculate convergence, there was a significant correlation between convergence and evidence for ASD risk based on the significance of each of the remaining genes to ASD (Spearman’s Rank Correlation, rho = −0.14, p = 4.75×10⁻⁶³) (Figure 3B). There was a significant enrichment of ASD genes implicated through rare variant analyses (q-value <0.05) among the positively convergent genes (Z > 2) (OR: 4.65 [2.99 - 7.10], p = 5.74×10⁻¹¹). This effect size was larger when only including Bonferroni-significant (p < 0.01) positively convergent genes (OR: 9.56 [3.85 - 28.41], p = 6.56×10⁻⁹). That is, genes most positively convergent were also most likely to be associated to ASD through rare variant analyses.

ASD risk genes are biologically more likely to be intolerant. Convergent genes also are significantly more likely to be intolerant (mean LOEUF for convergent genes = 0.74, mean LOEUF for other genes 0.99, Wilcoxon p-value = 8 × 10⁻¹²⁷). However, after splitting genes based on intolerance, we found a significant correlation between convergence and ASD association for both tolerant (LOEUF > 0.35, Spearman’s Rank Correlation, rho=0.17, P=5.51×10⁻¹⁶) and intolerant genes (LOEUF < 0.35, Spearman correlation, rho=0.15, P=2.07×10⁻⁵) (Supplementary Figure 5). A significant relationship also existed between convergence and differential expression of ASD post-mortem brain tissue compared to controls with the positive convergence being correlated with downregulation in ASD (Spearman’s R = −0.22, p = 3.62×10⁻⁴¹, Figure 3C).

Potential to identify novel ASD genes missing from current genetic analyses

Genes with shorter coding sequences will have reduced power in genetic analyses that require multiple deleterious variants among cases to quantify risk. Similar issues could exist for genes where deleterious variation is inviable but less deleterious modulation could contribute to ASD risk. Given the significant relationship between convergence and ASD association, we asked whether our set of convergent genes included potential ASD risk genes with properties that current genetic studies might be underpowered to identify. Among the most significantly associated 71 ASD risk genes, there was a highly significant skew towards longer coding sequence (median = 3,642bp) compared to 1,293bp among all other genes. After splitting our positive convergent genes into those with even a weak association to ASD (BF > 2) or not, we identified significantly increased median coding sequence lengths in the ASD associated genes compared to the rest (median ASD-associated coding sequence length = 2,702 bp; median non-ASD -associated coding sequence length = 1,809 bp, p = 1.94×10⁻²⁶ Wilcoxon Signed Rank Test, Figure 4A). The most positively convergent genes not associated to ASD had similar intolerance scores to those that were linked to ASD (Figure 4B). We show that positive convergence is more significantly correlated with the ASD effect size driven by missense than protein-truncating variation (Spearman’s Rank Correlation, missense type A, missense variants with missense badness, Polyphen-2, and constraint [MPC]⁴¹ scores between 1 and 2 : rho = 0.177, p = 2.57×10⁻¹³; missense type B, missense variants with MPC scores greater than or equal to 2: rho = 0.113, p = 3.66×10⁻⁶; protein-truncating: rho = 0.109, p = 6.99×10⁻⁶).

• Download figure
• Open in new tab

Figure 4. Convergence captures novel intolerant and small genes implicated in ASD risk.

(A) Distributions and medians (colored vertical lines) of coding sequence lengths for significant positively convergent genes stratified by their association with ASD. High confidence ASD exome genes (purple) consist of the 71 genes with a q-value < 0.001. A Bayes factor (BF) > 2 in the ASD exome data is considered ASD-associated. The positive convergent genes are split by those with BF > 2 (green) and those with BF < 2 (yellow). (B) Relationship between convergent coexpression and intolerance (LOEUF scores) stratified by association with ASD amongst significantly convergent genes. Significant convergence defined as P_bonf < 0.01.

Finally, we sought to assess whether convergence may increase power to detect genes that have some level of evidence linking them to ASD risk but have not been implicated in exome-consortia analyses. To address this question, we leveraged the manually curated database of genes implicated in ASD susceptibility by the Simons Foundation Autism Research Initiative (SFARI, https://gene.sfari.org/). Genes in this database come from large-scale sequencing studies but also from functional studies, clinical reports and genetic studies of common variation. After excluding all genes associated with ASD using the same criteria as before (BF > 2), there was a significant enrichment of the remaining positively convergent genes among this broader ASD geneset (Fisher’s Exact Test, OR: 3.31 [2.48 – 4.35], p = 3.65×10⁻¹⁴) (Supplementary Table 4). This finding suggests that convergence may provide a useful supplement to sequencing studies in searching for additional ASD risk genes.