Toward a Kinh Vietnamese Reference Genome: Constructing a De Novo Genome Assembly Using Long-Read Sequencing and Optical Mapping

TL;DR

This study creates a high-quality reference genome for the Kinh Vietnamese population using advanced sequencing techniques, improving genomic accuracy for this group.

Contribution

The paper introduces a new Kinh Vietnamese reference genome (VHG1.2) using long-read and optical mapping technologies.

Findings

The VHG1.2 assembly has high accuracy (QV: 48), completeness (BUSCO: 92%), and continuity (N50: 50 Kbp).

Using VHG1.2 revealed significant genetic variants compared to the standard hg38 reference genome.

The hybrid sequencing approach proved effective for de novo assembly of population-specific genomes.

Abstract

Background: Population-specific reference genomes are essential for improving the accuracy and reliability of genomic analyses across diverse human populations. Although Vietnam ranks as the 16th most populous country in the world, with more than 86% of its population identifying as Kinh, studies specifically focusing on the Kinh Vietnamese reference genome remain scarce. Therefore, constructing a Kinh Vietnamese reference genome is valuable in the genetic research of Vietnamese. Methods: In this study, we combined PacBio long-read sequencing and Bionano optical mapping data to generate a de novo assembly of a Kinh Vietnamese genome (VHG), which was subsequently polished using multiple Kinh Vietnamese short-read whole-genome sequences (WGSs). Results: The final assembly, named VHG1.2, comprised 3.22 gigabase pairs of high-quality sequence data, demonstrating high accuracy (QV: 48), completeness (BUSCO: 92%), and continuity (295 super scaffolds, super scaffold N50: 50 Kbp). Using multiple bioinformatic tools for variant calling, we observed significant variants when the population-specific reference VHG1.2 was used compared to the standard reference genome hg38. Conclusions: Overall, our genome assembly demonstrates the advantages of a long-read hybrid sequencing approach for de novo assembly and highlights the benefit of using population-specific reference genomes in population genomic analysis.

Figures

Full text

1. Introduction

Over the past few decades, human genomics and precision medicine have undergone tremendous advancements with the proliferation of multiple next-generation sequencing technologies. Short-read sequencing technologies in particular have been instrumental; they are routinely used in large-scale population analyses and molecular diagnostic applications owing to their high throughput, low cost, and high-accuracy output [1,2,3]. To this day, Illumina remains the most dominant sequencing platform on the market, responsible for generating more than 90% of the world’s sequencing data and serving as the primary sequencing platform for diagnostic genome sequencing globally [4,5,6]. However, these technologies have inherent limitations. Short-read sequencing requires alignment to a well-defined reference; if a short read is aligned to the wrong position, any difference between the reads and the reference can be misinterpreted as a mutation [3,7,8]. This is especially true regarding genomic variations between individuals or regions that short reads can only map to ambiguously [9,10,11].

Since the completion of the Human Genome Project and the release of the finished sequence in 2004, there remains a notable bias in genomic research, with a predominant focus on populations of European descent [12,13,14]. The most widely used human reference genomes, GRCh37 (hg19) and GRCh38 (hg38) are primarily constructed from individuals of European and African ancestry, resulting in the underrepresentation of other populations [15,16,17]. This imbalance introduces reference bias, where reads containing alternative alleles common in underrepresented populations are less likely to align correctly to the reference genome. This misalignment can lead to inaccurate variant calling and potentially misleading interpretation in downstream analyses [18,19,20].

To improve the reference accuracy for short-read sequencing analysis and to address the gap in representation, numerous groups across the world have constructed population-specific reference genomes from high-coverage, de novo assembled individual human genomes, such as the Chinese genomes [21,22,23], Korean [24,25], Danish [26], Swedish [27], and Japanese [8,28]. All of these genomes are continuously being upgraded to improve their completeness and specificity [29,30,31,32]. Subsequent analyses using these population-specific genomes revealed substantial differences in mapping and variant calling results [20,28,33]. These findings highlighted the necessity of developing population-specific reference genomes to enhance the reliability of genomic analyses across diverse human populations.

Long-read sequencing technologies have now achieved a level of accuracy and throughput that significantly expands their applications in genomics. Technologies such as single-molecule real-time (SMRT) sequencing by Pacific Biosciences and nanopore sequencing by Oxford Nanopore Technologies generate reads that can span tens of kilobases (Kbp) to megabases (Mbp) [1,34]. These long reads enable the resolution of complex and repetitive regions, facilitating more accurate and contiguous genome assemblies [35,36,37].

In addition to long-read sequencing, optical mapping has been intensively used to provide genome scaffolding and detect large-scale structural variations. Optical maps have an average molecule length of 225 Kbp and can easily span genomic regions that are difficult for DNA sequencing while not having the limitations of linked-read sequencing technologies [37]. Recent studies have demonstrated the success of combining long-read sequencing with optical mapping to produce high-quality de novo [27,38,39]. These combined approaches are particularly effective in creating de novo assemblies comparable to hg38 in terms of completeness and contiguity.

Previous studies have shown that the population structure of the Kinh Vietnamese (KHV), a population of 100 million people, has a main component that is most prevalent in the Chinese population and significantly diverges from not only European and African populations but also South Asian neighboring populations [10,40]. Thus, there have been calls to construct a population-specific reference genome to address the inherent biases in using hg38 reference for the Vietnamese population. A previous attempt to construct a KHV reference genome was limited by using exclusively Illumina’s short-read sequencing technology [41], which is not a viable platform to use in the context of creating a low-error assembly without a reference guide [42]. In this study, we combined PacBio HiFi sequencing and Bionano optical mapping to generate a high-quality de novo assembly of a Kinh Vietnamese genome and polish this construct using other Vietnamese short-read whole-genome sequences (WGSs) to create a population-specific genome. Though still in need of further work, our genome assembly represents a significant step toward a new population-specific reference for genomic analysis using short-read sequencing technologies in Vietnam.

2. Materials and Methods

2.1. Ethics Declarations

This study and sample collection were approved by the Scientific Council of the Institute of Biotechnology, Vietnam Academy of Science and Technology, on 1 June 2021.

2.2. Donor Selection

Three adult male Vietnamese volunteers, aged 30–59 years, were recruited for the study. This age range was chosen based on previous studies indicating it is the most biologically stable period with minimal genetic and health-related variability, which could influence genetic results [43]. The donors self-reported being healthy, with no inherited conditions, and were of Kinh ancestry, with their families having lived in Vietnam for more than eight generations. The selected donor has a healthy son without any inherited conditions or medical history.

2.3. Genomic DNA Preparation

Genomic DNA was extracted from fresh blood using the Blood & Cell Culture DNA Midi Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s protocol. The resulting DNA molecules were 80–100 Kbp in size. DNA quality was assessed through gel electrophoresis, a Qubit 3 Fluorometer (Invitrogen, Waltham, MA, USA), and a 4200 TapeStation System (Agilent, Santa Clara, CA, USA) (Figure S1).

2.4. Pacbio HiFi Library Preparation and Sequencing

PacBio HiFi library preparation and sequencing were performed following the manufacturer’s standard protocol with minor modifications. Genomic DNA was sheared into fragments of 10–15 Kbp using g-tubes (Covaris, Brighton, UK) and processed using the SMRTbell Express Template Prep Kit 2.0 (PacBio). Each library preparation step was followed by cleanup using AMPure PB beads (PacBio, Menlo Park, CA, USA) and quality control using the Bioanalyzer 2100 system (Agilent). The final DNA library was loaded onto SMRT Cells (PacBio, Menlo Park, CA, USA) and sequenced using the PacBio Sequel and Sequel II systems, with sample setup in the SMRT Link portal version 9.0.

2.5. Bionano Optical Map Generation

Genome optical mapping was performed using the Saphyr system (Bionano Genomics) following the manufacturer’s protocol. A total of 1.5 mL of deep-frozen whole blood was shipped on dry ice to Bionano Genomics (San Diego, CA, USA) and processed according to the Bionano SP Blood and Cell Culture DNA Isolation protocol for 100X Human Genome Sample Analysis. Genome assembly, scaffolding, and structural variant detection were performed using Bionano Solve software version 3.2.

2.6. Genome Assembly of SMRT Sequencing Reads

The raw data were initially processed and evaluated through a primary analysis on a host computer before being transferred to a secondary analysis server via a local network. The data were analyzed using SMRTLink Portal v9.0 with a minimum predicted accuracy of 0.99. HiFi contig assembly was generated from HiFi reads using the HiFiasm [44,45] version 0.16.1-r37.

2.7. Hybrid Scaffold Assembly

Scaffolds constructed from Bionano data were combined with the HiFiasm contigs through Bionano Access version 1.7 on the Bionano platform to create super scaffolds assembly with parameters “Resolve Conflicts” and “Trim Overlapping Sequence = On” selected [38]. Conflicts between the assembly and optical maps were automatically resolved in favor of the optical maps. Later, a second round of super scaffolding was performed between the output of the hybrid process and the consensus optical map from the DLE-1 restriction enzyme [25,27,46]. The resulting hybrid scaffold assembly was termed super scaffolds.

2.8. Assembly Correction and Improvement with Pacbio HiFi Reads

The super scaffolds underwent polishing and correction with HiFi reads using Inspector [47] version 1.0.1, where it was applied with the error correction module and datatype PacBio HiFi to generate a draft assembly. The corrected and polished assembly was named VHG1.1.

2.9. Variant Calling for VHG

A small variant calling for VHG was performed in parallel using GATK version 4.4.0.0, according to a previous publication [38]. First, HiFi reads were mapped to reference genome hg38 using Minimap2 [48] version 2.24 (r1155), and the resulting data were converted to FASTQ format via Samtools to be used as input data for variant calling in the next steps. After removing duplicates, the alignment data were used to identify single-nucleotide polymorphisms (SNPs) and short insertions–deletions (indels) with default parameters. Then, filtered small variants were subjected to variant calling.

Structure variants (SVs) from HiFi reads were detected using Pbsv version 2.9.0 from PacBio (https://github.com/PacificBiosciences/pbsv (accessed on 6 May 2024)). The HiFi reads were aligned to the hg38 reference genome using Pbmm2 version 1.16.0 with default parameters and processed to identify SVs [46]. SVpack was used to filter the called SVs in VCF format using a series of filtering operations (https://github.com/PacificBiosciences/svpack (accessed on 7 May 2024)). Initially, only SV calls with a “PASS” flag in the “FILTER” column were filtered to obtain confident SV calls, and duplicates (if any) were removed for the following stage of the analysis. Next, the SVs were filtered based on the SV types (insertions, duplications, and deletions) to obtain each call set. Finally, the SVs were filtered for length ≥ 50 bp [27,46].

2.10. Polishing Assembly Using Illumina Short-Read Whole-Genome Sequences

To refine VHG for the Kinh Vietnamese population, the VHG1.1 assembly underwent further polishing using Jellyfish-based Assembly Sequence Polisher for Error Reduction (JASPER) tool with optimized parameters (-k 37, -t 24, and -p 4) to handle the inherent complexity of genomic assembly. Given the critical importance of representing the genetic diversity of the Kinh population, seven male Kinh Vietnamese (KHV) genomes from Ho Chi Minh City, Vietnam (HG01842, HG01846, HG02064, HG02076, HG02122, HG02058, HG01872), were sourced from the 1000 Genomes Project (1KGP) database (http://www.internationalgenome.org (accessed on 14 October 2024)). These genomes, totaling approximately 415 Gbp of Illumina sequencing data, were used to refine and enhance the VHG1.1 assembly toward population-representative reference, resulting in the VHG1.2 assembly.

2.11. Ordering and Aligning of Super Scaffolds to Chromosomes

The assembly was positioned onto the chromosomes of the hg38 reference genome using RagTag [49] version 2.1.0 with scaffold command [1,37]. Each group of assigned super scaffolds was double-checked through pairwise genome alignments to hg38 through Lastz [50] version 1.04.15 with the following parameters: --gfextend –chain –gapped. The data were visualized using R packages karyoploteR and GenomicRanges version 4.4.1.

Assemblies were mapped with hg38 and T2T references using Minimap2 version 2.24 (r1155) [48]. D-GENIES dotplot was used to quickly compare and visualize the alignment between the VHG1.2 assembly and hg38 or T2T references [51].

2.12. Genome Assembly Evaluation

Genome assemblies were evaluated by Inspector [47] version 1.0.1 with the datatype HiFi option for evaluation assembly using HiFi reads and the hg38 reference genome. Assembly completeness was assessed with Busco [52] version 5.5.0 with the -l parameter and primates_odb10 dataset. Statistical results were exported and compared to identify the highest-quality assembly. The final assembly was also evaluated with the QUAST toolkit [53] version 5.0.2.

2.13. Reference Genome Set

The hg38 (accession GCA_000001405.15) reference genome used in this study was described in a previous study [54]. It includes primary hg38 sequences (autosomes, chromosomes X and Y), mitochondrial genome, and unlocalized scaffolds that belong to a chromosome without a definitive location and order.

2.14. SNP and Indel Calling of KHV Genome

Three low-coverage WGS KHV genomes (HG01852, HG02070, and HG02079) were sourced from the 1KGP to be utilized for variant calling analysis. Common mapped reads between VGH1.1, VHG1.2, and hg38 were identified and extracted according to the method by Takayama et al. in 2021 [38]. In detail, the Illumina reads of each genome were mapped to each reference genome using bwa-mem [55] version 0.7.17-r1188 and processed with Samtools version 1.3.1. Then, mapped reads are extracted from bam files (-F 4 option), and their read names are isolated. The read names are then sorted and filtered for uniqueness using sort and uniq. Finally, common read names are identified using the awk command, and the corresponding reads are filtered from the bam files with the -N option in Samtools.

The output files of the above step were used for small variant calling with bioinformatic tools, have been evaluated, tested, and benchmarked, and are widely used in clinical genetics. First, we employed GATK version 4.4.0.0, the “industry standard” for short SNPs and indels calling [38,56]. The process starts by removing duplicate reads with MarkDuplicatesSpark and generates raw variant calls using HaplotypeCaller. Variants are categorized into SNPs and indels using SelectVariants, and both types were filtered based on quality metrics like Quality-by-Depth (QD < 2.0), low quality scores (QUAL < 30.0), high strand bias (FS > 200.0), or positional bias in reads (ReadPosRankSum < −20.0) with VariantFiltration. Finally, the remaining high-confidence variants were counted with VCF tools.

To prove that VHG1.2 is the first version of the Kinh Vietnamese population-specific reference genome, the same output files from Samtools were used to recall SNPs and indels with DeepVariant via Docker, a deep learning-based variant caller, following these default parameters: --model_type = WGS for whole-genome sequencing, --ref for the reference genome (hg38.fa, VHG1.1.fa, or VHG1.2.fa), and --reads for the BAM file containing mapped reads [57]. The variant calling data were exported to output_vcf and --output_gvcf files. The variants were counted using Excel for Microsoft 365.

2.15. Structural Variant Calling of KHV Genome

Cue, a deep learning SV calling framework [58], was utilized to detect SVs larger than 5 Kbp in KHV genomes (HG02134, HG01840, HG01852) from the 1KGP. The FASTQ sequencing reads were mapped with the hg38 reference genome and VHG1.2 by bwa-mem [50] version 0.7.18. The resulting BAM files and index files were processed using Cue with model cue.v2.pt.

2.16. Comparison VHG1.2 with East Asian Assemblies

To evaluate the genetic similarity of VHG1.2 to other East Asian assemblies, genome-to-genome alignments were performed using Minimap2 against hg38 and three publicly available East Asian reference genomes: AK1 (Korean, GCA_001750385.2), HX1 (Chinese, GCA_001708065.2), and JG1 (Japanese, GCA_014905855.1).

2.17. Data Visualization and Summary Statistics

Statistical analysis and data visualization were performed using Excel, R version 4.4.1, and Python version 3.9.16.

3. Results and Discussion

3.1. De Novo Assembly of a Kinh Vietnamese Genome

The study commenced with HiFi sequencing and optical mapping of DNA extracted from the blood of an anonymous Kinh Vietnamese donor with a normal karyotype who passed a general health examination (Figure S2). In total, 20X whole-genome coverage was achieved from SMRT and 1600 Gbp long-read sequencing with a mean subread length of 10.29 Kbp and an N50 subread length of 12.19 Kbp. Quality control measures using SMRTLink Portal v9.0 revealed 5.28 million subreads with a mean read length of 11.73 Kbp and read quality of Q36 (Tables S1 and S2).

A de novo assembly was performed from PacBio CCS reads, using HiFiasm software version 0.16.1-r37, generating 2198 HiFi contigs (N50 = 8.78 Mbp), with the longest contig at 78.82 Mbp and a total size of 3.01 Gbp (Table 1). In addition, a VHG mitochondrial sequence was assembled using a similar approach (Figure S3).

At the same time, an optical map with a depth coverage of 129.57 X was also generated (Table S3), as previously described, with a label density of 15,77 per 100 Kbp [38]. After filtering, this optical map was de novo assembled to produce a genome map. The size of the optical genome map was 5.8 Gbp, containing 1042 optical scaffolds. The N50 of the scaffolds was approximately 39.14 Mbp, and the longest scaffold was 108 Mbp (Table S4).

To obtain whole-genome assemblies, the HiFi contigs were aligned to the optical mapping scaffolds using the Bionano Access platform, forming super scaffolds [11,38]. This process resolved many gaps in the HiFi reads (Figure S4). The hybridization of 2198 HiFi contigs and 1042 Bionano scaffolds produced 295 super scaffolds, with an N50 of 50.6 Mbp (Table S4). When compared to other assemblies using the same method, our Kinh Vietnamese assembly (VHG) showed significant improvements in both sequencing and assembly metrics (Table 1). Notably, VHG had fewer super scaffolds than the most up-to-date JG1 hybrid assembly (295 and 624, respectively), suggesting that VHG is a more optimized scaffold-level assembly.

3.2. SNPs and Indels Detection of VHG’s HiFi Reads

The study identified 3,905,585 SNPs and 2,095,662 indels in the VHG by comparing HiFi reads with hg38 (Table S5). These numbers can be compared with 4,823,475 SNPs and 974,100 indels (length ≤ 100 bp) detected in the KHV [41] and 3,518,309 SNPs and 378,750 indels detected in the HX1 genome [21]. The differences in the number of small variants can be attributed to variations in sequencing technologies, variant calling approaches, and population characteristics.

Structural variants (SVs) in the VHG were identified through a comparative analysis with the GRCh38 genome assembly. The analysis found a total of 21,760 SVs, with sizes starting at 50 bp, comprising 9514 insertions, 8871 deletions, 60 inversions, and 3315 duplications (Table S6). The total length of the SVs spanned approximately 12.9 Mbp, accounting for roughly 0.4% of the genome. Notably, most SVs (90.53%) are smaller than 1 Kbp, reflecting a similar trend observed with ZF1 SVs (He et al. 2020) (Table S6). Among the 21,760 SVs in VHG, 20,400 SVs, accounting for 93.75%, were assigned to 22 chromosomes (including X and Y chromosomes), with no inversion found on chromosomes 13, 15, 18, 19, and 22 (Figure 1a). The maximum length of insertion, deletion, inversion, and duplication SVs were 15,835, 76,914, 90,932, and 681,903 bp, respectively (Figure 1b, Table S6). The median lengths of insertion, deletion, inversion, and duplication SVs, excluding outliers, were 244, 128, 1567, and 103 bp, respectively (Table S6, Figure 1c,d). The differences in the numbers of VHG SVs also can be compared with the 18,210, 17,900, 20,175, and 17,936/17,687 SVs detected in the AK1, ZF1, HX1, and Swe1/2 genomes, respectively [21,26,28,47]. Our findings indicated an equivalent analysis in detecting structural variants compared to previous studies on similar PacBio long-read datasets, which reported structural variant variation across various populations and highlighted the sensitivity of long-read sequencing technologies employed in research [21]. Additionally, the diversity and abundance of SVs observed in the VHG Kinh genome emphasize the genetic uniqueness of the Vietnamese population. This is in line with other recent genomic studies that underscore the importance of including diverse ethnic groups to capture the full spectrum of human genetic variation. This comprehensive catalog of SVs not only enhances our understanding of the Kinh Vietnamese genome but also contributes to the broader field of human genomics, providing valuable insights into the genetic diversity and structural complexity present within Southeast Asian populations.

Furthermore, an analysis of the distribution of SVs across the 24 chromosomes revealed a strong correlation between the length of our assembled chromosome and the number of SVs (R2 = 0.83, P = 4.59 × 10^−5^) (Figure 1e). This pattern has also been observed in a previous study [47].

3.3. Polishing and Quality Evaluation of the De Novo Assembly

The hybrid super scaffold assembly was evaluated and corrected using Inspector to generate higher-quality contigs, removing misassembles and false structural variants (Figure 2). This process improved completeness, as measured by BUSCO score from 91.5% to 92, and reduced the number of small-scale assembly errors per Mbp decreased from 18.2 to 2.44. The corrected assembly, named VHG1.1, had a quality value (QV) evaluated by Inspector was 48.45%, corresponding to an accuracy of 99.9% and high contiguity with an N50 of 50.627 Mbp (>23 Mbp) [47,59] (Table S7).

To avoid the idiosyncratic nature of using the sequence of a single individual as a population-specific representative, we aimed to generate a more population-reflective reference. A typical KHV genome can differ from the reference human genome by 4 to 4.2 million sites; however, among these variants, around 40,000 are common (>5%) within the Kinh Vietnamese population [10]. In this study, the VHG1.1 assembly was “corrected” to VHG1.2 by removing rare and uninformative homozygous variants where other KHV genomes had consensus, thus adjusting these positions to the major, “common” allele in the population. JASPER polishing did not alter the assembly’s contiguity values; however, it increased the number of misaligned sequences to hg38 in VHG1.2 compared to VHG1.1.

Throughout the polishing process, all assemblies maintained a total genome length of 3.22 Gbp, closely approximating the expected size of the hg38 reference genome (3.21 Gbp). The longest scaffold length remained consistent across the assemblies, at approximately 145.9 Mbp, reflecting stability in long scaffold generation. The N50 values, a measure of assembly continuity, ranged from 50.63 to 50.64 Mbp across the assemblies, with minimal variation between them (Table 2, Table S8). 99.8% of contigs were mapped to scaffolds, and 95% of super scaffolds (283/295) were mapped to the 24 chromosomes. The size of each chromosome super scaffold and the unknown sequence are summarized in Table S9.

Comparing the final assembly VHG1.2 to the hg38 and T2T references showed high (>75%) identity scores of 86.3% and 92.1%, respectively (Table 2, Table S10, Figure S5). Notably, 12 super scaffolds totaling 8.78 Mbp failed to align with all references. Similar results were observed in previous publications, where the unaligned sequences ranged from 10.3 Mbp to 13.8 Mbp [1,28], and even in hg38 with 11.5 Mbp of unplaced and unlocalized sequences [16]. This improved genomic alignment when the T2T reference indicated that the VHG1.2 assembly, despite being incomplete, contains sequences not resolved in hg38. Interestingly, the VHG1.1 assembly showed no significant difference in matching with the T2T compared to hg38, suggesting that the extra polishing step not only corrected the assembly to a KHV consensus but also improved its completeness. This improvement was not reflected through the BUSCO score as this method only tested completeness using a set of conserved genes. We hypothesized that these changes in alignment were in genomic regions where the Illumina reads combined with PacBio HiFi reads reached a reliable threshold for JASPER to correct, and these regions were covered in T2T but not hg38. It was also possible that in hg38, these regions contained large-scale polymorphisms from the various contributing donors, whereas both T2T and VHG comprised only one donor’s genome or one fairly homogeneous population’s genome like the Kinh Vietnamese [16,60].

3.4. Toward the First Version of the Kinh Vietnamese Population-Specific Reference Genome

To validate VHG1.2 as a population-specific assembly for the Kinh Vietnamese, we used three additional KHV genomes from 1KGP we selected for variant calling analysis. The variant callers used in this study have been previously evaluated, tested, and benchmarked and are widely adopted in clinical genetics [61,62,63].

We observed a similar trend across the output of different methods, where the number of SNPs and indels called was reduced when a VHG assembly was used as the reference compared to the hg38 reference (Figure 3a,b, Table S11). This result aligned with previous comparisons of the Han Chinese CN1 and HX1 reference genomes with T2T-CHM13 and hg38, respectively, in variant calling for samples of Chinese descent [20,21], as well as variant calling results using the Korean genome reference KOREF [24]. Polishing our assembly with other KHV genomes also affected the number of called variants, potentially reducing false negatives by using a more population-specific reference. In particular, as the reference instead of hg38, the de novo assembly VHG1.1 reduced the amount of SNP and indel calls by 800,000, and the population-specific VHG1.2 assembly reduced this number by a further 100,000 (approximately 3% of the total called SNPs and indels).

Interestingly, there was a significant difference in SNV calling results between GATK and DeepVariant, with GATK detecting up to 30% fewer SNPs and 15% fewer indels than DeepVariant. This discrepancy is likely due to the limitations of GATK, which make it more sensitive to artifacts and data quality issues, particularly sequencing depth and read mapping ambiguity [17,61].

DeepVariant, on the other hand, has been shown to perform better with low-coverage Illumina datasets of 12X or less [7,61], which applied to our test set. Despite these differences, both methods still showed a statistically significant reduction in the number of called SNPs and indels when VHG1.2 was used as a reference instead of hg38 (p < 0.05, see Table S12).

In addition, we utilized the test set for structural variant (SV) detection using Cue, a deep-learning-based variant caller. SV discovery, using only short-read Illumina data, remains highly challenging due to limited long-range contiguity [64]. However, despite being restricted to detecting only three types of SVs (deletions, inversion, and duplication) of 5 + Kbp or larger, Cue has demonstrated state-of-the-art performance among short-read callers, achieving high F1 scores even with low-coverage datasets [55]. By calling SVs in our test set using Cue, we were able to obtain highly confident estimates of the number of SVs per sample, as well as the number of SVs shared between individuals, enabling us to assess the specificity of our assembly in comparison with hg38 as a reference (Figure 3c, Table S13).

For large SVs greater than 2 Kbp, approximately 70% were shared across most super-populations, with 31% identified as unique to East Asian populations [65]. Previous studies also reported thousands of superfluous SVs in worldwide samples that most likely come from the lack of population-specific sequences in the hg19 and hg38 [34,66], as East Asian haplotypes only account for 4.3% of the Genome Reference Consortium’s reference genomes. Furthermore, Takayama et al. suggested that most of the Japanese population shared many of the identified SVs [34]. Based on this, we expected a significantly smaller number of SVs to be detected when using VHG1.2 as a reference compared to hg38. Indeed, Cue identified 3 to 4 times fewer SVs using VHG1.2, compared to hg38, mostly of deletion and inversion types (Figure 3c). The amount of shared SVs between all three individuals with hg38 as reference was 11.5% of the total SVs detected, but between each pair of samples, the amount of shared SVs reached up to 45.5%. In contrast, with VHG1.2 as a reference, only one SV (2.2%) was shared between all samples, and the amount of shared SVs between each pair of samples was approximately 25% (Table S13). Unfortunately, as we currently not yet be able to do this analysis with matching coordination between the two references, we could not ascertain how many of the called SVs were shared between hg38 and VHG1.2. This will be another aspect to look into when VHG1.2 can be completed and annotated by using higher coverage Illumina KHV datasets and Oxford Nanopore Technology (ONT) reads in the future, as it is another indicator of how accurate our assembly was compared to hg38.

3.5. Representativeness of VHG1.2 for the Vietnamese Population

Genome-to-genome alignments revealed that VHG1.2 shares a high degree of similarity with East Asian assemblies, with over 75% identity observed in AK1, HX1, and JG1 as 86.95%, 84.58%, and 82.29%, respectively. The mapping data show the greatest divergence between VHG1.2 and hg38, with only 79.25% of the sequence aligning at >75% identity (Figure 4 and Figure S6, Table S14).

Interestingly, among the East Asian references, JG1 exhibited the most distant alignment with VHG1.2. This is consistent with previous findings that, although JG1 is of East Asian origin, it tends to localize outside the main East Asian cluster in population structure analyses [38]. Although VHG1.2 is still in its early version and has several limitations, these findings collectively support its representativeness for the Vietnamese population and its close genetic relationship with other East Asian populations.

4. Conclusions

VHG1.2 can be considered the highest-quality Kinh Vietnamese genome assembly to date with high completeness, accuracy, and contiguity, comparable to other high-quality population-specific de novo genome assemblies. This marks an important first step toward constructing a Vietnamese reference genome, which will facilitate the development of personalized medicine approaches in Vietnam. However, the VHG1.2 assembly is not entirely complete. This can be explained by the limitation of the sequencing data we produced and the capacity of our computational power.