An Icelandic pangenome reference | Nature

Pangenome reference

The pangenome reference used in this study is a sequence graph augmented with phased small variants. It is stored in two data files: a sequence graph in reference rGFA format¹⁶ containing the structural variants and a supplementary phased VCF³³ file containing small variants. The rGFA sequence graph is constructed first from haplotype assemblies, then the small variants are extracted from the assemblies and added to the VCF file³⁰.

The sequence graph G(V,E) has vertices v ϵ V and edges e ϵ E. A segment is a sub-sequence s of one or more input haplotype assemblies. Each segment is made of two vertices, having respectively the forward and reverse directions. Traversing the forward vertex spells out s while traversing the reverse vertex spells out its reverse complement. An edge e connects two segments if their respective sub-sequences are adjacent in one or more input assemblies. The edges are bidirected such that each end of an edge is either pointing forward or reverse. The direction determines the traversal directions of the connected segments. Thus, there are four forms of edges: forward–forward, forward–reverse, reverse–forward and reverse–reverse.

The rGFA graph is incrementally constructed from a set of input haplotype assemblies (Supplementary Methods). The first input assembly serves as the reference assembly, which is the basis of the coordinate system. When another haplotype assembly is inserted in the graph, existing segments might be split up and additional vertices with new sequences introduced. The sequences spelled out by paths in the graph from a previous insertion can always be spelled out in the updated graph. Since these sequences do not change, they are denoted stable. To keep the sequences stable across insertions, we give each edge and vertex a rank corresponding to the input haplotype assembly that it is originally from. Traversal of a path of rank r spells out a sub-sequence in the r-th haplotype assembly sequence.

Weaver

Weaver maps paired-end short reads in FASTQ format to a pangenome reference, such as HPRC-ICE. The reference can either be linear in FASTA format or a graph in rGFA format with an optional phased VCF file containing small variants. It outputs alignments in Sequence Alignment/Map (SAM) format⁵³, typically compressed as a BAM or CRAM file. Before mapping, Weaver preprocesses the pangenome reference once in an indexing step.

Indexing

The Weaver index has two components. The first component is a seed index that stores the minimizers found in the pangenome reference alongside their locations. A minimizer is a sub-sequence of length k (23 by default), called a k-mer, with the minimum hash value in a sliding window containing w overlapping k-mers (11 by default). We use the minimap2 (ref. ⁵⁴) hash function, which returns the same hash value for a k-mer and its reverse complement. As a result, the same set of minimizers is extracted regardless of the orientation of sequences in the pangenome reference. The minimizers are extracted from all haplotype sequences represented in the graph and in the VCF. They are then added to a seed index along with their graph locations (Fig. 3a). The minimizers may therefore contain bases from alternative alleles in the VCF file. The second component is the ICU index, which stores all pairs of vertices that see each other. We define a vertex u as seeing vertex v if there exists a path from u to v in the graph, which is d bp (d = 1,500 by default) or less.

Mapping

During Weaver read mapping, minimizers are extracted from the read sequence in the same way as during indexing. The minimizers are then anchored to the graph using the seed index (Fig. 3b). Weaver links together seed index hits to form chains (Fig. 3c) when the seeds anchor onto vertices of same rank and are interspersed by the same distances in the read and the graph. If both reads in a pair have long chains that see each other, other small chains that do not see chains on the other read are discarded (Supplementary Methods 2). The chains are then extended along the stable sequence by traversing paths of the same rank in the graph. Then we perform a pairwise alignment between the stable sequence and the read sequence (Fig. 3d). An alignment score S is calculated as

where h is the haplotype sequence, and A = 1, B = 4, C = 6, O = 7 and E = 1 by default. n_match and n_mismatch are the numbers of sequence matches and mismatches in the alignment, respectively. n_clip is the number of soft clips at the beginning or end of the read sequence. Gaps in the alignment are penalized using an affine cost, where n_gapOpen and n_gapExtend are the numbers of gap openings and extensions, respectively. We choose the alignment maximizing this score with h as the stable haplotype sequence and, in case of a tie, we use the one that has gaps as far left as possible.

Variation-aware scoring

After all chains have been aligned, Weaver estimates which of them are most likely at their correct genomic location. Comparing the previously calculated alignment scores would be biased towards the arbitrary selected stable sequence. Instead, Weaver calculates a weighted average alignment score, wS, across all haplotype sequences h ϵ H

where the weight P(h|S_h) is the probability that the read was sequenced from haplotype sequence h. Common variants are observed in many haplotypes and thus impact the alignment score more than rare variants. The haplotypes h ϵ H include the stable sequence and all the sequences represented in the VCF file. Weaver selects the alignment with the maximum wS as the primary alignment.

Emblask

Emblask is a diploid genome assembly pipeline that produces a set of two haplotype assemblies, also called dual assembly, from the offspring of a parent–offspring trio (Supplementary Methods 3). The method is a hybrid approach using both noisy long reads and accurate short reads, specifically ONT R9.4 and paired-end Illumina reads in this study. Emblask takes as input long reads for the genome of the offspring as well as short reads for the three members of the trio.

In the following, we refer to long reads as LRs and short reads as SRs. We define cov(A) and cov(A,s) as the LR coverage of assembly A and the coverage of sequence s ∈ A, respectively. Phasing refers to mapping LRs to a sequence or a set of sequences, calling small variants with PMDV²⁷ and assigning the alleles of the called heterozygous variants to a haplotype H₁ or H₂ using the LR overlap between adjacent variants. A phase set delineates a region in which two or more heterozygous variants are phased. In a phase set, the haplotype for which the phased variants have the most reference alleles is the reference haplotype H_r while the other haplotype is the alternate haplotype denoted H_a. Haplotagging refers to assigning a haplotype tag H₁ or H₂ to LR alignments in a phase set based on the phased variants they overlap.

Global error correction

The first step in Emblask is to decrease the LR error rate using the SRs with Ratatosk²⁵. After correcting the LRs, Emblask trims sub-sequences with low correction scores reflecting uncorrected or low-quality corrected bases.

Collapsed assembly

Emblask assembles the corrected LRs into a collapsed assembly A_c with Flye²⁶. Each assembled sequence, called a contig, contains a combination of alleles from the paternal and maternal haplotypes. Each locus is therefore represented in at most one contig.

Local error correction

Emblask improves the corrected LR error rate by mapping SRs and corrected LRs to the contigs of A_c to perform local corrections. Non-overlapping segments of 200 kb are defined on A_c to split the read alignments into different windows that are corrected separately. A paired-end SR can occur in multiple windows if it cannot be mapped uniquely on A_c.

Haplotype-resolved assembly

Emblask assembles the corrected LRs into a haplotype-resolved assembly A_r with Flye. The output haplotigs represent sub-sequences of the paternal or maternal haplotype but without distinction to which of the two haplotypes each sequence is from. Furthermore, the assembly is still fragmented and incomplete because for any locus, only one of the two parental haplotypes might have been assembled into a haplotig.

Haplotig cleaning

Each haplotig must be assigned to either the paternal or maternal haplotype to create the final assembly. Errors in the haplotigs can lead to an incorrect parental assignment that would result in fragmentation, gaps and false duplications in the final assembly. Therefore, haplotigs must be filtered, split and polished before performing the parental assignment. The coverage of each haplotype assembly is initially expected to be half the coverage of the haplotigs in A_c:

Corrected LRs are mapped to A_r and only haplotigs h′ with coverage within the range \(\frac{c}{2} < \mathrm{cov}({A}_{r},{h}^{{\prime} }) < 3c\) are kept to eliminate haplotigs that are the result of erroneous duplications or collapsing during assembly. Corrected LRs mapping to the remaining haplotigs are then phased and haplotagged with Margin⁵⁵. Haplotigs with coverage greater than 50% of c are annotated as collapsed coverage (CC) candidates. CC candidates containing multiple phase sets are split between phase sets to ensure haplotype phasing integrity. CC candidates are also split by removing the sub-sequences of phase sets for which at least 25% of the heterozygous SNPs have a phase inconsistent with the reference and alternate haplotypes. The resulting haplotigs are then polished with Flye using only the untagged and reference-tagged alignments. Additional fine-grained haplotig cleaning takes place by refining the expected haplotype coverage with the coverage of high-quality phased SNPs in the polished haplotigs. The final set of cleaned and polished haplotigs is denoted A_r′.

Haplotig trio binning

A set of haplotigs B_r′ that closely approximates the missing haplotigs of A_r′ is produced by haplotagging LRs with respect to A_r′ and polishing haplotigs of A_r′ with the alternate-tagged and untagged alignments. Haplotigs of A_r′ and B_r′ are then assigned to the paternal or maternal haplotype with which they share the most sub-sequences. For any trio-binned haplotig h_a in A_r′, if there exists an approximated alternate haplotig h_b in B_r′ assigned to the same parent as h_a, the parental assignment with the lower confidence is flipped (Supplementary Methods 3).

LR trio binning

Haplotigs in A_r′ have been assigned to a parental haplotype but A_r′ is still incomplete and fragmented. To resolve both issues, LRs are mapped to A_r′, phased and haplotagged. For each haplotig assigned to a parental haplotype H, LRs from the reference-tagged primary alignments are assigned to H while LRs from alternate-tagged primary alignments are assigned to the other parental haplotype. LRs from untagged primary alignments outside of phase sets are assigned to H or evenly distributed between the two parental haplotypes if the local coverage indicates the presence of two haplotypes (Supplementary Methods 3).

Dual assembly

Each group of LRs assigned to either parental haplotype is assembled independently with Flye, resulting in two preliminary haplotype assemblies, which are filtered further with the aim of removing erroneous duplications caused by incorrect parental assignment. Each haplotype assembly is then polished with the LRs and the SRs.

HPRC-ICE

The HPRC-ICE pangenome reference is first composed of the HPRC year 1 pangenome graph⁴ built with Minigraph¹⁶ using GRCh38 as the reference assembly. Non-reference vertices present in fewer than nine assemblies were removed from the graph. The HPRC-ICE is also composed of a VCF file representing phased small variants called from the ICE, HPRC Y1 and T2T-CHM13 assemblies. In the following, we refer to the Icelandic assemblies created from PacBio HiFi reads as ICE-PB and the assemblies created from ONT-Illumina reads as ICE-ONT/ILMN (Supplementary Methods 5).

Genome assembly

All PacBio HiFi samples were assembled using hifiasm or hifiasm-trio if the parental short reads were available. All ONT R9.4 samples were assembled with trio short reads using the Emblask pipeline (Supplementary Methods 3, 5 and 6). ONT datasets were automatically downsampled to 50× by Emblask prior to each assembly step in the pipeline to keep running time and memory usage tractable.

Variant calling

Phased small variants were called from all the ICE dual assemblies with a modified version of dipcall⁵⁶ using minimap2 (ref. ⁵⁴) v2.24 and wider z-drop score parameters to improve the contiguity of the assembly alignments⁵⁷. The output variant calls were left-aligned and normalized, and multi-allelic variants were split into bi-allelic. Furthermore, all structural variants, variants with a star allele and variants with a missing genotype or a half-missing half-reference were filtered out. Variant calls from the ICE-ONT/ILMN assemblies were then filtered and polished using the trio short reads (Supplementary Methods 8). Finally, each diploid genotype was split into two haploid genotypes, one for each haplotype assembly.

Quality control

Small variants from the ICE haplotype assemblies fulfilling all the following criteria were merged into the HPRC-ICE set: minimum 95% k-mer completeness, 90% BUSCO single-copy completeness, QV45, 0.8 Mb haplotig N50, in addition to maximum 2% phase switch error rate and 2.5% false duplication rate. In the ICE-ONT/ILMN set, 636 haplotype assemblies (96.65%) passed all the quality control criteria and 22 assemblies (3.35%) failed at least one quality control requirement. In the ICE-PB set, 44 assemblies passed all quality control requirements and 54 passed all quality control requirements except the phase switch error rate, which cannot be computed without parental Illumina reads. Among the 17 individuals with retained ICE-PB and ICE-ONT/ILMN assemblies, variant calls from the ICE-PB assemblies of 16 individuals were merged and the remaining individual was set aside for validation purposes.

Small variants merging

We used the HPRC Y1 Minigraph-Cactus v1.1 VCF file, which contains variants converted to VCF format from the sequence graph produced by the Minigraph-Cactus pipeline for the haploid T2T-CHM13 assembly and the HPRC dual assemblies with many variants nested in a snarl. To obtain non-overlapping sites, bubbles were popped with vcfbub⁵⁸ and only sites with alleles shorter than 100 kb were initially kept. All structural variants and variants with a star allele were then removed, followed by a left-alignment and normalization of all remaining variants. Diploid genotypes were then split into haploid genotypes. The HPRC haploid variant calls were merged with their ICE haploid counterparts using bcftools⁵⁹. The resulting multi-sample VCF contains 787 haplotype samples: 602 samples from the ICE-ONT/ILMN haplotype assemblies, 96 samples from the ICE-PB haplotype assemblies, 88 samples from the HPRC haplotype assemblies (HG002, HG005 and NA19240 were set aside for internal validation) and the haploid T2T-CHM13 assembly. Complex variants were decomposed into simpler primitives with vcfwave⁶⁰ and duplicated primitives were merged.

Icelandic DNA data

Whole-genome sequencing of 57,630 Icelanders followed standard Illumina TruSeq PCR-Free methodology using HiSeqX, NovaSeq and NovaSeqX machines. All the samples were sequenced with minimum 20× coverage. Illumina SNP chip-typing was performed on 182,338 Icelanders for long-range phasing⁶¹ and imputation, as described previously⁶². Among the Illumina-sequenced Icelanders, 312 were sequenced with ONT R9.4 long reads and 49 were sequenced with PacBio HiFi (Supplementary Methods 5).

All participating subjects signed informed consent. The personal identities of the participants and biological samples were encrypted by a third-party system approved and monitored by the Data Protection Authority. The National Bioethics Committee and the Data Protection Authority in Iceland approved these studies.

Statistical analyses

We used logistic regression with an additive model to test for the association between sequence variants and binary traits. The reported P values are based on two-sided tests with age and sex as covariates. No statistical methods were used to predetermine sample size for association testing. Reported correlations are Pearson’s correlation coefficients (r).

In Figs. 2 and 4, the box plots show the distributions of data points with the interquartile range (IQR) from the 25th to the 75th percentiles represented as a box, the median value represented as a line within the box, the whiskers as lines extending the box to the minimum and maximum values within 1.5 times of the IQR and the outlier values as data points extending beyond the whiskers.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.