Isolation of neuronal nucleiHuman prenatal brain tissue was obtained from the NIH NeuroBioBank. Tissue from the frontoparietal lobes was sectioned into 40 μm sections with a cryostat. The cortical plate was dissected using a microscalpel under a light microscope and dissociated into nuclei using gentle dounce homogenization in 1.5 ml lysis buffer (10 mM Tris-HCl pH 8.0, 250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 0.1% Triton-X, 1× cOmplete Mini Protease, 1 mM DTT). Nuclei were pelleted by centrifugation for 10 min at 900g. The supernatant was removed, and the nuclei were resuspended in sorting buffer (1× PBS pH 7.4, 0.8% BSA, 0.1 μg ml−1 DAPI). Flow cytometry was performed, gating for FSC, SSC and DAPI content to isolate single nuclei. For postnatal brain nuclei, staining was performed with both DAPI and Alexa Fluor 488 conjugated anti-NeuN at 1:1,000 for 1 h at 4 °C (Sigma Aldrich MAB377X) as previously described50.Single-nuclei whole genome library preparationTbA has been used in multiple contexts for scWGA under different names, with a core workflow comprising (1) nucleosome depletion; (2) Tn5 tagmentation and direct adapter addition; and (3) library amplification. Single nuclei were lysed in microwells in 2 μl lysis buffer (20 mM NaCl, 20 mM Tris-HCl pH 8.0, 0.15% Triton-X, 1 mM EDTA, 250 mM DTT, 0.05 μl Thermolabile Proteinase K (NEB P8111S)) at 30 °C for 1 h, followed by proteinase inactivation at 55 °C for 10 min. Proteinase K digests histones within nucleosomes and improves accessibility of the genome by eliminating chromatin structure. Tn5 was loaded according to the manufacturer’s instructions (Diagenode C01070010-20), and 1 μl of Tn5 diluted 1,000-fold was used per single nucleus. Tagmentation was performed using manufacturer’s buffers and instructions (Diagenode). Single-cell genome libraries were generated by 16 cycles of PCR using NexteraXT v.2 indexes. Rows of 12 single cells were pooled for cleanup using the Zymo Clean and Concentrator Kit. Each row of single cells was eluted in 42 μl Zymo elution buffer for subsequent AMPure XP size selection at 0.8×. Four rows (48 cells) were pooled for sequencing in equimolar concentration, and sequencing was performed on NovaSeq 6000 (Novogene) for a target library size of 0.5 Gb per single-cell library. Illumina reads were aligned to the human reference with decoy sequence GRCh37d5 (hs37d5) using bwa-mem51, followed by duplicate marking and base quality score recalibration according to GATK best practices52.Single-cell library quality assessmentTo ensure analysis of only high-quality single-cell genomes, we first removed cells with fewer than 10 million total covered bases (the ‘total’ column from GATK DepthOfCoverage, that is, summed per-base depth across all positions), and then used an NPC line dataset as the basis for transfer learning to distinguish high-quality from poor-quality sequencing profiles. The TbA libraries from single cells of the NPC line allowed us to establish a benchmark of sample quality, that is, on staurosporine treatment 24 h to induce apoptosis, we created a labelled TbA dataset comprising (1) high-quality, recently viable cells; (2) cells marked for apoptosis by Annexin V staining with nuclear envelope still intact; and (3) cells with compromised nuclear integrity marked by propidium iodide (Extended Data Fig. 1c). The NPC line was differentiated as previously described53 using an induced pluripotent cell line (iPSC) derived from well-characterized human fibroblasts54. NPC identity was confirmed (21 March 2021) by immunohistochemistry using antibodies against PAX6 (BD Pharmingen 561462; diluted 1:300) and Nestin (Millipore Sigma HPA026111; diluted 1:500).To categorize cells computationally, we first normalized read counts for each cell in the NPC dataset using the read-depth normalization function of HiScanner34 with bin size set to 500 kb, resulting in a cell-by-bin matrix. Next, we used Scanpy55 to scale the read counts, standardizing variance across cells, followed by principal component analysis for dimensionality reduction. Cell-to-cell similarity was quantified by constructing a neighbourhood graph with scanpy.pp.neighbors. Clustering was performed using the Leiden algorithm with a resolution parameter of 0.2. After visual inspection and correlating cluster-wise read count profiles with existing replication timing tracks, library sizes and MAPD scores, we determined that the three Leiden clusters (QC Leiden 0, 1 and 2) in NPC data are high quality, atypical and poor quality, respectively. We subsequently applied scanpy.tl.ingest to project these labels and embeddings from NPC to human brain single-cell TbA data with default parameters. Only cells labelled as ‘high quality’ (Leiden 0) and MAPD Single-cell CNV calling pipelineAfter filtering at the cellular and sample level as described above, we called CNVs on single-cell genomes. For TbA data, CNV calling was based on read depth followed by validation using FOD (Supplementary Methods). For PTA data, CNV calling jointly leveraged read depth and B-allele frequency (BAF).TbA: CNV calling based on read depthHiScanner34 was applied in rdr_only mode with specific parameters optimized for our analysis. The bin size was set to 500 kb mappable positions. The hg19/GRCh37 genome was used as the reference. λ was set to 64, and ploidy restriction was enabled (max_wgd=1 and restrict_gamma=true). The multisample option was disabled (multisample=false) to process each cell individually. HiScanner v.1.6.1 was installed through pip. To ensure the reliability of detected CNVs, we implemented a stringent filtering approach: CNV candidates with P ≥ 0.05 or those spanning less than 5 Mb were excluded. We also removed candidates within the pseudoautosomal regions (PARs) of chromosomes X and Y, as defined in the hg19/GRCh37 genome assembly.TbA: CNV validation based on FODTo refine CNV calls, we first defined a cell-specific diploid FOD distribution using the five longest segments with predicted copy number 2 (assuming no whole-genome doubling), from which the mean and standard deviation were derived. FOD z-scores for each CNV candidate were then computed against this baseline (Supplementary Information). Calls were retained if the FOD z-score exceeded 1.96 for gains (CN > 2) or fell below −1.96 for losses (CN PTA: CNV calling based on read depth and BAFMatched bulk whole-genome sequencing data were processed using the same alignment and preprocessing pipeline. For PTA data, HiScanner34 was applied in BAF+RDR mode (rdr_only=false) with the matched bulk as reference, with a bin size of 500 kb mappable positions, λ set to 64, max_wgd=1, ADO threshold of 0.2, and restrict_gamma=true. CNVs were reported as allelic copy numbers (CNA|CNB), from which total copy number was derived as the sum. Post-processing filters were applied sequentially: (1) bins called as non-diploid across cells from three or more donors were treated as recurrent artefacts and reset to CN = 2; (2) CNV segments spanning less than 5 Mb were excluded; (3) putative homozygous deletion calls (allelic pattern X∣0) inconsistent with the observed VAF (VAF > 0.1) were removed; (4) extreme copy number states (CN > 4) were excluded; and (5) gains (CN > 2) with VAFestimate Copy number signature analysisTo identify copy number signature exposures at the single-cell level, we used SigProfilerAssignment (v.0.1.4) in conjunction with SigProfilerMatrixGenerator (v.1.2.26; refs. 56,57). First, we generated a 48-channel copy number summary using the channels defined by COSMIC. Given the lack of allele specificity in our single-cell copy number calls, we developed a conversion scheme to transform our integer copy number states into pseudo-allele-specific states. Specifically, we mapped copy numbers into allelic states as follows: 0 to 0|0, 1 to 1|0, 2 to 1|1, 3 to 2|1, 4 to 3|1, and so on, up to a maximum copy number of 11 (10|1). This approach enabled construction of 48-channel copy-number summary vectors compatible with existing signature-analysis frameworks. Next, we linearly decomposed the copy number summary into copy number signature exposures. SigProfilerAssignment was run using default parameters, using the COSMIC v.3.4 pan-cancer signature panel (CN1–22, also referred to in this study as SigCN1–22; ref. 36) as the reference. Before signature clustering, we excluded three cells for which CNV calling produced no callable output and cells for which SigProfilerAssignment could not produce a valid signature decomposition.To cluster cells by signature exposure profile, we first normalized the signature exposure vector of each cell to sum to 1, then added a small amount of Gaussian noise (σ = 10−4) to break degeneracies in the sparse exposure space. For TbA cells passing quality filters (MAPD k-nearest-neighbour graph (k = 30, cosine distance) on these perturbed exposures using Scanpy (v.1.11.0; ref. 55), computed a two-dimensional UMAP embedding, and applied the Leiden community detection algorithm (resolution = 0.1) to identify four discrete groups of cells: cluster 0 (diploid predominant), cluster 1 (gain predominant), cluster 2 (focal loss predominant) and cluster 3 (broad loss predominant). PTA cells (after the same noise perturbation) were subsequently assigned to clusters by mapping onto the TbA reference embedding using the ingest function of Scanpy, which propagates cluster labels through k-nearest-neighbour label transfer in the reference UMAP space. To assess the robustness of the signature decomposition, we computed for each cell the cosine similarity between its observed copy number profile and the profile reconstructed from its assigned signature weights. A cosine similarity close to 1 indicates that the assigned signatures faithfully recapitulate the copy number landscape of the cell, whereas lower values flag cells whose profiles are poorly explained by the reference panel. Cosine similarities were aggregated and compared across clusters, serving as a per-cluster measure of decomposition confidence.Breakpoint distribution analysisTo investigate the relationship between CNV breakpoints and genomic features, we analysed replication timing and DNA double-strand break (DSB) distributions. For replication timing analysis, we used the replication timing track for the BG02ES (human embryonic stem cell) line generated by the ENCODE project58,59 downloaded from the UCSC Genome Browser60. The average replication timing signal was calculated for each 500 kb mappable bin and normalized between 0 and 1. For DSB analysis, BLISS coordinates were downloaded from FigShare (https://doi.org/10.6084/m9.figshare.18530531.v2)61, and DSB BED files for each replicate were intersected with the 500 kb genomic bins to quantify DSB signals, which were then normalized to the per-track maximum.To establish a background distribution for CNV enrichment analysis, 10,000 random CNVs were generated for each original CNV, matching their length distribution within autosomes (chr1–22). Peri-centromeric and telomeric regions were excluded to minimize biases, with centromeric regions defined by extending UCSC cytoband coordinates by ±5 Mb and telomeric regions defined as the first and last 5 Mb of each chromosome. Replication timing data were divided into deciles representing quantiles across the genome. For each decile, the observed-to-expected (O/E) ratio was calculated by normalizing observed breakpoints (that is, start and end sites of a CNV) to the expected number based on random simulations. The correlation between replication timing O/E ratios for CNVs from different categories (age group, signature group) was calculated, and statistical significance was evaluated using Pearson correlation tests.To analyse the relationship between CNV breakpoints and DSBs, we compared CNV breakpoint distributions to BLISS data from six neural cell preparations: two biological replicates each from human neuroepithelial stem cells, NPCs and postmitotic neurons (NEUs). BLISS read counts were summed into 500 kb genomic bins and normalized to the per-track maximum. DSB hotspots were defined independently for each track as bins exceeding the 99th percentile of normalized signal; a bin was classified as a hotspot if it met this threshold in any of the six tracks. For each CNV, 10 length-matched random CNVs were generated by placing segments of identical length at uniformly sampled positions across autosomes, and breakpoint positions (start and end of each CNV) were extracted for both observed and simulated sets. O/E ratios were calculated per individual as the count of breakpoints overlapping hotspot or non-hotspot bins divided by the mean count across the 10 simulated replicates. Neuroepithelial stem cell replicates were excluded from the primary analysis because of their strong correlation with replication timing (Pearson r = 0.50–0.53, FDR-corrected P r = 0.17–0.23) and NEU (r = 0.15–0.24); NPC and NEU tracks were therefore used for the aggregated hotspot analysis. O/E ratios were capped at 4 for visualization. Statistical significance was assessed using paired t-tests comparing the hotspot O/E ratio of each individual with their non-hotspot O/E ratio (n = 15 pairs each for gains and losses), with P values adjusted for multiple testing using the Benjamini–Hochberg FDR procedure.CNV breakpoints were further stratified by recurrence. A breakpoint was classified as recurrent if the same genomic position was observed in more than one donor; breakpoints detected in more than two donors were used for full-CNV recurrence analysis. A breakpoint was classified as private if it was detected in exactly one cell from a single donor. Within-donor clonal breakpoints—shared by two or more cells from the same individual—were analysed separately. BLISS hotspot enrichment analysis was performed independently for each recurrence class.CNV gene set enrichment analysisTo identify hot and cold spots in diploid-predominant neurons, we adapted an approach similar to that in ref. 18. We defined genomic bins of 5 Mb by aggregating the existing 500 kb bin coordinate set that was used in CNV calling, with any remaining small segments at the end of the chromosome discarded. Each 5 Mb bin was assigned an empirical overlap count, representing the number of overlaps with observed CNVs. We then performed 10,000 permutations, as outlined in the previous section, preserving the CNV sizes. For each permutation, the number of overlaps with random CNVs was also calculated, allowing us to generate a null distribution for each bin. Bins were classified as hotspots or coldspots using unadjusted one-sided empirical P values, with upper-tail P P 62 (pLI > 0.9), SFARI63 autism genes, neurodevelopmental disorder genes64, schizophrenia64 genes, postsynaptic genes65 and cancer driver genes from COSMIC66, and epilepsy genes (OMIM67 symptom search ‘epilepsy’). To assess gene set enrichment within hotspots and coldspots, hypergeometric tests were performed to evaluate over-enrichment of each gene set within these regions compared with all genomic bins.Single-nuclei multiomic ATAC/RNA library preparationNuclei were isolated without dissection as described in section ‘Isolation of neuronal nuclei’ and permeabilized after flow cytometry according to the Nuclei Isolation from Complex Tissues for Single Cell Multiome ATAC + Gene Expression Sequencing protocol from 10× Genomics (Demonstrated Protocol, CG000375). Chromium Next GEM Single Cell Multiome ATAC + Gene Expression libraries were prepared according to the manual with adaptations to also enable detection of complementary DNA strands (duplexes)68. Duplex information was not used for this study.Cell type and CNV annotation of multiomic snATAC and snRNA dataThe multiomic dataset consists of snRNA-seq and snATAC-seq data from prenatal and postnatal brain samples, produced following the 10× Genomics protocol. Data were processed using the 10× Genomics Cell Ranger ARC pipeline (v.2.0.1). Subsequently, the dataset was analysed with Seurat (v.5.0.1; ref. 69) and Signac (v.1.12.0; ref. 70) R packages.QC filtering of nuclei from snATAC-seqFor both prenatal and postnatal brain samples, nuclei were filtered on several metrics. First, multiplet detection was performed on snATAC-seq data using AMULET, after which cells classified as multiplets were filtered out. Next, filtering thresholds for several metrics were determined by visual inspection of rank plots of each metric for all samples. For total RNA molecule count, RNA feature count (number of genes), and total ATAC fragment count, a lower threshold of 1,000 was used to remove empty background nuclei, and no upper threshold was used given multiplet removal had already been performed. Finally, for both nucleosome signal and transcription start site (TSS) enrichment values from snATAC-seq, nuclei were removed with lower (0.32 for nucleosome signal and 1.75 for TSS enrichment) and upper (1.59 for nucleosome signal and 6.80 for TSS enrichment) thresholds determined by the rank plots. Moreover, a lower threshold of 0.5 was determined from the rank plots for the fraction of reads in peaks.Cell-type annotation based on snRNA-seqCell-type annotation was performed by mapping each sample to a reference sample. We chose separate references for prenatal and postnatal because of the differences in cell-type composition between prenatal and postnatal cells. For the prenatal samples, H29180_F4 was chosen as reference, and for the postnatal samples, the downsampled reference was used from a previously published dataset71. The prenatal reference was manually annotated by first performing unsupervised clustering on normalized, principal component analysis-transformed RNA count data. The resulting clusters were then manually annotated using known marker genes for different cell types. A UMAP (uniform manifold approximation and projection) representation was generated for the prenatal reference. For the postnatal reference, the pre-existing annotations from ref. 71, based on clustering performed in the paper, were used. The remaining non-reference ‘query’ samples were then projected onto the reference UMAPs using the MapQuery function of Seurat. The MapQuery function transferred the manually annotated cell-type labels from reference to query samples, resulting in complete cell-type annotation for all samples in the multiome dataset.CNV inference from snATAC-seqCNV calling was performed using the R package epiAneufinder (v.1.1.3; ref. 72), with a window size of 100 kb and a minimum fragment number of 20,000 per sample. The parameter minCNVsize was set to 0, meaning that the minimum callable CNV size was one bin, that is, 100 kb. After CNV calling was complete, the CNV calls were combined across all samples by taking the intersection of bins in all samples. For the resulting bins, for each cell across all samples, the percentage of bins with a CNV (either gain or loss) was calculated. The resulting distribution of the percentage of bins affected was used to determine a cutoff for our definition of CNV Neurons or CNV Cells. Tukey’s rule was applied for identifying outliers, such that any cell with a percentage of bins greater than the third quartile plus 1.5 times the interquartile range (Q3 + 1.5IQR) was classified as an outlier. These outliers were then classified as CNV neurons, and non-outliers as non-CNV neurons.DGM analysisWe analysed the snRNA-seq part of single-cell multiome data separately for prenatal samples (n = 3; UMB4827, UMB1932 and 29180_F4) and postnatal samples (n = 5; APGEQ, UMB1790, UMB4643, UMB1278 and UMB1864) because of the differences in cell-type composition. For each age group, sample-specific Seurat objects were merged into a combined dataset while preserving sample identifiers. Gene expression values were log-normalized with a scale factor of 10,000. To assess pathway activation, we calculated module scores using gene sets from the Molecular Signatures Database73, including Gene Ontology Cellular Component terms and curated gene sets (C2). For C2, we only included gene sets that matched the keywords ‘senescence’, ‘apoptosis’, ‘autophagy’, and ‘caspase’—hereafter referred to as ‘cell elimination C2 pathways’. Module scores were computed using the AddModuleScore function of Seurat with default parameters, which calculates the average expression of each gene set in a given cell after controlling for gene expression program complexity by subtracting the aggregated expression of randomly selected control gene sets with similar average expression levels.To test for associations between gene module scores and CNV neuron status, we implemented a stratified permutation testing approach. For each gene module, we calculated the observed difference in mean module scores between CNV and non-CNV cells (hereafter referred to as ‘absolute observed difference’). Statistical significance was assessed through 10,000 permutations, in which CNV neuron labels were randomly shuffled within each stratum defined by sample and cell type, preserving the underlying data structure. P values were calculated as the proportion of permuted differences that exceeded the observed difference in absolute value and were adjusted for multiple testing using Bonferroni correction. To avoid confounding by cell-type-specific gene expression differences, this analysis was performed separately for each major cell type (excitatory neurons, inhibitory neurons, progenitors and other cell types) and developmental stage (prenatal and postnatal).Clustering of CNV neurons based on DGM and cell death pathwaysTo identify distinct subpopulations of CNV neurons based on their pathway activation patterns, we performed hierarchical clustering on the normalized gene module scores from a selected set of pathways. As clustering on all available pathways could introduce noise from non-informative or weakly differential pathways, we focused on the top 10 most differentially activated gene modules from Gene Ontology Cellular Component and cell elimination C2 pathways, ranked by the absolute observed difference (see the previous section for definition of ‘absolute observed difference’ and the selection of cell elimination C2 pathways). The module score matrix was normalized by subtracting the mean and dividing by the standard deviation across cells. Hierarchical clustering was performed using average linkage with the cosine distance metric. We determined the optimal cluster number (n = 2) through visual inspection of the dendrogram and cluster stability.To visualize pathway activation patterns, we generated heatmaps using the seaborn clustermap function with row clustering enabled for gene sets and column clustering for cells. The resulting visualization showed two main clusters of cells with distinct pathway activation patterns. Cell cluster assignments were validated through multiple approaches, including assessment of cluster sizes (cluster 1: 1,667 cells; cluster 2: 2,487 cells) and evaluation of enrichment for CNV neurons using Fisher’s exact test. The association between cluster membership, and CNV neuron status was quantified using odds ratios with 95% confidence intervals. Multiple testing correction was performed using the Bonferroni procedure.Ethics statementHuman materials were used in this study according to state and federal regulations with approval from the Boston Children’s Hospital and Mass General Brigham Institutional Review Boards. Postmortem human brain tissue was obtained from the NIH NeuroBioBank (University of Maryland Brain and Tissue Bank), the Autism BrainNet Collection, the University of Washington Birth Defects Research Laboratory (BDRL) and the Mayo Clinic Brain Bank. Informed consent for research was obtained from donors or their next of kin by the contributing biobanks as approved by their respective Institutional Review Boards. The NPC line was derived from a previously described control iPSC line generated from commercially available human fibroblasts54.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
