Using the linear references from the pangenome to discover missing autism variants

Domain: neuroscience
Published: 2026-01-23
Editorial outcome: minor_revision

Abstract
To better understand large-effect pathogenic variation associated with autism, we generated long-read sequencing (LRS) data to construct phased and near-complete genome assemblies (average contig N50 = 43 Mbp, QV = 56) for 189 individuals from 51 families with unsolved cases. We applied read- and assembly-based strategies to facilitate comprehensive characterization of de novo mutations, structural variants (SVs), and DNA methylation. Using LRS pangenome controls, we efficiently filtered >97% of common SVs exclusive to 87 offspring. We find no evidence of increased autosomal SV burden for probands when compared to unaffected siblings yet observe a suggestive trend toward an increased SV burden on the X chromosome among affected females. We establish a workflow to prioritize potential pathogenic variants by integrating autism risk genes and putative noncoding regulatory elements defined from ATAC-seq and CUT&Tag data from the developing cortex. In total, we identified three pathogenic variants in TBL1XR1, MECP2, and SYNGAP1, as well as nine candidate de novo and biallelic inherited homozygous SVs, most of which were missed by short-read sequencing. Our work highlights the potential of phased genomes to discover complex more pathogenic mutations and the power of the pangenome to restrict the focus on an increasingly smaller number of SVs for clinical evaluation.

Introduction
Autism is a class of neurodevelopmental disorders (NDDs) characterized by challenges in social interaction, communication, and repetitive behaviors, with symptoms and severity varying widely among individuals. Currently, the median prevalence of autism globally is approximately 1%, with a male-to-female ratio of about 4 to 1. Genetic investigations into rare variants have typically focused on gene-disruptive de novo mutations (DNMs), rare inherited variants, and copy number variants (CNVs) discovered through meta-analyses of short-read sequencing (SRS) and earlier array CGH. These variants account for an estimated 20% of autism cases and have led to the discovery of ~1000 risk genes and numerous CNVs associated with NDDs. The underlying mutations for the remaining cases, including structural variants (SVs) mapping to repetitive regions, remain poorly understood.

SVs, including deletions (DEL), insertions (INS), inversions (INV), translocations, large-scale CNVs, and other complex rearrangements, are defined as affecting ≥50 base pairs of DNA. They have been shown to have larger effects because they can disrupt coding or noncoding regulatory regions, alongside protein-coding genes, thereby playing a critical role in gene regulation and human disease. However, many SVs occur in technically and methodologically challenging regions, particularly repetitive sequences, making them difficult to detect and completely characterize using conventional SRS approaches associated with whole-genome sequencing (WGS) or whole-exome sequencing (WES). Similarly, a subset of smaller single-nucleotide variants (SNVs) and indels can be missed by SRS because of their association with low-complexity and repetitive DNA.

Long-read sequencing (LRS) data (15-30 kbp) significantly enhances the sensitivity of variant detection, especially in repetitive DNA regions. Recent studies have revealed that LRS data provide access to ~91% of the human genome, substantially increasing DNM discovery by ~30% and SV discovery by over 47% when compared to SRS datasets. Consequently, LRS has been increasingly applied to a variety of unsolved patients and disorders to enhance pathogenic variant discovery, although most studies to date have involved relatively modest cohorts focused almost entirely on read-based discovery. For example, Hiatt et al. reanalyzed 10 NDD families and 86 probands using Pacific Biosciences (PacBio) high-fidelity (HiFi) LRS and found an additional yield of 7.3% beyond SRS, mainly in the coding regions. Sanchis-Juan et al. applied Oxford Nanopore Technologies (ONT) LRS to complement SRS in four probands. Moreover, the ability to accurately call methylated CpGs, especially from ONT LRS, has the added benefit of simultaneously discovering potential imprinting defects.

Beyond read-based variant discovery, the combination of LRS technologies (ONT and PacBio) has facilitated the construction of near-complete telomere-to-telomere (T2T) genome assemblies as part of the Human Genome Structural Variation Consortium (HGSVC) and the Human Pangenome Reference Consortium (HPRC). These consortia recently made hundreds of diverse human genomes publicly available. This resource is potentially valuable to the clinical genetics community because variant discovery is more complete, providing a control to assess the frequency of variants in regions typically unassayable by SRS and therefore absent or unreliable in associated databases such as gnomAD. Moreover, assembly-based comparisons between offspring and parental genomes have been shown to further increase the power to discover DNM by essentially eliminating reference biases. Notwithstanding, the number of samples still remains modest, and much larger sample sizes will be required to understand the full spectrum, especially for those with lower minor allele frequency (MAF).

Using LRS assembly approaches, we sought to construct reference-free genome assemblies comparable to the HGSVC and HPRC controls for all members of autism families. In this study, we present our initial LRS and assembly resource of 189 individuals from 51 unsolved autism families where no pathogenic variant was previously identified in the proband via conventional Illumina whole-genome, whole-exome, or gene panel testing with SRS. To build reference resources, we constructed near-complete genomes for each individual and assigned a workflow to systematically identify variants from high-quality assemblies; we then compared them with HPRC and HGSVC population controls to discover and validate de novo and rare variants (<0.5%) as further candidates for autism. Importantly, constructing genomes comparable to pangenome references allowed us to dramatically reduce the variant search space for SVs, highlighting the increasing utility of the pangenome for disease variant discovery.

Results
We focused on the sequence and assembly of genomes from 51 unsolved simplex autism families. The set included 46 families (174 individuals) with idiopathic autism from the Simons Simplex Collection (SSC) and Study of Autism Genetics Exploration (SAGE) and five families with a diagnosis of Rett syndrome (15 individuals). The Rett families had been previously screened using either gene panels and/or WES with no MECP2 pathogenic mutation reported by clinical testing labs after multiple attempts. Similarly, the 46 families were part of large-scale CNV and WGS initiatives over the last decade, where no pathogenic variant had been reported by multiple groups, including our own. For 36 families, there was also an unaffected sibling serving as a genetic control (17 of these quads were sex-matched).

We developed an LRS workflow to enhance variant discovery using reference-level quality genomes with a particular emphasis on characterizing previously undetected SVs and DNMs. First, we sequenced all 189 genomes using PacBio HiFi sequencing technology from peripheral blood (n = 139), cell lines (n = 16), or a mixture of both (n = 34) when DNA from peripheral lymphocytes was limited. Per sample, we generated an average of 36-fold sequence coverage with an overall N50 read length of 19 kbp after extensive quality control. Using parental Illumina reads and hifiasm, we generated haplotype-resolved genome assemblies (i.e., phased by parent-of-origin) for each of the 87 offspring. The resulting assemblies are highly contiguous (average contig N50 of 43 Mbp) and highly accurate (quality value (QV) = 56). The assembly quality is comparable to that of control pangenomes (with mean QV of 56) from the HPRC and HGSVC.

For each family, SNVs and small indel call sets were generated by the GATK and DeepVariant callers using T2T-CHM13v2.0-aligned HiFi reads mainly from blood samples (n = 172). Putative DNMs were further validated by confirming their presence in ONT and Illumina data, as previously described. We discovered on average 96 DNMs per child (n = 78), with approximately 83% predicted to be germline and the remainder arising postzygotically. SVs in each individual were identified from the phased assemblies using PAV with GRCh38 as the reference genome. GRCh38 was selected due to the broad availability of annotation resources mapped to this assembly, many of which were integrated into our analyses to facilitate comprehensive functional characterization of putative disruptive variants. Each SV was considered validated if supported by at least one of the alignment-based SV callers, either PacBio structural variant caller (PBSV), Sniffles, or both. We aggregated the validated SVs from all 189 study samples and compared them with SVs observed in the 108 population controls from the HPRC and HGSVC via Truvari, focusing on rare SVs exclusive to autism families. As a result of subsequent LRS of additional 1000 Genomes Project (1KGP) samples recently released as part of the HPRC and other efforts, we expanded our SV call set from controls in a staged manner, allowing us to access variants of reduced MAF.

Based on PAV analysis of each child’s assembly, we identified an average of 27,576 SVs per diploid genome of which 20,716 SVs (95% Confidence Intervals = 9) were also supported by alignment-based methods. We note that two Rett-like samples where DNA was limited (HYZ204_p1 and HYZ207_p1) have lower contig N50 and lower assembly QV, resulting in fewer validated SVs. The average Mendelian concordance rate for SVs across 87 trios is 90.4% with many of the discrepant alleles associated with multiallelic variants, such as variable number tandem repeats (VNTRs) where both sequence and alignment artefacts confound variant calling. This high-confidence SV call set, thus, affects 11,074,300 bp (95% Confidence Intervals = 4273) per sample (0.4% of the genome). We applied the same SV discovery approach to 108 pangenome controls resulting in on average 24,341 ± 36 validated SVs per diploid genome. This approach improved the sensitivity of control SV detection and facilitated more effective filtering of SVs in the sample set. Both pangenome controls and study samples are from diverse superpopulations. And, as expected, genomically diverse African samples had higher numbers of SVs than non-Africans. We integrated the two call sets for a total of 271,375 non-redundant SVs and observed an expected SV size distribution with modes at 300 bp and 6 kbp, corresponding to Alu and LINE retrotransposition events, respectively.

After filtering with the pangenome we identified a total of 33,548 nonredundant SVs (57,716 genotyped SVs) that were exclusively observed in the 87 children (51 probands and 36 unaffected siblings). At the individual sample level, we effectively filtered ~97% of SVs per child, resulting in approximately 663 rare SVs per child for further consideration (affected vs. unaffected, Z = -0.23, p = 0.82, two-sided Mann-Whitney U Test). After testing for Mendelian inheritance, the autism set was further reduced to 25,272 nonredundant SVs. We then classified rare SVs into six categories, including autosomal heterozygous (n = 33,701), autosomal homozygous (n = 6639), X chromosome heterozygous (n = 1129), X chromosome homozygous (n = 277), and hemizygous events from males on the X (n = 565) and Y (n = 141) chromosomes. We evaluated the remaining low-confidence SV calls through a suite of tools that test for transmission and de novo variants. We validated a total of 36 de novo SVs in 51 probands (n = 21) and 36 siblings (n = 15) that were absent from 108 controls (affected vs. unaffected, p = 1, χ² test).

Because the majority of rare SVs (96%) map to noncoding DNA, we annotated all SVs for regulatory potential as well as association with known NDD genes. To define putative regulatory elements (REG), we integrated published datasets from ENCODE with regulatory sequences predicted from single-cell and bulk ATAC-seq and CUT&Tag experiments performed on developmentally staged material (16–24 weeks) from the cerebral cortex. The additional annotation allowed us to identify another 6171 SVs associated with brain-derived regulatory regions (brainREG) beyond those defined by ENCODE intersection and consequently a 45% (6171/13,773) increase in potential SVs affecting noncoding regulatory DNA. For NDD candidate genes, we primarily focused on those previously reported. We considered three classes of SVs, namely: de novo, homozygous, and private inherited SVs defined as those observed only once in the parental population. On autosomes, we identified 9474 private SVs in 51 probands and 7034 private SVs in 36 unaffected siblings. The rates of autosomal SVs do not differ significantly between probands and unaffected siblings across various categories (nominal p > 0.05, OR < 1, χ² test).

For X chromosome SVs, we analyzed males and females separately. Female probands exhibited a slight excess of X chromosome private SVs compared to the unaffected female siblings (nominal p = 0.29, OR = 1.16, adjusted p = 1, χ² test), particularly for deletions (nominal p = 0.09, OR = 1.49, adjusted p = 1, χ² test). We also considered 1417 private homozygous SVs that were inherited from two heterozygous parents on autosomes in only one family but were never observed as homozygous in any controls. Once again, we observed a slight enrichment in the number of these private homozygous SVs in probands relative to unaffected siblings (p > 0.05, χ² test).

In addition to SVs, we constructed nearly complete X and Y chromosomes by leveraging LRS and parental sequenced data. Excluding the pseudoautosomal region (PAR), centromere, and highly repetitive Yq12 heterochromatic regions, we estimate that on average 95% of the X chromosome and 70% of the Y chromosome can be aligned. The assemblies served two purposes: they validated DNMs on the sex chromosomes and allowed parental transmission to be fully assessed (i.e., paternal Y and maternal vs. paternal X chromosomes) without the use of a reference genome.

The analysis highlighted three hemizygous tandem repeat (TR) expansion outliers in male probands, each inherited from the mother in the intron of IL1RAPL1 (Interleukin 1 Receptor Accessory Protein Like 1, SFARI score 2 gene), the intron of F9 (Coagulation Factor IX), or the intergenic region between MXRA5 (Matrix Remodeling Associated 5) and SNORA48B (Small Nucleolar RNA, H/ACA Box 48B). These longer TR noncoding variants have only been observed in heterozygous states in females and variants of such lengths have yet to be observed in controls. Notably, the corresponding probands exhibited features commonly associated with neurodevelopmental delay: IQ scores of 18, 13, and 20, and calibrated severity scores (CSS) of 9, 6, with the CSS data for the third not available.

Because LRS data allow CpG methylation to be robustly called, we used the 889 CpG islands across the X chromosome to assess X chromosome inactivation (XCI) skewing in the blood of female probands and their unaffected sisters. The analysis revealed extreme examples of XCI skewing, including preferential inactivation of the maternally inherited X chromosome (see 11071), potentially consistent with the mother carrying a damaged X chromosome.

Using the LRS and assembly data, we comprehensively searched for both pathogenic and potential candidate variants that were missed by SRS analyses. In total, we identified three DNMs classified as pathogenic. Among the 46 idiopathic autism samples, we discovered a de novo stop-gain variant (G1124*) in SYNGAP1, a well-known autism-associated gene encoding a Ras GTPase-activating protein essential for synaptic function and cognitive development. This pathogenic DNM in 12237_p1 is clearly supported by both HiFi and ONT reads from the proband but was not reported in three prior SRS analyses of this family. An analysis of the Illumina sequence data, however, confirms the presence of the variant in a GC-rich region of the genome where a cluster of rare and additional false calls were present, likely resulting in this region being subsequently filtered during QC.

Two additional pathogenic variants were discovered by LRS among the five autism-diagnosed females with features reminiscent of Rett syndrome. This included an 874bp de novo DEL in MECP2 in HYZ207_p1, which effectively disrupts the last exon of the gene and introduces a premature stop codon, truncating the protein by approximately 140 amino acids. Similar deletions involving exon 4 have been reported in Rett patients. This pathogenic variant was previously missed in three rounds of clinical testing, including two gene panel sequencing tests through ARUP Laboratories and Quest Diagnostics, and one test of WES through Ambry Genetics. It was confidently identified in our LRS analysis and subsequently validated using all three sequencing platforms. We also identified a de novo missense mutation within TBL1XR1 classified by ClinVar as pathogenic (rs1057517933). TBL1XR1 encodes a transducin (beta)-like 1 X-linked receptor 1–that directly interacts with MECP2. This de novo TBL1XR1 (D370N) missense variant was previously reported in the same patient and two other cases in DECIPHER and was recently classified as a pathogenic variant. We also observed a corresponding decrease in methylation at this CpG site within the exon of TBL1XR1. As part of this analysis, we also note two DNMs called by both SRS and LRS mapping to promoters of genes strongly implicated in neuronal development (POGZ and DDX3X). POGZ encodes a zinc finger protein involved in chromatin remodeling and transcriptional regulation, while DDX3X, an ATP-dependent RNA helicase, plays a crucial role in RNA metabolism, translation regulation, and neuronal development. Both genes have been implicated in NDDs, including autism. These two DNMs have never been observed in gnomAD and given their critical location are candidates for functional testing using massively parallel reporter assay (MPRA) to determine if they significantly reduce expression levels.

Among de novo SVs, we identified several candidates of potential regulatory consequence. For example, we identified a 71bp de novo TR INS within the intron of CNTN3, a SFARI gene encoding Contactin 3, which mediates cell surface interactions during nervous system development and the outgrowth and guidance of axons and dendrites. This INS is predicted to disrupt HNRNPK transcription factor (TF) binding sites in 11201_p1; notably, HNRNPK also functions as an RNA-binding protein. We attempted to recall this SV using SRS-based callers, including Manta, Smoove, CNVnator, and Canvas, and re-genotyped it using Paragraph based on Illumina alignments. All SRS tools failed to detect this variant, likely due to its mapping within a CT-rich TR. We also identified a 73bp de novo TR INS predicted to interrupt regulatory regions of CPT1C and TF binding clusters in 14455_p1. CPT1C encodes carnitine palmitoyltransferase 1C, a neuron-specific protein located in the endoplasmic reticulum, and has been associated with spastic paraplegia with an emerging role in neuropsychiatric conditions. This de novo INS was missed by all SRS-based callers, likely due to high-GC content (63% within 100bp).

Finally, we evaluated rare biallelic inherited homozygous SVs for potential pathogenicity because of the unique capability of LRS to phase almost all variants. We identified six candidates SVs inherited biparentally and associated with SFARI risk genes or cortex-specific regulatory regions. None of these have been reported as homozygous in SRS genomic controls from gnomAD or were identified in pangenome controls. This set includes: a homozygous DEL overlapping an enhancer located in the 3′ untranslated region (UTR) of CLN8 (a SFARI score 2 gene) in proband 11611_p1; a homozygous DEL disrupting the enhancer of ARHGEF10 (a SFARI score 2 gene), encoding a Rho guanine nucleotide exchange factor (GEF), in 13414_p1; a homozygous cortex-specific cis-regulatory element DEL in individual 12651_p1 mapping to the 3’ UTR of LRPAP1, which encodes LDL receptor-related protein associated protein 1 and has been linked to dementia and late-onset Alzheimer’s disease; a homozygous DEL within an intronic enhancer of PREX1 (a SFARI score 2 gene), in proband 12861_p1; a 332bp homozygous INS disrupt TF binding in TBC1D5 (a SFARI score 2 gene) in 11918_p1, and a 500bp homozygous DEL mapping to the promoter region of both LMF2 and NCAPH2 (a SFARI score 3 gene) in individual 14350_p1. This promoter region is predicted to harbor regulatory activity based on multiple datasets, including chromatin accessibility and enhancer marks in the developing brain, suggesting potential cis-regulatory effects on the expression of one or both genes.

To assess the potential impact of candidate variants, we performed differential expression analysis using long-read RNA sequencing of lymphoblastoid cell lines generated with the PacBio Kinnex platform. Gene-level counts were obtained with IsoQuant and analyzed with DESeq2 for probands carrying genes with two promoter mutations and five homozygous deletions that were expressed in blood. Across all seven candidate loci, no significant reduction in gene expression was detected in blood-derived RNA samples from the target proband compared to the remaining samples. The absence of detectable expression changes may reflect tissue-specific regulatory effects, as these variants are more likely to exert functional consequences in brain tissues, which are directly relevant to the disorder etiology.

To enhance the power of pangenome filtering to effectively exclude common SVs and focus on a high-confidence pool of rare variants, we expanded the control cohort size from 108 to 285 and then to 569 individuals. The 285-control set contains 177 newly LRS and assembled samples from the HPRC (HPRC Y2), while the 569-control set includes an additional 284 publicly available 1KGP samples sequenced using ONT. Increasing the number of controls, especially samples of African origin, nearly doubles the number of SVs from 271,375 (108 controls) to 445,142 (569 controls) nonredundant SVs. We generated a population-level SV reference and compared the MAF distribution between the 108 controls and the largest dataset comprising 569 controls and 102 unrelated parents (n = 1342). Based on this comparison, we estimated the fraction of variants captured at different MAF thresholds. Notably, only 16.2% of SVs with MAF < 0.1% in the larger dataset were detected in the 108-control set, illustrating the limited sensitivity of the smaller reference for identifying ultra-rare variants. In an ancestry- and sex-matched analysis restricted to individuals of European ancestry (38 autism families and 73 controls), each child carried an average of 754 rare SVs. Focusing on rare variants absent from the largest 569-control set, the number of rare SVs concomitantly drops to 202 events per sample, corresponding to ~74% of rare variants being filtered out by the more diverse population controls. This demonstrates the substantial gain in sensitivity for removing common and low-frequency variants when using a larger, heterogeneous reference panel. As a result, 99% of common SVs are excluded per individual forming a more tractable and potentially biomedically relevant set of rare variants for downstream interpretation and enrichment analyses. Although still not statistically significant (p = 0.23, two-sided Mann-Whitney U test), we note that applying 569 controls results in a larger difference with respect to SV burden between probands and unaffected siblings (Z = 1.20). Of note, only one of our previously 13 proposed candidates (73bp de novo TR INS in CPT1C) was excluded and none of our pathogenic variants were excluded.

We also compared our LRS SV datasets with a comprehensive SRS SV call set (INS and DEL) generated by GATK-SV from 63,046 unrelated genomes, as represented in gnomAD v4.1. To ensure a fair comparison between SRS and LRS datasets that differ in resolution and breakpoint precision, we used relaxed matching criteria in Truvari bench: a minimum of 50% reciprocal overlap, 50% size similarity, a maximum breakpoint distance of 500 bp, and no sequence similarity requirement. We added the gnomAD SV IDs to the supplementary datasets. Within the 108-control set, 162,916 of 271,375 total nonredundant SVs were observed in the autism families. Of these, 77.6% (126,377 SVs) were present in the 108 LRS controls, yielding 663 rare SVs per child. In contrast, only 24.1% (39,234 SVs) overlapped with the SRS-based gnomAD SVs from 63,046 unrelated controls, corresponding to 14,110 rare SVs per child when considering only short-read controls. Thus, despite its much larger sample size, the short-read dataset can be considered an auxiliary resource for further filtering rare SVs beyond those already excluded by the LRS controls, underscoring the substantially greater power of LRS data for rare variant filtering. Specifically, incorporating the SRS dataset would have further excluded 0.71% (147/20,716 per child) and 0.27% (55/20,716 per child) of SVs from the 108- and 569-control sets, respectively. Given the inherent differences in breakpoint resolution and variant representation between LRS and SRS platforms, these overlaps reflect potentially inflated estimates. A reassuring sign of the stringency and completeness of our filtering strategy is that all 12 candidates identified were rare variants absent from both the LRS and SRS control datasets.

Methods
Samples selection and previous sequence characterization
Illumina WGS was previously applied to the blood DNA of the 174 individuals corresponding to 46 families with idiopathic autism from the SSC and SAGE. Potential pathogenic variants were screened from SRS based on the studies by Wilfert et al. and Fu et al., including (1) SNVs: de novo likely gene-disruptive (LGD), de novo missense, and rare inherited LGD variants in NDD-related genes or genes with pLI scores ≥ 0.9; (2) CNVs in the morbidity map or span of a gene with pLI score ≥ 0.9. No known genetic cause was identified in the 46 probands from SRS data. In addition, the 46 probands selected in this study do not exhibit exceptional polygenic risk scores among the 13,989 individuals examined and have known IQs ranging from 13 to 91. Similarly, five girls diagnosed with Rett-like syndrome from Baylor College of Medicine had no causal variants identified in MECP2 by prior gene panel testing or WES. The selected families were predominantly those with female probands due to the interest in discovering X chromosome variants, and the large-effect variants are more likely to be discovered. As shown in Fig. 1a, 17 quads with sex-matched offspring (12 female-female and 5 male-male quads), 19 quads with sex-mismatched offspring (16 female-male and 3 male-female quads), and 17 trios (15 female and 2 male trios) were included in this study. Briefly, the 51 probands (41 females and 10 males) were selected to represent cases that pose difficulties in pinpointing the cause of autism using Illumina SRS data.

LRS data generation and QC
We generated PacBio HiFi and ONT sequencing data at the University of Washington (UW) for 189 individuals. Illumina WGS for Rett-like trios was generated from blood DNA using the TruSeq library kit and sequenced on a NovaSeq with paired-end 150 bp reads at the Northwest Genomics Center. For the five Rett-like trios, DNA was extracted from blood using the Monarch HMW DNA Extraction Kit for Cells & Blood from NEB (T3050L) (n = 6) or the Qiagen Puregene Blood Core Kit (158023) (n = 9), following the manufacturer’s specifications. The whole-blood DNA from the SSC was extracted previously as part of that biobank. The cell line DNA was extracted from lymphoblastoid cell lines with either a modified Gentra Puregene (Qiagen) protocol when used for ONT sequencing or with the Monarch HMW kit (NEB T3050L) when used for PacBio sequencing. ONT libraries were constructed using the Ligation Sequencing Kit (ONT, LSK110 and LSK114) with modifications to the manufacturer’s protocol. The library was loaded onto a primed R9.4.1 or R10.4.1 flow cell (FLO-PRO002 or FLO-PRO114M) for sequencing on the PromethION, with two nuclease washes and reloads after 24 and 48 hours of sequencing.

PacBio HiFi data from family 14455 (n = 4) were published in Noyes et al. Data from three families (n = 11) were graciously generated by PacBio. Remaining individuals’ HiFi data were generated from blood or cell line HMW DNA according to the manufacturer’s recommendations. At all steps, quantification was performed with Qubit dsDNA HS (Thermo Fisher Scientific, Q32854) measured on DS-11 FX (Denovix) with the size distribution checked using FEMTO Pulse (Agilent, M5330AA and FP-1002-0275.) The samples’ incoming size distribution determined shearing conditions, either no shear (n = 12), or sheared with the Megaruptor 3 (Hologic Diagenode, B06010003 & E07010003) system using one (n = 36) or two (n = 130) sequential runs to target a peak size of ~20 kbp. After shearing, the DNA were used to generate PacBio HiFi libraries using the Express Template Prep Kit v2 (n = 12, PacBio, 100-938-900) or SMRTbell prep kit 3.0 (n = 166, PacBio, 102-182-700). Size selection was performed with Pippin HT using a high-pass cut-off between 9-17 kbp based on shear size (Sage Science, HTP0001 and HPE7510). Libraries were sequenced either on the Sequel II platform on SMRT Cells 8M (PacBio, 101-389-001) using Sequel II sequencing chemistry 2.0 (n = 16, PacBio, 101-842-900), 2.2 (n = 4, PacBio, 102-089-000), or 3.2 (n = 32, PacBio, 102-333-300) with 2h pre-extension and 30h movies on SMRT Link v9-11.1, or on the Revio platform on Revio SMRT Cells (PacBio, 102-202-200) and Revio polymerase kit v1 (n = 126, PacBio, 102-817-600) with 2h pre-extension and 24 or 30h movies on SMRT Link v12.0-13.1.

To ensure the use of high-quality reads for constructing robust assemblies, we first filtered out nonhuman contamination reads from both HiFi and ONT data. We employed highly accurate Illumina reads and utilized yak (commit f389bad, https://github.com/lh3/yak.git) to calculate the QV for each read, and a two-sided z-test was conducted on the resulting QV values. Reads with a z-score less than −2, indicating a potential risk of contamination, were compared against the Kraken2 (v2.1.3) database, and those identified as nonhuman in origin were excluded. We masked the Y chromosome from GRCh38 to generate the GRCh38noY reference genome. HiFi reads were aligned to GRCh38 for males and to GRCh38noY for females using pbmm2 (v1.13.1, https://github.com/PacificBiosciences/pbmm2). To ensure the family pedigree, the relatedness between sample pairs and ancestry prediction were conducted using Somalier (v0.2.19) based on the alignment. ONT reads were aligned to the reference genome using minimap2 (v2.28.0), and the family pedigree was confirmed with VerifyBamID (v2.0.1). In addition, ntsm (v1.2.1) was applied to each HiFi and ONT fastq for sample swap detection.

Phased genome assembly construction
The HiFi assembly was constructed by hifiasm (v0.16.1). Parental short reads were processed with yak (v0.1, https://github.com/lh3/yak.git) and then hifiasm trio-binning mode was used for phasing child samples, while the parental assemblies were partially phased by default. The sex chromosome contigs from father samples were aligned to the T2T-CHM13v2.0 reference to reassign the Y chromosome contigs to hap1 (or paternal haplotype) and X chromosome contigs to hap2 (or maternal haplotype). Assembly QVs were evaluated by merqury (v1.3) with k-mers from Illumina data (meryl v1.4); next, the completeness of phased assemblies relative to the reference and contig N50 values were calculated.

Variant discovery
SNVs and small indels in 73 children, with available ONT data, were recalled with DeepVariant (v1.4.0) and GATK (v4.3.0.0) based on the HiFi alignments to the T2T-CHM13v2.0 reference. DNMs including both de novo and postzygotic mutations were further validated by ONT and/or Illumina reads using the method described in Noyes et al. We annotated DNMs using the Ensembl Variant Effect Predictor (VEP, v110.1) and referred to the predicted impact and scores from dbNSFP (v4.8a), CADD score (v1.3), and gnomAD genome allele frequency (v4.1.0) by lifting coordinates over to both GRCh38 and GRCh37 using UCSC LiftOver. SVs were called using the phased assembly variant caller (PAV, v2.3.4) by aligning the assembled genomes to a reference (GRCh38 or GRCh38noY). The alignment-based SVs were called by PacBio SV calling and analysis tool (PBSV, v2.9.0, https://github.com/PacificBiosciences/pbsv) and Sniffles (v2.2). SVs from the control population were additionally called by Delly (v1.2.6), Sawfish (v0.12.4), and cuteSV (v2.1.0). The ratio of insertions to deletions, the ratio of heterozygous to homozygous variants, and the size distribution of SVs were evaluated for each VCF file.

SV merging and filtering strategy
The strategy involved three major steps. Briefly,
1. Caller set validation. SVs detected in each sample from different SV callers were first normalized, sorted, and merged with BCFtools (v1.20) on the basis of PAV callset, and then collapsed by Truvari (v4.3.1). SVs located within know genomic gaps, telomeric regions, centromeres, and PARs on GRCh38 from UCSC Genome Browser tracks were excluded.
2. Inter-sample merge. We extracted SVs supported by PAV and at least one of the alignment-based callers for each individual and then merged SVs from both controls and autism families using a list of VCFs with Truvari. In this study, we have applied three control sets consisting of 108, 285, and 569 individuals.
3. Rare SV pool discovery. We developed a custom script to extract six categories of rare SVs as described in the main text (v1.0.0, https://github.com/EichlerLab/asap). For autosomal SVs, we retained heterozygous and homozygous SVs present only in the children but not in the controls. For SVs on the sex chromosomes, we performed sex-matched comparisons and filtered SVs seen in controls with the same sex. Separately, we compared our LRS SV datasets with the SRS control set from 63,046 unrelated genomes in gnomAD v4.1 to evaluate the overlap between platforms. SV comparison was conducted using Truvari bench with relaxed matching parameters to account for differences in breakpoint precision between LRS and SRS datasets.

Transmission curation of rare SVs
Parental genotypes corresponding to each rare SV in the children are provided. Variants that followed Mendelian inheritance patterns were designated as high-confidence SVs. The remaining SVs were subjected to transmission curation using the following toolchain:
1. Initial caller support using Truvari. To minimize the loss of inheritance information in parents lacking SV caller support, we collapsed SVs from parents with the child’s SVs using truvari bench.
2. Callable region evaluation using BoostSV (v1.0). To ensure the SVs fall within confidently callable regions across samples in a single family, we developed a tool, BoostSV (v1.0, https://github.com/jiadong324/BoostSV), leveraging a machine-learning approach trained on control samples. This tool evaluates read support, mapping quality, and data quality metrics from alignments surrounding the target SVs in each parent. A quality threshold of 0.5 was applied to obtain the transmission.
3. Genotyping support using kanpig. We applied the k-mer-based genotyper kanpig (v0.3.1) to parental HiFi alignments for each SV to assess allele presence and genotype consistency.
4. Rare TR expansions/contractions and multiple sequence alignment (MSA). SVs overlapped with TR catalogs derived from the four-generation control family in Porubsky et al. were genotyped by TRGT (v1.4.1) using HiFi alignments.
5. Read-based support validation using subseq. To further assess SV transmission, subseq (v1.0, https://github.com/EichlerLab/subseq-smk) was used to quantify read support for each SV in parental genomes.
6. Manual inspection using IGV. For SVs lacking sufficient support in Steps 1-5, we conducted visual inspection of supporting reads in the Integrated Genomics Viewer (IGV, v2.16.0).

SV annotation
A customized script was implemented to annotate previously published NDD candidate genes and regulatory elements, as well as integrate annotations from AnnotSV (v3.4). In terms of REG, we integrated published datasets from UCSC Genome Browser tracks, including candidate cis-regulatory elements (ENCODE Regulation, ENCODE cCREs, ORegAnno, GeneHancer) and ENCODE histone marks (H3K27Ac, H3K4Me1, H3K4Me3). Additionally, we incorporated epigenomic profiles from the cerebral cortex (brainREG), particularly a cis-regulatory element map generated from 27 male and 21 female prenatal human cortex samples by ATAC-seq and consensus maps for CTCF, H3K27ac, H3K27me3, and H3K4me3 generated from six male and five female prenatal human cortex samples by CUT&Tag. For NDD candidate genes, we primarily focused on those previously reported. Additional ENCODE TF Clusters, UCSC noncoding RNA (tRNA, snRNA, lincRNA, sno_miRNA), repetitive regions (UCSC SegDup, UCSC RepeatMasker, UCSC Simple Repeats, TRs), noncoding constraint Gnocchi score and CADD-SV score (v1.1.2) were annotated.

Sex chromosome assemblies and transmission
We partitioned the sex chromosomes (T2T-CHM13v2.0, excluding the PAR, centromere, and highly repetitive Yq12 heterochromatin regions) into 1 Mbp windows and identified those covered by contigs that aligned to ≥95% of the window sequence with no more than three overlapping contigs. The coverage percentage was calculated as the number of qualified windows divided by the total number of windows, representing coverage relative to the reference sex chromosomes. We designed a pipeline to assemble contiguous X and Y chromosomes and validate transmission patterns within families (https://github.com/projectoriented/contiguous-X). The pipeline consists of two main steps: (1) scaffolding contigs with at least 50% of their sequence aligning to chromosome X (minimap2) via RagTag (v2.1.0), and (2) visualizing alignments across haplotypes within the family using the SVbyEye R package.

Methylation analysis
We developed a pipeline to extract phased methylation signals from ONT alignments (https://github.com/projectoriented/continuous-methylation). Briefly, when parental Illumina reads were available, ONT reads from the offspring were phased using the Canu (v2.1.1) trio-binning method. The phased ONT reads were then aligned to the GRCh38 reference genome using minimap2 (v2.24.0) and haplotagged accordingly. For individuals (e.g., parents) without available parental Illumina data, SVs were identified using Sniffles (v2.2) and small variants were called with Clair3 (v1.0.2), followed by haplotagging with LongPhase (v1.7.2). Methylation tags from unmapped BAM files were linked to the phased alignments using methylink (v0.6.0, https://github.com/projectoriented/methylink). Finally, Modkit (v0.3.1, https://github.com/nanoporetech/modkit) was used to generate methylation calls in BED format via pileup function. We analyzed mean methylation differences between the two X chromosome haplotypes across 889 CpG islands, including ±5 kbp flanking regions. We then used deepTools (v3.5.5, https://github.com/deeptools/deepTools) to compute the methylation matrix and generate the plots.

Long-read Kinnex generation and analysis
Total RNA was extracted from 5 million lymphoblast cells using the RNeasy mini kit from Qiagen. After initial QC using DS-11 FX (Denovix) for UV-Vis quantification and Bioanalyzer RNA 6000 Nano (Agilent, P/N G2939A and 5067-1511) for quality score calculation, 150–250 ng per sample were processed into cDNA and sequencing library concatamers using the Iso-Seq Express 2.0 and Kinnex full-length RNA kits (PacBio, P/N 103-071-500 and 103-072-000) according to manufacturer’s protocols with the following modifications: amplified cDNAs were analyzed for average size using Femto Pulse (Agilent P/N M5330AA and FP-1003-0275) and quantified using fluorescence (Qubit HS DNA Fisher Scientific P/N Q32854) before equimolar pooling and an additional 1.1X SMRTbell Cleanup Bead wash to remove residual primer-dimer before proceeding to Kinnex PCR. Kinnex concatamers were checked for length using Femto Pulse (P/N FP-1002-0275) before loading on 1 SMRT Cell 25M on the PacBio Revio sequencing platform using SPRQ chemistry with Adaptive Loading and a 30-hour acquisition time. After sequencing, raw data were processed with the Read Segmentation and Iso-Seq workflow in SMRT Link v25.3 using default parameters to generate sample-demultiplexed full-length reads (flnc.bam) files. BAM files were converted to FASTA format using SAMtools (v1.21) and aligned to the GRCh38 reference genome with minimap2 (v2.28) using the parameters -ax splice:hq -uf --secondary=no --eqx -K 200M. Gene and transcript quantifications were performed with IsoQuant (v3.10.0) using GENCODE v49 annotations. Raw gene counts from the target proband and remaining samples were used as case and control groups, respectively, for differential expression analysis with DESeq2 (v1.50.0), incorporating sex and diagnosis as variables.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Sui et al. present a study using long read sequencing (LRS) data from 51 families of unsolved autism cases to perform deep characterization of de novo mutations, structural variants, and methylation. They use newly available LRS/assemblies from controls from HPRC/HGSVC to substantially decrease the search space and identify potential rare/denovo pathogenic variants at structural variants (SVs). Although the study did not find an overall increased SV burden in probands, this is likely due to insufficient power at this sample size.\nThis study is a valuable proof of principle in showing the promise of using LRS as a diagnostic tool. A major challenge in moving to LRS is a lack of large population control databases of SVs generated by these technologies (compared to up to at least hundreds of thousands of short read datasets that are used as controls for SNVs). Although the control sample set used here is still pretty small, it enables filtering 97% of variants (99% with the expanded control set) which is promising. It is also encouraging that the study found known things (e.g. SYNGAP1, MECP2, TBL1XR1 mutations) that were not found or filtered by short reads. Finally, the study generates high quality LRS data and assemblies that can presumably be used as a resource in future studies.\nOverall the paper is well written and easy to follow. I have several comments regarding clarifying some of the terminology used and toning down some of the claims that are made or strengthening them with some additional analysis:\n1. The manuscript oversells claims about the completeness of the control variant callsets. In particular, if I understood the sample size correctly (108 controls), most of the variants in controls are probably common and many ultra-rare things are missing. It would be helpful to provide an estimate of the % of variants in the population at different MAF thresholds that are estimated to be in the control set (or some other sort of metrics that could quantify the completeness of the reference. Fig. 5a shows the cumulative SV number but doesn’t assess how low of MAF we think that reference captures).\nSome specific places in the text that are relevant:\nLine 130: The claim that variant discovery from HPRC/HGSVC is complete is misleading. This is true in the sense that all of the variants in each genome are theoretically discovered. But not true in the sense that these resources still have limited sample size and thus are missing many rare variants. This nuance should be clarified.\nLine 144: what MAF threshold is being used to define ultra-rare? The “Variant discovery” section of the results mentions 108 population controls (line 189). Which would mean the lowest frequency that can be detected is >0.1%, and likely a large number of rare things are missing from that reference by chance.\nIn the last section of the paper, the control pool is expanded to 569 individuals. It wasn’t clear why not to just use that larger control set in the first place?\n2. How much improvement do the LRS control datasets give for filtering control SVs compared to other datasets based on short reads? The authors state that 97% of variants can be filtered using the population controls, which is exciting. Are there available short read SV catalogs that could have been used and would those add anything beyond the LRS controls? It is clear that per genome the number of SVs discovered is far smaller than that from LRS, but the available sample sizes are much larger which could help with rare variants.\n3. I find the usage of the word pangenome here and elsewhere to be confusing. This is largely because the use of the term in general in the field (not specifically in this paper) is unclear. Technically, a pangenome is a collection of genomes from individuals of the same species, which this paper uses. But the word has been commonly associated specifically with using that collection as a reference (rather than a single linear reference) and also to graph-based representations (e.g. as in minigraph, pggb, etc although other paradigms to represent pangenomes exist). This paper does neither of these, and is rather using SVs from these genomes as a set of population controls. If those control SVs had been derived from short reads, e.g. from 1000 Genomes, would that also be considered a pangenome analysis? The title, “Pangenome discovery of missing autism variants”, gave me the impression the pangenome would actually be used for genotyping e.g. as a reference rather than mostly for downstream filtering, but it appears genotyping was still done using GRCh38 as a reference. I would suggest revising the title to be less misleading about that and also clarifying in the intro what constitutes pangenome analysis.\nAdditional points:\nThe X chromosome result of a trend toward increased SV burden on chrX among affected females is highlighted in the abstract but is not actually statistically significant. While the result is promising, it should be clarified in the abstract that this is only a suggestive trend.\nLine 207: The average Mendelian Inheritance rate for SVs across the trios is 90.4%. How does this compare with other SV studies (seems maybe comparable to Liao et al Nature 2023 HPRC paper?). I assume the mutation rate of SVs is far less than 10% so most of the inconsistencies are driven by genotyping errors. What are the implications of that for the potential for false positives here?\nLine 255: “The brain-derived regulatory regions (brainREG) annotation resulted in a 31% increase in the total number of regulatory SVs” -> compared to what as a baseline?\nLine 250 section “SV burden analysis”: This paragraph starts off about regulatory annotations, but I think the second half of the paragraph is focusing on overall burden (starting on line 263). It would be good to clarify which results are overall vs. across the categories tested. If line 263 onwards is not about the regulatory annotations that might be better as a new paragraph.\nLine 341: MPRA is not defined\nLine 347: worth noting that HNRNPK is also an RNA binding protein", "reviewer_2": "The manuscript submitted by Evan Eichler, Huda Zoghbi and colleagues is first and foremost a novel conceptual paper in the implementation of near-complete T2T sequencing into the clinic. In this way it sets standards and guides future T2T diagnostic papers.\nThe authors analysed by near-complete genome sequencing 189 individuals from 51 families with unsolved autism. Specifically the authors claim that most of the variants identified could not have been detected by srGS.\nBioinformatic tools were more or less the todays standard tools in the field showing, that wetlabtechnology and bioinformatics are ready to enter the clinic. There will be further improvements over time as the current GRCh38 reference data set still has numerous gaps or errors, however, this manuscripts adds to additional samples to be used as reference data sets.\nThe authors state nicely why srGS might have failed to detect several clear causative variants. However, it seems to me that standard long read GS may have detected these variants as well and a de novo assembly might not have been necessary, please comment on this.\nA couple of potential novel finding ahve been made like identification of repeat expanions exclusively in males inherited from their mothers and several patients with skewing of X inactivation. As there are not suppl data on clinical presentation and as the cohort most likely is quite heterogeneous (already based on the molecular identification of Rett syndrome in the cohort) it may add to the manuscript to mention whether there are any specifics on these patients (such as developmental issues, IQ or severe forms of autism), please add.\nFinally, in several of the potentially identified mutations (Table 1) novel regulatory mutations may have been identified. Many of them are easily detected by RNAseq analysis. It would clearly add to the manuscript to have these verified whether or not the variants led to a reduced gene expression.\nMinor comments:\npage 4, line 102 \"conventional SRS approaches such as Illumina whole-genome sequencing...\"\nI guess this would have been similar with BGI or Ultima (or other srGS) based technologies, so I would deleted the word \"Illumina\"\npage 4, line 119.. it seems to me that the references do not fit entirely\npage 5, line 141 ...in the proband via conventional Illumina SRS... please specify, do you mean Illumina whole genome SRS? or was this to a certain percentage also WES?\npage 5, line 168, I did not understand why you had to use parental Illumina reads for hifiasm to generate the parental haplotypes as you had already each person sequenced with PacBio.\npage 10, lines 281 - 285, the authors identified three hemizygous repeat outliers in males, all inherited from the mother.\nBased on the limited reference data today, it can not be judged whether or not these are causal repeat expansions. In this chapter the authors refer only to the IL1RAPL1 but not to the F9 and intergenic repeat expansion (all provided correctly in suppl S8), however for easier readability it would be nice to list them here all. Any specifics on the country of origin?\npage 10, line 285: I guess between \"....SFARI score 2 gene). and .....Because LRS data allow...\" whould be a new paragraph", "reviewer_3": "The genetic basis of autism remains underexplored. In this study, the authors performed genome sequencing of 189 individuals from 51 families with unsolved cases and generated a set of high-quality assemblies. These assemblies were compared with publicly available high-quality assemblies from healthy donors, leading to the discovery of three pathogenic variants in TBL1XR1, MECP2, and SYNGAP1, as well as nine candidate variants, most of which were missed by short-read sequencing. The study is well designed and solid. I have a few concerns as outlined below.\nMajor comments\nMECP2 is structurally complex with diverse SVs and CNV events. I am wondering whether this case also harbors any other events besides the discovered deletion. Moreover, what is the actual deleted sequence in the last exon region, and what is the precise effect of this deletion on the last exon?\nThe paper includes a large control cohort, which is very good. However, controls appear underrepresented for EUR ancestry (e.g., 39 EUR among 569). Please report the ancestry and sex composition of controls vs. families, and perform ancestry- and sex-matched analyses.\nIn Fig. 1, the authors showed the N50 of contigs in autism family assemblies. However, providing NGx curves against GRCh38 and CHM13 would still be valuable. Comparing family assemblies with control assemblies would also be helpful.\nThe authors claim a “modest rate of pathogenic variant discovery (5.9%).” Did the authors consider evaluating compound heterozygous and oligogenic models (e.g., small variants, or small variants + SVs)? This may help recover additional pathogenic candidates.\nPlease describe unpublished tools (e.g., subseq, BoostSV) and provide versions, key parameters, and workflow availability to enhance reproducibility.\nMinor comments\nFig. 2a: Distinguish proband vs. sibling more clearly (shapes or colors).\nFig. 4c,d: Correct “Genehance” → “GeneHancer.”\nTypo: “TBXLR1” → “TBL1XR1.”\nI would like to thank the authors for their contribution to the field through this work."}, "author_response": "We are grateful to the reviewer for their careful reading and valuable suggestions.\n1. The manuscript oversells claims about the completeness of the control variant callsets. In particular, if I understood the sample size correctly (108 controls), most of the variants in controls are probably common and many ultra-rare things are missing. It would be helpful to provide an estimate of the % of variants in the population at different MAF thresholds that are estimated to be in the control set (or some other sort of metrics that could quantify the completeness of the reference. Fig. 5a shows the cumulative SV number but doesn’t assess how low of MAF we think that reference captures).\nThis is a good suggestion. We added some comments to the text regarding the “completeness” of the 108-control set and the possibility of missing variants based on sample size. To address the comment of how many SVs are likely missing in that limited set, we created a population SV set and plotted the MAF in the 108-control set to the larger sample set. From that, we computed the number of missing variants in the 108 controls as a function of MAF in 569 controls and 102 unrelated parents (n=1,342). We added the following statement and one SI figure (with underline indicating the modifications in response to Reviewer #3):\n“Increasing the number of controls, especially samples of African origin, nearly doubles the number of SVs from 271,375 (108 controls) to 445,142 (569 controls) nonredundant SVs (Supplementary Data 5-8). We generated a population-level SV reference and compared the MAF distribution between the 108 controls and the largest dataset comprising 569 controls and 102 unrelated parents (n = 1,342). Based on this comparison, we estimated the fraction of variants captured at different MAF thresholds (Supplementary Fig. 12). Notably, only 16.2% of SVs with MAF < 0.1% in the larger dataset were detected in the 108-control set, illustrating the limited sensitivity of the smaller reference for identifying ultra-rare variants. In an ancestry- and sex-matched analysis restricted to individuals of European ancestry (38 autism families and 73 controls), each child carried an average of 754 rare SVs. Focusing on rare variants absent from the largest 569-control set, the number of rare SVs concomitantly drops to 202 events per sample (Fig. 5b), corresponding to ~74% of rare variants being filtered out by the more diverse population controls. This demonstrates the substantial gain in sensitivity for removing common and low- frequency variants when using a larger, heterogeneous reference panel. As a result, 99% of common SVs are excluded per individual forming a more tractable and potentially biomedically relevant set of rare variants for downstream interpretation and enrichment analyses…. ”\nSupplementary Figure 12. Effect of control cohort size on SV discovery and minor allele frequency (MAF) spectrum. The plot shows the fraction of SVs captured (blue) in an initial control set of 108 genomes across different MAF thresholds derived from an expanded dataset comprising 569 controls and 102 unrelated parents (1,342 haplotypes in total). Only 16.2% of SVs with MAF < 0.1% in the larger dataset were detected in the smaller 108-control reference, underscoring the limited sensitivity of the smaller cohort for identifying ultra-rare variants.\nSome specific places in the text that are relevant:\nLine 130: The claim that variant discovery from HPRC/HGSVC is complete is misleading. This is true in the sense that all of the variants in each genome are theoretically discovered. But not true in the sense that these resources still have limited sample size and thus are missing many rare variants. This nuance should be clarified.\nWe revised the sentence to reflect this nuance as suggested:\n“This resource is potentially valuable to the clinical genetics community because variant discovery is more complete, providing a control set to assess the frequency of variants in regions typically unassayable by SRS and therefore absent or unreliable in associated databases such as gnomAD (Collins et al. 2020; S. Chen et al. 2024). Moreover, assembly-based comparisons between offspring and parental genomes have been shown to further increase the power to discover DNM by essentially eliminating reference biases (Porubsky et al, 2025). Notwithstanding, the number of samples still remains modest and much larger sample sizes will be required to understand the full spectrum especially for those with lower minor allele frequency (MAF).”\nLine 144: what MAF threshold is being used to define ultra-rare? The “Variant discovery” section of the results mentions 108 population controls (line 189). Which would mean the lowest frequency that can be detected is >0.1%, and likely a large number of rare things are missing from that reference by chance.\nThat is a good catch as they are hardly ultra-rare at a frequency of <1/216 from 108 population controls. We revised the sentence accordingly:\n“…we then compared with HPRC and HGSVC population controls to discover and validate de novo and rare variants (<0.5%) as new candidates for autism.”\nIn the last section of the paper, the control pool is expanded to 569 individuals. It wasn’t clear why not to just use that larger control set in the first place?\nOur experiments were staged over a year and reflect a continued expansion of the 1000 Genomes Project (1KGP) samples as part of ongoing efforts from our laboratory and others. We felt that showing this progression might be useful to the research community as the controls used have diverse ancestry. As described both in the first and last Results sections, we used the initial control set of 108 individuals from the HPRC and HGSVC, previously published in Liao et al. 2023, Logsdon, Rozanski, et al. 2024 and Ebert et al. 2021. We then analyzed 177 newly released HPRC samples ourselves, expanding the set to 285. To further increase the power, we subsequently included 284 additional ONT controls from the 1KGP, generating the larger set of 569 controls. These stepwise expansions reflect the availability of data at the time of analysis, while also highlighting that larger control sets substantially increase statistical power and robustness, allowing us to better distinguish rare signals from background common variants as now nicely illustrated by the suggested Supplementary Fig. 12. To clarify, we now include citations of published controls and revised the text to reflect this aspect of the experimental design:\n“We aggregated the validated SVs from all 189 study samples and compared them with SVs observed in the 108 population controls from the HPRC and HGSVC (Ebert et al. 2021; Logsdon, Rozanski, et al. 2024; Liao et al. 2023) via Truvari (English et al. 2022), focusing on rare SVs exclusive to autism families. As a result of subsequent LRS of additional 1000 Genomes Project (1KGP) samples recently released as part of the HPRC and other efforts (Gustafson et al. 2024), we expanded our SV callset from controls in a staged manner, allowing us to access variants of reduced MAF (see below Supplementary Fig. 12).”\n2. How much improvement do the LRS control datasets give for filtering control SVs compared to other datasets based on short reads? The authors state that 97% of variants can be filtered using the population controls, which is exciting. Are there available short read SV catalogs that could have been used and would those add anything beyond the LRS controls? It is clear that per genome the number of SVs discovered is far smaller than that from LRS, but the available sample sizes are much larger which could help with rare variants.\nWe performed the suggested comparison using GATK_SV calls from 63,000 genomes from gnomAD SRS data and the filtering gain was marginal (0.27%, 55/20,716 SVs). By contrast, if we limited to SRS, a large fraction of SVs (~68%, 14,110/20,716) remained undetected. Because we believe this analysis is of interest to the general readership, we added the following paragraph (in bold) to the last result section of the main text:\n“To enhance the power of pangenome filtering to effectively exclude common SVs and focus on a high-confidence pool of rare variants, we expanded the control cohort size from 108 to 285 and then to 569 individuals. …. \nWe also compared our LRS SV datasets with a comprehensive SRS SV callset (INS and DEL) generated by GATK-SV (Collins et al. 2020) from 63,046 unrelated genomes, as represented in gnomADv4.1 (Collins et al. 2020). To ensure a fair comparison between SRS and LRS datasets that differ in resolution and breakpoint precision, we used relaxed matching criteria in Truvari bench: a minimum of 50% reciprocal overlap, 50% size similarity, maximum breakpoint distance of 500 bp, and no sequence similarity requirement. We added the gnomAD SV IDs to the supplementary datasets (Supplementary Data 2, 5 and 7, Supplementary Fig. 14). Within the 108-control set, 162,916 of 271,375 total nonredundant SVs were observed in the autism families. Of these, 77.6% (126,377 SVs) were present in the 108 LRS controls, yielding 663 rare SVs per child. In contrast, only 24.1% (39,234 SVs) overlapped with the SRS-based gnomAD SVs from 63,046 unrelated controls (supplementary Fig. 14), corresponding to 14,110 rare SVs per child when considering only short-read controls. Thus, despite its much larger sample size, the short-read dataset can be considered an auxiliary resource for further filtering rare SVs beyond those already excluded by the LRS controls, underscoring the substantially greater power of LRS data for rare variant filtering. Specifically, incorporating the SRS dataset would have further excluded 0.71% (147/20,716 per child) and 0.27% (55/20,716 per child) of SVs from the 108- and 569-control sets, respectively. Given the inherent differences in breakpoint resolution and variant representation between LRS and SRS platforms, these overlaps reflect potentially inflated estimates. A reassuring sign of the stringency and completeness of our filtering strategy is that all 12 candidates identified were rare variants absent from both the LRS and SRS control datasets.”\nAlso, in the Methods:\n“SV merging and filtering strategy. The strategy involved three major steps. Briefly,\n1. Callerset validation. ...\n2. Inter-sample merge. ...\n3. Rare SV pool discovery. …\nSeparately, we compared our LRS SV datasets with the SRS control set from 63,046 unrelated genomes in gnomAD v4.1 to evaluate the overlap between platforms. SV comparison was conducted using Truvari bench with relaxed matching parameters to account for differences in breakpoint precision between LRS and SRS datasets:\ntruvari bench -c {gnomad.v4.1.INSDEL50.non_neuro_controls.sites.vcf.gz} -b {108/285/569ctr_189asd_collapsed.vcf.gz} --pctsize 0.5 --pctseq 0.0 --pctovl 0.5 --sizefilt 50 --sizemax 100000 -o {output}\nWe incorporated gnomAD SV IDs, overlap metrics (PctSizeSimilarity:PctRecOverlap:SizeDiff), and allele counts (heterozygous, homozygous, and sex-specific; overall and within the non-neuro subset) into the collapsed SV tables (Supplementary Data 2, 5 and 7). Similar to step 3 above, we applied sex- and zygosity-aware filtering to identify rare SVs absent from all SRS genomes.\nSupplementary Figure 14. Comparison of nonredundant SVs observed in LRS and SRS control sets. Collapsed nonredundant SVs from the three LRS control datasets were compared against SRS-based SVs (INS and DEL) generated by GATK-SV from 63,046 unrelated genomes in gnomAD v4.1 (corresponding to Supplementary Data 2, 5 and 7). Overlap between LRS and SRS SVs was defined using relaxed matching criteria in Truvari bench: a minimum of 50% reciprocal overlap, 50% size similarity, and maximum breakpoint distance of 500 bp, with no sequence similarity requirement. Cross-hatched bars indicate SVs that would have been detected by SRS.\n3. I find the usage of the word pangenome here and elsewhere to be confusing. This is largely because the use of the term in general in the field (not specifically in this paper) is unclear. Technically, a pangenome is a collection of genomes from individuals of the same species, which this paper uses. But the word has been commonly associated specifically with using that collection as a reference (rather than a single linear reference) and also to graph-based representations (e.g. as in minigraph, pggb, etc although other paradigms to represent pangenomes exist). This paper does neither of these, and is rather using SVs from these genomes as a set of population controls. If those control SVs had been derived from short reads, e.g. from 1000 Genomes, would that also be considered a pangenome analysis? The title, “Pangenome discovery of missing autism variants”, gave me the impression the pangenome would actually be used for genotyping e.g. as a reference rather than mostly for downstream filtering, but it appears genotyping was still done using GRCh38 as a reference. I would suggest revising the title to be less misleading about that and also clarifying in the intro what constitutes pangenome analysis.\nYes, there does seem to be confusion in the field in part because terms need to be clearly defined. There are three definitions/activities that are distinguishable: 1) the generation of the pangenome, which as indicated by the referee is simply the collection of nearly complete linear references; 2) construction of a graph-based pangenome, which combines these linear references to a singular reference where variants are captured and correctly represented by the underlying graph structure; and 3) the use of that graph-based reference or the collection of linear references themselves for (hopefully improved) genotyping into whole-genome sequencing data. Our paper focuses on #1 and #3 without the use of a graph, i.e., we only use the collection of linear references. While it is true that SRS datasets from many human genomes might technically be considered a pangenome if SRS were assembled, implicit in a pangenome design is the construction of an ab initio assembly first, prior to variant discovery as opposed to aligning reads to a singular reference, which is how SVs are typically discovered by SRS and would not qualify as a pangenomic approach.\nWhile we prefer the shorter title and believe it is valid, a more descriptive title that would make it clear that we are not employing a graph could be:\n“Using the linear references from the pangenome to discover missing autism variants”\nAdditional points:\nThe X chromosome result of a trend toward increased SV burden on chrX among affected females is highlighted in the abstract but is not actually statistically significant. While the result is promising, it should be clarified in the abstract that this is only a suggestive trend.\nWe revised the abstract to clarify that the observation represents a suggestive trend rather than a statistically significant result. The revised abstract now reads:\n\"…yet observe a suggestive trend toward an increased SV burden on the X chromosome among affected females.”\nWe also made this clearer in the legend:\n“No significant differences (χ² test p-values exceeding 0.05) in the number of SVs between probands and siblings were observed across these categories, with a still insignificant trend observed on the X chromosome for enrichment of SVs on affected females compared to unaffected sisters.”\nLine 207: The average Mendelian Inheritance rate for SVs across the trios is 90.4%. How does this compare with other SV studies (seems maybe comparable to Liao et al Nature 2023 HPRC paper?). I assume the mutation rate of SVs is far less than 10% so most of the inconsistencies are driven by genotyping errors. What are the implications of that for the potential for false positives here?\nOur Mendelian inheritance rate is comparable to other studies that have investigated using LRS and called SVs (Duan et al., BMC Genomics, 2002; Schloissnig et al., Nature, 2025). There are at least two factors: the individual callers used and the type of variants considered. Mendelian errors are highest for multi-allelic variants (e.g., VNTRs) where there have been reports as high as 15%. It should be noted however that such variants are particularly prone to alignment artefacts as well as sequencing artefacts–so there is also a reproducibility in calling. We have expanded upon this to make the implications with respect to the nature of the false positives clearer:\n“The average Mendelian concordance rate for SVs across 87 trios is 90.4% with many of the discrepant alleles associated with multiallelic variants, such as variable number tandem repeats (VNTRs) where both sequence and alignment artefacts confound variant calling (Schloissnig et al. 2025).”\nLine 255: “The brain-derived regulatory regions (brainREG) annotation resulted in a 31% increase in the total number of regulatory SVs” -> compared to what as a baseline?\nWe apologize for this confusing sentence and the 31% statistic was not correct. The additional annotation provided by ATAC-seq and CUT&Tag experiments allowed us to identify an additional 6,171 SVs that intersect with brain-specific regulatory regions that were not defined by ENCODE.\n\"The additional annotation allowed us to identify another 6,171 SVs associated with brain-derived regulatory regions (brainREG) beyond those defined by ENCODE intersection and consequently a 45% (6,171/13,773) increase in potential SVs affecting noncoding regulatory DNA (Fig. 2c-d).\"\nLine 250 section “SV burden analysis”: This paragraph starts off about regulatory annotations, but I think the second half of the paragraph is focusing on overall burden (starting on line 263). It would be good to clarify which results are overall vs. across the categories tested. If line 263 onwards is not about the regulatory annotations that might be better as a new paragraph.\nThis is a good suggestion so we split the paragraph as recommended:\n“.... The rates of autosomal SVs (Fig. 2c, Supplementary Fig. 6) do not differ significantly between probands and unaffected siblings across various categories (nominal p > 0.05, OR < 1, χ² test).\nFor X chromosome SVs, we analyzed males and females separately. …”\nLine 341: MPRA is not defined\nWe now define the MPRA acronym as “massively parallel reporter assay” in the text.\nLine 347: worth noting that HNRNPK is also an RNA binding protein\nWe revised the text to indicate that HNRNPK, in addition to acting as a transcription factor, also functions as an RNA-binding protein:\n\"This INS is predicted to disrupt HNRNPK transcription factor (TF) binding sites in 11201_p1; notably, HNRNPK also functions as an RNA-binding protein.\"\nReviewer #2 (Remarks to the Author):\nThe manuscript submitted by Evan Eichler, Huda Zoghbi and colleagues is first and foremost a novel conceptual paper in the implementation of near-complete T2T sequencing into the clinic. In this way it sets standards and guides future T2T diagnostic papers.\nThe authors analysed by near-complete genome sequencing 189 individuals from 51 families with unsolved autism. Specifically the authors claim that most of the variants identified could not have been detected by srGS.\nBioinformatic tools were more or less the todays standard tools in the field showing, that wetlabtechnology and bioinformatics are ready to enter the clinic. There will be further improvements over time as the current GRCh38 reference data set still has numerous gaps or errors, however, this manuscripts adds to additional samples to be used as reference data sets.\nWe thank the referee for their constructive and helpful comments.\nThe authors state nicely why srGS might have failed to detect several clear causative variants. However, it seems to me that standard long read GS may have detected these variants as well and a de novo assembly might not have been necessary, please comment on this.\nThe referee raises a good point. While we found local assembly helpful for confirming parent-of-origin, the missing SRS variants did not, in most cases, require an assembly. We added that point to the first paragraph of the discussion:\n“Our analysis of 51 families mostly with daughters affected with autism (often more severely) where we attempted to sequence and assemble the entire euchromatic portion of each genome suggests a more modest rate of pathogenic variant discovery (5.9%). We consider this yield of new pathogenic variant discovery low given that LRS application increased genome-wide sensitivity of DNM detection by 20-40% (Noyes et al. 2025) and we purposefully sequence affected females where the probability of discovery of a large effect mutation is expected to be higher (Iossifov et al. 2014). It should be noted that LRS alignment to a reference was sufficient to detect these variants previously missed by SRS. Genome assemblies, while helpful for parent-of-origin determination and phasing variants, were not critical to the discovery of most variants.”\nA couple of potential novel finding ahve been made like identification of repeat expanions exclusively in males inherited from their mothers and several patients with skewing of X inactivation. As there are not suppl data on clinical presentation and as the cohort most likely is quite heterogeneous (already based on the molecular identification of Rett syndrome in the cohort) it may add to the manuscript to mention whether there are any specifics on these patients (such as developmental issues, IQ or severe forms of autism), please add.\nThank you for pointing this out. We further added available Seizure and Calibrated Severity Score (CSS) from the ADOS to Supplementary Data 1 and evaluated the IQ, CSS and Seizure features in the patients with potential novel findings.\nWe revised the text as follows:\n“The analysis highlighted three hemizygous tandem repeat (TR) expansion outliers in male probands (Supplementary Fig. 9, Methods), each inherited from the mother in the intron of IL1RAPL1 (Interleukin 1 Receptor Accessory Protein Like 1, SFARI score 2 gene), the intron of F9 (Coagulation Factor IX), or the intergenic region between MXRA5 (Matrix Remodeling Associated 5) and SNORA48B (Small Nucleolar RNA, H/ACA Box 48B). These longer TR noncoding variants have only been observed in heterozygous states in females and variants of such lengths have yet to be observed in controls. Notably, the corresponding probands exhibited features commonly associated with neurodevelopmental delay (Supplementary Fig. 9): IQ scores of 18, 13, and 20, and calibrated severity scores (CSS) of 9, 6, with the CSS data for the third not available.”\nSupplementary Figure 9. Outlier TR expansions transmitted from mothers to male probands. a) Allele lengths (generated by TRGT) of three TR catalogs (Porubsky, Dashnow, et al. 2025) are shown in gray for unaffected individuals (285 controls and 138 unaffected samples) and in blue for affected probands. The three outlier expansions are highlighted in red and dark gray. b) Transmission patterns from maternal haplotypes to male probands are illustrated using multiple sequence alignments (MSAs), with comparisons to the GRCh38 reference. Note that females, with two X chromosomes, display bidirectional TR lengths (both alleles), whereas males, with only one X chromosome, show unidirectional TR length from a single allele. c) IQ and calibrated severity scores (CSS) of probands are shown (corresponding to Supplementary Data 1). Values for the three samples carrying outlier TR expansions on the X chromosome are indicated with red lines. Missing information denotes unavailable data.\nFinally, in several of the potentially identified mutations (Table 1) novel regulatory mutations may have been identified. Many of them are easily detected by RNAseq analysis. It would clearly add to the manuscript to have these verified whether or not the variants led to a reduced gene expression.\nThis was a very interesting suggestion. We attempted to address this comment by performing RNA-seq as well as Kinnex (long-read transcriptomic) sequencing on homozygous deletions (n=5) and promoter substitutions (n=2) where the associated genes would be expressed in the lymphoblastoid cell line. This limited the analysis to seven families, and we did not observe significantly reduced expression for these candidates, probably because they were mainly derived from blood samples.\nWe added these additional experimental details to the results and methods section:\nResults section: Pathogenic and autism candidate variant discovery.\n“To assess the potential impact of candidate variants, we performed differential expression analysis using long-read RNA sequencing of lymphoblastoid cell lines generated with the PacBio Kinnex platform. Gene-level counts were obtained with IsoQuant (Prjibelski et al. 2023) and analyzed with DESeq2 (Love et al. 2014) for probands carrying genes with two promoter mutations and five homozygous deletions that were expressed in blood. Across all seven candidate loci, no significant reduction in gene expression was detected in blood-derived RNA samples from the target proband compared to the remaining samples (Supplementary Fig. 11). The absence of detectable expression changes may reflect tissue-specific regulatory effects, as these variants are more likely to exert functional consequences in brain tissues, which are directly relevant to the disorder etiology.”\nMethods section:\n“Long-read Kinnex generation and analysis. Total RNA was extracted from 5 million lymphoblast cells using the RNeasy mini kit from Qiagen. After initial QC using DS-11 FX (Denovix) for UV-Vis quantification and Bioanalyzer RNA 6000 Nano (Agilent, P/N G2939A and 5067-1511) for quality score calculation, 150-250 ng per sample were processed into cDNA and sequencing library concatamers using the Iso-Seq Express 2.0 and Kinnex full-length RNA kits (PacBio, P/N 103-071-500 and 103-072-000) according to manufacturer’s protocols with the following modifications: amplified cDNAs were analyzed for average size using Femto Pulse (Agilent P/N M5330AA and FP-1003-0275) and quantified using fluorescence (Qubit HS DNA Fisher Scientific P/N Q32854) before equimolar pooling and an additional 1.1X SMRTbell Cleanup Bead wash to remove residual primer-dimer before proceeding to Kinnex PCR. Kinnex concatamers were checked for length using Femto Pulse (P/N FP-1002-0275) before loading on 1 SMRT Cell 25M on the PacBio Revio sequencing platform using SPRQ chemistry with Adaptive Loading and a 30-hour acquisition time. After sequencing, raw data were processed with the Read Segmentation and Iso-Seq workflow in SMRT Link v25.3 using default parameters to generate sample-demultiplexed full-length reads (flnc.bam) files.\nBAM files were converted to FASTA format using SAMtools (v1.21; Danecek et al. 2021) and aligned to the GRCh38 reference genome with minimap2 (v2.28; Li 2021) using the parameters -ax splice:hq -uf --secondary=no --eqx -K 200M. Gene and transcript quantifications were performed with IsoQuant (v3.10.0; Prjibelski et al. 2023) using GENCODE v49 annotations. Raw gene counts from the target proband and remaining samples were used as case and control groups, respectively, for differential expression analysis with DESeq2 (v1.50.0; Love et al. 2014), incorporating sex and diagnosis as variables.”\nSupplementary Figure 11. Blood-derived long-read RNA sequencing (PacBio Kinnex) of seven candidate genes. Normalized read counts were generated using the counts() function from DESeq2 for 11 samples (including seven target probands and randomly selected probands and siblings from families 11201 and 14232). Adjusted p-values were calculated by comparing each target proband to the remaining samples using the Wald test in DESeq2.\nMinor comments:\npage 4, line 102 \"conventional SRS approaches such as Illumina whole-genome sequencing...\"\nI guess this would have been similar with BGI or Ultima (or other srGS) based technologies, so I would deleted the word \"Illumina\"\nWe agree that the statement applies broadly to short-read genome sequencing technologies. We revised the text as suggested removing the reference to Illumina:\n“... using conventional SRS approaches associated with whole-genome sequencing…”\npage 4, line 119.. it seems to me that the references do not fit entirely\nThe cited studies are intended to support the general statement that long-read sequencing (LRS) is increasingly applied to unsolved patients and disorders, in addition to the specific examples mentioned. We revised the sentence to make this clearer:\n“Consequently, LRS has been increasingly applied to a variety of unsolved patients and disorders to enhance pathogenic variant discovery, although most studies to date have involved relatively modest cohorts focused almost entirely on read-based discovery (Hiatt et al. 2021, 2024; Sanchis-Juan et al. 2023; Pauper et al. 2021; Sinha et al. 2025; Steyaert et al. 2025). For example, Hiatt et al. reanalyzed 10 NDD families and 86 probands using PacBio high-fidelity (HiFi) LRS and found an additional yield of 7.3% beyond SRS mainly in coding regions (Hiatt et al. 2021; 2024). Sanchis-Juan et al. applied ONT LRS to complement SRS in four probands (Sanchis-Juan et al. 2023).”\npage 5, line 141 ...in the proband via conventional Illumina SRS... please specify, do you mean Illumina whole genome SRS? or was this to a certain percentage also WES?\nFor 46 of the families it was WGS but for the 5 Rett-like samples it involved WES. As described in detail in the Methods section, the unsolved probands were screened using Illumina WGS (n=46, SSC and SAGE cohorts) as well as WES or targeted gene panel testing (n=5, Rett cohort). We revised the text to:\n“In this study, we present our initial LRS and assembly resource of 189 individuals from 51 unsolved autism families where no pathogenic variant was previously identified in the proband via conventional Illumina whole-genome, whole-exome, or gene panel testing with SRS (Methods).”\npage 5, line 168, I did not understand why you had to use parental Illumina reads for hifiasm to generate the parental haplotypes as you had already each person sequenced with PacBio.\nAlthough all individuals were sequenced with PacBio HiFi, trio-binning in hifiasm requires parental Illumina WGS data to provide reliable parental-specific k-mers for phasing. These k-mers are then used to assign the child’s HiFi reads to maternal (hap2) and paternal (hap1) haplotypes (also called phased or haplotype-resolved assemblies). Importantly, by leveraging the family pedigree, we are able to generate more accurate and complete phased genomes for the children, making it straightforward to determine the parental origin of variants — a key strength of our study.\nTo clarify, we revised the sentence to:\n“Using parental Illumina reads and hifiasm (Cheng et al. 2021), we generated haplotype-resolved genome assemblies (i.e., phased by parent-of-origin) for each of the 87 offspring.”\nThese details were also described in the Methods section:\n\"Phased genome assembly construction. The HiFi assembly was constructed by hifiasm (v0.16.1, Cheng et al. 2024). Parental short reads were processed with yak (v0.1, https://github.com/lh3/yak.git) and then hifiasm trio-binning mode was used for phasing children samples, while the parental assemblies were partially phased by default.\"\npage 10, lines 281 - 285, the authors identified three hemizygous repeat outliers in males, all inherited from the mother. Based on the limited reference data today, it can not be judged whether or not these are causal repeat expansions. In this chapter the authors refer only to the IL1RAPL1 but not to the F9 and intergenic repeat expansion (all provided correctly in suppl S8), however for easier readability it would be nice to list them here all. Any specifics on the country of origin?\nWe agree with the reviewer and revised the text accordingly:\n\"The analysis highlighted three hemizygous tandem repeat (TR) expansion outliers in male probands, each inherited from the mother (Supplementary Fig. 9): the intron of IL1RAPL1 (Interleukin 1 Receptor Accessory Protein Like 1, SFARI score 2 gene), the intron of F9 (Coagulation Factor IX), and the intergenic region between MXRA5 (Matrix Remodeling Associated 5) and SNORA48B (Small Nucleolar RNA, H/ACA Box 48B). These longer TR noncoding variants have only been observed in heterozygous states in females, and variants of such lengths have yet to be observed in controls. Notably, the corresponding probands exhibited features commonly associated with neurodevelopmental delay (Supplementary Fig. 9): IQ scores of 18, 13, and 20, and calibrated severity scores (CSS) of 9, 6, with the CSSdata for the third not available.\"\nRegarding ancestry, the 13900 family has a 100% probability of AMR ancestry, while families 13700 and 14317 each have over 98% probability of EUR ancestry. These predictions were made using Somalier (v0.2.19 (Pedersen et al. 2020)) and were also reported in Supplementary Data 1.\npage 10, line 285: I guess between \"....SFARI score 2 gene). and .....Because LRS data allow...\" whould be a new paragraph\nWe agree with the reviewer and revised the manuscript to start a new paragraph at “Because LRS data allow...”, as suggested.\nReviewer #3 (Remarks to the Author):\nThe genetic basis of autism remains underexplored. In this study, the authors performed genome sequencing of 189 individuals from 51 families with unsolved cases and generated a set of high-quality assemblies. These assemblies were compared with publicly available high-quality assemblies from healthy donors, leading to the discovery of three pathogenic variants in TBL1XR1, MECP2, and SYNGAP1, as well as nine candidate variants, most of which were missed by short-read sequencing. The study is well designed and solid. I have a few concerns as outlined below.\nMajor comments\nMECP2 is structurally complex with diverse SVs and CNV events. I am wondering whether this case also harbors any other events besides the discovered deletion. Moreover, what is the actual deleted sequence in the last exon region, and what is the precise effect of this deletion on the last exon?\nBecause we have nearly complete sequence data, we revisited this patient’s genome and although there are 490 rare SVs that are not observed in controls, none are plausible candidates for Rett syndrome (Supplementary Data 4). The deletion we discovered disrupts the open-reading frame of MECP2, leading to a premature stop codon that truncates the protein by approximately 140 amino acids (as shown below). Similar deletions involving exon 4 have been reported in Rett patients (Hardwick et al. 2007, Poleg et al. 2025); after consultation with our colleague Dr. Zoghbi, we are quite confident that this is the causative variant for this patient.\nWe revised the text to:\n“Two additional pathogenic variants were discovered by LRS among the five autism-diagnosed females with features reminiscent of Rett syndrome. This included an 874 bp de novo DEL in MECP2 in HYZ207_p1, which effectively disrupts the last exon of the gene and introduces a premature stop codon, truncating the protein by approximately 140 amino acids (Fig. 4b). Similar deletions involving exon 4 have been reported in Rett patients (Hardwick et al. 2007, Poleg et al. 2025). This pathogenic variant was previously missed in three rounds of clinical testing, including two gene panel sequencing tests through ARUP Laboratories and Quest Diagnostics, and one test of WES through Ambry Genetics. It was confidently identified in our LRS analysis and subsequently validated using all three sequencing platforms.”\nThe paper includes a large control cohort, which is very good. However, controls appear underrepresented for EUR ancestry (e.g., 39 EUR among 569). Please report the ancestry and sex composition of controls vs. families, and perform ancestry- and sex-matched analyses.\nThe reviewer is correct—the pangenome has generally been biased against Europeans in favor of samples of African origin in order to maximize diversity. In Supplementary Data 1.2, there are a total of 73 Europeans in the full set of 569 controls. For each control, we indicated which nested control sets it belongs to, with the smallest set listed first. For example, a sample present in all three sets is labeled “108controls, 285controls, 569controls” while a sample in the medium and largest sets is labeled “108controls, 285controls” and a sample present only in the largest set is labeled “569controls.” This format allows readers to easily track the inclusion of each sample across the nested control sets. We added a detailed breakdown of ancestry and sex for both controls and autism families as part of our supplement (Supplementary Data 1.3), as shown below.\nAncestry F M total\n108controls 60 48 108\nAFR 30 22 52\nAMR 14 10 24\nEAS 8 7 15\nEUR 3 5 8\nSAS 5 4 9\n285controls 144 141 285\nAFR 56 39 95\nAMR 28 24 52\nEAS 21 38 59\nEUR 19 15 34\nSAS 20 25 45\n569controls 320 249 569\nAFR 115 92 207\nAMR 48 31 79\nEAS 51 52 103\nEUR 43 30 73\nSAS 63 44 107\n189samples 107 82 189\nAFR 2 3 5\nAMR 12 12 24\nEAS 4 1 5\nEUR 87 64 151\nSAS 2 2 4\n87children 56 31 87\nAFR 1 1 2\nAMR 8 4 12\nEAS 2 1 3\nEUR 44 24 68\nSAS 1 1 2\nWe have already performed a sex-matched analysis. To clarify, we moved the statement, \"For autosomal SVs, we retained heterozygous and homozygous SVs present only in the children but not in the controls. For SVs on the sex chromosomes, we performed sex-matched comparisons and filtered SVs seen in controls with the same sex” to “3. Rare SV pool discovery.”\nPer the referee’s request, we performed an ancestry-matched analysis where we selected 38 families of primarily European ancestry (140 individuals, ~74%) and repeated the rare SV discovery using the 73 European controls. From the largest merged set with 569 controls, we extracted 177,522 nonredundant SVs (~40% of the total 445,142 SVs). The ancestry- and sex-matched analysis resulted in an average of 754 rare SVs per child and ~74% of which were filtered out by diverse populations. We incorporated both of the reviewer’s suggestions into the same paragraph and added the requested analysis, with the revised text underlined to indicate changes made in response to Reviewer #1:\n“Increasing the number of controls, especially samples of African origin, nearly doubles the number of SVs from 271,375 (108 controls) to 445,142 (569 controls) nonredundant SVs (Supplementary Data 5-8). We generated a population-level SV reference and compared the MAF distribution between the 108 controls and the largest dataset comprising 569 controls and 102 unrelated parents (n = 1,342). Based on this comparison, we estimated the fraction of variants captured at different MAF thresholds (Supplementary Fig. 12). Notably, only 16.2% of SVs with MAF < 0.1% in the larger dataset were detected in the 108-control set, illustrating the limited sensitivity of the smaller reference for identifying ultra-rare variants. In an ancestry- and sex-matched analysis restricted to individuals of European ancestry (38 autism families and 73 controls), each child carried an average of 754 rare SVs. Focusing on rare variants absent from the largest 569-control set, the number of rare SVs concomitantly drops to 202 events per sample (Fig. 5b), corresponding to ~74% of rare variants being filtered out by the more diverse population controls. This demonstrates the substantial gain in sensitivity for removing common and low-frequency variants when using a larger, heterogeneous reference panel. As a result, 99% of common SVs are excluded per individual forming a more tractable and potentially biomedically relevant set of rare variants for downstream interpretation and enrichment analyses. ….”\nIn Fig. 1, the authors showed the N50 of contigs in autism family assemblies. However, providing NGx curves against GRCh38 and CHM13 would still be valuable. Comparing family assemblies with control assemblies would also be helpful.\nWe added this NGx curve to the Supplementary Figure 2 and mention it in the main text:\n“The resulting assemblies are highly contiguous (average contig N50 of 43 Mbp (Fig. 1c and Supplementary Fig. 2)) and highly accurate (QV=56).”\nSupplementary Figure 2. NG(x) curves show cumulative assembly contiguity across sequencing batches. The NG(x) statistic represents the contig length at which x% of the estimated genome size is covered. Each color corresponds to a distinct sequencing batch: reference genome CHM13 (purple), GRCh38 (black), autism families (yellow), HGSVC controls (blue), HPRC controls (green), and 1KGP ONT controls (pink).\nThe authors claim a “modest rate of pathogenic variant discovery (5.9%).” Did the authors consider evaluating compound heterozygous and oligogenic models (e.g., small variants, or small variants + SVs)? This may help recover additional pathogenic candidates.\nWe did indeed consider both of these as we thought similarly. We found no evidence of compound heterozygotes and the sample size is unfortunately far too small to observe any oligogenic SV burden. Our previous analysis suggests that at least 2000 families (Wilfert et al., 2021) are required to observe a signal in aggregate for protein-disruptive mutations. All the promising anecdotal examples we identified were unfortunately observed in siblings.\nPlease describe unpublished tools (e.g., subseq, BoostSV) and provide versions, key parameters, and workflow availability to enhance reproducibility.\nThe GitHub documents of subseq (https://github.com/EichlerLab/subseq-smk, Ebert et al. 2021) and BoostSV (v1.0, https://github.com/jiadong324/BoostSV) have been updated with more details. The overall workflow, including all the steps and tools used in this study, can be found on GitHub: https://github.com/EichlerLab/asap.\nTo make the parameters and version clearer, we added them to the Methods:\n“2. Callable region evaluation using BoostSV (v1.0). To ensure the SVs fall within confidently callable regions across samples in a single family, we developed a tool, BoostSV (https://github.com/jiadong324/BoostSV), leveraging a machine-learning approach trained on control samples (Porubsky, Dashnow, et al. 2025). This tool evaluates read support, mapping quality, and data quality metrics from alignments surrounding the target SVs in each parent. A quality threshold of 0.5 was applied to obtain the transmission.\n…\n5. Read-based support validation using subseq. To further assess SV transmission, subseq (https://github.com/EichlerLab/subseq-smk, (Ebert et al. 2021)) was used to quantify read support for each SV in parental genomes. A dynamic window size was determined based on the SV size, and the number of reads traversing the window were counted. The {size50_1_1} parameter was used, requiring a minimum of one read supporting the SV while allowing a 50% size deviation penalty.”\nMinor comments\nFig. 2a: Distinguish proband vs. sibling more clearly (shapes or colors).\nThe colors distinguishing probands and siblings have been updated in Fig. 2a to improve clarity.\nFig. 4c,d: Correct “Genehance” → “GeneHancer.”\nThe label has been corrected from “Genehance” to “GeneHancer” in Fig. 4c,d.\nTypo: “TBXLR1” → “TBL1XR1.”\nThe text has been corrected to TBL1XR1. Thank you for catching this.\nI would like to thank the authors for their contribution to the field through this work.\nWe greatly appreciate the reviewer’s detailed and helpful comments.\nReviewer #3 (Remarks on code availability):\nI would appreciate for providing the documents for some unpublished tools.\nUnpublished tools with descriptions have been added to GitHub."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have thoroughly addressed all of my comments. The manuscript has improved substantially and is an important contribution to the field.", "reviewer_3": "I have carefully reviewed the authors’ revisions, and I think all of my previous concerns have been fully addressed. I have no further comments at this time."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have thoroughly addressed all of my comments. The manuscript has improved substantially and is an important contribution to the field.\nWe thank the reviewer for their time and efforts and recognizing the improvements to the manuscript.\nReviewer #1 (Remarks on code availability):\nThe code and documentation for the analysis is very detailed and the authors have done a great job.\nI noticed many links say \"Internal path\" and link to an Eichler Lab site which requires login info, and therefore is not publicly accessible or fully reproducible.\nWe removed any internal use paths/links in the GitHub (https://github.com/EichlerLab/asap) and all code is now publicly accessible to allow for full reproducibility.\nReviewer #3 (Remarks to the Author):\nI have carefully reviewed the authors’ revisions, and I think all of my previous concerns have been fully addressed. I have no further comments at this time.\nWe thank the reviewer for their time and efforts reviewing this manuscript.\nReviewer #3 (Remarks on code availability):\nSource code has been properly deposited in several public repositories. I thank the authors for their efforts.\nWe confirm code is available on GitHub (https://github.com/EichlerLab/asap)."}]

Decision letter