Isolation of neuronal nuclei
Human prenatal brain tissue was sourced from the NIH NeuroBioBank. Sections, about 40 μm thick, were cut from the frontoparietal lobes using a cryostat. The cortical layer was carefully dissected with a microscalpel under a microscope, and the nuclei were separated through gentle dounce homogenization in a lysis buffer comprised of 10 mM Tris-HCl (pH 8.0), 250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 0.1% Triton-X, 1× cOmplete Mini Protease, and 1 mM DTT. After centrifuging for 10 minutes at 900g, the nuclei were pelleted. The supernatant was discarded, and the nuclei were resuspended in sorting buffer containing 1× PBS (pH 7.4), 0.8% BSA, and 0.1 μg/ml DAPI. Flow cytometry was conducted, focusing on FSC, SSC, and DAPI content to sort single nuclei. For postnatal brain nuclei, staining involved both DAPI and Alexa Fluor 488 labeled anti-NeuN at a dilution of 1:1,000 for one hour at 4°C, using a product from Sigma Aldrich.
Single-nuclei whole genome library preparation
TbA has been utilized for scWGA in various contexts under different labels, mainly involving three steps: (1) nucleosome depletion; (2) Tn5 tagmentation and direct adapter addition; and (3) library amplification. Single nuclei were lysed in microwells with 2 μl of lysis buffer (20 mM NaCl, 20 mM Tris-HCl (pH 8.0), 0.15% Triton-X, 1 mM EDTA, 250 mM DTT, and 0.05 μl Thermolabile Proteinase K) at 30°C for one hour, followed by a 10-minute inactivation at 55°C. This process improved genome access by digesting histones in nucleosomes. Following manufacturer’s instructions, Tn5 was prepared, using 1 μl diluted 1,000-fold for each nucleus. Tagmentation was then performed, leading to single-cell genome libraries through 16 cycles of PCR with NexteraXT v.2 indexes. Rows of 12 single cells were combined for cleanup using the Zymo Clean and Concentrator Kit. Each row was eluted in 42 μl of Zymo elution buffer for AMPure XP size selection at 0.8×. Four rows, totaling 48 cells, were pooled for sequencing at an equimolar concentration, executed on a NovaSeq 6000 platform, targeting a library size of 0.5 Gb per single-cell library. Illumina reads were aligned to the human reference genome with decoy sequence GRCh37d5 using bwa-mem, followed by duplicate marking and base quality recalibration per GATK best practices.
Single-cell library quality assessment
To analyze only high-quality single-cell genomes, we initially discarded cells with fewer than 10 million total covered bases, indicated by the total column from GATK DepthOfCoverage. We then used a dataset from an NPC line for transfer learning to differentiate between high- and low-quality sequencing profiles. Constructing a benchmark of sample quality was based on staurosporine treatment for 24 hours to induce apoptosis, creating a labelled TbA dataset that included (1) high-quality viable cells; (2) cells marked for apoptosis by Annexin V staining while maintaining intact nuclear envelopes; and (3) cells with compromised nuclear integrity indicated by propidium iodide. The NPC line was differentiated as previously noted with an iPSC derived from well-characterized human fibroblasts. NPC identity was confirmed through immunohistochemistry using PAX6 and Nestin antibodies.
To categorize cells computationally, we normalized read counts within the NPC dataset via the HiScanner’s read-depth normalization function, creating a cell-by-bin matrix. Then we scaled the read counts using Scanpy, followed by principal component analysis for dimensionality reduction. Similarities between cells were assessed by constructing a neighbourhood graph, and clustering was completed through the Leiden algorithm with a resolution parameter set to 0.2. After visual inspections and correlating cluster read count profiles with known replication timing tracks, we identified the Leiden clusters within NPC data. Thus, we labeled clusters as high quality, atypical, and poor quality. The scanpy.tl.ingest function was subsequently utilized to project these labels and embeddings from NPC to human brain single-cell TbA data with default parameters, focusing on cells identified as high quality.
Single-cell CNV calling pipeline
Post filtering at both cellular and sample levels, CNVs were identified in single-cell genomes. For TbA data, CNV calling relied on read depth and was validated using FOD. In PT data, CNV calling combined read depth with B-allele frequency (BAF).
TbA: CNV calling based on read depth
HiScanner was employed in rdr_only mode with tailored parameters for our analysis. A bin size of 500 kb was used, referencing hg19/GRCh37. Key settings included λ at 64 and enabling a ploidy restriction for individual cells. The multisample option was disabled to analyze each cell separately. To ensure credibility in detected CNVs, a rigorous filtering approach was adopted: CNV candidates with P ≥ 0.05 or spanning less than 5 Mb were excluded, alongside candidates within pseudoautosomal regions of chromosomes X and Y as defined in the referenced genome assembly.
TbA: CNV validation based on FOD
For refining CNV calls, a cell-specific diploid FOD distribution was defined using the five longest segments with a predicted copy number of 2. From this, the mean and standard deviation were calculated. FOD z-scores for CNV candidates were computed relative to this baseline. Retaining calls required the FOD z-score to exceed 1.96 for gains or fall below -1.96 for losses.
PTA: CNV calling based on read depth and BAF
Matched bulk whole-genome sequencing data underwent the same processing pipeline. PTA utilized HiScanner in BAF+RDR mode with matching bulk as reference, setting a bin size of 500 kb, λ at 64, max_wgd at 1, and an ADO threshold of 0.2. CNVs were reported in terms of allelic copy numbers, from which total copy number was inferred. A series of post-processing filters was sequentially applied.
Copy number signature analysis
To pinpoint copy number signature exposures at the single-cell level, we employed SigProfilerAssignment alongside SigProfilerMatrixGenerator. We created a 48-channel copy number summary, mapping integer copy number states into pseudo-allele-specific states. This allowed us to form vectors that align with existing signature-analysis frameworks. Next, we decomposed the copy number summary into exposures using SigProfilerAssignment with default parameters from the COSMIC v.3.4 panel. Before signature clustering, certain cells were excluded based on invalid outputs from CNV calling and signature decomposition.
To cluster cells by their signature exposure profiles, we normalized each cell’s signature exposure vector before adding a small amount of Gaussian noise. For TbA cells meeting quality filters, we constructed a k-nearest-neighbour graph and identified clusters based on pathway activation patterns.
Breakpoint distribution analysis
We examined the connection between CNV breakpoints and genomic features by analyzing replication timing and DNA double-strand break distributions. We used data from the ENCODE project for replication timing, with average signals calculated for mappable bins. For DSB analysis, we downloaded BLISS coordinates and quantified DSB signals within those bins.
To create a background distribution for CNV enrichment analysis, we generated 10,000 random CNVs matching the length distribution of original CNVs. The observed-to-expected ratios were calculated, and statistical significance was determined through correlation tests.
When analyzing the relationship between CNV breakpoints and DSBs, we compared distributions across various neural cell preparations, generating random CNVs to extract breakpoint positions. Observed and simulated breakpoint data were compared to identify hotspots.
CNV gene set enrichment analysis
To determine hot and cold spots in diploid-dominant neurons, we implemented a method similar to previous research, defining genomic bins of 5 Mb and assigning empirical overlap counts related to observed CNVs. A series of permutations was conducted, allowing us to categorize bins as hotspots or coldspots based on statistical testing.
Single-nuclei multiomic ATAC/RNA library preparation
Nuclei were isolated without dissection, following the protocol outlined for multiome ATAC and RNA sequencing from 10× Genomics. Libraries were prepared per the manual, with adaptations included to detect complementary DNA strands, although duplex information wasn’t utilized in this study.
Cell type and CNV annotation of multiomic snATAC and snRNA data
The resultant multiomic dataset integrated snRNA-seq and snATAC-seq data from both prenatal and postnatal samples, processed using the Cell Ranger ARC pipeline. Analysis followed through Seurat and Signac R packages.
QC filtering of nuclei from snATAC-seq
Nuclei from both sample types were filtered based on several metrics, employing AMULET for multiplet detection and defining filtering thresholds through rank plots. Empty background nuclei were removed by setting a lower threshold, while upper thresholds for nucleosome signal and TSS enrichment were determined similarly.
Cell-type annotation based on snRNA-seq
Cell-type annotation involved mapping each sample to a reference, with distinct references utilized for prenatal and postnatal samples. Manual curation based on clustering was applied for prenatal samples, while pre-existing annotations served for postnatal. Non-reference samples were labeled accordingly, achieving comprehensive cell-type annotation.
CNV inference from snATAC-seq
CNV calls were executed using the epiAneufinder package for single-cell ATAC-seq data, and the calls across all samples were aggregated by taking the intersections. A percentage calculation of bins with CNVs was utilized to define CNV neurons versus non-CNV neurons.
DGM analysis
We separately analyzed snRNA-seq data from prenatal and postnatal samples due to differing cell-type compositions. Sample-specific Seurat objects merged into a comprehensive dataset while maintaining identifiers. Log-normalization was applied before calculating module scores based on gene sets relevant to cellular pathways.
For testing associations between gene module scores and CNV neuron statuses, a stratified permutation testing method was used. Each gene module’s observed mean differences were compared against permuted samples to evaluate their significance, while controlling for potential confounding factors.
Clustering of CNV neurons based on DGM and cell death pathways
We sought to classify CNV neurons based on pathway activation patterns through hierarchical clustering on normalized gene module scores, focusing on the top differential pathways. Heatmaps were generated to visualize these patterns, validating cluster assignments through various assessments.
Ethics statement
This study’s human materials were utilized aligning with state and federal regulations, having approval from both Boston Children’s Hospital and Mass General Brigham Institutional Review Boards. Postmortem brain tissue was acquired ethically from various sources, with informed consent obtained as approved by their respective boards.
Reporting summary
Additional information regarding the research design can be found in the linked summary.
