Benchmarking nanopore sequencing and rapid genomics feasibility: validation at a quaternary hospital in New Zealand

The rise of clinical genomic testing, utilizing exome and whole genome sequencing, has enabled the detection of genomic changes (i.e. single nucleotide variants [SNVs], small insertions and deletions [indels], copy number variants [CNVs], structural variants [SVs]), and elucidated the underlying genetic basis of rare Mendelian disorders and cancers^{1,2,3,4,5,6,7}. Over the past decade, the increasing adoption of genomic testing has generated substantial evidence supporting precision medicine^2,4,6. Such approaches have enabled molecular diagnoses for genetic disorders, guiding tailored medical interventions^6,7,8,9.

Genomic testing has been established using cost-effective, short-read (150 bp in fragments) sequencing platforms from Illumina (i.e. NovaSeq, HiSeq, MiSeq, NextSeq), Thermo Fisher (Ion Torrent sequencer), and BGI (i.e. BGISEQ and MGISEQ). However, short-read sequencing has well recognised limitations^10,11. Firstly, it is difficult to uniquely align short reads to complex repetitive genomic regions involved in short tandem repeat (STR) expansion disorders (e.g. Fragile X syndrome and Huntington’s disease)^10,12. Secondly, the requirement for PCR amplification in short-read sequencing may contribute artefacts and hinder the identification of native base modifications¹⁰. Thirdly, short read lengths hinder the identification and precise phasing of alleles in large SVs¹⁰.

The recently developed long-read sequencing (LRS) platforms (i.e. Oxford Nanopore Technologies [ONT] and Pacific Biosciences [PacBio]) employ direct inspection of single molecules during DNA synthesis, yielding long phaseable reads (>10 kb) in real-time^12,13,14,15. Consequently, long reads generate highly reliable complete genome assemblies¹⁵, which can serve as benchmarks for short-read data. The utilization of ONT long reads as a standalone sequencing platform in clinical diagnosis has been demonstrated^16,17,18. In the research setting, LRS has been used to: a) identify and fine-map structural variations at single-nucleotide resolution; and b) resolve the haplotypes of heterozygous SVs^13,14,19. Novel pathogenic variants have been uncovered by LRS technology in human diseases with a previously unknown underlying genetic cause²⁰. Additionally, long reads facilitate the characterization of pathogenic repeat expansions in genomic regions that are challenging to sequence using short-read sequencing technology^21,22.

The clinical application of LRS^{13,14,15,16,17,18,19,20,21,22} requires confidence in the accuracy of variant calling for SVs, CNVs, STRs, and SNVs/indels. However, high per-base error rates in low-complexity and homopolymer sequences^12,23, and other issues have led to concerns about the application of ONT in clinical settings. Thus, there is a need for comprehensive benchmarking to: 1) confirm precision relative to the routinely used short-read technologies; and 2) illustrate the benefits and limitations of LRS technology for application to CGS.

The Global Alliance for Genomics and Health (GA4GH) developed benchmarking protocols to evaluate the performance of sequencing platforms and variant-calling methods before their integration into clinical practice^24,25,26. Benchmarking is essential to ensure adherence to standards and relies upon datasets where the relationship between input and output is known. This facilitates testing of consistency between the expected and observed outcomes (true positives)²⁴. The Genome in a Bottle (GIAB) consortium offers reference samples (e.g. HG001 from the HapMap project and trios of Ashkenazi Jewish and Han Chinese ancestry from the Personal Genome Project) with ground-truth calls for SNVs, small indels, and SVs^27,28. Notably, GIAB recently provided a curated benchmark of challenging medically relevant genes through haplotype-resolved whole-genome assembly²⁹. The GA4GH resources enable performance assessment, optimization, and analytical validation of CGS assays and workflows for detecting genomic variations^24,25. Indeed, GIAB datasets and benchmarks are considered the gold standard for evaluating sequencing technologies and variant calling pipelines²⁷.

The ONT platform generates sequencing data in real-time, allowing samples to be distributed across flow cells to reduce the sequencing time, where each additional flow cell reduces the sequencing time needed on a sample by 1/n (where n is the number of flow cells). Notably, this has been demonstrated in a recent study that sequenced a single human genome across 48 flow cells, generating high-depth genome-wide data (200 Gigabases) and candidate variant identification in less than eight hours³⁰. The ONT platform is also capable of targeted sequencing through adaptive sampling, which removes the need to design custom probes to capture genes or regions of interest through a dynamic and modifiable process during the sequencing run³¹. DNA and RNA base modifications, including 5-methylcytosine (5mC), 5-hydroxymethylcytosine (5hmC), and N⁶-methyladenine (6 mA), can also be detected computationally on raw ONT data without the need to perform special library preparations such as bisulphite conversion³², which is known to cause DNA damage and can lead to overestimation of the 5mC level³³.

We have established an expandable rapid genomic testing pipeline based around the ONT PromethION2 (P2) solo system connected to AI-driven genomics analysis and interpretation software (i.e. Fabric GEM™ software) for tertiary analysis. In the establishment phase, we benchmarked our pipeline using GAG4H tools and GIAB reference cell lines HG002 – HG007 for SNVs and small indels analysis. In addition, we used CNVPANEL01 (Coriell Institute) to measure our ability to detect large-scale chromosomal abnormalities. Finally, we present the results of the pipeline validation phase, performed in parallel with a clinically accredited short-read rapid genomic testing service.

Benchmarking samples and truth sets

We acquired CNVPANEL01 as 3 µg genomic DNA (at 100 µg/ml) per sample and GIAB reference samples (i.e. HG002 – HG007) with available truth sets, from the Coriell repository (Coriell Institute for Medical Research, 403 Haddon Avenue Camden, NJ 08103, USA).

Library preparation and nanopore sequencing

DNA samples (1500 ng) were sheared to 10–15 kb using Covaris g-TUBES (Covaris) in a bench-top centrifuge for 1 min at 2000 RCF (room temperature). Nanopore sequencing libraries were prepared according to the genomic DNA Ligation Sequencing Kit V14 (SQK-LSK114) protocol (ONT, Oxford Science Park, OX4 4DQ, UK). Prepared libraries were loaded on PromethION flow cells (R10.4) and sequenced (i.e. depth of between 24-42X) with the PromethION 2 (P2) solo device using Kit 14 chemistry and MinKNOW v23.07.8 (Oxford Nanopore Technologies [ONT], Oxford Science Park, OX4 4DQ, UK).

Base calling of nanopore reads and variant calling

Base calling of raw ONT signal data was completed using Dorado v0.3.3 (https://github.com/nanoporetech/dorado) with the high accuracy (hac) model ([email protected]). In addition, base calling of the HG002 sample was also completed with the super accuracy (sup) model ([email protected]). The resulting FASTQ files, with a Phred quality score (Q score) > 9, in the fastq_pass folder, were processed with EPI2ME Labs’ wf-alignment pipeline (https://github.com/epi2me-labs/wf-alignment; v0.5.2). Briefly, FASTQ files were aligned to the GRCh38 reference genome using minimap2 (v2.26)³⁴. EPI2ME Labs’ wf-human-variation pipeline (https://github.com/epi2me-labs/wf-human-variation; v1.7.0) was subsequently employed for genomic variant processing, including SNV and small indel calling with Clair3 (v1.0.4)³⁵, SV calling with Sniffles2 (v2.2)³⁶, and CNV calling with QDNAseq (v1.38)³⁷ using default parameters, with a VNTR annotation file provided for accurate SV identification. Repeat expansions were genotyped using Straglr (https://github.com/philres/straglr)³⁸ as implemented in EPI2ME Labs’ wf-human-variation pipeline v1.7.0.

Benchmarking of variant calling

Variant comparison tools (https://github.com/ga4gh/benchmarking-tools)²⁴ are integral to genomic benchmarking as they identify shared variations between ground-truth calls and comparison results (i.e., true positives [TP]), along with variants unique to each set (i.e., false negatives [FN]), and additional variants (i.e., false positives [FP]). We compared called SNVs and small indels with GIAB ground-truth variants (benchmark version v4.2.1)²⁴ using hap.py v0.3.15 (https://github.com/Illumina/hap.py), and each variant was labelled as TP, FP, or FN. Hap.py also provides precision (positive predictive value [PPV]), recall (sensitivity) and F1 scores (harmonic mean of precision and recall) calculated as follows:

For SVs, we employed Truvari v4.1.0 (https://github.com/ACEnglish/truvari)³⁹ to benchmark variants with GIAB ground-truth SVs. Each variant was categorized as TP, FP, or FN based on this comparison.

Rarefaction benchmarking analysis

Rarefaction was performed to evaluate the sensitivity and reliability of long-read variant calling across different sequencing depths. Subsampling of the Binary Alignment Map (BAM) files was performed using Samtools⁴⁰, by randomly selecting subsets of reads from the original alignment files. The subsampled BAMs were then subjected to variant calling analysis as described in the variant calling section. Benchmarking for SNVs and small indels was conducted as detailed in the benchmarking section. Rarefaction curves were generated using python v3.10.8 and the seaborn v0.12.2 library to illustrate the relationship between sequencing depth and the called variants, enabling the evaluation of variant calling performance and reliability across varied sequencing depths.

Benchmarking analysis for challenging clinically relevant genes

We called SNVs and small indels across genomic regions overlapping challenging clinically relevant genes²⁹ using the original BAM files and pipeline outlined in the variant calling section. Benchmarking for SNVs and small indels was conducted as detailed the benchmarking section.

Methylation analysis

Raw ONT signal data in POD5 files (https://github.com/nanoporetech/pod5-file-format) was base called (Dorado v0.5.0) using the high accuracy (hac) DNA base modification model ([email protected]_5mCG_5hmCG@v2) to detect modified bases (i.e. 5-methylcytosine [5mC], 5-hydroxymethylcytosine [5hmC]). The modified BAM files (modBAMs) were aligned to the GRCh38 reference genome and modkit v0.2.3 (https://github.com/nanoporetech/modkit) employed to generate genome-wide summary counts of modified and unmodified bases into bedMethyl files. Haplotype-specific 5mC differentially methylated regions (DMRs) in the HG002 genome were identified using ont-methylDMR-kit (https://github.com/NyagaM/ont-methylDMR-kit), which utilizes the Bioconductor DSS (Dispersion Shrinkage for Sequencing) package⁴¹. We used the following DSS⁴¹ parameters for calling DMRs: delta (threshold for defining DMRs) at 10%, p-value at < 0.01; minimum DMR length of 100 bps, and at least 10 CpG sites per DMR (https://github.com/NyagaM/ont-methylDMR-kit). The pipeline supports haplotype-specific analysis, DMR detection between two samples, group methyl analysis, as well as genomic annotation of significant DMRs (https://github.com/NyagaM/ont-methylDMR-kit).

Benchmarking results visualization

Plots were generated using the seaborn v0.12.2 and matplotlib v3.7.1, and python v3.10.8.

Newborn Genomics Programme (NBG) study design

This is a research study to determine the medical and economic impacts of rapid whole genome sequencing (rWGS) within the New Zealand health care landscape. Ethics approval was obtained from the Northern B Health and Disability Ethics Committee for the study entitled: Newborn Genomics – Te Ira oo Te Arai (Ethics reference: 2023 FULL 15542). Locality approval was obtained from the Research Review Committee Te Toka Tumai Auckland for the project entitled: Newborn Genomics – Te Ira oo Te Arai (Reference A + 9855 [FULL 15542]). This study is registered in ClinicalTrials.gov (Newborn Genomics Programme; NCT06081075; 2023-10-12). The clinical protocol was adopted and modified as per Lunke et al., 2023⁴².

NBG study participants and recruitment criteria

Children with suspected genetic conditions and their families were recruited into the study from the neonatal and paediatric intensive care units (i.e. NICU and PICU, respectively) and the National Metabolic Service at Te Toka Tumai | Auckland City Hospital (New Zealand) between November 2023 and August 2024. Within NICU, participation was limited to proband-parent trios of critically sick neonates with evidence of a suspected genetic condition, without a clear non-genetic aetiology, or who developed an abnormal response to standard therapy for an underlying condition within the preceding seven days. For infants within PICU or under the care of Metabolic Services, participation was limited to proband/parent trios of children with an acute or chronic illness with evidence of a suspected genetic condition without a clear non-genetic aetiology.

All participants continued to receive the standard of care, irrespective of whether they were included in the study.

Potential participants were referred to the geneticist on-call (by telephone) for a formal genetic review, mainly by a neonatologist or a paediatric intensivist or the lead paediatric subspecialist for the patient when a genetic condition was suspected, or when the aetiology of a condition was unclear and a genetic cause needed to be ruled out to guide further clinical management.

Inclusion and exclusion criteria were modified from Dimmock et al.⁴³ and McKeown et al.⁴⁴.

The inclusion criteria was:

• acutely ill inpatient

• admitted to NICU or PICU between April 2023 – March 2026

• under the care of the National Metabolic Service between April 2023 and March 2026

• within 1 week of hospitalization or within 1 week of developing abnormal response to standard therapy for an underlying condition

• suspected genetic condition, without a clear non-genetic aetiology

The exclusion criteria was:

• patients whose clinical course is entirely explained by

  1. ▪

     isolated prematurity

  2. ▪

     isolated unconjugated hyperbilirubinemia

  3. ▪

     infection or sepsis with expected response to therapy

  4. ▪

     a previously confirmed genetic diagnosis that explains the clinical condition

  5. ▪

     isolated transient neonatal tachypnoea

  6. ▪

     meconium aspiration syndrome

  7. ▪

     trauma

• inability to source blood or buccal samples for DNA extraction from at least the mother and child

Of note, participants were only considered for the study if they were referred to clinical genetics as a part of their standard of care workup.

Following the referral of potential participants to the study, a multidisciplinary meeting (MDM) was convened via video conference to evaluate the eligibility of the referral based on the study’s inclusion and exclusion criteria. At a minimum, the MDM was comprised of a clinical geneticist, genetic counsellor, principal investigator, project manager, representative from the genomic analytical team (bioinformatician, variant curator) and the referring clinician. After agreeing to participate, the patients were registered on RedCap, and a study reference was generated. Subsequently, the clinical geneticists and genetic counsellors completed the clinical information, including phenotypic characterization using HPO terms, and facilitated informed consent.

Parents or guardians of the proposed probands were informed of the details of the study using the HDEC-approved Newborn Genomics Programme Participant Information Sheet (Supplementary Note 1), and had the opportunity to ask questions to an on-call geneticist and genetic counsellor from the Genetic Health Service New Zealand. Written informed consent was obtained from parents or guardians before any study-specific processes were undertaken (Newborn Genomics Programme consent form; Supplementary Note 2).

Clinical geneticists and subspecialists performed clinical phenotyping, which was recorded on RedCap using the Human Phenotype Ontology (HPO) terms (https://hpo.jax.org/app/) to optimize phenotypic data exchange during the curation stages of the analysis. At the same time, a phenotype-focus gene list was generated using PanelApp (Australia [https://panelapp.agha.umccr.org/] and the UK [https://panelapp.genomicsengland.co.uk/]) and shared with the genomic analytical team for inclusion in the Bayesian AI-based clinical decision support tool (Fabric GEM™ software).

NBG sample collection, DNA extraction, library preparation, sequencing, and variant calling

After obtaining consent, duplicate blood samples were collected: 4 mL EDTA blood samples from the mother and father, and 500 µL EDTA blood samples from the child. One set of samples was sent to the Liggins Institute newborn genomics laboratory for sequencing and variant analysis, while the second set was sent to the clinical laboratory, Victorian Clinical Genetics Services (VCGS) in Melbourne, Australia, for a concurrent, independent, short-read-based analysis as described in Lunke et al. 2023⁴².

High molecular weight DNA was extracted from 300 µl of the whole blood using the Puregene DNA extraction Kit (Qiagen) following the manufacturer’s protocol, and the extracted DNA eluted in nuclease-free water (Thermo Fisher Scientific). The quantification and purity assessment of the DNA samples were performed using the Qubit system (Thermo Fisher Scientific) and a spectrophotometer (Implen NanoPhotometer). The library preparation and sequencing procedures were carried out as detailed in the library preparation section. Finally, the base calling of sequenced reads and variant calling analysis was conducted following the methods described in the base calling and variant calling sections. Of note, we have developed a Standard Operating Procedure (SOP) to lock pipelines to specific versions and outline procedures for updates and upgrades since the software for base calling and variant analysis is frequently updated (Supplementary Note 3).

NBG candidate variant prioritization and genomic results reporting

For variant filtering and prioritisation, we used Fabric GEM^TM as the primary interpretation platform and QCI Clinical Insights Interpret-translational (QIAGEN Inc., https://digitalinsights.qiagen.com/) as the secondary ‘research-confirmatory’ platform. We retained variants with ≥ 10 reads, a Variant Allele Frequency (VAF) ≥ 20% in the proband, a frequency ≤ 1% in gnomAD v3.1 and those located in the exonic regions or within +/−20 bases of exon/intron boundaries. Additionally, intronic variants beyond ±20 bases from exon start/end predicted to affect splicing by MaxEntScan were also retained.

In the initial analysis, we focused on variants in candidate genes, HPOs and panels (PanelApp Australia⁴⁵) suggested by clinicians. If no pathogenic or likely pathogenic variants were identified (based on the American College of Medical Genetics and Genomics [ACMG] guidelines⁴⁶), we expanded the analysis to all genes. Genes in the incidentalome (PanelApp Australia Version 0.308⁴⁵) were excluded, unless they were relevant to the patient’s phenotype as indicated by the clinician. In so doing, we adopt a judicious approach to the reporting of variants of uncertain significance (VUS) in the acute care setting, only including those that are deemed related to the patient’s phenotype and are typically very close to being classified as ‘likely pathogenic’ according to ACMG criteria⁴⁶. This stratification of VUS is recommended by the Association for Clinical Genomic Science (ACGS) in the UK.

A multidisciplinary review meeting (MDM) was then held to evaluate the results. The review MDM comprised the same clinicians and study representatives who attended the recruitment MDM. During the meeting, the genomic data analysts presented the quality control report and discussed the prioritized variants, and the evidence for pathogenic or likely pathogenic variants, for genotype-phenotype correlation. The VCGS results were not shared with the NBG team, ensuring they remained blinded to the clinically validated results until the variant review MDM. Finally, the clinical geneticist and genetic counsellor disclosed and discussed the molecular diagnosis based on the accredited acute care genomics service (VCGS) results.

Upon completing the variant review meeting, the NBG study team generated a research report. Simultaneously, the clinical laboratory (VCGS) produced its validated clinical report, which was directly returned to the clinical team for disclosure to the families. Finally, the genetic counsellor communicated the study report findings to the study participants and addressed any discrepancies identified in the reports.

Ethics statement

This study was performed in line with the principles of the Declaration of Helsinki. Ethics approval was obtained from the Northern B Health and Disability Ethics Committee for the study entitled: Newborn Genomics – Te Ira oo Te Arai (Ethics reference: 2023 FULL 15542). Locality approval was obtained from the Research Review Committee Te Toka Tumai Auckland for the project entitled: Newborn Genomics – Te Ira oo Te Arai (Reference A + 9855 [FULL 15542]).

Patient consent statement

Parents of the participating newborns provided written informed consent.

Overview of genomic benchmarking workflow for acute care clinical pipeline

We have performed genomic variant benchmarking of an expandable acute care clinical pipeline, using the set standards and guidelines provided by GAG4H²⁴. DNA samples from the Coriell repository, including the characterized GIAB reference samples (HG002 – HG007) with available truth sets, were used for benchmarking variant calling of SNVs (i.e. single base substitutions) and small indels (i.e., insertions and deletions < 50 bps) (Fig. 1). Sample HG002 was used to benchmark SVs (i.e., genomic alteration >50 bps encompassing insertions, deletions, duplications, inversions, and translocations). Coriell samples carrying pathogenic variants (i.e., GM06936, GM06870, GM01416, GM20556, GM09367, GM05966, GM05067, GM09216) were used to evaluate the performance of long reads in the identification of large-scale chromosomal abnormalities (Fig. 1).

Genomic DNA sequencing and variant calling was performed using Genome in a Bottle cell lines (HG002 – HG007) and Coriell Institute for Medical Research CNVPANEL01 cell lines. Nanopore sequencing libraries were prepared and sequenced using PromethION flow cells (R10.4). Variant calling was performed using EPI2ME pipeline that includes: a clair3 for single nucleotide variants and small indels analysis, b sniffles2 for structural variant calling; and (c) QDNAseq for copy number variant analysis. Single nucleotide variants and small indels were benchmarked against Genome in a Bottle ground-truth variants using hap.py v0.3.15 (NIST v4.2.1), and Truvari v4.1.0 was employed to benchmark structural variants (NIST v0.6). Modified base calling and alignment to GRCh38 reference genome was performed using dorado v.5.0 and genome-wide summary counts of 5-methylcytosine (5mC) and 5-hydroxymethylcytosine (5hmC) were generated using modkit v0.2.3. The figure was created using Adobe Illustrator.

Full size image

ONT long reads provide high precision for small variant calls

Following sequencing using the R10.4 pore, base calling, and variant calling, we obtained mean read lengths of 4.2 kb, mean alignment accuracy of 97.2%, read N50 of 7.6 kb, and average depth of coverage of 37.9X for six samples (Fig. 2a). Benchmarking was performed to evaluate the general performance of ONT reads on SNVs and small indels (up to 50 bps) calling from GIAB samples HG002 – HG007 (see Methods). Precision, recall, and F1 scores were computed against truth sets (National Institute of Standards and Technology [NIST)] benchmark v4.2.1; https://github.com/ga4gh/benchmarking-tools/blob/master/resources/high-confidence-sets/giab.md), for: 1) high-confidence regions excluding homopolymers, defined as four or more consecutive identical nucleotides ±1 base pair on each side; 2) genome-wide coding regions; and 3) 273 challenging medically relevant genes for the HG002 genome (CMRG v1.0; https://data.nist.gov/od/id/mds2-2475).

a Violin plots showing mean read length (kb), mean alignment accuracy (%), N50 (kb), and mean depth of coverage for Genome in a Bottle samples. Plots showing the relationship between read length and average per sequence quality scores (Q scores) for these samples are available as Supplementary Figs. 1–12. b Comparisons of precision (positive predictive value), recall (sensitivity) and F1 scores (harmonic mean of precision and recall) for single nucleotide variants and small indels called from nanopore sequencing data compared to Genome in a Bottle high-confidence regions (left) and coding regions excluding homopolymers and difficult-to-map genomic regions (right) for Genome in a Bottle samples HG002-HG007. A summary of benchmarking metrics across two separate sequencing runs for HG002-HG007 samples is available in Supplementary Table 1. CDS: coding sequence; N50: the length of the shortest read among the longest sequences, encompassing ~50% of the total nucleotides in a set of sequences.

Full size image

Across two separate sequencing runs, the average SNV precision and recall were 0.998 and 0.992), respectively, while small indel precision and recall were 0.922 and 0.831, respectively, within GIAB high-confidence regions (Fig. 2b; Supplementary Table 1; Supplementary Table 2). When assessing variants exclusively within coding regions and regions excluding homopolymers and difficult-to-map genomic regions, SNVs achieved an F1 score of 0.994, and small indels reached 0.968 (Fig. 2b; Supplementary Table 3; Supplementary Table 4). Furthermore, we assessed the performance of long reads on identifying variants in CMRG. ONT LRS demonstrated precision and recall scores > 0.967 and > 0.978, respectively, for SNVs and >0.836 and >0.701, respectively, for small indels within the 273 genes in the CMRG set (Table 1). We observed slightly improved F1 scores for small indels called from the super accuracy (sup) base called HG002 genome (i.e., F1 score = +0.01) for high-confidence regions; and an F1 score of 0.776 for the CMRG set (Supplementary Table 5; Supplementary Table 6). Overall, these results are consistent with previous benchmarking reports⁴⁷ and validate the efficacy of the EPI2ME Labs’ implementation of Clair3³⁵ in generating high-quality small variant calls comparable to gold-standard results.

Increasing sequencing depth beyond 40X does not improve small variant detection

Sequencing depth has been identified as being critical for the accurate identification of variants for the diagnosis of genetic diseases⁴⁸. We determined the optimal genomic depth for precise small variant identification from ONT reads. We downsampled the HG005 alignment file ( ~ 40X) by randomly extracting sets of reads (i.e. at proportions of 0.12, 0.25, 0.3, 0.4, 0.5, 0.6, and 0.75 of the total set) to simulate sequence data of the same sample at sequencing depths of 4.8X, 10X, 12X, 16X, 20X, 24X, and 30X. Upsampling was performed at proportions of 1.25, 1.5, 1.75, and 2 to mimic depths of 50X, 60X, 70X, and 80X. Our analysis revealed that beyond a sequencing depth of 40X, there were no significant improvements in the detection of SNVs and small indels (Fig. 3). These findings indicate that ~40X is the optimal depth for accurate small variant discovery from ONT reads, and additional depth beyond this threshold does not enhance the accuracy or sensitivity of small variant detection.

Rarefaction analysis of subsampled HG005 BAM files reveal consistent F1 scores for single nucleotide variant at 40X and 80X, with scores at 0.995 and 0.992, respectively (left panel). Similarly, F1 scores for small indels (<50 bp) vary by <2% (i.e. 0.925 and 0.941) at depths of 40X and 80X, respectively (right panel).

Full size image

ONT reads accurately identify structural variations and haplotype-specific tandem repeats

SVs were identified from the HG002 genome using Sniffles2 with a tandem repeat bed file provided to improve SV calling in repetitive regions (Methods)³⁶. SV calling performance from two separate sequencing runs, together with publicly available ONT SVs (https://labs.epi2me.io/giab-2023.05/), were benchmarked against the publicly available high-confidence GIAB ground-truth SVs (from GRCh37 reference genome) and SVs called within the CMRG genes using Truvari v4.1.0. ONT reads demonstrated high precision ( > 0.951) and recall ( > 0.943) for high-confident regions (Table 2). Additionally, we observed precision and recall metrics of > 0.889 and 0.946, respectively, for SVs within CMRG genes (Table 3), consistent with published benchmarks on ONT long reads⁴⁷.

We profiled genomic regions associated with repeat expansions since these regions contribute to the development of numerous neurodevelopmental disorders (NDDs) (e.g. congenital and childhood-onset myotonic dystrophy type 1⁴⁹). Using Straglr³⁸, we genotyped and quantified 37 clinically relevant tandem repeat regions (including DMPK and NOTCH2NLC) in HG001, HG002, HG003, and HG004 genomes. Our findings indicate that long reads enable the genotyping of these regions (Supplementary Figs. 13–16), while also providing haplotype-specific information for the repeat elements (Supplementary Fig. 17).

Accurate identification of copy number variation at 2X sequencing depth

CNVs occur in neonatal disorders (e.g., 22q11.2 deletion syndrome⁵⁰). In clinical testing, whole genome sequencing is poised to replace chromosomal microarray for the detection of CNVs⁴⁸. As such, the clinical utility of long-read sequencing in identifying CNVs has been demonstrated¹⁷. We assessed whether our variant calling workflow accurately detects clinically relevant pathogenic large chromosomal aberrations from sixteen Coriell samples (Table 4; Supplementary Figs. 18–33). Benchmarking results confirmed the reliable and accurate detection of these pathogenic CNVs using long-read sequencing (Table 4). Notably, we successfully identified two CNVs at ~2X coverage by downsampling the BAM files to ~2.6 M reads: 1) the isodicentric chromosome CNV (i.e., 47,XY,+idic(15)(q13).ish idic(15)(q13)(D15Z1 + + ,D15S11 + + ,GABRB3 + +).arr Yq11.223q11.23(23920264-27079691)x2,15q11.1q13.3(18276329-30557740)x4); and 2) the XXXX syndrome CNV (i.e., 48,XXXX) (Supplementary Fig. 34). Collectively, these results underscore the potential for long-read genomic testing in neonatal intensive care for the diagnosis of suspected genetic conditions resulting from large chromosomal events.

High concordance between ONT methylation calls and bisulphite sequencing

DNA methylation (5mC) is implicated in the pathogenic mechanism of FMR1-related disorders (e.g. fragile X syndrome), with expanded alleles typically exhibiting promoter hypermethylation and silencing of FMR1⁵¹. As such, DNA methylation profiling is essential for complete genetic diagnosis of FMR1-related disorders⁵¹. Notably, ONT sequencing facilitates concurrent profiling and quantification of DNA methylation (5-methylcytosine [5mC], and 5-hydroxymethylcytosine [5hmC]). Genome-wide 5mC characterization of the HG002 genome using ONT sequencing identified 28.8 million CpG sites (98% of total GRCh38 CpG sites). We did not detect any high levels of 5hmC methylation across the genome. Comparing ONT 5mC calls to standard whole-genome bisulphite sequencing (WGBS) of the HG002 genome, acquired from the ONT open data repository (https://labs.epi2me.io/gm24385-5mc), identified a strong correlation (r = 0.949; Fig. 4a), consistent with highly accurate methylation calling. Notably, using the ont-methylDMR-kit pipeline (https://github.com/NyagaM/ont-methylDMR-kit), we identified haplotype-specific 5mC DMRs within gene promoters for imprinted genes (Fig. 4b, c), as well as in novel DMRs (Supplementary Fig. 35; Supplementary Data 1). This illustrates the potential utility of the haplotype-level resolution offered by ONT-based sequencing reads.

a Density plots of the level of 5mC base modification across CpG sites detected through nanopore sequencing (green) and bisulphite sequencing (orange) from the HG002 sample. Pearson correlation analysis identified a strong correlation for methylation levels across CpG sites between the nanopore sequencing and whole genome bisulphite sequencing technologies. b SNURF-SNRPN and (c) PEG3 loci exhibits haplotype-specific differential methylation, linked to maternal imprinting (paternal expressed allele) as detected by ont-methylDMR-kit pipeline (https://github.com/NyagaM/ont-methylDMR-kit) and visualized using modbamtools (https://github.com/rrazaghi/modbamtools). Methylated CpG sites are denoted in red, while unmethylated sites are represented in blue. A methylation frequency plot and gene loci are visualized above the haplotype reads. A summary of genome-wide haplotype-specific 5mC differentially methylated regions and their genomic annotations is available in Supplementary Data 1.

Full size image

Application of the long-read pipeline in acute care genomic diagnosis

Our goal was to develop a scalable acute care genetic diagnostic pipeline by harnessing the capabilities of the ONT PromethION 2 solo system integrated with Fabric GEM™ (an AI-driven genomics analysis and interpretation software; https://fabricgenomics.com/). This integration was designed to facilitate precise genome annotation with rapid variant interpretation and prioritisation. The ultimate goal was to provide clinicians with access to actionable information pertaining to SNVs, small indels, and SVs (Fig. 5).

The process begins with the recruitment of critically ill newborns, infants, and children suspected to have genetic disorders admitted to the neonatal and paediatric intensive care (NICU and PICU) (in this instance at Starship Child Health) into the Newborn Genomics programme (NBG) as described in the Methods. Once consented, blood samples (from mother, father) and a buccal or blood sample (from the child) are collected in duplicate. One set of samples is sent to the NBG laboratory for sequencing, and the second to the clinical laboratory (in this instance at the Victorian Clinical Genetics Services [VCGS] in Melbourne, Australia) for variant confirmation. Samples sent to the NBG lab are sequenced using the nanopore PromethION 2 solo system and analysed as described in the Methods. Subsequently, precise genome annotation and rapid variant interpretation is conducted using Fabric GEM™. A further variant review multi-disciplinary meeting is convened as described in the Methods, and orthogonal confirmation (i.e. Sanger Sequencing) of the candidate genetic variants performed by accredited genetic testing facilities (in this instance at VCGS). Finally, clinical (from accredited genetic testing facilities) and research (from NBG) reports are generated, summarizing the evidence for the identified variants, with the clinical report forwarded to the genetic counsellor for communication of the genetic diagnosis to the participants. The figure was created using Adobe Illustrator.

Full size image

During the establishment phase of this pipeline, ten critically sick children in the neonatal and paediatric intensive care (NICU and PICU) at Te Toka Tumai/Starship Child Health were referred for rapid long-read genomic sequencing (Ethics: Approval from Health and Disability Ethics Committee [Reference 2023 FULL 15542]; Locality approval: A + 9855 [FULL15542]). In parallel, samples were provided to a clinically accredited genomics laboratory (VCGS; Melbourne, Australia) for rapid short-read Illumina genomic sequencing as described in Lunke et al. 2023⁴². Sequencing, genome variant curation and analysis were independently undertaken at each site. For the 30 samples, we obtained a mean read length of 8.2 kb, mean alignment accuracy of 96.8%, read N50 of 7.6 kb, and an average depth of coverage of 36X, with probands reaching 43.9X (Supplementary Fig. 36). Plots showing the relationship between read length and average per sequence quality scores (Q scores) for the 10 proband samples are available in Supplementary Figs. 37–46).

Genomic results from the accredited laboratory were independently provided to the clinician so that the long-read provider was unaware of the accredited results until after their results were presented to the MDM. Identical results were obtained across the ten proband-parent trios recruited into the programme. Six received a genomic diagnosis, while four did not (Table 5). The concordance of genomic findings in both positive and negative cases demonstrates the applicability and reliability of our pipeline for acute clinical care. Finally, four additional trios were sequenced (Supplementary Table 7). In one of these trios, we identified complex heterozygous genotypes in the GBA gene, consisting of a 55 bp deletion on one haplotype and a missense variant on the other (Supplementary Fig. 47). This highlights the potential of the long-read system to accurately assign variants to the correct haplotype and improve diagnostic yield, a benefit we anticipate will become apparent with larger sample sizes.

Rapid genomic diagnosis offers potential benefits for critically sick patients that include guiding clinical management and improving prognosis^16,30,52. Although long-read sequencing technology is opening new opportunities for rapid diagnosis and treatment of rare genetic disorders, its adoption in clinical settings has been limited. This is despite evidence suggesting that integrating long reads increases genetic testing capabilities beyond SNVs and small indels to include SVs, CNVs, and STRs^{13,14,15,16,17,18,19,20,21,22}, as well as the possibility of moving testing closer to point-of care⁵³. The underutilization of this technology in the acute care of patients admitted in neonatal and paediatric intensive care units represents a gap, given the potential benefits it could offer for patient care and disease management.

The incorporation of nanopore sequencing into clinical practice presents challenges given perceptions about the stability and performance of LRS technology^12,23,54,55. For example, the 97% alignment accuracy can lead to false positives. We contend that these issues are partly due to the inherent limitations of the GRCh38 reference genome. While accurate and comprehensive, it contains gaps and repetitive sequences, such as transposable elements, homopolymer regions, satellite DNA, and segmental duplications, making unique alignment difficult⁵⁶. Additionally, since the reference genome is a composite of multiple individuals, it may not perfectly match any single genome, leading to alignment discrepancies⁵⁶. To address these concerns, we have extensively benchmarked and validated our pipeline, focusing on ensuring consistent and stable performance. Furthermore, we contend that a sequencing depth of 40X, particularly for the sick probands, is sufficient to enhance statistical confidence in genotyping variant calls, even in homopolymer regions, as previously reported⁵⁷. Moreover, in line with clinical and technical recommendations^45,58,59, we have developed standardized protocols and quality control measures for long-read data to facilitate the integration of nanopore sequencing into clinical practice. These protocols ensure data consistency and reliability across different laboratory personnel, analytical teams, and platforms. Implementing such standardization measures is crucial for enabling the broader adoption of LRS technology in clinical settings.

Our established acute care clinical workflow integrates a highly scalable genome sequencing platform (i.e. PromethION 2 solo) with an AI-driven genomics analysis and interpretation software (Fabric GEM™; https://fabricgenomics.com/) to facilitate precise genome annotation, rapid variant interpretation, and prioritization. This pipeline has been established to provide clinicians with actionable information regarding SNVs, small indels, and SVs. In the initial benchmarking phase, we employed GA4GH benchmarking tools and GIAB samples sequenced to >30X coverage, demonstrating highly accurate variant calling metrics, particularly in high-confidence and coding regions of the human genome. Additionally, our workflow accurately detected large-scale CNVs in all sixteen samples from the Coriell Copy Number Variation Reference Panel 1 (CNVPANEL01), showcasing the capability of genotyping STRs and reliably profiling DNA methylation genome-wide. We believe that future updates to the latest ONT library preparation Kit 14 chemistry, along with R10.4.1 flow cells and continuous updates to the base calling models, will further improve precision and recall of SNVs, small indels, and SVs⁶⁰.

Finally, in the validation phase, the clinical utility of long-read trio sequencing was demonstrated through the concordance of genomic findings with an established Acute Care Genomics service (VCGS; Melbourne, Australia), which utilized Illumina short-read sequencing technology⁴². Notably, this concordance was achieved in all ten acute cases examined. However, we did not validate the limit of detection of mosaicism for SNV, small indels, CNVs and SVs. Nonetheless, our findings reinforce the potential of long-read sequencing for comprehensive genomic analysis and its applicability in clinical diagnostics.

The implementation of nanopore sequencing in the clinical care pathway is crucial^{16,17,18,42,48,53}. This study demonstrates the feasibility of such a pathway with the availability of funding, technology, skilled laboratory personnel, and researchers supporting rapid genomic testing for critically ill patients. The specialist multidisciplinary team (MDT) model is ideal for complex cases and provides clinicians with input for rare diagnoses that often lack established clinical management guidelines^61,62. In our opinion, acute care genetics requires a turnaround-driven approach to optimize care and offer economic benefits. In our experience to date, the 5-day average turnaround time from our study seems to meet clinical requirements⁶³. In the evaluation of the cohort reported in this study, some cases remained undiagnosed, with at least one being non-genetic. In one case, the diagnosis of biallelic PCSK1 congenital diarrhoea⁶⁴ informed management and surveillance, with the family being reassured about the self-limiting nature of the diarrhoea and surveillance initiated through paediatric endocrine services. Three families used the results for reproductive risk assessment, while reanalysis of the negative trios will continue in the research setting. As expected in critically ill patients, many babies died. However, obtaining a diagnosis offers closure, enabling families to plan subsequent pregnancies.

Deployment of a system like this in New Zealand, a country that arguably has a challenging geography and population distribution, can be managed in two ways. One option is to use the Oxford Nanopore P2 machines as part of a distributed system. In theory, the distributed system can comprise of fully independent, or interdependent nodes that incorporate machines in different locations to optimize return times to a central analytical facility. However, both models could place significant strain on the workforce, requiring highly skilled workers who are already in short supply. The alternative is to centralize facilities in the main population centres in the North and South of the country. Despite New Zealand’s challenging geography, we have systems in place to transport samples to these centralized location within hours of collection. These rapid transit programs ensure that turnaround time is not significantly impacted, and the model reduces workforce impacts. We believe this approach will be effective in practice.

In conclusion, we have successfully implemented a scalable clinical pipeline for rapid trio long-read whole-genome sequencing in an acute care setting, aiming to provide prompt and actionable genomic information to clinicians. Through this effort, we have demonstrated the feasibility of achieving rapid precision medicine for critically sick children on a national scale using long-read technology.