Population-scale Long-read Sequencing in the All of Us Research Program

Participant selection and LRS

We selected 1,027 self-identified Black or African American participants over 20 years of age who met all requirements for sequencing library construction (Methods) from the AoU CDR version 7 release of 245,388 participants recruited from sites across the U.S. (Fig. 1a-d). This focus was motivated by three considerations: (1) Black and African American self-identity reflects a socio-political population construct that captures the lived experiences, social context, and structural factors that shape health outcomes. However, lower representation of individuals who self-identify as Black and African American in genomics research has made it challenging to characterize the full impact of how these factors intersect with genetic variation to influence health outcomes and disease risk in this population; (2) many, but not all, individuals who self-identify as Black or African American have substantial African genetic ancestry, which may enrich certain variants; and (3) populations of African genetic ancestry harbor the highest levels of genetic diversity across geographic populations, which increases the power of variant discovery, particularly for SVs. Approximately 96% of sequenced participants had associated phenotypic data as part of their EHR (general phenotype distribution in Fig. 1e). Importantly, this information was not used to inform selection criteria.

• Download figure
• Open in new tab

Fig. 1: LRS of self-identified Black or African American participants from All of Us.

a, 1,027 participants selected for PacBio HiFi sequencing (black vertical lines) from 245,388 AoU participants (top multicolored bar) with existing short-read data in CDRv7 releases and beyond. Computed genetic ancestry from srWGS data shown for: AFR–African, AMR–Admixed American, EAS– East Asian, EUR–European, MID–Middle Eastern, SAS–South Asian, OTH–other (typically admixed). b, PCA of SNV genetic data from participants selected for PacBio LRS (gray diamonds), overlapping participants with Oxford Nanopore Technologies (ONT) data (red, filled diamonds), and participants from the Human Pangenome Reference Consortium (HPRC) for comparison and evaluation (blue triangles), projected onto the AoU srWGS data. c, Geographic sampling locations for 1,027 long-read participants collected across the continental U.S., shown at the US Census Bureau Census Division level (no LRS participants sampled from Alaska, Hawaii, or other territories). d, Coverage distributions of selected LRS participants, stratified by sequencing technology. e, Mean read lengths of selected LRS participants, stratified by sequencing technology. f, Clinical phenotype distribution in AoU-LR participants based on EHR data. Diagnoses were mapped using the Observational Medical Outcomes Partnership (OMOP) concept IDs and interpreted via standard Observational Health Data Sciences and Informatics (OHDSI) concept mappings. For each disease category, all descendant concept IDs were identified using the OHDSI CONCEPT_ANCESTOR table, and counts reflect the number of unique participants with a diagnosis matching any related concept.

Participant selection for LRS occurred in 2020 before either genotyping array or srWGS data generation for All of Us participants began, precluding the use of either for selection or ancestry confirmation prior to sequencing. Nevertheless, after srWGS data became available, a principal component analysis (PCA) with ancestry labels computationally inferred using a random forest classifier trained on 1000 Genomes Project (1KGP) and the Human Genome Diversity Project (HGDP) autosomal variants from gnomAD v3.1³ confirmed that the selected participants were of majority genetically similar to African reference populations or admixed across reference (Fig. 1b).

For LRS data generation, we used the PacBio (Sequel IIe) platform, which provides highly accurate reads through circular consensus sequencing (CCS) and outperforms other LRS technologies in base-level accuracy²¹. To enable scalable population-level variant discovery, we adopted a mid-pass sequencing strategy targeting ∼8x coverage. Our groups have previously shown that this strategy increases SV discovery sensitivity at substantially lower cost than deep LRS using multiple technologies to completely assemble diploid genomes^20,22. Compared to Illumina srWGS performed on the same individuals, the average length of LRS reads was 115-fold longer (mean=17.5 kbp, Table 1). As a control for variant discovery, a subset (n=50) of participants were also sequenced at 35x sequence coverage using Oxford Nanopore Technologies (ONT) (Table 1, Fig. 1). The combined dataset (termed the “CDRv7 AoU Long-Read dataset”), along with variant calls and local assemblies, is already available to approved researchers via the AoU Researcher Workbench (RW).

Variant discovery and harmonization

We identified small genomic variants (SNVs and indels <50 bp) using the DeepVariant deep learning-based variant caller²³. All SVs (≥50 bp) were detected using an ensemble approach that combined de novo assembly-based discovery via Phased Assembly Variant Caller (PAV)¹² with two read-based algorithms: Sniffles2 and PBSV^21,24.

The ensemble algorithm approach improved sensitivity to diverse variant types, including those detectable by only a single method. However, integrating results from multiple tools, particularly at moderate (∼8x) coverage, poses major challenges for variant filtering and merging while attempting to distinguish true variants from false positives from LRS²⁵. To address this and construct a unified SV callset across the cohort, we developed a custom harmonization pipeline, including filtering, normalization, and merging (Supplementary Fig. 1). Briefly, we clustered and normalized raw calls within each participant using Truvari²⁶, selected for its high recall and minimal over- and under-merging when benchmarked against high-quality diploid assemblies from the Human Pangenome Reference Consortium (HPRC)¹⁴. We then trained a per-participant XGBoost²⁷ classifier to filter likely false positives using genotyping features extracted with Kanpig²⁸ (i.e., allele depth, allele balance, variant length, and multi-caller support). Training labels were derived from Truvari matches to dipcall²⁹ variants on HPRC assemblies, enabling receiver operating characteristic (ROC) curve-based thresholding to approximate target true positive rates (TPRs). We applied thresholds to generate both lenient (TPR ≈ 0.9) and stringent (TPR ≈ 0.7) filtered callsets. Finally, filtered per-participant callsets were merged across the cohort using Truvari, and the resulting cohort-level VCF was regenotyped with Kanpig to maximize sensitivity and recover shared variants across participants (Supplementary Methods).

In our final variant discovery callsets, the number of SNVs was comparable between the PacBio LRS and Illumina srWGS platforms analyzed using 990 overlapping participants (Table 2). As expected, the number of indels was substantially reduced in the PacBio LRS as a consequence of the challenges of sequencing short tandem repeats on different sequencing platforms^17,21. By contrast, there was a dramatic increase in SV discovery ranging from 60% to 290% depending on the stringency level applied to the putative SV callset before comparisons. This result is consistent with prior studies demonstrating the superior sensitivity of long-read technologies for SV detection^12,15,30. We have therefore released two SV callsets—a lenient callset of 1,213,876 SVs maximized for sensitivity and a strict callset of 665,869 SVs optimized for high specificity. Both SV datasets include the HPRC controls and were subjected to similarly rigorous QC steps with all metrics provided (Fig. 2).

• Download figure
• Open in new tab

Fig. 2: Structural variant calling and quality control.

a, Receiver operating characteristic (ROC) curves for per-participant SV filtering performance, evaluated against 47 high-quality diploid assemblies from the HPRC. True positive rate (TPR) and false positive rate (FPR) are shown under strict and lenient filtering thresholds, stratified by variant type and genomic context (tandem repeat [TR] vs. non-TR regions). One participant is highlighted for clarity. b, De Finetti diagram depicting ratios of homozygous reference, heterozygous, and homozygous alternate genotypes at each SV site in strict (left, n=616,411 SVs) and lenient (right, n=1,157,726 SVs) callsets, as well as the fraction of SV sites in Hardy-Weinberg equilibrium. Every site is assumed to be biallelic. c, Distribution of SV lengths for insertions and deletions, stratified by callset. Peaks corresponding to known mobile element families (Alu, SVA, and L1) are labeled. d, Per-participant SV counts across different groups: 47 HPRC samples used for comparison, 1,027 AoU participants with long reads (subsets with other data types indicated by horizontal bars; PacBio: purple; ONT: blue; Illumina: red). Purple and blue curves represent expected variant counts based on raw Sniffles2 calls on ∼30x PacBio data from HPRC samples and ∼35x ONT data from AoU, respectively, illustrating expected SV yield at higher coverage. e, Histograms of unfiltered allele frequencies from the AoU+HPRC panel (x-axis) and 2,504 unrelated 1KGP short-read samples imputed against it (y-axis), for SV-length insertion/deletion (left/right, Pearson correlation coefficient 0.84/0.94) bubbles. f, Number of imputed AoU LRS SVs per participant across five continental population groups in 2,504 unrelated 1KGP short-read samples. g, UpSet plots show the number of protein-coding genes (including intronic loci) intersected by an SV across the five continental groups in 1KGP.

Evaluation and quality control

We applied a series of quality controls to assess the accuracy of both the lenient and strict SV callset. First, to evaluate Mendelian consistency as a function of sequence coverage, we selected three trios with high-coverage HiFi sequencing data¹¹. Each participant’s LRS data were downsampled to 8x, 16x, and 25x coverage using GATK DownsampleSam³¹. For each coverage level for each participant, we applied our full pipeline, merging and regenotyping separately for each coverage level. Using the GRCh38 reference, we found a mean Mendelian discordance rate of 6.4% for the TPR=0.9 leniently filtered callset at 8x coverage, with the rate dropping to 3% at 25x coverage (Supplementary Fig. 2). Notably, this discordance rate from our analysis of 8x coverage data was consistent to that reported cohort-wide by the Human Genome Structural Variation Consortium (HGSVC) using high-coverage LRS (6%), further validating the utility of the mid-pass approach¹². The strictly filtered TPR=0.7 callset had a similar discordance rate of 6.8% at 8x coverage, dropping to 2.9% at 25x coverage. Discordance rates were similar but slightly elevated on the T2T-CHM13 reference (Supplementary Fig. 3). Both references exhibited substantially improved concordance rates for variants outside of tandem repeats (Supplementary Fig. 4).

Both datasets show the expected SV size distribution with notable peaks at 300 bp and 6 kbp, consistent with SINE and LINE mobile element insertions (MEIs)¹² as well as an increasing discovery curve for insertions and deletions (Fig. 2c,d) consistent with expectations. After collapsing MEIs at the same locus, we detected a total of 28,101 (24,663) Alu, 3,760 (2,873) LINE-1, and 1,750 (1,253) SVA insertions in the lenient (and strict) callsets. Per participant, we detected a median of 1,620 (1,500) Alu, 166 (140) LINE-1, and 83 (68) SVA insertions in the lenient (and strict) callsets (Supplementary Fig. 5a), with expected length distributions for each class (Supplementary Fig. 5b). For comparison, using high-coverage LRS, the HGSVC¹¹ detected 1,585 Alu, 212 LINE-1, and 98 SVA insertions per participant of African genetic ancestry. We observed extensive allelic heterogeneity in SVA length (Supplementary Fig. 5c,d), with the top-most variable SVAs having length differences exceeding 2 kbp across participants. These length differences were predominantly driven by copy number changes in the central CG-rich variable number tandem repeat (VNTR) region (Supplementary Fig. 5d).

As a final validation approach to assess precision, we selected 50 AoU participants that overlapped with the PacBio data analysis and sequenced them with ONT (Fig. 1, Table 1). We utilized Kanpig²⁸ to genotype the HiFi-based SV calls on GRCh38 and T2T-CHM13. We observe high support for the lenient SV callset with 95.5% precision for GRCh38 and 90.5% precision for T2T-CHM13. Subsequent genotyping (see below) against the 1KGP samples confirmed that the majority of common SVs are observed among all population groups, but that SVs specific to African genetic ancestry are the second most abundant (Fig. 2e-g).

Creation of an LRS reference panel for srWGS genotyping, phasing, and imputation

To build a phased, haplotype-resolved reference panel optimized for SV imputation into srWGS data, we integrated the AoU and HPRC LRS callsets (spanning SNVs, small indels, and SVs) using the stringent (TPR=0.7) SV callset to minimize false positives. We first performed read-based physical phasing of both callsets jointly using HiPhase³², followed by postprocessing to remove (1) singleton variants, (2) duplicate variants between SV and SNV callsets, (3) short variants with AF < 0.5% and >25 kbp from any SV, and (4) positionally conflicting variant calls yielding inconsistent haplotypes. We then applied SHAPEIT4 for statistical phasing and imputation and postprocessed the phased haplotypes by concatenating adjacent alleles to define non-overlapping “bubble” variant sites, generating a pangenome graph representation (see Supplementary Methods for further details)³³. Finally, using the KAGE³⁴ genotyping and GLIMPSE³⁵ phasing and imputation methods, we called variants in this LRS-derived reference panel from srWGS case participants.

The resulting phased, imputed, and haplotype-resolved SNV/indel/SV reference panel of 1,074 AoU+HPRC participants covers the autosomal chromosomes and consists of 19,942,647 variant sites in bubble representation, which contain 30,918,204 variant bubble alleles; 189,382 of these bubbles are SV-length and 74,550 are further multiallelic. For comparison, the 44-sample HPRC panel contains 21,304,582 (15,684,910) bubble sites in the autosomes, of which 71,873 are SV-length and 55,156 are further multiallelic. After decomposing bubble alleles back to their original representations, the panel contains 27,143,444 distinct variant alleles, including 770,473 SVs (653,740 insertions and 117,698 deletions).

Both LRS-panel and srWGS-case pipelines were validated using HPRC assembly-based dipcall truth with the vcfdist alignment-based benchmarking method³⁶. Performing a leave-one-out experiment with 40 HPRC samples on a single evaluation chromosome, the LRS panel (unfiltered srWGS cases) achieved 94.0 ± 1.1% (95.5 ± 1.0%) precision and 91.4 ± 1.3% (87.9 ± 2.0%) recall for SVs with length <5 kbp in dipcall-confident, non-repeat/homopolymer regions. This accuracy is comparable across both pipelines, which demonstrates the additional power that prior knowledge derived from LRS resources can unlock in readily available srWGS data. Further details of the panel creation and case pipelines, as well as our validation procedures, are provided in the Supplementary Methods.

Characterization of SVs enriched in African genetic ancestry

Using our AoU+HPRC reference panel, we genotyped and imputed variants in srWGS data using all 3,202 individuals from the 1KGP³⁷, generating a comprehensive SV callset across five continental populations (Fig. 2f). 1KGP samples with African genetic ancestry showed the highest number of SVs per individual compared to other groups, consistent with previous findings of greater genetic diversity among individuals of African genetic ancestry³⁸. A PCA of the SV callset achieved clear population stratification, closely matching patterns observed in the published 1KGP SNV dataset³⁹ and supporting the reliability of our SV genotypes for population genetic analyses (Supplementary Fig. 6a,b).

We next intersected SVs with coding regions of genes and regulatory elements (Figs. 2g, Supplementary Fig. 6c, Supplementary Table 1). The majority of SV-overlapping elements (40.3% of coding; 43.7% of regulatory) were shared across all populations, indicating that our panel captures globally common variation. The second-largest category included elements specific to African genetic ancestry (15.1% coding; 17.1% regulatory), reflecting both the elevated diversity within genomes of African genetic ancestry and their strong representation in the AoU cohort.

The majority of the imputed SVs (315,307 of 541,049; 58.3%) had higher allele frequencies in African genetic ancestry than in non-African genetic ancestry. To further assess population differentiation, we calculated the fixation index (F_st) between African and non-African groups and identified 8,249 SVs with F_st ≥ 0.15, indicating strong enrichment in African genetic ancestry (Supplementary Fig. 6d; Supplementary Table 1). These variants were distributed across the genome and intersected with 776 high-priority disease-associated genes (Supplementary Fig. 6e). One notable example is a 61 bp deletion in BARD1 (African ancestry AF = 0.12; non-African ancestry AF = 0.002; F_st = 0.19; Supplementary Fig. 6f), a gene linked to elevated breast cancer risk⁴⁰. We also observed SVs specific to African genetic ancestry in genes associated with conditions disproportionately affecting individuals of African ancestry, such as APOB and MYH7, both of which are implicated in cardiovascular disease⁴¹ (Supplementary Table 1).

Characterization of clinically relevant and population-enriched SVs

Because this dataset represents the largest LRS survey to date of participants self-identifying as Black or African American, we conducted an in-depth analysis of SVs, which have greater potential for novel allelic variant discovery and outsized impact on health⁴². We first investigated known disease-associated loci (such as those involving triplet repeat instability and complex pharmacogenomic regions like CYP2D6) to assess the added value of LRS in identifying novel variants of potential clinical significance. Second, we focused on SVs that are enriched in groups self-identifying as Black or African American and absent from existing variant catalogs. These variants, more likely to have larger effect sizes, were tested for association with gene expression and disease. Together, these analyses provide a valuable resource for understanding genetic diversity and improving health outcomes in populations that continue to face structural barriers to receiving high-quality tailored healthcare.

Tandem repeat loci and disease

Pathogenic expansions of short tandem repeats (STRs) have been implicated in over 50 inherited diseases⁴³, where both tract length and purity of the repeat influence instability. However, many of these disease-associated loci have been difficult to fully characterize by srWGS or PCR due to their size or GC-rich content. To address this, we performed detailed characterization of two well-known triplet repeat loci: the CGG repeat in the 5’ untranslated region (UTR) of FMR1, associated with Fragile X syndrome (FXS) (Fig. 3a,b), and the CAG repeat in the coding region of HTT, associated with Huntington’s disease (Fig. 3d,e). Repeat sequences were derived with TRGT⁴⁴ for a total of 1,753 (sex-filtered) haplotypes at FMR1 and 1,944 at HTT. In both cases, we observed the most common haplotypes reported in the literature, including characteristic AGG interruptions occurring every 9-10 CGG units within the FMR1 repeat (Fig. 3a,d)^45,46.

• Download figure
• Open in new tab

Fig. 3: Expanded triplet repeat loci and structure of disease-associated loci.

a, Most common and b, longest FMR1 repeat alleles. Sequences were extracted from 1,753 sex-filtered long-read haplotypes with TRGT and plotted according to repeat and interruption patterns. AGG interruptions are included in repeat counts. Arrows indicate unstable transmission risk alleles with 25 to 33 (light blue) or at least 34 (dark blue) uninterrupted CGG repeats, premutation alleles with at least 55 repeats (red), and an allele with five interruptions (yellow). c, Length distributions indicating the number of repeat units in extracted TRGT sequences for FMR1 and HTT. d, Most common and e, longest HTT repeat alleles. Sequences were extracted from 1,944 long-read haplotypes with TRGT and plotted. Blue arrows indicate intermediate-length unstable transmission risk alleles with 27 to 35 CAG repeats, and red arrows indicate potential low-penetrance Huntington’s disease alleles with at least 36 repeats. f, Boxplot of the distribution of tandem repeat length variation values (measured by median absolute deviation) for the entire catalog (first box), a set of loci specifically expanded in humans relative to nonhuman primates (second box), the CODIS loci (third box), and the known pathogenic loci (fourth box). g, 5’ UTR CGG repeats and h, protein-coding CAG repeats with median (50^th percentile) repeat length for each locus plotted against the 99^th percentile repeat length, both measured in trinucleotide units. At each locus, the longest pure (uninterrupted) repeat was segmented. A y=x trend line is included for reference, indicating where the median and 99^th percentile repeat lengths are equal. Out of 4,115 CGG loci and 1,483 CAG loci, 11 and 12 disease-associated loci are highlighted in red/green with annotated genes, while non-disease-associated loci are shown in purple/orange.

Because expanded alleles are associated with disease/disorder risk, we specifically focused on the tail-end of each length distribution (Fig. 3b,e). In FMR1, we predicted ten of the longest alleles to be potentially unstable due to the loss of the 3’ AGG interruption, extending the length of the longest pure CGG repeat tract^47,48. We also identified four premutation-sized alleles (≥55 repeat units), of which the longest occurred in a participant for whom we estimate an ∼83-93% probability of expansion to a full mutation in offspring on maternal transmission. This estimate was phenotypically supported by the participant’s EHR, which describes a child with symptoms commonly associated with FXS.^49,50. In HTT, we observed 12 potentially unstable intermediate-length alleles and two reduced-penetrance alleles (≥36 CAG repeat units), corresponding to a carrier allele frequency of ∼0.1% (2/1,944 haplotypes). We identified several noncanonical and rare interruptions in the FMR1 and HTT repeats (Supplementary Note 1). In addition, we identified a participant with an ATXN3 repeat expansion exceeding the pathogenic, high-penetrance threshold of >55 CAGs for spinocerebellar ataxia type 3 (Supplementary Fig. 7) with a diagnosis of cerebellar ataxia recorded in their EHR, supporting this genotype⁵¹.

To identify potentially novel loci associated with triplet repeat instability, we applied a catalog-driven approach⁵² to GRCh38-aligned participants and identified several classes of tandem repeats exhibiting significantly greater length variability than the catalog as a whole (Fig. 3f). These included loci known to be expanded in humans relative to nonhuman primates, autosomal CODIS loci used in forensic analyses, and most recently described pathogenic tandem repeats^53–56. We further characterized this set by plotting the 99^th percentile length of the longest pure repeat segment against the median (50^th percentile) longest pure repeat segment length for every 5’-UTR CGG repeat (Fig. 3g) and every coding-region CAG repeat (Fig. 3h) in the catalog. Consistent with the initial observation, known pathogenic repeats for both groups clearly segregate from the majority of loci, occupying the right tail of the length variability distribution. In addition to these known loci, we identified 13 CAG repeat loci and 15 CGG repeat loci that exhibit similar distributions of variation to the known pathogenic loci. We hypothesize that these loci represent strong candidates for further study as potentially pathogenic STRs. A companion manuscript further explores these observations and candidates⁵⁷.

CYP2D6 variation

Given the crucial role of CYP2D6 in the metabolism of commonly prescribed medications and the known difficulty of resolving this locus using srWGS due to complex structural variation^58,59, we evaluated whether LRS assemblies could provide higher-resolution genotyping. We first developed a k-mer-based method to determine the identity and copy number of CYP2D6 and its paralog CYP2D7 (Methods) and applied it to 1,266 AoU haplotypes (62%) contiguously assembled over the full CYP2D6-7 region, following HPRC validation. This analysis revealed 13 distinct structural configurations (Fig. 4a), including 38 CYP2D6 duplications, 87 full-gene deletions, and 18 hybrid alleles corresponding to four structures: CYP2D6-CYP2D7::6 (n=6), CYP2D6-CYP2D7::6-CYP2D7 (n=5), CYP2D6::7 (n=4), and CYP2D7::6-CYP2D7 (n=3). We next identified star alleles, defined as haplotypes associated with enzymatic activity levels. Consensus calls from four validated short-read genotypers supported duplications for 28/33 participants and deletions for 56/75 participants, and classified seven hybrid alleles as *68+*4. Among the five duplications with discordant calls, three had indeterminate short-read consensus, and the Pangu genotyper from PacBio-supported duplications in three cases⁶⁰. Among 19 deletion discrepancies, one was indeterminate by srWGS, and Pangu called the *5 full-gene deletion allele in one participant and reported only a single haplotype in 14 others. srWGS calls identified four of ten hybrid discrepancies as potential *4 hybrids, while Pangu called three as the *4.013 exon 9 conversion.

• Download figure
• Open in new tab

Fig. 4: CYP2D6-7 haplotype structure and content.

a, Copy number variation and hybrid structures for 1,266 AoU-LR Phase 1 haplotypes. Hifiasm assemblies with a single contig over the full CYP2D6-7 locus were divided into 100 bp k-mers and colored by the best-mapping reference gene. An example plot for each structure is shown along with the number of haplotypes the structure was seen in. b, All (left) and novel (right) missense variants. Predicted haplotype coding sequences were compared with the GRCh38 coding sequence, which was substituted with changes for the assigned (sub)star allele to find novel variants. Variants are colored by AlphaMissense classification, and asterisks indicate presence in the srWGS and LRS AoU Hail matrix tables.

We developed an LRS star-allele inference pipeline that selects the closest-matching reference allele by maximizing variant matches and alignment score. Applying this to the AoU LRS assemblies, we observed 91% concordance with srWGS star allele calls. The most frequent unambiguous disagreements involved differences in copy number or hybrid structure (Supplementary Table 2). We restricted our analysis to 255 participants with a canonical CYP2D6-7 structure, for which diplotypes (and thus metabolizer status) could be fully inferred. Among these, 172 were categorized as CYP2D6 normal metabolizers, 77 as CYP2D6 intermediate metabolizers, and 6 as CYP2D6 poor metabolizers. To examine novel mutations beyond the SNVs and indels found in the closest-matching star allele, we predicted the coding and protein sequence for each haplotype and compared them to the reference translation. Among 861 haplotypes with unambiguous assignments and a single predicted transcript, 13 carried missense variants not present in the corresponding suballele reference. Of these, two were classified by AlphaMissense as likely pathogenic and two as having ambiguous pathogenicity (Fig. 4b)⁶¹. Only one of these four variants was also supported by the srWGS data.

Functional and clinical interpretation of AoU SVs

Using the stringent SV callset, we evaluated the ability of the LRS dataset to detect SVs in other biologically and clinically significant regions. We identified 251,723 SVs mapped to protein-coding loci, including intronic regions. Among these, 68.2% of deletions and 63.4% of insertions were not detected by srWGS, yet intersected 8,219 OMIM genes⁶², 2,801 genes prioritized for inherited diseases⁶³, and 49 from the American College of Medical Genetics and Genomics (ACMG) dataset⁴¹, suggesting the considerable value of LRS data for improving the discovery of medically relevant variants (Supplementary Table 3). LRS also enabled the resolution of SVs in complex, highly polymorphic regions. For example, 723 LRS-unique SVs were identified within the major histocompatibility complex, including an 816 bp deletion within HLA-A (Supplementary Fig. 8a).

We then prioritized the stringent SVs using CADD-SV⁶⁴, which classified over 99% (574,531 SVs) as benign (Fig. 5a). However, 473 SVs received PHRED-scaled scores ≥20, placing them among the top 1% of predicted deleterious variants relative to a healthy population reference (Supplementary Table 3). As expected, most putatively deleterious SVs were rare within the AoU cohort (Fig. 5a, Supplementary Fig. 8b), particularly those with scores ≥30, consistent with purifying selection acting against highly damaging variants^5,65,66. We further compared the AoU long-read SVs with recently published LRS datasets, including the 1KGP ONT callset, HSGVC, and HPRC^12,14,67. SVs absent from these external resources generally had a lower allele frequencies in the AoU cohort than shared SVs (Fig. 5a, Supplementary Fig. 8b), suggesting they are both rare and potentially population-specific. Among 273 high-priority, previously unreported SVs, 172 overlapped 170 medically relevant genes⁶², 73 intersected 69 high-priority disease-related genes⁶³, and 15 affected 14 genes associated with cancer risk⁶⁸ (Supplementary Table 3). These SVs represent compelling candidates for studying disease risk, especially in individuals of African genetic ancestry.

• Download figure
• Open in new tab

Fig. 5: Functional impact of AoU SVs.

a, Distribution of CADD-SV score for shared SVs in the AoU strict cohort. Each data point represents an SV. The x-axis indicates the number of participant containing the SV, and the y-axis shows the PHRED-scaled CADD score. “Known” SVs are those identified in at least one of the 1KGP-ONT, HGSVC, or HPRC datasets^12,14,67, while “AoU detected” SVs are those absent in those datasets. The dashed gray horizontal line denotes the score threshold above which SVs are considered likely pathogenic. b, Relationship between the best-fit slope (β) derived from OLS regression and gene-level q-values. The eGenes shown are medically significant and listed in OMIM⁶², with high-confidence q-values (<1 × 10^-60). c, Distribution of genotypes and gene expression values in 731 participants for the BID-associated deletion, with a q-value of 5.92 × 10^-13. d, Manhattan plot of the 322 bp deletion and nearby SNVs, with log₁₀ P values. The deletion is the top variant, and points are colored by their LD (r²) with this SV. The enhancer overlapping the SV is shown in green. e, Histogram showing the number of SVs in LD (r² ≥ 0.5) with SNVs from the GWAS catalog. SVs located within trait-associated genes are shown in orange, and those outside these genes are shown in blue. Bars with diagonal hatching (///) indicate SVs linked to disease- or disorder-related traits. f, LD heatmap of an insertion (marked in blue) along with nearby GWAS variants and SVs located within ±100 kbp of the insertion. GWAS variants in LD are labeled with their IDs. GWAS-associated genes are shown below the heatmap.

We then leveraged matched EHR data to explore the clinical outcomes of participants carrying rare, high-impact SVs not found in previous LRS datasets. For instance, this highlighted a 52 bp insertion in SLC2A12 (CADD-SV score: 51) found in a participant diagnosed with Type 2 diabetes and hyperglycemia (Supplementary Fig. 8c, Supplementary Tables 3 and 4). As SLC2A12 encodes GLUT12, a glucose transporter implicated in metabolic regulation^69,70, this suggests a potential link to disease. In addition, a 536 bp insertion in TOX (CADD-SV score: 43), a gene involved in asthma pathogenesis⁷¹, was found in a participant with asthma (Supplementary Fig. 8d, Supplementary Tables 3 and 4). Several high-impact SVs were predicted to disrupt protein structure. One example is a 226 bp insertion in PSIP1 (CADD-SV score: 31), associated with hearing loss⁷², which introduced a premature stop codon and was predicted by AlphaFold to truncate the protein (Supplementary Fig. 8e,f, Supplementary Tables 3 and 4)⁷³. Notably, the participant carrying this variant presented with sensorineural hearing loss. In another case, a 4,025 bp deletion (CADD-SV score: 21.3) removed three exons of TARS2 in a participant with chronic kidney disease, consistent with prior studies linking TARS2 mutations to this condition⁷⁴ (Supplementary Fig. 8g,h, Supplementary Tables 3 and 4).

Measuring effects of SVs on gene expression and complex traits

To evaluate the regulatory impact of SVs, we performed expression quantitative trait loci (eQTL) mapping using the imputed SV callset from 1KGP and paired gene expression data from the MAGE dataset⁷⁵. This analysis identified 3,835 significant SV-eQTLs (q-value < 0.05), involving 2,511 medically relevant OMIM genes, 747 high-priority disease genes, and 117 cancer-associated genes, highlighting potential mechanisms through which SVs contribute to disease risk (Fig. 5b, Supplementary Table 1). Notably, 889 of the eGenes were absent from two earlier long-read SV-eQTL analyses involving 31 diverse long-read genomes and 65 samples from the 1KGP-ONT dataset^67,76. The majority (84.6%) of SVs in SV-eQTL pairs were common across all populations (AF ≥ 0.05), allowing robust cross-population analysis of regulatory effects (Supplementary Fig. 9a). Consistent with the AoU LRS reference panel was enriched for participants genetically similar to African reference populations, we further identified 34 statistically significant SV-eQTLs in which the SVs were common in individuals of African genetic ancestry (AF ≥ 0.05) but rare in others (AF < 0.05), indicating newly identified regulatory factors associated with African genetic ancestry (Supplementary Table 1). For example, a 512 bp insertion (AFR AF=0.56; Non-AFR AF=0.02, F_st=0.65) was associated with expression changes in SHPK, a gene involved in sedoheptulokinase deficiency and broader metabolic regulation (Supplementary Fig. 9b-c, Supplementary Tables 1 and 5).^77,78

To identify SVs with likely causal effects on gene expression, we applied the fine-mapping tool CAVIAR⁷⁹ to quantify the posterior probability of causality while accounting for linkage disequilibrium (LD) with SNVs within a 1 Mbp window. This analysis identified 144 SV-eQTLs in which the SV had the highest CAVIAR score compared to the flanking SNVs for the eGene. One notable example is a 322 bp deletion within an enhancer upstream of BID, a pro-apoptotic BCL-2 family member involved in cancer progression and immune regulation^80,81. SV carriers showed decreased BID expression, consistent with the expectation that disruption of enhancer activity would alter gene expression (Fig. 5c, Supplementary Fig. 9d, Supplementary Tables 1 and 5).

This deletion is the lead variant, with a 386-fold stronger statistical significance than the strongest nearby SNP with a concordant direction of effect, indicating the importance of SVs in gene regulation (Fig. 5d, Supplementary Fig. 9e).

SVs in LD with GWAS-significant SNVs

To explore the relationship between SVs and complex traits, we intersected imputed SVs with the NHGRI-EBI GWAS Catalog⁸², a comprehensive repository of genome-wide association study (GWAS) summary statistics. We identified SNV-SV pairs in moderate LD (r² ≥ 0.5) where the SNV had genome-wide significance. This analysis revealed 6,019 traits associated with 5,549 SVs, including 1,103 traits related to diseases or disorders (Fig. 5e; Supplementary Fig. 10a, Supplementary Table 1). For instance, a 321 bp insertion in TSPAN8, a negative regulator of airway mucin secretion⁸³, showed high LD with the asthma-associated rs11178649 (β = −0.04, p = 1 × 10^-10, r² = 0.69; Supplementary Fig. 10b). Using matched EHRs from our AoU cohort, carriers of this SV had reduced asthma risk (odds ratio [OR] = 0.46, 95% confidence interval [CI]: 0.25−0.87, p = 0.01; Supplementary Tables 1 and 5), supporting a potential protective role.

We next evaluated 677 SVs in high LD (r² ≥ 0.8) with SNVs implicated in 537 medically relevant traits. One example is a 1,375 bp deletion in LD with rs16856925, an SNV linked to breast cancer risk (OR = 1.19, p = 7 × 10^-10, r² = 0.93; Supplementary Fig. 10c). This SV was also significantly associated with increased breast cancer risk in the AoU cohort (OR = 2.91, 95% CI: 1.21−7.01, p = 0.02; Supplementary Tables 1 and 5). Furthermore, some SVs exhibited potential pleiotropic effects, a phenomenon previously observed in a multi-trait SNP-based PheWAS⁸⁴. For instance, a 63 bp insertion at position 75,493,495 on chromosome 7 was in strong LD (maximum r² = 0.88) with six GWAS SNVs (rs10953513, rs17207196, rs58862095, rs62477684, rs62477735, rs794377) associated with cardiovascular and metabolic traits, including atrial fibrillation, hypertension, obesity, and metabolic syndrome (Fig. 5f, Supplementary Fig. 10d). Interestingly, the EHR data from SV carriers revealed elevated risks of hypercholesterolemia (OR = 4.73, 95% CI: 1.48−15.07, p = 0.009) and ventricular tachycardia (OR = 3.35, 95% CI: 1.07−10.52, p = 0.039; Supplementary Tables 1 and 5), which are closely related to the GWAS-identified traits, suggesting this SV may contribute broadly to cardiometabolic disease.

Phenome-wide SV analysis identifies disease risk factors

To further explore the clinical impact of SVs, we expanded our analysis to the srWGS and EHRs of 10,000 self-identified Black or African American participants in the All of Us cohort (CDRV7) (Supplementary Fig. 11a). For this, we performed genotyping and imputation based on the long-read-derived AoU reference panel, along with additional quality controls to exclude low confidence variants (Methods, Supplementary Fig. 11b-e). Using these results, we performed a phenome-wide association study (PheWAS) and identified 291 SV-disease associations (p < 1 × 10^-5) across the genome spanning 226 distinct conditions. (Fig. 6a-c; Supplementary Fig. 12; Supplementary Table 6). Of these associations, 148 (50.9%) involved SVs absent from the matched short-read callset derived from 990/1,027 of the same participants, indicating the importance of LRS in capturing clinically relevant variants.

• Download figure
• Open in new tab

Fig. 6: Genome-wide linkage and trait associations of SVs.

a, Genotyping and imputation with the long-read panel. Small variants and SVs are identified from 1,074 long-read participants (AoU + HPRC) and then genotyped and imputed in 10,000 short-read self-identified Black or African American participants from All of Us biobank, with the total number and allele frequency distributions of SNVs and SVs shown. b, Disease phenotypes are extracted from EHRs of the same participants and grouped into 11 categories. Conditions belonging to multiple categories are assigned to all relevant groups. c, Comparison of association significance between each SV and the strongest nearby SNV (within ±100 kbp) for the same phenotype (central panel). Points are colored by LD (r²) between the SV and SNP. Circles indicate genic SVs (overlapping gene bodies), and triangles indicate intergenic SVs. The top and right panels display stacked histograms of −log₁₀(p-values) for SNP- and SV-disease associations, respectively. Genic SVs are shown in purple; intergenic SVs are shown in light blue. d, Distribution of genome-wide SV-disease associations with large effect sizes (odds ratio [OR] ≥ 2.5 or ≤ 0.4) in which the SV is the lead variant. Orange indicates associations with SVs discovered in 1,027 AoU and 47 HPRC long-read samples and subsequently imputed in 10,000 AoU short-read participants, while blue represents associations with SVs identified in the short-read callset derived from the same 1,027 AoU participants. e, Relationship between allele frequency (log scale) and relative risk (OR > 1) for genome-wide disease-associated SVs with stronger signals than nearby SNPs. Circles indicate SVs within gene bodies, and triangles indicate intergenic SVs. Associations specific to African genetic ancestry are shown in yellow and all others in gray. Point size reflects association significance. f, Manhattan plots showing a 50bp deletion linked to hypertensive heart failure. The insertion is the lead variant, and surrounding points are colored based on their LD (r²) with the SV. g, Relationship between allele frequency (log scale) and relative risk (OR > 1) for coding-region disease-associated SVs with stronger signals than nearby SNPs. Associations specific to African genetic ancestry are shown in yellow and all others in gray. Point size reflects association significance. h, Manhattan plots showing a 200 bp insertion associated with atelectasis. The insertion is the lead variant, and surrounding points are colored based on their LD (r²) with the SV.

Our initial set of associations used a screening threshold of significance so that we could compare the relative significance of SNVs and SV across a broad range of traits. To prioritize likely causal variants, we conducted fine-mapping to estimate the posterior probability that each SV contributed to the disease risk, accounting for LD with nearby SNVs within ±100 kbp. This analysis identified 191 (65.6%) SV-disease pairs across 160 traits where the SV was the most likely causal variant within the locus (Fig. 6c, Supplementary Table 6). On average, each participant carried 9 strong associations (OR ≥ 2.5 or ≤ 0.4) in which the SV was the lead variant, compared with 4 in matched short-read SV callset, indicating the important contribution of SVs to participant disease risk, particularly observable with long-read sequencing (Fig. 6d). We further observed that common, low-penetrance SVs were generally associated with more modest relative risks than rare SVs, consistent with evolutionary constraints on high-impact variants. Of note, associations specific to the African genetic ancestry group tended to exhibit larger effect sizes (median AFR-specific = 5.85, median cross-ancestry = 2.72; p = 1.02 × 10^-8) and lower allele frequencies (median AFR-specific = 0.0004, median cross-ancestry= 0.028; p = 8.68 × 10^-11), which may represent high-risk variants enriched in individuals of African genetic ancestry (Fig. 6e). In addition, 98 (51.3%) of these lead SV associations contained variants not detected in the paired short-read callset. For example, a 50bp deletion in PARVB, which encodes β-parvin, a protein that regulates cardiomyocyte shape and sarcomere assembly, was associated with hypertensive heart failure (OR = 1.87; 95% CI: 1.52–2.31; p = 2.32 × 10^-8)(Fig. 6f, Supplementary Fig. 13, Supplementary Tables 6 and 7)⁸⁵.

Of particular interest are SV-disease associations where SVs intersect and potentially disrupt the coding sequence of genes. Given the potentially severe functional consequences of such disruptions, we applied a relatively permissive threshold (p < 1 × 10^-3) and identified 114 associations. Among them, 67 were linked to SVs absent from the short-read callset, further demonstrating the potential of LRS in PheWAS. In 62 of the 114 associations (54.3%), SVs showed stronger association signals than nearby SNVs within ±100 kbp, including 40 involving SVs not detected in short-read data (Fig. 6g, Supplementary Fig. 14, Supplementary Table 6). For example, a 200 bp insertion in CLCN7 was associated with atelectasis (OR = 2.60; 95% CI: 1.62−4.18; p = 3.4 × 10^-4) (Fig. 6h, Supplementary Fig. 15a-d, Supplementary Tables 6 and 7). This variant introduces a premature stop codon, resulting in a truncated protein lacking key functional domains, as predicted by AlphaFold (Supplementary Fig. 15e). CLCN7 mutations have been reported to cause severe alveolar atelectasis in mice, providing mechanistic support for this finding in humans⁸⁶. Additionally, we found a 190 bp insertion in COL14A1 associated with periodontal disease (OR = 4.97; 95% CI: 2.27−10.88; p = 8.9 × 10^-4) (Supplementary Fig. 16a-e, Supplementary Tables 6 and 7). Prior SNP-based GWAS have also identified an association between COL14A1 and periodontal disease⁸⁷, although chiefly through intronic variants potentially altering splicing rather than a direct protein coding change. Here the insertion disrupted the coding sequence and caused a truncated protein lacking a segment of the collagenous triple-helical domain and the entire C-terminal tail, likely disrupting protein function (Supplementary Fig. 16f).

Population and participants

Library preparation

Read mapping and assembly

Alignment and metrics calculation

We aligned reads to the GRCh38_no_alt and CHM13v2.0 reference sequences. For PacBio HiFi reads, we applied the following pbmm2 v1.12.0 command:

• Download figure
• Open in new tab

For ONT data, we applied the following minimap2^93,94 v2.24 command:

• Download figure
• Open in new tab

The MD tag encoding mismatched and deleted reference bases were computed post-alignment using the samtools calmd command. We computed genome-wide and per-chromosome coverage metrics using MosDepth, and read length/sequence identity metrics using NanoPlot and custom scripts.

De novo assembly

Diploid genome assemblies were generated using hifiasm v0.16.1, which constructs haplotype-resolved contig graphs from PacBio HiFi reads. Assemblies were run per sample using 32 threads and default parameters optimized for CCS reads. The command used was:

• Download figure
• Open in new tab

This produced assembly graphs in GFA format: a primary contig graph (<prefix>.bp.p_ctg.gfa) and per-haplotype contig graphs (<prefix>.bp.hap*.p_ctg.gfa). GFA files were converted to FASTA format using the following awk command:

• Download figure
• Open in new tab

For each haplotype-specific contig set, individual FASTA files were generated:

• Download figure
• Open in new tab

Variant detection

Small variant detection

Small variants, defined as single-nucleotide variants (SNVs) and insertions or deletions (indels) shorter than 50 bp, were detected using DeepVariant v1.4.0. BAM files containing aligned PacBio HiFi reads were used as input, along with the GRCh38 reference genome. DeepVariant was run with the PacBio model type and the following parameters:

• Download figure
• Open in new tab

The output VCF and gVCF files were retained for downstream genotyping and joint variant calling. Genomic positions with insufficient coverage or low confidence were annotated accordingly. Quality metrics were computed using standard tools (bcftools, vcftools) to assess depth, transition/transversion ratio, heterozygosity rate, and call concordance in samples with matched short-read data.

Joint small variant calling with GLnexus

Joint genotyping of small variants across samples was performed using GLnexus v1.4.1. DeepVariant gVCF files for each sample were first compressed and indexed with bgzip and tabix to prepare them for cohort-level processing. We used the --config DeepVariantWGS preset, which is optimized for whole-genome sequencing data generated by DeepVariant, and modified parameters for compatibility with PacBio HiFi reads.

The following command was used for joint genotyping:

• Download figure
• Open in new tab

Structural variants (SVs)

SVs ≥ 50 bp were identified using an ensemble of two reference-based algorithms, PBSV and Sniffles2, and one assembly-based algorithm, PAV. Each caller was run independently per chromosome to parallelize processing and manage computational resources efficiently. Per-chromosome outputs were later merged to generate per-sample VCFs for downstream integration and harmonization.

PBSV (v2.6.0) was run in discovery and calling mode using the --tandem-repeats flag to improve sensitivity for repeat expansions and contractions. Tandem repeat loci were derived from Tandem Repeat Finder annotations available from PacBio. The following commands were applied:

• Download figure
• Open in new tab

Sniffles2 (v2.0.6) was used to detect SVs with low coverage support using the following command:

• Download figure
• Open in new tab

The PAV v1.2.1 pipeline was used for assembly-based SV detection. Primary and haplotype-resolved contigs from hifiasm were aligned to the GRCh38_noalt or CHM13v2.0 reference using minimap2 with -x asm20 parameters. SVs were then called by comparing the assembled contigs to the reference. PAV identifies complex rearrangements, insertions, and deletions not easily captured by alignment-based methods.

The resulting VCFs from PBSV, Sniffles2, and PAV were retained separately and labeled with their source caller for traceability. These per-caller, per-chromosome VCFs were later normalized and harmonized as part of the SV integration and filtering pipeline (see Supplemental Materials: Building a cohort-level SV callset).

Cohort-level Callset Integration

Per-sample SV callsets from PBSV, Sniffles2, and PAV were normalized to a single callset, clustered and merged within each sample, filtered into stringent and lenient tiers, and then merged across all samples to produce cohort-level callsets. Tool selection and parameterization were guided by benchmarking against assembly-based truth sets. Finally, we re-genotyped the merged callset across the cohort to recover missed variants. Full implementation details, including software versions, parameters, and benchmarking comparisons, are provided in the Supplementary Methods.

Reference panel construction

We combined SVs from our strict callset on 1,027 AoU participants with 47 HPRC samples into a phased reference panel. Variants were first physically phased within samples using HiPhase³², filtered by allele frequency and proximity to SVs, statistically phased across the cohort using SHAPEIT4⁹⁵, and postprocessed to create a pangenome bubble graph Full implementation details are provided in the Supplementary Methods.

Imputation

Using the reference panel described above, we imputed SVs into (1) all 3,202 samples from the 1000 Genomes Project, and (2) 10,000 self-identified Black or African American AoU srWGS datasets. Genotyping was performed with KAGE to estimate genotype likelihoods from k-mer counts, and phasing and imputation was performed with GLIMPSE^34,35. Accuracy was evaluated through allele-frequency correlations, Hardy-Weinberg equilibrium, leave-out validation experiments, comparison against assembly-based truth, and Mendelian error rate in 1kGP trios. Full implementation and benchmarking details are described in the Supplementary Methods.

Repeat loci analysis

CYP2D6 variation

Hifiasm^98,99 haplotype-resolved assemblies from the 1,027 AoU participants were aligned to the T2T-CHM13v2.0 reference genome¹⁷. Haplotigs spanning the full CYP2D6-7 region were extracted, including up to 20 kbp of flanking sequence. Haplotypes with no or multiple haplotigs over the full locus were excluded.

Extracted haplotig sequences were divided into 100 bp k-mers and aligned to CYP2D6, CYP2D7, and CYP2D7-spacer T2T-CHM13v2.0 reference sequences. K-mers were colored according to best match by MAPQ score, and k-merized sequences were plotted (Code Availability). This pipeline was tested on 47 high-quality phased HPRC assemblies, for which the copy number of CYP2D6 and CYP2D7 showed 100% concordance to reported HPRC haplotypes, disregarding hybrid alleles¹⁴.

Short-read haplotyping from CRAM files for the 1,027 long read samples was completed using validated haplotype callers Aldy (v4.4)¹⁰⁰, Cyrius (v1.1.1)¹⁰¹, PyPGx (v0.20.0)¹⁰², and StellarPGx (v1.2.7)¹⁰³. A consensus call was calculated as a composite value of the haplotype call from all callers based on majority agreement between callers (minimum 2). An independent short-read analysis in which PyPGx was run equivalently except using control gene VDR produced the same PyPGx calls.

Star alleles were assigned by aligning PharmVar reference sequences to each non-CNV haplotype and selecting a best-matching star allele by CIGAR scoring, requiring the presence of all defining variants¹⁰⁴. Diplotypes and phenotypes were derived according to CPIC guidelines for samples for which both haplotypes were fully resolved¹⁰⁵.

To assess the accuracy of our star allele calling pipeline, star allele calls on 40 HPRC assemblies in the high-coverage 1KGP set were validated against the short-read consensus calls, achieving 95% concordance for unambiguously assigned non-SV haplotypes^14,39. For all disagreements our call was supported by at least one caller.

Missense variants were identified by predicting haplotype coding sequences using Augustus (v3.5.0)¹⁰⁶ with gene structure hints, aligning these to the GRCh38 coding sequence, and comparing the resulting amino acid sequences to those for reference star alleles. Lollipop plots were generated with ProteinPaint^107,108. AlphaMissense hg38 pathogenicity classifications⁶¹ were obtained for each base change from the CYP2D6 (P10635) AlphaFold heatmap.

Impact prediction of SVs

From the All of Us SV callset, we filtered out SVs with genotypes labeled as "./." or "0/0" across all samples and those ≤50 bp in length. We assessed the functional impact of the remaining SVs using CADD-SV⁶⁴ and analyzed the relationship between the number of samples sharing each SV and its predicted pathogenicity score. For SVs predicted to be likely pathogenic or absent from the matched short-read SV dataset, we examined their overlap with medically relevant genes, including those listed in OMIM⁶², high-priority genes associated with inherited diseases⁶³, oncology/tumor suppressor genes from OncoKB⁶⁸, and genes in the American College of Medical Genetics and Genomics (ACMG) list⁴¹.

To investigate the potential impact of these SVs on human health, we then accessed EHRs for each individual from the All of Us Researcher Workbench and mapped the clinical conditions to Systematized Nomenclature of Medicine - Clinical Terms (SNOMED CT). For high-scoring SVs intersecting coding regions of medically significant genes, we used AlphaFold3 to predict the resulting protein structures and compared them to the reference structures to assess potential functional disruption at the protein level.

Cross-population SV interpretation

We analyzed SV counts per sample and their allele frequencies across 2,504 unrelated 1KGP samples from five continental groups. To assess population stratification, we performed a PCA on genotyped SVs using PLINK¹⁰⁹, excluding variants within tandem repeat regions¹¹⁰. For comparison, we conducted a PCA using chromosome 1 SNVs from the 1KGP, applying the same filtering criteria. To investigate the functional impact of SVs, we mapped them to gene regions using GENCODE annotations, regulatory elements from Ensembl and analyzed their distribution across populations. F_st values were calculated for each SV between African and non-African populations with BCFtools¹¹¹.

SV-eQTL analysis

We performed eQTL calling following the method described by Kirsche et al⁷⁶ to investigate the impact of SVs on gene expression in 731 samples from the 1KGP with the matched RNA-seq data. Briefly, SV genotypes were categorized as 0 (homozygous reference), 1 (heterozygous alternative), and 2 (homozygous alternative), and only SVs with at least 5% of samples in each genotype class were selected to ensure statistical robustness. Using an ordinary least squares (OLS) linear regression method, we evaluated associations between SV genotypes and normalized gene expression levels, including sex as a covariate. Then we computed gene-level eQTL p-values by applying Bonferroni correction to the minimum eQTL p-value for each gene, adjusting by the number of eQTLs associated with that gene. These p-values were further corrected for multiple testing using the Benjamini-Hochberg method at a false discovery rate of ≤5%. All reported SVs had a mean genotype posterior (GP) ≥ 0.7 across the cohort. To validate genotyping accuracy, we selected 40 samples with both long-read and short-read data from the 1KGP cohort and computed per-genotype concordance for each sample (Supplementary Fig. 9f).

To compare the effects of SVs and SNVs on genes associated with both variant types, we used the previously generated SNV data from (https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased/), restricted to the same 731 samples. For each gene, we evaluated SNV-expression associations within a 1 Mbp window of the gene using the same OLS regression model applied to SVs. We then selected the top SV-eQTL and the 1,000 most strongest SNV-eQTLs per gene to calculate the z-scores and LD between variants. These values were used as input for the CAVIAR⁷⁹ fine-mapping analysis to estimate the posterior probability of each variant being causal.

Overlap between SVs and GWAS-significant traits

To investigate potential associations between SVs and phenotypes, we downloaded the GWAS Catalog v1.0 database and identified genome-wide significant variants (p < 5 × 10^-8) located within ±100,000 bp of our SVs. Genotype data for these GWAS variants was obtained from the 1KGP Phase 3 dataset³⁷, focusing on the same set of 2,504 unrelated samples used for SV genotyping. For each SV-GWAS variant pair, we calculated LD and retained those with r² ≥ 0.5. To characterize disease- and disorder-related traits, we used the ontologies from the Ontology Lookup Service (http://www.ebi.ac.uk/ols) and condition-domain concepts from the SNOMED vocabulary to classify traits, enabling further exploration of the potential medical implications of associated variants.

To further evaluate the medical relevance of SVs linked to GWAS variants, we identified SV carriers among 1,027 AoU participants and retained SVs present in at least 10 individuals. We then extracted EHR data from AoU Researcher Workbench and included 848 participants with available condition records. Clinical conditions with at least 10 cases and 10 controls were selected for analysis. We applied Fisher’s exact test to assess associations between SV presence and clinical diagnoses. Analyses of sex-specific conditions were limited to individuals of the corresponding sex. All reported SVs had a mean GP ≥ 0.7 and a minor allele frequency (MAF) ≥ 0.05 among the 2,504 1KGP samples. Genotyping accuracy was validated using the same 40 long-read/short-read 1KGP samples (Supplementary Fig. 10e).

Phenome-wide analysis of SVs in 10,000 self-identified Black or African American participants

To explore the clinical impact of identified SVs, we integrated genomic data with EHRs from an expanded AoU cohort. We firstly genotyped and imputed variants in 10,000 short-read samples from self-identified Black or African American participants, retaining SVs present in at least 10 individuals. Clinical condition data were obtained from AoU Researcher Workbench, and only conditions with at least 100 cases and 100 controls were included to ensure statistical power. To improve the accuracy of the phenome-wide association study, we excluded conditions not known to be influenced by genetic factors and those with overly broad or ambiguous names.

We applied Fisher’s exact test to identify conditions significantly associated with SV carriers. Analyses of sex-specific conditions were conducted only in individuals of the relevant sex. For SVs associated with disease, we extracted the GP values across all SV carriers and calculated the mean GP. SVs with a mean GP < 0.7 or deviating from Hardy-Weinberg equilibrium (HWE < 1×10^-5) were excluded to reduce false-positive associations. We next compared SV AFs between the 1,027 AoU long-read samples and the 10,000 short-read samples, excluding SVs with AF differences ≥ 1.5-fold.

Genotype confidence was further assessed by re-genotyping SVs with Locityper¹¹² and comparing concordance with imputed allele calls across carriers. This was performed in two steps: (1) building a variant database, for example with:

and (2) re-genotyping each variant against individual samples with:

SVs were removed if fewer than 70% of carriers showed allele differences < 20 between the two methods (Supplementary Fig. 11).

To evaluate the robustness of identified associations, we excluded samples with kinship > 0.177 or with fewer than 10 SV carriers and applied Firth’s penalized logistic regression, while adjusting for age, sex, sequencing center, and the first ten principal components. Associations were retained only if the effect direction was consistent across both methods and the OR estimates differed by no more than 1.5-fold (Supplementary Fig. 12).

We used a moderate significance threshold (p < 1×10⁻⁵) for the initial association scan to facilitate direct comparison of association strength between SNVs and SVs across traits, rather than to declare genome-wide significance. This screening threshold was chosen to capture a broader set of candidate variants for downstream evaluation. To determine whether the SV was the most significantly associated variant for the condition from this screening set, we identified SNVs within ±100 kbp and assessed their association with the corresponding phenotype using the same statistical method. Then, we applied CAVIAR for fine-mapping, incorporating each significant SV and the 1,000 nearby SNVs with the most significant association signals to estimate their posterior probabilities of causality. Finally, for SVs intersecting coding regions of medically significant genes with known disease relevance, we used AlphaFold3 to predict altered protein structures and compared them to reference models to evaluate potential functional consequences at the protein level.

We then assessed whether the disease-associated SVs were unique to the long-read dataset by comparing them with the short-read SV callset derived from 990 of the 1,027 samples used to construct the long-read reference panel (for comparison methods, see Supplementary Methods: Variant annotation and comparison to external datasets). We quantified the number of associations in which the SV was absent from the matched short-read callset. In addition, we identified associations specific to individuals with African genetic ancestry by selecting SVs that were only imputed among the 2,504 unrelated 1KGP samples with African genetic ancestry.

Furthermore, we performed a PCA using PLINK on 2,504 unrelated 1KGP individuals and the 10,000 AoU samples to assess population structure. For the reported associations in the results, AoU samples were further stratified into: (i) cases vs. controls, (ii) SV carriers among controls vs. all other individuals, and (iii) SV carriers among cases vs. all other individuals. The diffuse distribution of these subgroups across the PCA plots suggests the observed associations reflect real biological signals rather than batch effects or other artifacts in phenotypes coded at a single hospital or collection site.