Large mosaic copy number variations confer autism risk

Detection of mosaic copy number variants in ASD cohorts

We sought to characterize the contribution of mCNVs arising during early development to ASD risk. We analyzed blood-derived genotype array intensity data from 2,591 autism-affected families in the Simons Simplex Collection (SSC) cohort²⁶ and saliva-derived genotype intensity data from 8,866 autism-affected families in the Simons Powering Autism Research for Knowledge (SPARK) cohort²⁷. All SSC probands and siblings were 3-18 years old at enrollment; most SPARK probands and siblings were in or near the same age range, with a small fraction of older probands (1.2% between age 30-40 and 0.3% over age 40; Supplementary Fig. 1a). After data quality control (Methods), 12,077 probands and 5,500 siblings remained (Table 1). On average 900,935 genotyped variants remained in SSC samples and 579,300 in SPARK samples due to differences in genotyping density between arrays.

We performed haplotype phasing using both a population reference panel and the pedigree structure of the data to obtain near-perfect long-range phase information in offspring. We leveraged the phase information to sensitively detect mCNVs in autosomes of probands and siblings using MoChA³⁰ and checked parental genotypes to ensure events were not germline (Methods; see URLs); we excluded sex chromosomes to avoid confounding from the imbalanced sex ratio between probands and siblings (9,776:2,301 males:females in probands versus 2,718:2,782 in siblings). Following previous studies^29,31, we filtered mCNV calls that exhibited evidence of DNA contamination, and we restricted our analysis to events for which copy number state could be confidently determined (Methods, Supplementary Fig. 2). We further excluded mCNVs frequently observed in age-related clonal hematopoiesis (specifically, focal deletions at IGH and IGL and low-cell-fraction CNN-LOH events^29,31–33), which we expected to be present in a very small fraction of samples (<1%, given the young ages of participants) and unrelated to ASD status. We verified that genotyping intensity deviations within the remaining mCNVs were consistent with estimated mosaic cell fraction and copy number state (Supplementary Fig. 2).

We detected 64 mCNVs in 59 individuals (35 gains, 24 losses, and 5 CNN-LOH in 0.34% of SSC and SPARK samples; Table 1 and Supplementary Table 1) ranging in cell fraction from 2.8% to 73.8% (median = 27.1%) and in size from 49.3 kb to 249.2 Mb (median = 2.5 Mb) (Fig. 1a). All but one carrier was younger than 28 years old (oldest: 47 y.o.; median: 12 y.o.). Of the 64 detected mCNVs, 45 events were present in 44 unique probands (0.36%) and 19 events were present in 15 unique siblings (0.27%), with one sibling carrying five events on a single chromosome, reminiscent of chromothripsis³⁴ (Supplementary Fig. 3, Supplementary Note 1). We did not observe a significant increase in mCNV detection rate with increasing age (Supplementary Fig. 1b), consistent with our filtering of age-related clonal hematopoiesis events, nor did we observe a bias in the parental haplotype on which mCNVs were located (Supplementary Table 1, Supplementary Fig. 4, Methods).

• Download figure
• Open in new tab

Figure 1: ASD probands carry a burden of large mosaic CNVs.

a, Histogram of mosaic CNV sizes in probands (gold) and siblings (purple). b, Box-and-whisker plots of mCNV sizes in probands versus siblings across all events and stratified by copy-number state (Gain, Loss, or CNN-LOH); P-values, Mann-Whitney U-test. No CNN-LOH events were detected in siblings. c, Percent of probands and siblings carrying a mCNV >4 Mb in size combined across cohorts (filled diamonds) and stratified by cohort (unfilled circles); error bars, 95% CI (Wilson score interval). d, Percent of probands and siblings carrying a mCNV of length at least L, with L varying from 0-8 Mb; shaded regions, 95% CI. The burden in probands is robust to the choice of size threshold.

ASD probands carry a burden of large mosaic CNVs

We investigated whether mCNVs in probands had properties distinguishing them from mCNVs in siblings. The size distribution of mCNVs was markedly different between the two groups (Fig. 1a): probands carried mCNVs that were an order of magnitude longer on average than those in siblings (median length = 7.8 Mb vs. 0.59 Mb, P = 1.6×10⁻³ Mann-Whitney U-test, Fig. 1a,b), a trend that was consistent across copy number states and cohorts (Fig. 1b, Supplementary Fig. 5). We did not observe a significant difference between mosaic cell fractions of mCNVs in probands and siblings (Supplementary Fig. 6), although this may reflect our limited power to detect mCNVs present in small proportions of cells (Supplementary Note 2, Supplementary Fig. 7)

Of mCNVs >4 Mb long, 25 were carried by probands and only 1 was found in a sibling. This significant burden in probands of mCNVs >4 Mb (odds ratio (OR) = 11.4, 95% confidence interval (CI) = 1.5-84.2, one-sided Fisher’s exact P = 7.4×10⁻⁴) was consistent across cohorts (P = 0.043 in SSC and P = 6.6×10⁻³ in SPARK; Fig. 1c, Supplementary Fig. 8); was robust to multiple hypothesis correction for the choice of the 4Mb length threshold (P = 2.3×10⁻³ after adjusting for considering all possible thresholds; Methods; Fig. 1d); and was robust to the exclusion of carriers >20 y.o. (P = 1.7×10⁻³). These results suggest that large mCNVs – which appear to be extremely rare in unaffected individuals – explain ∼0.2% of ASD cases (95% CI=0.08-0.29%; Methods).

While we did not observe a significant increase in the frequency of smaller (<4 Mb) mCNVs in probands versus siblings (P = 0.99), we wondered whether some smaller mCNVs in probands might contribute to ASD by altering dosages of specific genes previously implicated in autism susceptibility (“ASD genes”). We analyzed overlap of mCNVs with a curated set of 222 high-confidence ASD genes from the SFARI Gene database (Methods). Smaller (<4 Mb) mCNVs in probands overlapped ASD genes more often than expected by chance (Expected = 1.42, Observed = 4; P = 0.044), in contrast to smaller mCNVs in unaffected siblings (Expected = 1.69, Observed = 1; P = 0.84), suggesting that smaller mCNVs also contribute to the etiology of ASD. (This analysis was uninformative for large mCNVs, most of which are expected to overlap at least one ASD gene by chance.)

When possible, we verified that probands carrying an mCNV did not carry other high-risk germline genetic mutations. Of 15 SSC probands with mosaic CNVs, four also carried previously reported dnCNVs¹¹; only one was >1 Mb in size and none overlapped ASD genes. One proband with a mCNV also carried a previously reported de novo loss-of-function (dnLoF) variant in AFM⁷, a gene with no known connection to ASD (Supplementary Table 2). (We were unable to perform an equivalent analysis for SPARK probands as curated sets of de novo germline CNVs and LoF variants are not yet available for this cohort.) These results indicate that mCNVs comprise orthogonal genetic aberrations that independently contribute ASD risk.

Differences between germline and mosaic CNVs

Interestingly, mCNVs in probands had characteristics different from germline dnCNVs previously reported in SSC probands. Mosaic CNVs were significantly larger than dnCNVs (median length = 7.8 Mb vs. 0.92 Mb, P = 7.3×10⁻⁵; Fig. 2a) and did not exhibit focal recurrence in any genomic location (Supplementary Fig. 9 and Supplementary Note 3). Furthermore, ASD-associated dnCNVs that are observed recurrently in ASD probands¹¹ (ASD-dnCNVs; e.g., 16p11.2 deletion/duplication, 22q11.2 deletion/duplication) were notably absent from ASD probands as mCNVs (Fig 2b).

• Download figure
• Open in new tab

Figure 2: Mosaic and germline CNVs have different properties and effects.

a, Sizes of mCNVs compared to sizes of de novo CNVs (dnCNVs) identified by Ref. 11 in SSC probands. De novo CNVs <100 kb in size were removed to account for our limited sensitivity to detect mosaic CNVs <100 kb in size. b, Percent of samples carrying a germline or mosaic CNV (gain or loss) in each of eight ASD-dnCNV regions in ASD cohorts (SSC + Autism Genome Project for germline; SSC + SPARK for mosaic) or the UK Biobank. Each marker indicates the percent of carriers of a specific ASD-dnCNV; markers corresponding to 16p11.2 CNVs are indicated with callouts. c, Effects of germline and mosaic 16p11.2 deletions on phenotypes previously associated with 16p11.2 deletions (units, s.d.).

We hypothesized that such mosaic analogues of ASD-dnCNVs 1) may be very rare or 2) may confer little or no ASD risk. To obtain further insight into both questions, we examined mosaic events previously detected in a population sample of 454,993 individuals of European ancestry in the UK Biobank (UKB)³⁰. Mosaic analogues of ASD-dnCNVs occurred much more rarely than their germline counterparts (Fig 2b, Supplementary Table 3); among eight previously-reported ASD-dnCNVs¹¹, only 16p11.2 deletions were detected recurrently in the mosaic state (in 73 UKB samples comprising 0.016% of the cohort; Supplementary Note 4). Mosaic status was not associated with mental health conditions (Supplementary Table 4), although our power was very limited by the sparsity of reported mental health diagnoses.

To better understand the phenotypic relationship between germline ASD-dnCNVs and mosaic analogues, we identified carriers of germline 16p11.2 deletions in the UK Biobank (Supplementary Fig. 10, Methods) and compared their phenotypes to those of mosaic 16p11.2 deletion carriers. While we were underpowered to directly measure ASD risk conferred by 16p11.2 deletions, we could compare the effects of germline and mosaic 16p11.2 deletions on quantitative traits measured in UKB. Consistent with previous reports^35–38, germline 16p11.2 deletions were strongly associated with several traits including fewer years of education, increased BMI, and decreased height. However, mosaic 16p11.2 deletions were not associated with any of these traits (Fig 2c) even when restricting to events at high cell fractions (Supplementary Table 5). These data reinforce our observation that the burden of mCNVs in ASD probands was driven by large mCNVs that disrupted large swaths of the genome; smaller mosaic CNVs may have limited phenotypic consequences, even when disrupting ASD-associated regions.

Mosaic CNV length associates with ASD phenotype severity

We next determined whether properties of mCNVs carried by probands were associated with ASD severity in these probands. ASD phenotypes were assessed with three measures common to both the SSC and SPARK cohorts, of which one measure – the Social Communication Questionnaire (SCQ) – was available for the majority of proband mCNV carriers in both cohorts (13 of 17 SSC carriers and 20 of 29 SPARK carriers; Supplementary Table 1). The SCQ is a standardized evaluation form completed by a parent rating an individual’s symptomatic severity throughout his or her developmental history; higher scores reflect a more severe ASD phenotype. Larger mCNV size significantly correlated with increased ASD severity as quantified by SCQ score (Fig. 3; Pearson correlation R = 0.43, P = 0.016). The longest mCNVs were CNN-LOH events; while such events do not modify gene dosage, they could produce phenotypic consequences by converting heterozygous gene-disrupting variants to the homozygous state (Supplementary Table 6, Supplementary Note 5). These results further highlight the important role of size when considering the potential pathogenicity of a mosaic event: larger mCNVs appear to be both more likely to result in ASD and to produce more severe phenotypes. We did not observe an association between mCNV cell fraction and phenotypic severity (Fig. 3, Supplementary Fig. 11).

• Download figure
• Open in new tab

Figure 3: Mosaic CNV size positively correlates with ASD severity.

ASD severity (quantified by the Social Communication Questionnaire (SCQ) summary score) versus mCNV size. For probands with more than one mCNV, the longest event size is used. Marker color indicates mosaic copy number state; marker size indicates mosaic cell fraction. Events discussed in the main text are labeled with black text; events discussed in Supplementary Notes are labeled with grey text. R, Pearson correlation coefficient. Shaded region, 95% CI. The association was robust to the scale used for CNV size (Spearman rank correlation Rs = 0.42, P = 0.019).

Mosaic CNVs correlate with individual-level clinical observations

We looked for evidence directly linking specific mCNVs to reported clinical symptoms. By cross-referencing the genes disrupted by individual mCNVs with known syndromes associated with those genes, we identified four probands in which the mosaic mutation appeared directly linked to the proband’s symptoms (Fig. 4, Supplementary Fig. 12).

• Download figure
• Open in new tab

Figure 4: Carriers of mosaic CNVs exhibit symptoms commonly observed with analogous germline mutations.

a, 18q distal deletion. b, NRXN1 deletion. c, 11p CNN-LOH. Top, diagrams of mosaic mutations altering inherited chromosomes in a fraction of cells. Paternal and maternal haplotypes are colored dark and light blue, respectively, with genes or regions of interest labeled below. Middle, description of mutations and observed clinical phenotypes. Bottom, maternal allele fraction at heterozygous SNPs (binned into groups of two adjacent SNPs) and total genotyping intensity (log R-ratio; LRR) at all SNPs genotyped on the chromosome (binned into groups of four adjacent SNPs); SNPs within the mCNV are highlighted.

Two probands carried mosaic 18q distal deletions removing 29.2 and 17.9 Mb of sequence from chromosome 18 in 27% and 16% of cells, respectively (Fig. 4a, Supplementary Fig 12). Germline 18q distal deletions (OMIM: 601808) are well-characterized causes of intellectual disability and ASD³⁹. While causal genes in this region have not been fully characterized, larger deletions which encompass the gene TCF4 and all distal genes produce profound intellectual disability and little to no verbal communication (Pitt-Hopkins Syndrome, OMIM:610954), while smaller deletions encompassing fewer genes result in relatively mild cognitive impairment^40,41. Of the two mosaic 18q distal deletions, the larger one extended beyond TCF4 and the smaller one did not. Consistent with their germline analogues, the proband with the TCF4 deletion was non-verbal, while the proband with the smaller deletion had an IQ in the normal range (full-scale IQ = 97, NVIQ = 98) and mild adaptive impairment by the Vineland Adaptive Composite Standard Score (VSS=66).

One proband carried a mosaic deletion of NRXN1 (OMIM:614332), which has a well-documented (but incompletely penetrant) association with ASD, intellectual disability and speech delay^42,43. The proband’s mosaic deletion encompassed the entirety of NRXN1 on his paternal haplotype in 8% of cells. Furthermore, on the maternal haplotype, the individual carried an inherited, rare start-loss variant of the beta isoform of NRXN1 (Fig. 4b, Supplementary Table 6, Supplementary Note 5). At age 13, the proband was reported to have an IQ in the range of 55-69 and slight language delay consistent with at least mild intellectual disability. Furthermore, the proband had ADHD, a condition also often associated with NRXN1 deletions⁴⁴. The mother exhibited no evidence of intellectual disability despite carrying the start-loss variant in NRXN1. This finding is consistent with previous reports that NRXN1 LoF variants are incompletely penetrant^45,46 and suggests that the observed germline-mosaic compound heterozygosity contributes to the proband’s clinical symptoms.

Another proband carried an acquired paternal uniparental disomy (UPD) of nearly the entirety of 11p in 32% of cells (Fig. 4c). The 11p15.5 region contains numerous paternally and maternally imprinted genes, and germline disruption of this region is known to produce syndromic growth disorders: Beckwith-Wiedemann syndrome (BWS, OMIM:130650; an overgrowth condition associated with hypermethylation) and Silver-Russell syndrome (SRS, OMIM:180860; an undergrowth condition associated with hypomethylation)⁴⁷. The proband exhibited abnormally slow growth, macrocephaly, and feeding difficulties, all of which are common symptoms of SRS. SRS has also been associated with increased risk of ASD and intellectual disability⁴⁸. While paternal 11p UPD is usually associated with BWS and maternal 11p UPD is usually associated with SRS, cases in which imprinting disruption led to the opposite phenotype have been reported⁴⁹. Interestingly, we observed one other case of a mosaic 11p UPD impacting 11p15.5 in a sibling with a reported genetic condition (and therefore excluded from our main analyses) (Supplementary Table 7). This individual also had a reported (unspecified) growth disorder.

These case studies reinforce our observation of an overall burden of large mCNVs in ASD probands with concrete examples in which specific mCNVs potentially underlie the disorder via a variety of plausible mechanisms. We also explored six other cases in which mCNVs deleted ASD genes but a direct connection to reported phenotypes was less clear due to the phenotypic heterogeneity of ASD⁵⁰ and the limited phenotype data provided for each proband (Supplementary Notes 6 and 7, Supplementary Fig. 13 and 14).

Mosaic CNVs are present in neurons

Our analyses of SSC and SPARK considered mCNVs ascertained from blood- and saliva-derived DNA. We reasoned that most of these mCNVs were likely present throughout the body given their moderate-to-high cell fractions⁵¹ and the young ages of carriers. To confirm that mCNVs are present in brain tissue, we performed whole-genome sequencing of post-mortem brain tissue from an additional 60 probands obtained through the NIH Neurobiobank and Autism BrainNet (Supplementary Table 8). We genotyped germline variants using GATK HaplotypeCaller best practices⁵² and identified mCNVs using MoChA (Methods).

We found two mosaic events (Supplementary Table 9): a mosaic 10.3 Mb gain of 2pcen-2q11.2 in sample AN09412 (Fig. 5a) and a mosaic loss of Y in BEAR12. We also discovered 15 germline CNVs overlapping ASD genes in other individuals, revealing potential causes of disease in several previously unresolved cases (Supplementary Table 10, Supplementary Fig. 15, Supplementary Note 8).

• Download figure
• Open in new tab

Figure 5: A complex mosaic chromosomal rearrangement present in neurons.

a, Phased allele fraction at heterozygous SNPs on chromosome 2, binned into groups of four adjacent SNPs. SNPs within the mCNV are highlighted, with distinct copy-number states indicated in different colors. Assembly gaps >1 Mb are shaded. b, Estimated mean copy number in each mCNV region as inferred from phased allele fractions (left) and sequencing read depths (right); error bars, 95% CI. Confidence intervals on allele fraction-based estimates are very narrow. c, Inferred structure of a complex duplication consistent with the observed data. Arcs on the ideogram indicate fusions supported by breakpoint analysis. Arrows are a schematic reconstruction of the event (not to scale); each arrow points in the 3’ direction relative to the GRCh37 reference genome. Black arrows indicate genomic regions with a single copy in the proper orientation within the duplicated region. The left breakpoint of the inverted duplication is approximate. d, Experimental validation of the three breakpoints, labeled according to their corresponding segment (inversion, INV; tandem duplication, TD; inverted duplication, ID). Left, fractions of cells containing each breakpoint estimated using digital droplet PCR (ddPCR) on DNA extracted from bulk brain tissue; error bars, approximate 95% CI. Right, validation of co-occurrence of breakpoints in single neurons. Observation of some but not all breakpoints in some neurons is probably explained by locus dropout, a common feature of single cell whole-genome amplification⁸⁰.

The gain event on chromosome 2 in AN09412 was unique in that it appeared to exhibit three segments with varying degrees of mosaicism (Fig. 5a). Using phased allele fractions of germline heterozygous SNPs and depth-of-coverage of sequencing reads, we estimated that the three segments were present in a ratio of 1:3:2 (Fig. 5b). Breakpoint analysis using split reads and discordantly mapped reads revealed three breakpoints (Supplementary Table 11): a tail-to-tail (T2T) inversion of 92.03-99.78 Mb, a tandem duplication (TD) of 99.87-101.94 Mb, and a head-to-head inversion (H2H) located at 102.38 Mb, each of which corresponded to one of the three segments. Using this information, we reconstructed a parsimonious linear structure of the event (Fig. 5c, Methods) consistent with gain of a single complex rearrangement present in 26% of cells (Fig. 5b).

Using quantitative digital droplet PCR (ddPCR), we confirmed that the three breakpoints were present in both neurons and non-neurons at a 26-36% mosaic cell fraction (Fig. 5d), indicating the mCNV arose in a fetal progenitor that gave rise to both neurons and glial cells. (Non-brain tissue was not available for this sample, so we could not investigate the presence of the CNV elsewhere.) We further confirmed that all three breakpoints occur within individual neurons using single cell ddPCR (Fig. 5d) and that none of the breakpoints were present in DNA from a control brain using gel electrophoresis (Supplementary Fig. 16), suggesting that the CNV arose from a single event, likely at a very early stage of development. The CNV spans two genes associated with ASD (MAL and NPAS2), and copy-number alterations of 2q11.2 have been associated with developmental delay^53,54, suggesting that the event may contribute to the proband’s disease.

We also validated the mosaic loss of Y in BEAR12 (Supplementary Fig. 16) and determined that the loss was limited to non-neuronal cell populations. This finding was unsurprising given that the BEAR12 donor was 74 y.o. (the oldest in the cohort) and age-related loss of Y has been reported extensively in blood^55,56 and more recently in aging brain tissue^57,58.

These results complement our analyses of mCNVs in large ASD cohorts, in which we analyzed DNA derived from blood and saliva under the assumption that mCNVs detected at moderate-to-high cell fractions were likely present throughout the body. Our validation of a mCNV in post-mitotic neurons of AN09412 indicates that mCNVs can arise during early development and propagate to multiple cell lineages in the adult body.

Genotyping intensity data

Genotyping intensity data for probands, siblings and parents in SSC and SPARK were obtained from SFARI Base. For each genotyped position, the data included the genotype call, the B allele frequency (BAF; proportion of B allele), and Log-R ratio (LRR; total genotyping intensity of A and B alleles) as estimated by Illumina’s GenomeStudio software.

Three types of genotyping arrays were used for SSC samples: Illumina 1Mv1 (n=1,354 individuals), Illumina 1Mv3 (n=4,626 individuals), and Illumina Omni2.5 (n=4,240 individuals). Details of data generation have been previously described in Sanders et al. 2015 (ref. 11). SPARK samples (n=27,376 individuals) were genotyped on the Illumina Infinium Global Screening Array-24 v.1.0. Details were previously described in Feliciano et al. 2018 (ref. 65). We did not analyze SPARK samples which had been previously genotyped on a different array as part of a pilot study (n=1,361 individuals).

We defined probands to be individuals with a diagnosis of ASD. We defined “unaffected siblings” as family members without an ASD diagnosis in the same generation as a proband (most of which were siblings). We defined parents as unaffected individuals with a proband as a biological child.

Converting Illumina Final Reports to BCF format

Genotyping intensity data for SSC were distributed in the Illumina Final Report format with genotyped positions reported with respect to the hg18 human reference genome. Positions were lifted over to hg19 coordinates based on rsID number. Positions without a rsID were discarded. Final Reports were converted to the BCF format and genotypes were converted from Illumina TOP-BOT format to dbSNP REF-ALT format using custom in-house scripts (positions for which TOP-BOT format could not be unambiguously converted to REF-ALT format were discarded). Samples from each of the three arrays were processed as separate batches.

Genotyping intensity data for SPARK were converted from PLINK PED format to BCF format using the recode option in plink1.9. Genotypes were converted from Illumina TOP-BOT format to dbSNP REF-ALT format using a modified version of the bcftools plugin fixref (URLs). Only single-nucleotide variants were retained for analysis.

LRR denoising for SPARK samples

We observed genome-wide spatial autocorrelation “wave” patterns⁶⁶ in numerous SPARK samples. Since the wave pattern was consistent across samples for each chromosome, we corrected the bias using the following algorithm based on principal components analysis:

1. Determine the mean LRR per chromosome per sample. For each sample, mean-shift the LRR signal genome-wide by the median of chromosome means for that sample.

2. For chromosome i:

   1. Determine the cohort-wide LRR deviation for the chromosome i as the median of mean chromosome i LRR signal across samples. Mean-shift each sample’s chromosome i LRR signal by the cohort-wide LRR deviation.

      1. To prevent confounding due to sex, this correction is performed independently for males and females.

3. For each chromosome i:

   1. Project the LRR matrix (number of samples by number of genotyped positions on chromosome i) onto the space spanned by its top k principal components. Subtract the projected matrix from the full LRR matrix.

Steps 1-2 of the algorithm mean-center the LRR signal genome-wide across an individual and per chromosome across the cohort. This is necessary to prevent PCA from projecting away mean-shifts due to large mosaic CNVs. Step 3 removes the variance explained by the top k principal components. In practice, we found k = 10 effectively removed the wave pattern (Supplemental Fig. 16).

PCA analysis was performed using the PCA method from the python package sklearn⁶⁷, which implements efficient PCA using randomized singular-value decomposition. LRR values were extracted from BCF files using ‘bcftools query’ and corrected values were incorporated into BCF files using ‘bcftools annotate’. One sample with >5% genotype missingness was excluded from the correction procedure. On average across autosomes, the top 10 PCs explained 57.1% of LRR variance in the SPARK cohort.

Variant-level quality control

We excluded genotyped variants with high levels of genotype missingness (>2%), evidence of excess heterozygosity (P < 1e-6, one-sided Hardy-Weinberg equilibrium test), and unexpected genotype correlation with sex (P < 1e-6, Fisher exact test comparing number of 0/0 genotypes vs. number of 1/1 genotypes in males and in females). We also exclude genotyped variants falling within segmental duplications with low divergence (<2%). Variant-level QC was performed for each array independently. The number of genotyped variants and number of variants excluded by QC are:

Sample-level quality control

We calculated two statistics to detect sample contamination: BAF concordance and BAF autocorrelation. Given a heterozygous SNP has BAF >0.5 (<0.5), BAF concordance is the probability that the following heterozygous SNP is >0.5 (<0.5). BAF autocorrelation is the correlation of the BAF at a heterozygous SNP with BAF at the neighboring (downstream) heterozygous SNP. For each sample, we calculated the statistic for each chromosome independently and took the median across all chromosomes as the sample value.

Neighboring positions with heterozygous genotypes in the genome are expected to have uncorrelated genotype intensity measures on an array. BAF concordance and BAF autocorrelation significantly higher than, respectively, 0.5 and 0, could reflect sample contamination with DNA from another individual because allelic intensities will be correlated at variants within haplotypes shared between the sample DNA and contaminating DNA. In SSC, we removed samples with BAF concordance > 0.51 or BAF autocorrelation >0.03, resulting in the exclusion of 10 probands and 9 siblings. We also excluded an additional proband (array ID: 7306256088_R02C01) with evidence of a large amplitude LRR wave pattern. In total, 2,594 probands and 2,424 siblings from SSC passed quality control.

In SPARK, we observed genome-wide evidence of BAF correlation between contiguous genotyped positions in high-quality samples. Thus, BAF concordance and BAF autocorrelation were not informative measures of contamination. Instead, we excluded samples with evidence of multiple very low cell-fraction CNN-LOH events (<10 % of cells and LRR deviation from zero < 0.2) because the probability of observing two or more true CNN-LOH events in a sample is exceedingly small given the young age of the individuals. We further removed any samples from individuals that had also participated in SSC (n=352) and one additional proband (SP0072755) that had an uncorrected LRR wave pattern after LRR denoising, resulting in exclusion of 622 probands and 54 siblings. Finally, we removed 37 siblings with a reported genetic diagnosis (of which one carried a mCNV; see main text). In total, 9,483 probands and 3,076 siblings from SPARK passed quality control.

Haplotype phasing

We used Eagle2 (ref. 28) (default settings) and the Haplotype Reference Consortium⁶⁸ phasing panel to perform statistical haplotype phasing of SSC samples. For probands and siblings we additionally used parental genotypes to correct phase-switch errors using the bcftools plugin trio-phase included with MoChA. Given the size of the SPARK cohort (>27,000 samples), we performed within-cohort statistical phasing using Eagle2. We additionally corrected proband and sibling phase estimates using parental genotyping data when available (at least one parent was also genotyped for the vast majority of probands and siblings). The combination of statistical haplotype phasing and pedigree-based phasing resulted in nearly perfect long-range phase information without phase-switch errors.

Discovery of mosaic CNVs

We applied MoChA to detect mosaic CNVs. The general statistical approach implemented in MoChA has been previously described²⁹. In brief, mCNVs result in allelic imbalance between the maternal and paternal haplotypes. Thus, the BAF of heterozygous SNPs within an mCNV will consistently deviate from the expected value of 0.5 toward either the paternal allele or the maternal allele. Such deviations can be sensitively detected even at low cell-fractions using long-range phase information provided the event is long enough to contain multiple genotyped heterozygous SNPs. Formally, MoChA uses a hidden Markov model (HMM) to search for consistent deviations. Gains (losses) also result in an increase (decrease) of total LRR signal with magnitude proportional to the cell fraction of the event; an HMM can also be used to detect LRR deviations from zero. Incorporation of phase information particularly increases sensitivity to detect large, low-cell fraction CNVs relative to previous models²⁹.

The details of MoChA differ from the previously described approach in two ways. First, MoChA uses two independent models to search for mCNVs: a haplotype-phase model (BAF+phase) as described in Loh et al. 2018 (ref. 29) and an LRR and (unphased) BAF model (LRR+BAF) similar to previous models for the detection of germline CNVs⁶⁹. A CNV is reported if it is discovered by either model. The introduction of the LRR+BAF model enables detection of germline (or very high cell fraction mosaic) losses and germline duplications including more than two haplotypes. Second, MoChA uses the Viterbi algorithm to search for deviations in either the phased BAF signal or the LRR signal instead of computing total likelihoods and applying a likelihood ratio test. The Viterbi algorithm is more direct, but its calibration is less precise when detecting very low cell fraction events. However, since we were interested in higher-cell-fraction mCNVs arising during early embryogenesis, such sensitivity was not necessary for this study.

Central to the sensitivity of MoChA is the quality of the long-range phase information. As discussed above, the combination of statistical haplotype phasing and pedigree phasing using parental genotypes resulted in near-perfect long-range phase information without phase-switch errors.

Classification of mosaic copy-number state

We needed to sensitively distinguish age-related and early-developmental mCNVs in a way that was robust to LRR signal noise due to e.g. GC content. Previous work on mosaic CNVs have not typically distinguished between age-related and early-developmental events. Thus we developed a new statistical method to classify events as gains, losses, CNN-LOH, or unknown using an Expectation-Maximization (EM) algorithm similar to K-means clustering where each cluster is defined by a line instead of a centroid. Let X = |ΔBAF| be the absolute deviation from 0.5 of phased BAF estimated across an event; let Y = |ΔLRR| be the absolute deviation from 0 for LRR estimated across an event, and let C ∈ {Gain, Loss, CNN-LOH} denote the copy-number state of the mosaic mutation. Then for gains, X and Y will linearly increase according to Y = Xβ_Gain + ∈, where β_Gain > 0; for losses, Y will linearly decrease as X increases according to Y = Xβ_Loss + ∈, where β_Loss < 0; and for CNN-LOH, Y = ∈, where ∈ ∼ N(0, σ²) is Gaussian noise in the estimation of X and Y.

Given a set of events, the parameters of the linear models and the copy-number state C_i for each event i are unknown. We thus iteratively apply the following EM algorithm:

1. Randomly initialize β_Gain ∈ (0, 3] and β_Loss ∈ [-3, 0) and set β_CNN-LOH = 0

2. Assign each event i a copy-number state C_i using least-squares classification:

3. Estimate the linear model parameters β_C for C ∈ {Gain, Loss} using univariate linear regression without an intercept term applied to all events assigned to class C in step 2:

   1. Since β_CNN-LOH = 0 is known, it is not re-estimated.

4. Repeat 2-3 until convergence.

5. Estimate σ_Gain and σ_Loss using univariate linear regression on the events classified as gains and losses, respectively.

To classify mCNV copy-number states in probands and siblings, the model was first trained on mCNVs in parents (after removal of germline CNVs). Events in probands and siblings were then classified using the linear model parameters estimated from the parents. The method implicitly accounts for error in LRR and BAF measures and thus is robust to noise in these signals.

We applied an additional step to improve classification of events extending to telomeres, given that CNN-LOH events generally arise due to mitotic recombination and therefore terminate at a telomere. To ensure apparent gains and losses terminating at a telomere were not misclassifications, we calculated the Bayes factor to compare the likelihood the event arose under the Gain or Loss model against the likelihood under the CNN-LOH model: where C ∈ {Gain, Loss} and σ_C is the standard deviation estimated from fitting the model on parental data. If B < 10 for a putative gain or loss terminating at a telomere, the copy-number state was reclassified as unknown.

Filtration of mosaic CNV calls

Determination of haplotype-of-origin

For mosaic gains and losses, the parental haplotype-of-origin was defined to be the haplotype carrying the mCNV. For CNN-LOH the parental haplotype-of-origin was defined to be the haplotype that was duplicated. To assign haplotype-of-origin, we calculated the average ALT allele frequency of heterozygous SNPs at which the ALT allele was unambiguously inherited from the father and the average ALT allele frequency of heterozygous SNPs at which the ALT allele was unambiguously inherited from the mother. For losses, the haplotype-of-origin was paternal if the average allele fraction of paternal SNPs was less than that of maternal SNPs; otherwise the haplotype-of-origin was maternal. For gains and CNN-LOH, the haplotype of origin was paternal if the average allele fraction of paternal SNPs was greater than that of maternal SNPs; otherwise the haplotype-of-origin was maternal.

Burden analysis

The statistical significance of the hypothesis that probands carry more mCNVs > 4 Mb than their unaffected siblings was quantified using a one-sided Fisher’s Exact Test. 95% confidence intervals for the percent of samples carrying an mCNV were calculated using Wilson’s score interval. To adjust the burden p-value for multiple possible choices of the size threshold for defining “large mCNVs,” we performed the following permutation analysis: proband and sibling labels of mosaic CNVs were randomly permuted based on the total number of probands and siblings in our study. We then determined the p-value of the most significant burden across all size thresholds for the permutation. This procedure was repeated 100,000 times. We calculated the threshold-adjusted p-value as where P_obs is the uncorrected p-value from the observed data, P_i is the maximum burden p-value from permutation i, and 1 is the indicator function.

The percent of ASD cases explained by large (> 4 Mb) mosaic CNVs was estimated as the difference between the percent of probands carrying a large mCNV and the percent of siblings carrying a large mCNV. The 95% CI between proportions was estimated using Wilson’s score interval as modified by Newcombe⁷¹.

Overlap of mCNVs with ASD genes

We downloaded all genes included in the SFARI Gene database of genes implicated in ASD. We restricted the list to the 222 genes that are classified as “Category 1” (high confidence), “Category 2” (strong candidate), or “S” (syndromic). We refer to this restricted list of genes as “ASD genes”. We determined whether mosaic CNVs overlapped ASD genes by annotating their overlap with all genes in the RefSeq database and intersecting the name of the RefSeq genes with the ASD gene list.

To determine whether a set of mCNVs overlapped ASD genes more often than expected by chance, we randomly permuted the mCNVs in probands around the genome K times, excluding assembly gaps >1 Mb in size in the hg19 reference. After each permutation we determined the number of segments overlapping an ASD gene. Let N_obs be the number of mCNVs overlapping ASD genes in the observed data. Let N_i be the number of permuted segments overlapping ASD genes in permutation i. The P-value of observing N_obs or more overlaps by chance is where 1 is the indicator function. When testing ASD gene overlap for short events (<4 Mb), we used K=10,000. For long events, we used K=1,000 for computational efficiency. We excluded CNN-LOH events when testing long events because they were too large to be randomly permuted.

Counts of germline ASD-associated CNVs

Counts of germline ASD-associated CNVs in ASD cohorts were obtained from Sanders et al.¹¹, Table 2 (which included samples from SSC and the Autism Genome Project, AGP). Counts of germline ASD-associated CNVs in UK Biobank individuals were obtained from Crawford et al.³⁷.

Identification of 16p11.2 germline deletion carriers in the UK Biobank

We extracted LRR and genotype calls from the 16p11.2 ASD-associated region listed in Table 2 of Sanders et al.¹¹ for individuals in the UK Biobank. Germline carriers of 16p11.2 deletions were defined as individuals with average LRR < −0.5 and <5 heterozygous SNP calls across the region (Supplementary Fig. 10).

Phenotype associations of germline and mosaic CNVs in ASD-associated regions

We defined high-confidence ASD-associated CNV regions as those listed in Table 1 and 2 in Sanders et al.¹¹ expanded by ∼1.5 Mb on either side (Supplementary Table 4 lists the exact expanded regions). We identified carriers of mosaic CNVs in the UK Biobank reported by Loh et al.³⁰ falling within the ASD regions. We refer to these individuals as ASD-dnCNV-analogue carriers. We used self-reported responses to the UK Biobank Mental Health Questionnaire to count the number of ASD-dnCNV-analogue carriers with a diagnosis of ASD, SCZ, BP, depression, or anxiety.

Following Owens et al.³⁶, we quantified the association between carrier status of germline or mosaic 16p11.2 deletions and phenotypes using the following linear regression model for continuous phenotypes: where y_i is the phenotype of individual i, x_c,i is the 16p11.2 CNV carrier status of individual i,x_age,i is the age of individual i, x_sex,i is the sex of individual i, x_array,i is the array used to genotype individual i, is the j^th genetic principal component of individual i, βs are the corresponding effect sizes and ∈_i ∼ N(0, σ²) is the remaining phenotypic variance. For binary phenotypes, we applied logistic regression with the same covariates. Continuous phenotypes were inverse-normal transformed within sex strata after adjusting for relevant covariates prior to analysis⁷². We restricted to individuals passing quality control filters from ref. 30 and of self-reported European ancestry.

We identified a set of quantitative traits and medical outcomes previously associated with 16p11.2 germline deletions^35–38. The association results for mosaic 16p11.2 deletions, high cell-fraction mosaic 16p11.2 deletions (CF > 0.3), and germline 16p11.2 deletions for all tested traits are reported in Supplementary Table 5. Medical phenotypes were coded as binarized versions of the following data fields from the UK Biobank Data Showcase: Renal failure: 132030, 132032, and 132034; Obesity: 130792; Heart failure: 131354.

Determining carriers of high-risk germline de novo variants

Curated germline de novo CNVs and LoF variants in SSC individuals^7,11,73 were obtained from ref. 11. We cross-referenced our list of mCNV carriers with carriers of de novo CNVs and LoF variants. For any mCNV carriers that also carried a de novo CNV, we determined whether the dnCNV overlapped an ASD gene as described above. The list of high confidence germline de novo CNVs was also used to estimate the size distribution of de novo CNVs in Fig. 2a; we removed de novo CNVs <100 kb in size to account for our limited sensitivity to detect mosaic CNVs below that size threshold.

Genotype-phenotype associations

We obtained phenotype data for individuals in SSC and SPARK from SFARI Base (SSC version 15 and SPARK version 2). Of the three ASD severity measures shared between SSC and SPARK (Development Coordination Disorder Questionnaire, DCDQ; Repetitive Behavior Scale-Revised, RBS-R; and Social Communication Questionnaire, SCQ) only SCQ was missing in less than 50% of SSC and SPARK samples. We measured association between SCQ score and mosaic CNV properties (size and cell fraction) using both Pearson and Spearman rank correlation. Z-normalizing SCQ scores independently in SSC and SPARK prior to association did not qualitatively change the results.

Identification of putative damaging variants within mCNVs in SPARK individuals

We obtained from SFARI Base exonic SNPs and indels detected in whole-exome sequencing data of SPARK individuals. In carriers of mosaic losses and CNN-LOH, we identified rare, putative damaging variants within the mCNV, defined as 1) variants with cohort variant allele frequency <1%; and 2) annotated as “High Impact” (start-lost, stop-lost, stop-gain, frameshift, splice-acceptor, splice-donor) or annotated as missense with CADD >20 (ref. 74) by Variant Effect Predictor⁷⁵.

Analysis of brain tissue

Inferred structure of a complex duplication

We inferred a linear structure of the complex duplication consistent with the following observations: three segments with relative abundance of +1 copy, +3 copies, and +2 copies; a tail-to-tail (T2T) inversion fusing 92.04 Mb to 98.78 Mb; a tandem duplication (TD) of 99.87-101.94 Mb; and a head-to-head (H2H) inversion fusing 102.382 Mb to 102.383 Mb. We first observed that each breakpoint corresponded to a segment with unique copy state: T2T inversion corresponded to +1 copy state; TD to +3 copy state, and H2H to a +2 copy state. We thus concluded that the tandem duplication must result in an additional three copies of 99.87-101.94 Mb and the H2H inversion is likely the result of an inverted duplication resulting in two copies of ∼102.0-102.382 Mb separated by a 1 kb segment (102.382-102.383 Mb) in the proper orientation (where the left breakpoint at ∼102.0 Mb is approximate because it is estimated based on discontinuity in allele fraction and read depth estimates rather than direct observation); we estimated via read depth that the segment 102.382-102.383 Mb is present in a +1 copy state. We further concluded that the duplication carries one copy of 92.04-98.78 in an inverted 3’-5’ orientation and one copy of 99.78-99.87 Mb in the proper 5’-3’ orientation.

Plotting mosaic CNV events

Mosaic CNV events with ideograms and gene / region annotations were plotted using a modified version of pyGenomeTracks⁷⁹.

Description of box plots

All box plots have the following properties: center line is the median, box limits are upper and lower quartile, and whiskers are 1.5x interquartile range. Outliers are not included in Figure 2a for clarity.