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

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

Introduction

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 cancers1,2,3,4,5,6,7. Over the past decade, the increasing adoption of genomic testing has generated substantial evidence supporting precision medicine2,4,6. Such approaches have enabled molecular diagnoses for genetic disorders, guiding tailored medical interventions6,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 limitations10,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 modifications10. Thirdly, short read lengths hinder the identification and precise phasing of alleles in large SVs10.

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-time12,13,14,15. Consequently, long reads generate highly reliable complete genome assemblies15, 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 demonstrated16,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 SVs13,14,19. Novel pathogenic variants have been uncovered by LRS technology in human diseases with a previously unknown underlying genetic cause20. Additionally, long reads facilitate the characterization of pathogenic repeat expansions in genomic regions that are challenging to sequence using short-read sequencing technology21,22.

The clinical application of LRS13,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 sequences12,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 practice24,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)24. 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 SVs27,28. Notably, GIAB recently provided a curated benchmark of challenging medically relevant genes through haplotype-resolved whole-genome assembly29. The GA4GH resources enable performance assessment, optimization, and analytical validation of CGS assays and workflows for detecting genomic variations24,25. Indeed, GIAB datasets and benchmarks are considered the gold standard for evaluating sequencing technologies and variant calling pipelines27.

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 hours30. 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 run31. DNA and RNA base modifications, including 5-methylcytosine (5mC), 5-hydroxymethylcytosine (5hmC), and N6-methyladenine (6 mA), can also be detected computationally on raw ONT data without the need to perform special library preparations such as bisulphite conversion32, which is known to cause DNA damage and can lead to overestimation of the 5mC level33.

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.

Methods

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)34. 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)35, SV calling with Sniffles2 (v2.2)36, and CNV calling with QDNAseq (v1.38)37 using default parameters, with a VNTR annotation file provided for accurate SV identification. Repeat expansions were genotyped using Straglr (https://github.com/philres/straglr)38 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)24 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)24 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:

$${Precision}={True; Positives}/({True; Positives}+{False; Positives})$$
(1)
$${Recall}={True; Positives}/({True; Positives}+{False; Negatives})$$
(2)
$$F1,{score}=frac{2,xleft({Precisionx; Recall}right)}{left({Precision}+{Recall}right)}$$
(3)

For SVs, we employed Truvari v4.1.0 (https://github.com/ACEnglish/truvari)39 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 Samtools40, 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 genes29 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) package41. We used the following DSS41 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., 202342.

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.43 and McKeown et al.44.

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. 202342.

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 GEMTM 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 Australia45) suggested by clinicians. If no pathogenic or likely pathogenic variants were identified (based on the American College of Medical Genetics and Genomics [ACMG] guidelines46), we expanded the analysis to all genes. Genes in the incidentalome (PanelApp Australia Version 0.30845) 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 criteria46. 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.

Results

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 GAG4H24. 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).

Fig. 1: An overview of the genomic benchmarking for the acute care clinical pipeline.
Benchmarking nanopore sequencing and rapid genomics feasibility: validation at a quaternary hospital in New Zealand

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).

Fig. 2: Sequencing metrics and variant calling accuracy ONT summary statistics.
figure 2

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 reports47 and validate the efficacy of the EPI2ME Labs’ implementation of Clair335 in generating high-quality small variant calls comparable to gold-standard results.

Table 1 Benchmarking metrics for SNVs and small indels within CMRG
Full size table

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 diseases48. 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.

Fig. 3: Sequencing depth exceeding 40X does not improve single nucleotide variant or small indels detection.
figure 3

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)36. 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 reads47.

Table 2 Performance metrics for identification and genotyping of high-confident SVs
Full size table
Table 3 Performance metrics for identification and genotyping SVs within CMRG genes
Full size table

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 149). Using Straglr38, 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 syndrome50). In clinical testing, whole genome sequencing is poised to replace chromosomal microarray for the detection of CNVs48. As such, the clinical utility of long-read sequencing in identifying CNVs has been demonstrated17. 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.

Table 4 Long reads accurately detect pathogenic CNVs from CNVPANEL01 samples
Full size table

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 FMR151. As such, DNA methylation profiling is essential for complete genetic diagnosis of FMR1-related disorders51. 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.

Fig. 4: Strong correlation between 5mC methylation calls from nanopore sequencing and whole-genome bisulphite sequencing.
figure 4

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).

Fig. 5: Flowchart for a scalable acute care clinical pipeline for rapid genome sequencing.
figure 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. 202342. 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.

Table 5 Clinical characterisation and genomic results from ten proband-parent trio genome sequencing using Illumina short reads and nanopore long reads
Full size table

Discussion

Rapid genomic diagnosis offers potential benefits for critically sick patients that include guiding clinical management and improving prognosis16,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 STRs13,14,15,16,17,18,19,20,21,22, as well as the possibility of moving testing closer to point-of care53. 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 technology12,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 difficult56. Additionally, since the reference genome is a composite of multiple individuals, it may not perfectly match any single genome, leading to alignment discrepancies56. 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 reported57. Moreover, in line with clinical and technical recommendations45,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 SVs60.

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 technology42. 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 crucial16,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 guidelines61,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 requirements63. 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 diarrhoea64 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.

Related Articles

SVLearn: a dual-reference machine learning approach enables accurate cross-species genotyping of structural variants

Structural variations (SVs) are diverse forms of genetic alterations and drive a wide range of human diseases. Accurately genotyping SVs, particularly occurring at repetitive genomic regions, from short-read sequencing data remains challenging. Here, we introduce SVLearn, a machine-learning approach for genotyping bi-allelic SVs. It exploits a dual-reference strategy to engineer a curated set of genomic, alignment, and genotyping features based on a reference genome in concert with an allele-based alternative genome. Using 38,613 human-derived SVs, we show that SVLearn significantly outperforms four state-of-the-art tools, with precision improvements of up to 15.61% for insertions and 13.75% for deletions in repetitive regions. On two additional sets of 121,435 cattle SVs and 113,042 sheep SVs, SVLearn demonstrates a strong generalizability to cross-species genotype SVs with a weighted genotype concordance score of up to 90%. Notably, SVLearn enables accurate genotyping of SVs at low sequencing coverage, which is comparable to the accuracy at 30× coverage. Our studies suggest that SVLearn can accelerate the understanding of associations between the genome-scale, high-quality genotyped SVs and diseases across multiple species.

Rapid brain tumor classification from sparse epigenomic data

Although the intraoperative molecular diagnosis of the approximately 100 known brain tumor entities described to date has been a goal of neuropathology for the past decade, achieving this within a clinically relevant timeframe of under 1 h after biopsy collection remains elusive. Advances in third-generation sequencing have brought this goal closer, but established machine learning techniques rely on computationally intensive methods, making them impractical for live diagnostic workflows in clinical applications. Here we present MethyLYZR, a naive Bayesian framework enabling fully tractable, live classification of cancer epigenomes. For evaluation, we used nanopore sequencing to classify over 200 brain tumor samples, including 10 sequenced in a clinical setting next to the operating room, achieving highly accurate results within 15 min of sequencing. MethyLYZR can be run in parallel with an ongoing nanopore experiment with negligible computational overhead. Therefore, the only limiting factors for even faster time to results are DNA extraction time and the nanopore sequencer’s maximum parallel throughput. Although more evidence from prospective studies is needed, our study suggests the potential applicability of MethyLYZR for live molecular classification of nervous system malignancies using nanopore sequencing not only for the neurosurgical intraoperative use case but also for other oncologic indications and the classification of tumors from cell-free DNA in liquid biopsies.

Comparative evaluation of SNVs, indels, and structural variations detected with short- and long-read sequencing data

Short- and long-read sequencing technologies are routinely used to detect DNA variants, including SNVs, indels, and structural variations (SVs). However, the differences in the quality and quantity of variants detected between short- and long-read data are not fully understood. In this study, we comprehensively evaluated the variant calling performance of short- and long-read-based SNV, indel, and SV detection algorithms (6 for SNVs, 12 for indels, and 13 for SVs) using a novel evaluation framework incorporating manual visual inspection. The results showed that indel-insertion calls greater than 10 bp were poorly detected by short-read-based detection algorithms compared to long-read-based algorithms; however, the recall and precision of SNV and indel-deletion detection were similar between short- and long-read data. The recall of SV detection with short-read-based algorithms was significantly lower in repetitive regions, especially for small- to intermediate-sized SVs, than that detected with long-read-based algorithms. In contrast, the recall and precision of SV detection in nonrepetitive regions were similar between short- and long-read data. These findings suggest the need for refined strategies, such as incorporating multiple variant detection algorithms, to generate a more complete set of variants using short-read data.

The role of gene copy number variation in antimicrobial resistance in human fungal pathogens

Faced with the burden of increasing resistance to antifungals in many fungal pathogens and the constant emergence of new drug-resistant strains, it is essential to assess the importance of various resistance mechanisms. Fungi have relatively plastic genomes and can tolerate genomic copy number variation (CNV) caused by aneuploidy and gene amplification or deletion. In many cases, these genomic changes lead to adaptation to stressful conditions, including those caused by antifungal drugs. Here, we specifically examine the contribution of CNVs to antifungal resistance. We undertook a thorough literature search, collecting reports of antifungal resistance caused by a CNV, and classifying the examples of CNV-conferred resistance into four main mechanisms. We find that in human fungal pathogens, there is little evidence that gene copy number plays a major role in the emergence of antifungal resistance compared to other types of mutations. We discuss why we might be underestimating their importance and new approaches being used to study them.

Integrated proteogenomic characterization of ampullary adenocarcinoma

Ampullary adenocarcinoma (AMPAC) is a rare and heterogeneous malignancy. Here we performed a comprehensive proteogenomic analysis of 198 samples from Chinese AMPAC patients and duodenum patients. Genomic data illustrate that 4q loss causes fatty acid accumulation and cell proliferation. Proteomic analysis has revealed three distinct clusters (C-FAM, C-AD, C-CC), among which the most aggressive cluster, C-AD, is associated with the poorest prognosis and is characterized by focal adhesion. Immune clustering identifies three immune clusters and reveals that immune cluster M1 (macrophage infiltration cluster) and M3 (DC cell infiltration cluster), which exhibit a higher immune score compared to cluster M2 (CD4+ T-cell infiltration cluster), are associated with a poor prognosis due to the potential secretion of IL-6 by tumor cells and its consequential influence. This study provides a comprehensive proteogenomic analysis for seeking for better understanding and potential treatment of AMPAC.

Responses

Your email address will not be published. Required fields are marked *