Search
Close this search box.

Utility of long-read sequencing for All of Us – Nature Communications

Optimizing variant detection in cell lines

We first focused on four established HapMap cell lines to assess the long-read ability to cover and comprehensively identify variants across medically relevant genes. We used the widely studied sample NA24385 (Caucasian male), as well as HG00514 (female Han Chinese), HG00733 (female Puerto Rican), and NA19240 (Yoruba male), for which assembly-based analyses are available33. The latter samples (HG00514, HG00733, NA19240) are not as well-curated as the NA2438520,38,39. Nevertheless, they provide valuable information, since variant calling tools have generally not been previously optimized across these samples to minimize overfitting.

First, we assessed the genome-wide coverage and read length capabilities across these control samples to quantify the ability of the different sequencing technologies at this basic level. Because of the randomized nature of whole genome sequencing, coverage and read length establish the fundamental limitations of variant identification. For example, in an extreme scenario, if there are zero reads spanning a given region of the genome, clearly no variants will be detected there. Moreover, if variants are detected by only a single read they are generally less trustworthy and are often filtered out. One simple but effective strategy is to require at least two or more supporting reads to identify a variant since it is substantially less likely two reads will have the same error at the same position (e.g., assuming a 1% sequencing error, there is a 0.01% chance of two reads having the same error at the same position at random)12. In low coverage situations, however, requiring two reads spanning a position to identify a variant will limit the recall of variants, especially for heterozygous variants since the haplotype-specific coverage is half of the total coverage. Nevertheless, an idealized analysis assuming a Poisson coverage distribution shows 6× coverage is sufficient to recall over 98% of homozygous variants and 8x is sufficient to recall over 90% of heterozygous variants (Fig. 1a, Methods). One notable exception to this model is capturing insertion variants when the length of the variant approaches or exceeds the length of the read, since the amount of coverage that spans the variant will be proportionally reduced by the length of the insertions (Fig. 1b, Methods). This idealized analysis establishes the fundamental limits for variation detection. In real-world analysis, there are additional considerations, especially non-random sequencing errors, repeats, alignment errors, overdispersion in coverage, and other biases that further complicate variant calling. This requires an empirical approach to measure and resolve that we discuss in later sections of the paper.

Fig. 1: Theoretical assessment of variant calling as a function of coverage and insert size.
figure 1

a Idealized recall of homozygous and heterozygous variants as a function of overall coverage, assuming at least two spanning reads are required to recall a given variant in the genome. Here, we assume the sequencing coverage follows the Poisson distribution centered on a given overall coverage level, and the coverage will be evenly distributed across maternal and paternal haplotypes following the binomial distribution with p = 0.5. b Idealized recall of insertion variants of different lengths assuming 20 kbp reads, 5×, 10× or 20× coverage overall, and at least 2 reads spanning the insertion. As plotted, this represents the recall of homozygous variants, although heterozygous variants follow the same distribution as a function of the haplotype specific coverage.

The short-read Illumina coverage was between 29.76× (NA19240) and 32.50× (NA24385) (see Methods), and across the long-read technologies (Fig. 2a), we obtained an average coverage of 45.29× for ONT and 35.70× for PacBio HiFi. It is important to highlight that ONT required only 1–2 PromethION flow cells per sample to produce this coverage, while PacBio HiFi required several times more using Sequel IIe (HG00514, HG00733, & NA19240 each required four flow cells, NA24385 six flow cells). ONT also generated longer aligned N50 lengths compared to PacBio HiFi, with an aligned N50 length of 20 kbp compared to an aligned N50 length of 11 kbp (Supplementary Fig. 1) (Fig. 2a).

Fig. 2: Summary comparison across HapMap samples.
figure 2

a Coverage on X-axis vs N50 on the Y axis, for samples NA24385, HG00514, HG00733, and NA19240 (HapMap). b Indels distribution for HapMap samples, indels length on the X-axis, and count and density (in red) on Y-axis (using a merge of DeepVariant and Clair3 for HiFi). c SV size distribution for HapMap samples, SVs length on the X-axis, and count and density (in red) on the Y-axis (using a merge of Sniffles and pbsv HiFi). d Circular plot for HapMap samples from outside in genes density (Brown), SNVs (substitutions) average density between HapMap data (Orange), SNVs (indels) average density between HapMap data (Aquamarine), deletions average density between HapMap data (Blue), insertions average count between HapMap data (Green).

Next, we assessed the ability of the technologies for the identification of SNVs and indels across these samples using state-of-the-art methods. Starting with Illumina, we utilized the Dragen pipeline (v3.4.12)40 with the exact specification as AoU, resulting in a high recall of 99.32% and precision of 99.63%, leading to a total F-score of 99.47% across GIAB HG002 benchmark v4.2.139 (see Methods).

For the long-read technologies, we evaluated several best practice small variant callers (Longshot41, Pepper/DeepVariant42, and Clair343) as well as various combinations of them. Figure 2b shows the indel size ranges that were detected, and Fig. 2d shows the comparison of SNV and indel density across samples. Importantly, Longshot can only identify substitutions (Supplementary Figs. 2 and 3), which makes it complicated to compare this method to the others, which can also identify short indels. DeepVariant and Clair3 both showed a similar number of base substitutions when applied to data generated by PacBio and ONT sequencers. Specifically, DeepVariant identified 4,335,053 substitutions in the PacBio data and 4,709,454 substitutions in the ONT data. Similarly, Clair3 identified 4,597,773 substitutions in the PacBio data and 4,578,168 substitutions in the ONT data. On the other hand, Longshot resulted in an approximately 1.57-fold increase in the number of substitutions for data generated by ONT, although the excess is almost entirely false positives (see below). Overall, a combination of Clair3 and DeepVariant using PacBio HiFi data achieved the best F-score (99.87%) at this coverage level (as shown in Supplementary Fig. 4). For ONT, the F-score is 98.74% from merging the results from Clair3 and DeepVariant (Supplementary Fig. 5). It is worth noting that the improvement from merging Clair3 and DeepVariant for PacBio is marginal compared to using only DeepVariant or Clair3 individually, with a gain of only 0.01% to 0.04% in precision, respectively. Conversely, for ONT, the merge F-score enhancement is 0.70% (DeepVariant) and 2.50% (Clair3) improved over using the individual methods. Furthermore, we assessed the correctness of genotypes across NA24385 and found that Clair3 achieved higher accuracy than DeepVariant. Given this result, we adjusted our merging strategy to utilize the genotype values (GT flags) from Clair3 when available. With the pipeline optimized for one sample, we extended the analysis to HG00514, HG00733, and NA19240 and benchmarked them against their respective gold-standard variant call sets44. Clair3 merged with DeepVariant achieved the best results. Regardless, the overall performance is reduced compared to NA24385 (HG002). The difference in the accuracy of variant calling between HG002 and other samples is attributed to the benchmark set regions for HG002 being filtered to exclude repetitiveness or polymorphic complexity regions, as mentioned in Wagner et al.38 and Li et al. studies45. Additionally, the variants in this sample have been extensively validated, and their accuracy is well-established.

We assessed the technologies’ abilities for SV calling in medically relevant genes using Manta starting with NA24385 (see methods). Using the state-of-the-art algorithm Manta developed by Illumina46, Illumina’s F-score is 0.45, due to its inability to identify large insertions. In contrast, we found that a combination of pbsv and Sniffles for ONT and PacBio achieved the best results, with an F-score of 0.93 using HiFi and 0.91 for ONT (see Supplementary Fig. 6). Furthermore, we conducted a titrated coverage comparison for sample HG002, focusing on the two long-read technologies, PacBio and ONT. We evaluated the SV detection F1-score for this comparison (Supplementary Fig. 18). Overall, the trends follow the idealized results presented in Fig. 1, with a sharp increase in accuracy between 5× to 15× coverage and align with other ongoing research47. We also observed that below 8× coverage, ONT demonstrated a higher F1 score for SVs compared to PacBio HiFi. However, beyond 8× coverage, there was a significant improvement in the F1-score for PacBio HiFi, exceeding 90%, while the F1-score for ONT remained around approximately 87%. When applied to the other samples, we found that using long-reads improved the accuracy of SV detection, although the F-score was slightly lower at 0.77 (as shown in Supplementary Fig. 7). When comparing the detection of SVs between the samples (as depicted in Fig. 2c), we observed that the frequency of detected SVs was similar.

Overall, the results of this study reproduce several previous findings that long-read sequencing enhances SV calling and achieves a high accuracy of genome-wide SNV and indel calling. We further established an improved variant calling pipeline for SNVs and SVs for long reads that work similarly well for PacBio HiFi as well as ONT and thus makes the comparison across the sequencing technologies faster and easier.

Performance assessment for All of Us samples across medically relevant genes

To provide a more realistic assessment of the value of long-read sequencing in a clinical research setting, we next utilized long-reads to benchmark the analysis of control samples commercially sourced by the AoU. These samples were sequenced multiple times across the individual AoU genome centers to establish and assess the variability of each center for different tissue sources and preparations. Specifically, we sequenced two anonymized AoU samples T662828295 and T668639440 using ONT, HiFi, and Illumina technologies for different tissue sources (white blood cells and whole blood cells, henceforth WBC and whole blood cells, respectively) and extraction methods (Chemagen and Autogen).

First, when assessing coverage, we observed that in this limited sample set, WBC achieved a better coverage with either Autogen or Chemagen extraction for ONT and HiFi, while the opposite was found for Illumina (whole blood cells achieve a better coverage). The minimum and maximum coverage for AoU samples T662828295 and T668639440 are shown in Supplementary Data S1: Illumina (31.95, 41.68, 34.19, 39.30), HiFi (6.46, 10.50, 5.11, 15.76) and ONT (28.92, 29.26, 28.09, 29.46). This coverage resulted from using one flow cell for ONT and maximum two flow cells for HiFi (Supplementary Fig. 8 and Supplementary Data S1). Furthermore, we observed that the extraction method impacted the N50 alignment length (Supplementary Fig. 9). The average N50 of aligned reads for HiFi (17,612) is greater than that of ONT (15,726 bp). However, ONT produced the longest single read alignment by a large margin (758,354 bp).

For small variant calls (substitutions, insertions, and deletions <50 bp), we assessed the concordance between technologies and found that long- and short-read, agree on approximately 67.55% (median 70.47%). We can explain this reduced concordance level by the lower HiFi coverage and the difficulty in accurately calling insertions and deletions between technologies. Nevertheless, when we compare substitutions only, the concordance reaches 79.00% on average (Supplementary Data S2). Unsurprisingly, among variants not found by all three technologies, HiFi (5.44%) shows greater SNV concordance with Illumina when compared to ONT (8.46%); except for sample T668639440 whole blood Chemagen, where ONT agreed two folds higher than HiFi with Illumina 15.49% and 6.46%, respectively. This is likely due to the much lower coverage from HiFi (5.11x) compared to ONT (28.09×) data set, as mentioned earlier. Additionally, we noticed an enrichment for Illumina-only identified variants (mean 12.10% and median 8.98%) compared to a lower unique identification on each of the long reads (ONT or HiFi ~2.29%). Of note, the higher Illumina average of 12.10% uniquely identified SNVs was chiefly due to the T668639440 WBC Autogen sample. This outlier substantially skewed the mean, causing 33.54% of variants to be identified uniquely by Illumina. Correspondingly, when we focused only on exon and intron regions (point mutations and indels), the concordance between long- and short-read increased to 81.33% and 77.46%, respectively. Additionally, HiFi showed higher concordance in exons and introns with Illumina (Supplementary Data S3).

We utilized SnpEff48 across the 24 merged data sets to annotate the merged variants obtained from the three technologies. A total of 10,939,305 variants were processed and annotated. Among them, we identified 7622 variants (0.02%) across the 24 data sets with high-impact annotations, including stop-gained or frameshift variants, which were of particular interest for further analysis.

When we compared the high-impact variants identified by ONT, HiFi, and Illumina, we discovered that the percentage of high-impact variants detected exclusively by Illumina was identical to those detected by both HiFi and ONT, which was just 0.04%. The same percentage was also observed for high-impact variants detected by the union of HiFi and Illumina. However, considering the high-impact variants detected by either ONT or HiFi alone, this percentage increased to 0.05%. Additionally, when evaluating all long-read technologies collectively, the percentage of high-impact variants identified rose to 0.09%, as depicted in Supplementary Fig. 19.

We next analyzed SVs using the three technologies (HiFi, ONT, and Illumina). We identified an average of 24,235 SVs per sample (for each tissue and extraction method), which aligns with previous studies44. The percentage of SVs agreed upon by all three technologies is approximately 22.00%; ONT and HiFi agreed on 53.86% of all SVs (31.90% identified only by long read plus 22.00% identified by all reads). However, 22.40% remains that is either detected exclusively by ONT or HiFi, meaning that 32% of the SVs identified by long read were not detected when using short read sequencing. Moreover, approximately 15% of SVs were found exclusively by Illumina.

We studied in depth each of these distinctive variants per technology. We found that 950 SVs were identified uniquely by Illumina in all datasets, and the majority are translocations (68.63%), followed by duplications (11.26%). For ONT, we identified 57 unique SVs, with the majority being deletions (56.14%), followed by insertions (42.11%), and duplications (1.75%). Furthermore, we found that PacBio had the least number of unique SV (54), with the majority being insertions (40.74%), followed by deletions (37.04%), duplications (16.67%), inversions (3.70%), and translocations (1.85%). Based on these results, it is likely that Illumina reports a higher number of false SVs, particularly translocations in healthy individuals, as seen in prior studies12,32. Meanwhile, deletions dominated the uniquely-identified SVs from ONT, while insertions dominated for HiFi. Additionally, we utilized vcfanno49 to annotate the SVs that we identified. We calculated the percentage of SVs that overlap genes (entire gene body) and presented these findings in Supplementary Fig. 19. This might be inflated as larger inversions, for example, might not have a direct impact on genes. Our analysis revealed that the percentage of SVs overlapping genes was lower in Illumina sequencing (~43%) compared to HiFi and ONT sequencing (~54%), as well as all SVs identified by long-read sequencing (~52%). Furthermore, when we filtered out inversions and recalculated the percentages of SVs overlapping genes, we noted a decrease of less than one percentage, as illustrated in Supplementary Fig. 19.

To investigate the clinical utility of long-reads for a program like All of Us, we focused on a set of 4641 genes that are reported to be medically relevant38. Most notably, these genes represent non-repetitive or generally non-complex genes of the human genome (see Fig. 3a). The coverage across these genes was similar to the genome-wide coverage across all technologies (see Fig. 3b). In Fig. 3b, we compared normalized gene coverage (average gene coverage divided by sample average coverage) between HiFi and ONT and found larger coverage variability for HiFi likely due to the overall lower coverage.

Fig. 3: Comparison across medically relevant genes.
figure 3

a Percentage of gene bases intersecting with hard-to-map regions in yellow is the 5207 genes and in blue 386 medical gene sets. b Normalized gene coverage (normalized by sample coverage) for HiFi (yellow) and ONT (blue) for T668639440 and T662828295. c Log percentage of SNVs per gene for the 4641 medical genes for samples T668639440 and T662828295.

Similarly, we compared the number of genes with average coverage less than one (henceforth, “uncovered genes”) across the sequencing technologies. For Illumina, we identified only three uncovered genes for T668639440 Chemagen (C4A, C4B, and OR2T5), while Autogen only had two uncovered genes (C4A and C4B) for the same sample and only one gene for sample T662828295 Autogen (C4A). PacBio HiFi covered all genes, except when using whole blood Chemagen. In sample T662828295 using whole blood Chemagen the gene PDE6G for sample T662828295 and 14 genes in T668639440 are uncovered. In contrast, we observed that all genes were covered using ONT. Overall, ONT and HiFi (regardless of the one sample that shows low coverage T668639440 whole blood Chemagen) showed the least amount of uncovered medically relevant genes.

To further assess the coverage from a medical perspective, we downloaded SNPs and indels from ClinVar that were reported pathogenic and checked whether the individual sequencing technologies adequately covered these variant locations. We used variants that are not conflicting in reporting their pathogenicity and submitted by multiple clinics (see Methods). This resulted in a set of 10,368 variant sites across the 4641 medically relevant genes.

We calculated the number of uncovered (coverage <1) variants per sample and condition: For the PacBio HiFi sample, on average, we have six variants without coverage for sample T662828295 (only two out of four datasets have uncovered variants) and 24 for T668639440 (three datasets have uncovered variants). ONT covered all variant sites, while Illumina, on average, did not cover six variants across the individual runs for sample T662828295 and seven for T668639440 (Supplementary Data S4).

Next, we analyzed the list of 73 American College of Medical Genetics and Genomics (ACMG) genes, in which mutations are commonly recommended to be reported to patients. These genes span both gene groups defined in this study, including mostly easily accessible genes (68) but also some challenging genes (5). Overall, the genes are well covered by each of the technologies and the normalized average gene coverage is one or more (Supplementary Fig. 10 and Supplementary Data S5). For the five challenging genes, we observed that the normalized average coverage is slightly higher for ONT than for HiFi, indicating a better mapping overall for the ONT data (Supplementary Fig. 11 and Supplementary Data S5). Additionally, for the pilot data T668639440 and T662828295, the normalized mean coverage for accessible genes are similar, HiFi (1.06) and ONT (1.04). Nevertheless, coverage differs in challenging genes, where HiFi is 0.88 and ONT is 1.03. For the benchmark of variants within genes available in the GIAB, the results from HiFi and Illumina are similar compared to ONT. However, the average F-score across these genes is higher for Illumina (93.64%) compared to HiFi (85.24%) and ONT (73.98%) that is due to low F-score for gene TNNI3 (Supplementary Data S5 and Supplementary Fig. 12). Based on the previous finding, we can conclude that the ACMG list is well covered by all technologies; Likewise, we can call variants with high accuracy using either Illumina or HiFi.

Resolution of highly challenging medically relevant genes

We assessed the utility of long-reads specifically within highly complex and repetitive medically relevant genes. In principle, this is where the long-read technology offers its greatest advantages, but it remains to be shown across primary tissues from patient donors. For this, we focused on a gene set from a recent GIAB publication, which proposed 386 genes that were found to be highly challenging for mapping and variant calling20.

We first assessed the normalized gene coverage (Fig. 4a) between HiFi and ONT. Here, HiFi has a lower median gene coverage than ONT similar to the genome-wide results. Further, normalized coverage distribution is centralized around median 1.05× for ONT compared to 0.97 for HiFi (Supplementary Fig. 13). Additionally, we counted the number of uncovered genes (i.e., genes that do not contain a single read).

Fig. 4: Comparison across 386 challenging medically relevant genes.
figure 4

a Normalized gene coverage distribution for the 386 genes in different datasets for ONT and HiFi. b The percentage of SNVs (substitutions and indels) in the data sets for the 386 genes. c Percentage of zero-base coverage per gene, where 50% or more of the gene body is not covered using HiFi and ONT data. d SVs breakpoints that intersect with medical genes across samples for different tissue and extraction methods T40 stands for sample T668639440 and T95 for T662828295, W for Whole Blood cells, B for WBC, and finally A for Autogen and C for Chemagen.

For Illumina, nine genes (CCL3L1, CRYAA, DGCR6, DUX4L1, H19, NAIP, PRODH, SMN1, and U2AF1) are uncovered in sample T668639440 and eight genes (CCL3L1, CFC1B, DUX4L1, H19, HLA-DRB1, SMN1, TAS2R45, and U2AF1) in sample T662828295. Thus, five genes are uncovered across both samples (Supplementary Data S6). Interestingly, in our analysis, we found that ONT and HiFi did not cover genes H19 and U2AF1. Nevertheless, previous studies found that these genes are incorrectly duplicated in GRCh38, which makes it hard, if not possible, to call variants in these genes19,50. However, genes CCL3L1 and DUX4L1 are covered, making them only challenging for short reads. Interestingly, the SMN1 gene differentiates between the ability of long-read to untangle this complex repetitive gene. While HiFi could not support coverage for the gene in all samples, ONT managed to cover it in ~50% bp of sample T662828295. Nonetheless, ONT failed to do the same in sample T668639440.

We next compared the percentage of genes where 50% or more of the gene body lacks coverage (Fig. 4c). In the majority of samples, we saw that the 386 genes group has a higher percentage of gene bodies that are not fully covered. However, for ONT, the percentage of uncovered gene bodies is always lower than HiFi for 386 and 5027 groups alike. Additionally, we saw only that the difference between the two gene sets is in guanine (G) and cytosine (C) GC content percentage (Supplementary Fig. 14). Thus, we conclude that the difference in the percentage of uncovered genes is due to the sample coverage.

We next assessed variant calling ability for long reads, starting with the GIAB sample that has a gold standard benchmark to compare. Specifically, we employed sample NA24385 to rank and characterize 273 genes (70.75% of our 386 gene panel) for which a GIAB benchmark is available. We excluded seven genes that do not report a variant in this benchmark. For each technology, we investigated the top and lowest ten genes ranked by F-score (Supplementary Data S7) and compared these genes across the other technologies. Importantly, the top ten genes that achieved the highest F-scores with Illumina, had the same or better F-score with HiFi. Meanwhile, for these same top ranked genes, ONT had a lower recall for three genes by failing to identify an insertion in each gene (PIGV and MYOT) and two substitutions in MYOT.

For the ten lowest performing genes using HiFi (CBS, CRYAA, GTF2IRD2, H19, KCNE1, KMT2C, MDK, MUC1, SMN1, and TERT), in six genes (CBS, CRYAA, GTF2IRD2, KCNE1, MUC1, and SMN1) HiFi still showed a better performance than Illumina. Likewise, ONT achieved a better F-score in these genes than Illumina. Moreover, in SMN1, KCNE1, and CBS genes, the ONT F-scores are better than HiFi and Illumina. However, in KMT2C and TERT genes, Illumina outperformed both HiFi and ONT F-scores (Illumina 67.47% and 59.74%), (HiFi 44.76% and 54.32%), and (ONT 32.47% and 47.78%), respectively.

For AoU samples T668639440 and T662828295, because we do not have an established benchmark for these samples, first, we compared the percentage of variants per gene to identify any abnormalities (Fig. 4b). When we consider the distribution of the variants in the 386 genes across datasets (tissue source and extractions), we found that the variant distribution is similar and the tissue source or extraction method did not substantially affect the variant distribution. We then analyzed the concordance of substitutions and indels across the 386 genes, and found the lowest concurrence in introns between the three technologies (approximately 60.47% for exons and 57.67% for introns).

For these 386 genes, we observed that the variant calling capability of Illumina seems to be highly impacted, even in exonic regions, compared to the easier set of 5027 medically relevant genes. This is most clear in the overall concordance between short and long reads: For the 5027 genes, we had a high concordance (81.33%) for all three platforms, while for the 386 the concordance drops to 60.74%. The latter is impacted by a reduced concordance of Illumina and showing also higher Illumina-only variants that are likely falsely identified (Supplementary Data S8). Maybe not surprisingly, this trend is amplified for the intronic regions (Supplementary Data S8)

We next investigated the coverage of the ClinVar pathogenic variants that intersect with these genes (see Methods). For HiFi sequenced samples, we achieved on average ~8× coverage for T662828295 and ~9× for the T668639440 sample per variant, both are lower than what the 5027 genes group achieved (9× and 10× respectively) (Supplementary Fig. 15). Likewise, not surprisingly, Illumina achieved the highest average coverage per variant, followed by ONT and HiFi respectively (Supplementary Data S4). Interestingly, when we compare the number of uncovered variants for each sample per technology, we can see a systematic distribution of variants without coverage in Illumina (SD: 1.39) compared to the more variable ONT results (SD: 5.93) (Supplementary Fig. 16 and Supplementary Data S4); However, we could not find a correlation between the gene composition (GC, GT, or AT rich) and the coverage (Supplementary Fig. 17).

In summary, HiFi outperformed the other technologies in both precision and recall (Supplementary Data S7). Furthermore, there are a few genes that are particularly challenging where all technologies missed variants like H19, which has one substitution (C/G) (Chr11 at position 1,996,209) that none of the three technologies managed to detect; KRTAP1-1, wherein both long-read technologies (HiFi and ONT) did not call any variants, while Illumina called two false-positive substitutions; and MDK, where ONT did correctly recall variants, but the other technologies called false-positive variants.

Interestingly, all the genes that ONT detected variants with low F-score were caused by uncalled indels (max three indels) like FLAD1 and PIGV (Supplementary Data S7 shows in more details genes names and F-score and the number of identified variants per technology). Finally, we compared the identified SVs between tissues and extraction in these samples (Fig. 4d) in the medically relevant genes. As we can see, bars 8 and 9 show the unique variants in AoU samples T662828295 and T668639440 supported by different tissue sources and extraction methods; likewise, bar 10 shows the effect of low coverage on sample T668639440 whole blood Chemagen as it lost around 4.4% of SVs that should share among all tissues and extraction methods.