Article A draft human pangenome reference ­­ https://doi.org/10.1038/s41586-023-05896-x Wen-Wei Liao1,2,3,60, Mobin Asri4,60, Jana Ebler5,6,60, Daniel Doerr5,6, Marina Haukness4, Glenn Hickey4, Shuangjia Lu1,2, Julian K. Lucas4, Jean Monlong4, Haley J. Abel7, Received: 9 July 2022 Silvia Buonaiuto8, Xian H. Chang4, Haoyu Cheng9,10, Justin Chu9, Vincenza Colonna8,11, Accepted: 28 February 2023 Jordan M. Eizenga4, Xiaowen Feng9,10, Christian Fischer11, Robert S. Fulton12,13, Shilpa Garg14, Cristian Groza15, Andrea Guarracino11,16, William T. Harvey17, Simon Heumos18,19, Published online: 10 May 2023 Kerstin Howe20, Miten Jain21, Tsung-Yu Lu22, Charles Markello4, Fergal J. Martin23, Open access Matthew W. Mitchell24, Katherine M. Munson17, Moses Njagi Mwaniki25, Adam M. Novak4, Hugh E. Olsen4, Trevor Pesout4, David Porubsky17, Pjotr Prins11, Jonas A. Sibbesen26, Check for updates Jouni Sirén4, Chad Tomlinson12, Flavia Villani11, Mitchell R. Vollger17,27, Lucinda L. Antonacci-Fulton12, Gunjan Baid28, Carl A. Baker17, Anastasiya Belyaeva28, Konstantinos Billis23, Andrew Carroll28, Pi-Chuan Chang28, Sarah Cody12, Daniel E. Cook28, Robert M. Cook-Deegan29, Omar E. Cornejo30, Mark Diekhans4, Peter Ebert5,6,31, Susan Fairley23, Olivier Fedrigo32, Adam L. Felsenfeld33, Giulio Formenti32, Adam Frankish23, Yan Gao34, Nanibaa’ A. Garrison35,36,37, Carlos Garcia Giron23, Richard E. Green38,39, Leanne Haggerty23, Kendra Hoekzema17, Thibaut Hourlier23, Hanlee P. Ji40, Eimear E. Kenny41, Barbara A. Koenig42, Alexey Kolesnikov28, Jan O. Korbel23,43, Jennifer Kordosky17, Sergey Koren44, HoJoon Lee40, Alexandra P. Lewis17, Hugo Magalhães5,6, Santiago Marco-Sola45,46, Pierre Marijon5,6, Ann McCartney44, Jennifer McDaniel47, Jacquelyn Mountcastle32, Maria Nattestad28, Sergey Nurk44, Nathan D. Olson47, Alice B. Popejoy48, Daniela Puiu49, Mikko Rautiainen44, Allison A. Regier12, Arang Rhie44, Samuel Sacco30, Ashley D. Sanders50, Valerie A. Schneider51, Baergen I. Schultz33, Kishwar Shafin28, Michael W. Smith33, Heidi J. Sofia33, Ahmad N. Abou Tayoun52,53, Françoise Thibaud-Nissen51, Francesca Floriana Tricomi23, Justin Wagner47, Brian Walenz44, Jonathan M. D. Wood20, Aleksey V. Zimin49,54, Guillaume Bourque55,56,57, Mark J. P. Chaisson22, Paul Flicek23, Adam M. Phillippy44, Justin M. Zook47, Evan E. Eichler17,58, David Haussler4,58, Ting Wang12,13, Erich D. Jarvis32,58,59, Karen H. Miga4, Erik Garrison11 ✉, Tobias Marschall5,6 ✉, Ira M. Hall1,2 ✉, Heng Li9,10 ✉ & Benedict Paten4 ✉ Here the Human Pangenome Reference Consortium presents a first draft of the human pangenome reference. The pangenome contains 47 phased, diploid assemblies from a cohort of genetically diverse individuals1. These assemblies cover more than 99% of the expected sequence in each genome and are more than 99% accurate at the structural and base pair levels. Based on alignments of the assemblies, we generate a draft pangenome that captures known variants and haplotypes and reveals new alleles at structurally complex loci. We also add 119 million base pairs of euchromatic polymorphic sequences and 1,115 gene duplications relative to the existing reference GRCh38. Roughly 90 million of the additional base pairs are derived from structural variation. Using our draft pangenome to analyse short-read data reduced small variant discovery errors by 34% and increased the number of structural variants detected per haplotype by 104% compared with GRCh38-based workflows, which enabled the typing of the vast majority of structural variant alleles per sample. The human reference genome has formed the backbone of human the Telomere-to-Telomere (T2T) consortium finished the first complete genomics since its initial draft release more than 20 years ago2. The pri- sequence of a haploid human genome, T2T-CHM13, which provides a mary sequences are a mosaic representation of individual haplotypes contiguous representation of each autosome and of chromosome X, containing one representative scaffold sequence for each chromo- with the exception of some ribosomal DNA arrays that remain to be some. There are 210 Mb of gap or unknown (151 Mb) or computationally fully resolved3. Using T2T-CHM13 directly improves genomic analy- simulated sequences (59 Mb) within the current GRCh38 release, con- ses; for example, discovering 3.7 million additional single-nucleotide stituting 6.7% of the primary chromosome scaffolds. Missing reference polymorphisms (SNPs) in regions non-syntenic to GRCh38 and better sequences create an observational bias, or streetlamp effect, which representing the true copy number variants (CNVs) of samples from limits studies to be within the boundaries of the reference. Recently, the 1000 Genomes Project (1KG) compared with GRCh38 (refs. 1,4). A list of affiliations appears at the end of the paper. 312 | Nature | Vol 617 | 11 May 2023 Although T2T-CHM13 represents a major achievement, no single one of these joins (in the paternal assembly of HG02080) was located genome can represent the genetic diversity of our species. Previous in a euchromatic, non-acrocentric region and was manually confirmed studies have identified tens of megabases of sequence contained within to be a misassembly. The remaining joins involved the short arms of structural variants (SVs) that are polymorphic within the population5. the acrocentric chromosomes (Fig. 1b and Supplementary Table 2). Owing to the absence of these alternative alleles from the reference This may be the result of misalignment, nonallelic gene conversion genome, more than two-thirds of SVs have been missed in studies that or other biological mechanisms that maintain large-scale homology used short-read data and the human reference assembly6–8, despite between the short arms of the acrocentrics—a phenomenon that we individual SVs being more likely to affect gene function than either have studied in an associated paper22. individual SNPs or short insertions and deletions (indels)9,10. To evaluate the resulting assemblies after manual correction of errors, To overcome reference bias, a transition to a pangenomic reference we developed an automated assembly quality control pipeline that has been envisioned11,12. Pangenomic methods have rapidly progressed combined methods to assess the completeness, contiguity, base level over the past few years13–15 such that it is now practical to propose that quality and phasing accuracy of each assembly (Supplementary Table 3 common genomic analyses use a pangenome. Here we sequence and and Methods). Haploid assemblies containing an X chromosome aver- assemble a set of diverse individual genomes and present a draft human aged a total length of 3.04 Gb, 99.3% the length of T2T-CHM13 (3.06 Gb), pangenome, the first release from the Human Pangenome Reference which also contains an X chromosome. Haploid assemblies containing Consortium (HPRC)15. These genomes represent a subset of the planned a Y chromosome averaged a total length of 2.93 Gb, which reflects the HPRC panel, which aims to better capture global genomic diversity size difference between the sex chromosomes (Fig. 1c). The average across the 700 haplotypes of 350 individuals. NG50 value, a widely used measure of contiguity, was comparable with the contig NG50 of GRCh38 (40 Mb compared with 56 Mb, respectively; Fig. 1d). Using short substrings (k-mer values of 31) derived from Illu- Assembling 47 diverse human genomes mina data, the Yak k-mer analyser20 estimated an average quality value We assembled 47 fully phased diploid assemblies from genomes selected (QV) of 53.57 for the assemblies, which corresponded to an average to represent global genetic diversity (Fig. 1a) and for which consent of 1 base error per 227,509 bases (Fig. 1e). To validate these QV esti- had been given for unrestricted access. All assemblies have been made mates, we benchmarked the HG002 and HG005 assembly-based variant publicly available, along with all data and analyses. These assemblies calls against the small variants called using Genome in a Bottle (GIAB; include 29 samples with long and linked read sequencing data generated v.4.2.1). We estimated QVs of 54 for HG002 and 55 for HG005, which were entirely by the HPRC and 18 samples sequenced by other efforts16–18. In highly similar to the k-mer QVs estimated using Yak. Consistent with some cases, we supplemented the 18 additional samples with further our manual observation that most errors were primarily small indels in sequencing. We selected the 29 HPRC samples from the 1KG lymphoblas- low-complexity regions, we found that approximately 32% of the indel toid cell lines, limiting selection to those lines classified as karyotypically errors were in homopolymers longer than 5 bp and an additional 48% normal and with low passage (to avoid artefacts from cell culture). We were in tandem repeats and low-complexity regions. Moreover, about also ensured that the cell lines were derived from participants for whom 42% of the indel errors were genotype errors, mostly heterozygous whole-genome sequencing (WGS) data were available for both parents variants incorrectly called as homozygous variants due to collapsed (for haplotype phasing). Cell lines meeting these criteria were prioritized haplotypes in the two assemblies of an individual (Supplementary by genetic and biogeographic diversity (Methods). Table 4). We next used Yak to analyse the phasing accuracy between We created a consistent set of deeply sequenced data types for every the maternal and paternal assemblies using k-mer values derived from sample (Supplementary Table 1). The data included Pacific Biosciences Illumina sequencing of the parents. An average haplotype switch error (PacBio) high-fidelity (HiFi) and Oxford Nanopore Technologies (ONT) rate of 0.67% was observed and a Hamming error rate of 0.79% (Fig. 1f). long-read sequencing, Bionano optical maps and high-coverage Hi-C We also calculated phase accuracy using Pstools23,24, which uses Hi-C Illumina short-read sequencing for all HPRC samples. We also gath- sequence data of the sample not used to create the assembly. Pstools ered previously generated high-coverage Illumina sequencing data reported slightly lower switch error rates than Yak but comparable for both parents of each participant19. We generated on average 39.7× Hamming error rates (Supplementary Fig. 2). Taken together, the above HiFi sequence depth of coverage for the 46 HPRC samples (excluding results indicate that the assemblies are highly contiguous and accurate. HG002, which had around 130× coverage). This depth of coverage is consistent with the requirements for high-quality, state-of-the-art Regional assembly reliability assemblies20 and facilitates comprehensive variant discovery irrespec- To determine which portions of the assemblies are reliable, we devel- tive of allele frequency (AF). The N50 value, which represents the short- oped a read-based pipeline, Flagger, that detects different types of est read length at which 50% of the total sequenced bases are covered misassemblies within a phased diploid assembly (Fig. 1g and Methods). by considering only equal or longer reads, was 19.6 kb on average for The pipeline works by mapping the HiFi reads to the combined maternal the HiFi reads (Supplementary Table 1; excluding HG002 because it was and paternal assembly in a haplotype-aware manner. It then identifies sequenced using a different library preparation protocol). coverage inconsistencies within these read mappings that are likely to For the core assembler, we chose Trio-Hifiasm20 after detailed bench- be due to assembly errors. This process is similar to likelihood-based marking of several alternatives21. Trio-Hifiasm uses PacBio HiFi long-read approaches, which assess the assembly given the reads25, but is adapted sequences and parental Illumina short-read sequences to produce near to work with long reads and diploid assemblies. Using Flagger, we iden- fully phased contig assemblies. The complete assembly pipeline (Sup- tified only 0.88% (26.4 Mb) of each assembly as unreliable (Fig. 1h and plementary Fig. 1 and Methods) included steps to remove adaptor and Supplementary Table 5). Using T2T-CHM13, we estimated that 0.09% of nonhuman sequence contamination and to ensure a single mitochon- reliably assembled blocks were falsely labelled as unreliable (Methods). drial assembly per maternal assembly. Compared with the distribution of contig sizes, the unreliable blocks were short (54.6 kb N50 average). We intersected the unreliable blocks Assembly assessment in the assemblies from Flagger with different repeat annotations (Fig. 1i We first searched for large-scale misassemblies, looking for gene duplica- and Supplementary Table 6). We estimated that the following percent- tion errors, phasing errors and interchromosomal misjoins (Methods). age of elements were correctly assembled: 95.4% of alpha satellites; We manually fixed three large duplication errors and one large phasing 91.5% of human satellites 2 and 3; 97.7% of segmental duplications (SDs); error, but left smaller errors, which are difficult to definitively distin- 94.3% of variable number tandem repeats (VNTRs); 94.2% of short guish from SVs. We found 217 putative interchromosomal joins. Only tandem repeats (STRs); and 98.8% of all human repeats26. Nature | Vol 617 | 11 May 2023 | 313 Article b c a Chr. 22 Africa 51% Chr. 13 3.06 15.711 Mb 17.509 Mb 0.10 3.04 Total lengh (×109) 3.02 Principal component 2 Americas 34% 0.05 3.00 Asia 13% Chinese-GIAB (1) 2.98 Europe 2% Ashkenazim Jewish-GIAB (1) 2.96 0 2.94 GWD (8) ASW (1) ACB (7) PJL (1) 2.92 CHS (3) Chr. 21 –0.05 PUR (8) ESN (1) 11.306 Mb Chr. 14 2.90 CLM (4) KHV (1) Maternal- Maternal- Paternal- Paternal- MKK-HapMap (1) 12.708 Mb female male female male PEL (4) –0.10 MSL (4) Chr. 15 Haplotype YRI (2) 17.694 Mb –0.10 –0.05 0 0.05 0.10 Principal component 1 Coverage distribution d e f g 10 Erroneous Frequency (Mb) CHM13 Duplicated 250 Average NG50: 40 Mb 8 Haploid Fit a Gaussian mixture model Collapsed 6 2.5 Maternal Actual 200 58 Paternal 4 Hamming error rate (%) Contig length (Mb) 2.0 2 150 GRCh38 (N-Split) 56 Paternal QV 0 1.5 0 20 40 60 80 100 54 Coverage 100 1.0 Erroneous False duplication 52 Haploid Collapsed 50 0.5 50 0 0 48 0 0.5 1.0 1.5 2.0 2.5 3.0 48 50 52 54 56 58 60 0 50 00 50 25 75 25 75 Window 0. 1. 1. 0. 0. 1. 1. Cumulative coverage (Gb) Maternal QV Switch error rate (%) Find reliable or unreliable blocks h 3,300 3,200 3,100 3,000 2,900 Total length (Mb) 70 60 50 Haploid 40 Unknown 30 Collapsed 20 Duplicated 10 Erroneous 0 CH ge 0 71 1 0 91 02 41 0 21 0 29 0 09 0 59 00 09 0 79 0 9 23 0 3 0 23 0 53 01 43 02 73 13 0 22 0 30 0 0 40 03 52 0 92 HG 1 2 2 8 0 48 0 40 0 8 8 H 098 0 58 0 58 0 28 0 86 0 86 0 06 0 6 0 86 6 02 75 0 35 0 55 0 45 0 17 0 57 0 05 NA 136 0 3 NA 208 00 57 NA 81 3 NA 197 1 10 a HG 18 HG 07 HG 10 HG 06 HG201 HG213 HG 25 HG 21 HG 35 HG 11 HG 07 HG 27 HG 34 HG 12 HG 06 M HG 26 HG 26 HG 35 HG 19 HG 34 HG 1 HG192 HG 04 HG 13 HG 12 HG 19 HG 34 HG 28 HG189 HG 35 HG 24 HG011 HG 7 HG 20 HG 21 HG 27 HG 22 HG G0 er 1 HG G 0 0 Av H HG Female Male i j 100 2,900 95 2,850 90 2,800 85 2,750 Total length (Mb) Reliability (%) 80 2,700 75 175 70 150 65 125 100 60 Only 1 alignment 75 Female 50 Not covered 55 Male 2 Alignments 25 0 >2 Alignments 50 Paternal Maternal Paternal Maternal al l al l al l al l al l al l na na na na na na rn rn rn rn rn rn er er er er er er te te te te te te Female Male at at at at at at Pa Pa Pa Pa Pa Pa M M M M M M AlphaSat HSat2/3 SD RepeatMasker VNTR STR Repeat annotation Fig. 1 | Presenting 47 accurate and near-complete diverse diploid human accuracy showing the switch error percentage versus Hamming error genome assemblies. a, Selecting the HPRC samples. Left, the first two percentage. g, Flagger read-based assembly evaluation pipeline. Coverage is principal components of 1KG samples showing HPRC (triangles) samples, calculated across the genome and a mixture model is fit to account for reliably excluding HG002, HG005 and NA21309. Right, summary of the HPRC sample assembled haploid sequence and various classes of unreliably assembled subpopulations (three letter abbreviations) on a map of Earth as defined by the sequence. For each coverage block, a label is assigned according to the most 1KG. ACB, African Caribbean in Barbados; ASW, African Ancestry in Southwest probable mixture component to which it belongs: erroneous, falsely duplicated, US; CHS, Han Chinese South; CLM, Colombian in Medellin, Colombia; ESN, Esan (reliable) haploid, collapsed, and unknown. h, Reliability of the 47 HPRC in Nigeria; GWD, Gambian in Western Division; KHV, Kinh in Ho Chi Minh City, assemblies using read mapping. For each sample, the left bar is the paternal Vietnam; MKK, Maasai in Kinyawa, Kenya; MSL, Mende in Sierra Leone; PEL, and the right bar is the maternal haplotype. Regions flagged as haploid are Peruvian in Lima, Peru; PJL, Punjabi in Lahore, Pakistan; PUR, Puerto Rican in reliable (green), constituting more than 99% on average of each assembly. Puerto Rico; YRI, Yoruba in Ibadan, Nigeria. b, Interchromosomal joins between The y axis is broken to show the dominance of the reliable haploid component acrocentric chromosome short arms. Red, the join is on the same strand; and the stratification of the unreliable blocks. i, Assembly reliability of six blue, otherwise. c, Total assembled sequence per haploid phased assembly. types of repeats. AlphaSat, alpha satellites; HSat2/3, human satellites 2 and 3. d, Assembly contiguity shown as a NGx plot. T2T-CHM13 and GRCh38 contigs j, Completeness of the HPRC assemblies relative to T2T-CHM13. The number of are included for comparison. e, Assembly QVs showing the base-level accuracy reference bases covered by none, by one, by two or by more than two alignments of the maternal and paternal assembly for each sample. f, Yak-reported phasing are included. (excluding chromosome X) on average with exactly one alignment. Completeness and CNV For all other assemblies (excluding chromosome Y), about 94.1% on To assess the completeness and copy number polymorphism of the average was single-copy covered (Fig. 1j and Supplementary Table 7). assemblies, we aligned them to T2T-CHM13 (Methods). The pater- On average, around 136 Mb (4.4%) of T2T-CHM13 was not covered by nal assemblies of male samples covered about 92.8% of T2T-CHM13 any alignment, which indicates that some parts of the genome are either 314 | Nature | Vol 617 | 11 May 2023 a b c 3,000 99.5 2,500 AMR AFR EAS Percentage annotated 99.0 2,000 SAS CHM13 Frequency Duplicated genes per genome 60 98.5 1,500 1,000 40 98.0 500 20 97.5 0 0 Protein- Noncoding Protein- Noncoding 0 5 10 15 20 Genome by number of duplications coding genes coding transcripts Gene copy number genes mapped transcripts mapped mapped mapped d e f g 25 90 25 Number of samples with additional copies 350 80 00 Chr. 6: 101295790−103818969 00 Additional copy number count Additional copy number count Chr. 10: 45541736−48056658 300 00 20 70 20 48 0 Total additional copies 00 00 60 250 SPDYE2B 30 10 15 15 SPDYE2 50 00 200 0 00 00 00 00 40 GPRIN2 25 47 10 10 150 10 30 SPDYE2 0 00 100 GPRIN2 GPRIN2 SPDYE2B 00 20 GPRIN2 SPDYE2 00 5 5 20 GPRIN2 SPDYE2 00 10 50 00 10 46 00 56 000 00 00 00 00 00 00 00 00 00 0 0 0 0 00 00 00 00 00 00 00 00 00 15 CLPSL1 OR9G1 CCL3L1 FAM90A1 GAGE1 CCL3L1 CLPSL1 OR9G1 DUSP22 ARL17A ARL17A TSPY4 BOLA2B BOLA2B SULT1A3 MAP2K3 FAM111B LINC02203 TSPY3 SULT1A3 LINC02203 MAP2K3 FAM111B PRAMEF18 NBPF9 PRAMEF18 ZNF595 GPRIN2 DUSP22 GPAT2 SPDYE2 CLPS NBPF9 RP11-577H5.5 ZNF595 SPDYE2 GPRIN2 GPAT2 CLPS RP11-577H5.5 CCDC200 CCDC200 KMT2C KMT2C PDPR PDPR RP11-514P8.6 RP11-514P8.6 TSPY9P RP11-514P8.7 RP11-514P8.7 GAGE2E TPTE GAGE2E TPTE 00 50 00 50 00 00 00 00 00 10 56 57 57 58 49 50 51 52 HG01361-2 HG00621-1 JAGYYW010000028.1:48955433−52399444 JAHBCD010000047.1:56010204−58018759 Gene family Gene family Fig. 2 | Transcriptome annotation of the assemblies. a, Ensembl mapping families in the HPRC assemblies out of all 1,115 duplicated genes, ordered by pipeline results. Percentages of protein-coding and noncoding genes and the number of samples with additional copies relative to GRCh38. Grey bars, transcripts annotated from the reference set in each of the HPRC assemblies. the number of samples with additional copies. Blue circles, the number of Orange points represent T2T-CHM13 for comparison. b, Frequency of gene additional copies per sample, with the size of the circle proportional to the copy number. Individual genes may have separate copy number states among number of samples. e, The top 30 most individually copied CNV genes or gene genomes, and the frequency reflects 3,210 observed copy number changes families in the HPRC assemblies, ordered by total number of additional copies among the HPRC genomes. c, Number of distinct duplicated genes or gene observed. Blue circles, the number of additional copies per sample. Grey bars, families per phased assembly relative to the number of duplicated genes the total number of additional copies summed over the samples. f, Dotplot annotated in GRCh38 (n = 152). The GRCh38 gene duplications reflect families illustrating haplotype-resolved GPRIN2 gains in the HG01361 assembly relative of duplicated genes, whereas the counts in other genomes reflect gene to GRCh38. g, Dotplot illustrating SPDYE2–SPDYE2B haplotype resolved gains duplication polymorphisms. The assemblies are colour coded according to within a tandem duplication cluster of the HG00621 assembly relative to their population of origin. d, The top 25 most commonly CNV genes or gene GRCh38. systematically unassembled or cannot be reliably aligned. About 90% assembled SDs and 1.8% more identical (98.9 compared with 97.1). This of these regions were centromeric or pericentromeric27 (Extended Data result reinforces earlier findings that the length and identity of SDs Fig. 1). Despite the majority of unaligned bases occurring within and play an important part in assembly accuracy28. around centromeres, on average, 90% of divergent and monomeric alpha satellites, gamma satellites and centromeric transition regions were covered by at least one alignment. Excluding the T2T-CHM13 Annotating 47 diverse genomes centromere and satellites3 and including only the expected sex chro- We developed a new Ensembl mapping pipeline to annotate GENCODE29 mosome for each haploid assembly, on average, around 99.12% of the genes and transcripts within each new haploid assembly (Methods). remaining reference was covered by exactly one alignment (Supple- A median of 99.07% of protein-coding genes (range = 98.08–99.40%) mentary Table 7). and 99.42% of protein-coding transcripts (range = 98.29–99.66%) The average number of T2T-CHM13 bases with two or with more than were identified in each of the HPRC assemblies (Fig. 2a and Sup- two alignments was about 32.4 Mb (around 1.0%) and about 20.0 Mb plementary Table 8). Similarly, a median of 98.16% of noncoding (around 0.6%), respectively. On average, per haploid assembly, these genes (range = 97.23–98.60%) and 98.96% of noncoding transcripts duplicated regions had about 82.20% and 39.82% overlap with the (range = 97.94–99.28%) were similarly annotated. Running this pipeline pericentomeric or centromeric satellites and SDs, respectively, and on T2T-CHM13 produced similar, slightly higher, results. Intersect- around 94.62% had overlap with either of them. We characterized the ing the HPRC annotations with the assembly reliability predictions, accuracy of regions aligned to SDs in T2T-CHM13 (excluding chromo- a median of 99.53% of gene and 99.79% of transcript annotations some Y) using a liftover of the assembly read-depth-based evaluation occurred wholly within reliable regions, which indicated that most (Extended Data Fig. 2). On average, we estimated that only 2.5% (4.99 of the annotated transcript haplotypes were structurally correct. To out of 199 Mb) of the SD sequence that could be lifted onto T2T-CHM13 examine transcriptome base accuracy, we looked for nonsense and was in error according to the read depth. To identify SDs associated with frameshift mutations in the set of canonical transcripts (one representa- these errors, we took all 5 kb windows across the unreliable regions tive transcript per gene; Supplementary Fig. 3, Supplementary Table 8 and intersected them with the longest and most identical overlapping and Methods). We found a median of 25 nonsense mutations and 72 SD. The median length of SDs overlapping sequences in error was 3.0 frameshifts per assembly. A median of 21 (84%) and 58 (80%) of these times longer (288 kb compared with 96.3 kb) than those in correctly nonsense mutations and frameshifts per assembly, respectively, were Nature | Vol 617 | 11 May 2023 | 315 Article supported by the independently generated Illumina variant call sets. and inversions. MC extends the Minigraph pangenome with a base These numbers were within the range of previously reported numbers level alignment of the homology relationships between the assemblies of loss-of-function mutations (between 10 and 150 per person, depend- using the Cactus genome aligner37 while retaining the structure of ing on the level of conservation of the mutation)1,30. Conservatively, if all the Minigraph pangenome. PGGB constructs a pangenome from an the non-confirmed frameshifts and nonsense mutations are assembly all-to-all alignment of the assemblies. Although both T2T-CHM13 and errors, this would predict 18 such transcript-altering errors per tran- GRCh38 are used to partition contigs into chromosomes, the PGGB scriptome (1 per 1.7 million assembled transcriptome bases). graph is otherwise reference free (that is, it does not base itself on a There were 1,115 protein-coding gene families within the Flagger- chosen reference assembly). predicted reliable regions of the full set of assemblies that had a gain in copy number in at least one genome (Fig. 2b). Each assembly had Measuring pangenome variation an average of 36 genes with a gain in copy number relative to GRCh38 The different algorithmic approaches used to construct a pangenome within its predicted reliable regions, with a bias towards rare, low-copy graph influence graph properties while representing the same under- CNVs (Fig. 2c). In detail, 71% of CNV genes appeared in a single haplo- lying sequences. The basic properties of the three graphs produced type. Previous studies using read depth found that rare CNVs generally with the different pangenome methods are shown in Supplementary occur outside regions annotated as being enriched in SDs31. The genome Table 10. The Minigraph chart, by virtue of being limited to structural assemblies confirmed this observation in sequence-resolved CNVs. variation, is the smallest, with more than two orders of magnitude fewer When stratifying duplicated genes on the basis of AF into singleton nodes and edges than the base level graphs. Its length (3.24 Gb), meas- (present in one haplotype), low frequency (<10%) and high frequency, ured as the total bases of all nodes, is similar to the MC graph (3.29 Gb) 15% (118 out of 771) of the singleton CNVs mapped to SDs as annotated despite the latter adding many small variants. This difference is due in GRCh38. Duplicated genes with a higher population frequency had to the MC graph also aligning a significant number of sequences left a greater fraction in SDs: 59% (83 out of 140) of low frequency and 81% unaligned by Minigraph. The PGGB graph contains roughly 5 Gb more (44 out of 54) of high frequency. Overall, 58 genes were CNVs in 10% or sequence because it includes highly structurally divergent satellite more of haploid assemblies, and 16 genes were amplified in the majority regions omitted from the other approaches (Methods) and does not of individuals relative to GRCh38 (Fig. 2d and Supplementary Table 9). implement any trimming or filtering of the input assembly contigs. Many of these genes were individually highly copy-number polymor- To characterize variants in the pangenome graphs, we used graph phic and part of complex tandem duplications (Fig. 2e). For example, decomposition to identify ‘bubble’ subgraphs that correspond to GPRIN2 is a copy-number polymorphic32 based on read depth and has non-overlapping variant sites. We then classified variant sites into small a sequence resolution of one to three additional copies duplicated variants (<50 bp) and SVs (≥50 bp) of different types (Methods). We in tandem in the pangenome (Fig. 2f). SPDYE2 is similarly resolved as found similar numbers of each variant type in each pangenome, with one to four additional copies duplicated in tandem (Fig. 2g). Other 22 million small variants in the MC graph (21 million in PGGB) (Fig. 3b), CNV genes were not contiguously resolved and reflect limitations of and 67,000 SVs in the MC graph (73,000 in PGGB, 75,000 in Minigraph) the current assemblies (see the associated article33). For example, the (Fig. 3c). We assessed variation in each individual assembly by tracing defensin gene DEFB107A has three to seven additional copies assembled their paths through the graphs and found similar numbers of small across all samples; however, this gene was assembled into three to seven variants and SVs within confident genomic regions defined by Dipcall38. separate contigs that do not reflect the global organization of this gene. Specifically, there were 5.34 million small variants per sample and 16,800 SVs per haplotype on average in the MC graph (5.35 million and 17,400, respectively in PGGB) (Fig. 3e,f). Differences in variant counts among Constructing a draft pangenome samples from different ancestry groups recapitulated previous observa- We used a sequence graph representation for pangenomes12,14 in which tions1. There was a total of 90 Mb of non-reference sequence in the SV nodes correspond to segments of DNA. Each node has two possible sites, excluding difficult-to-align centromeric repeats, in the MC graph orientations, forward and reverse, and there are four possible edges (55 Mb for PGGB, 86 Mb for Minigraph). Alu, L1 and ERV SVs appeared between any pair of nodes to reflect all combinations of orientations largely biallelic, whereas VNTRs frequently had three or more distinct (bidirected graph). The underlying haplotype sequences can be rep- alleles per site. The minor AF in the pangenomes of biallelic variants resented as walks in the graph. The model represents a generalized was similar for SNPs and for L1, Alu and VNTR variants, although VNTRs multiple alignment of the genome assemblies from which we built it, showed a slight shift towards more common alleles (Fig. 3d). whereby haplotypes are aligned where they co-occur on a given node We quantified the amount of euchromatic autosomal non-reference (Fig. 3a). (GRCh38) sequence that each of the 44 diploid genomes incremen- The process of generating a combined pangenome representation tally contributes to the pangenome (Fig. 3g and Methods) for both is an active research area. The problem is nontrivial both because of MC and PGGB graphs. We limited the analysis to the euchromatic computational challenges (there are hundreds of billions of bases sequence because we were generally confident in its assembly and of sequence to align) and because determining which alignments to alignment, and much of the heterochromatic sequence was omitted include is not always obvious, particularly for recently duplicated and from the MC graph (Methods). Overall, the euchromatic autosomal repetitive sequences. We applied three different graph construction non-reference sequence added up to about 175 Mb in the MC graph methods that have been under active development for this project: (around 190 Mb in PGGB), out of which about 55 Mb (around 105 Mb in Minigraph34, Minigraph-Cactus (MC)35 and PanGenome Graph Builder PGGB) was observed only on a single haplotype. Our analysis further (PGGB)36 (Extended Data Fig. 3 and Methods). The availability of these suggested that about 5 Mb and 70 Mb in the MC graph (around 10 Mb three models provided us with multiple views into homology rela- and 60 Mb in PCGB) could be attributed to core (present in ≥95% of all tionships in the pangenome while supporting validation of discov- haplotypes) and common genomes (present in ≥5% of all haplotypes), ered variation by independent methods. We included GRCh38 and respectively (Supplementary Table 11). We also estimated the growth T2T-CHM13 references within the pangenomes, and three samples of the euchromatic autosomal pangenome independent of the order (HG002, HG005 and NA19240) were held out to permit their use in of genomes by sampling 200 permutations (Supplementary Fig. 4) and benchmarking (hence 90 haplotypes total). In brief, Minigraph builds recording the median pangenome size across all samples in the MC a pangenome by starting from a reference assembly, here GRCh38, and graph. Our results indicated that the second genome added around iteratively and progressively adds in additional assemblies, recording 23 Mb of euchromatic autosomal sequence to the pangenome, whereas only SVs ≥ 50 bases. It admits complex variants, including duplications the last genome tended to add only about 0.64 Mb. These numbers are 316 | Nature | Vol 617 | 11 May 2023 a b ATTGGGCATCGGGTGAGAGTGACCCTTTAAGGCAGG 20 ATTGGGC-----GTGAGCGTGACCCCTTAAAGCAGG variant sites (×106) Number of small b d ATCGG A a c g 10 MC f h PGGB ATTGGGC GTGAG GTGACCC CTTAAA GCAGG e >3 alleles C 3 alleles 0 2 alleles SNP MNP INDEL c Minigraph e Alu MC MC PGGB L1 PGGB 6 small variants Number of SVA >3 alleles (×106) 3 alleles 3 ERV 2 alleles Other TE 0 Sample Sample VNTR STR PJL EAS AMR AFR Other LCR f DEL INS SNP Other repeat 20 Number of SVs SegDup (×103) 10 Low repeat 0 0 10 20 0 20 40 Haplotype Haplotype Number of SV sites (×103) Sum of non-reference node lengths (Mbp) d h Common Low frequency Rare Singleton Missing Fraction of biallelic sites 100 MC PGGB Alu MC L1 10–1 VNTR PGGB SNP Dipcall 10–2 HiFi–DeepVariant 10–3 0 50 100 0 0.5 0 0.5 Fraction of small variants (%) MAF MAF MC g PGGB 175 175 150 150 125 125 PJL EAS 100 100 AMR Mb Mb AFR Cumulative 75 75 all/depth ≥ 2 Common 50 50 Core 25 25 0 0 Sample Sample Fig. 3 | Pangenome graphs represent diverse variation. a, A pangenome regions without a clear VNTR or STR pattern. Other repeat, a site involving variation graph comprising two elements: a sequence graph, the nodes of mixed classes of repeats. SegDup, segmental duplication. Low repeat, a small which represent oriented DNA strings and bidirected edges represent the fraction of the longest allele in a site involving repeats. d, Pangenome minor AF connectivity relationships; and embedded haplotype paths (coloured lines) (MAF) spectrum for biallelic SNP, VNTR, L1 and Alu variants in the MC and PGGB that represent the individual assemblies. b, Small variant sites in pangenome graphs. e,f, Number of autosomal small variants per sample (e) and SVs per graphs stratified by the variant type and by the number of alleles at each site. haplotype (f) in the pangenome. Variants were restricted to the Dipcall-confident MNP, multinucleotide polymorphism. c, SV sites in the pangenome graphs regions. Samples are organized by 1KG populations. g, Pangenome growth stratified by repeat class and by the number of alleles at each site. Other TE, a curves for MC (left) and PGGB (right). Depth measures how often a segment is site involving mixed classes of transposable elements (TEs). VNTR, variable- contained in any haplotype sequence, whereby core is present in ≥95% of number tandem repeat, a tandem repeat with the unit motif length ≥7 bp. STR, haplotypes, common is ≥5%. h, Small variants in the GIAB (v.3.0) ‘easy’ regions short tandem repeat, a tandem repeat with the unit motif length ≤6 bp. Other annotated with AFs from gnomAD (v.3.1.2). LCR, low-complexity regions with mixed VNTR and STR and low-complexity conservative owing to additional highly polymorphic sequence resid- We annotated the small variants overlapping the GIAB (v.3.0) ‘easy’ ing in assembly gaps. Extrapolating under Heaps’ Law39 (Methods), we regions (covering 74.35% of GRCh38) with AFs from gnomAD (v.3.1.2) anticipate that at least an additional 150 Mb of euchromatic autosomal (Fig. 3h and Supplementary Table 12). In the MC graph, about 60.2% sequence will be added in the pangenome graph when HPRC produces (around 9.7 million variants) had an AF of 1% or greater. About 35.7% 700 haplotypes in the future. were rare, having an AF less than 1% but above zero. About 1.7% were Nature | Vol 617 | 11 May 2023 | 317 Article a b c MC MC PGGB MC PGGB 100 8,000 100 Recall (%) Variant calls Value (%) 90 50 4,000 80 70 0 0 100 PGGB Precision (%) 100 8,000 90 Variant calls Value (%) 80 50 4,000 70 Confident Easy Low map Tandem Other Confident Easy Low map Tandem Other and SegDup repeat difficult and SegDup repeat difficult 0 0 Average no. Average no. of variants 5.12 3.30 0.38 0.17 0.11 of variants 31.3 5.8 11.8 21.1 3.8 Recall 2, 0 –1 500 00 –7 00 –5 00 –2 00 0 0 5, 0 7, 0 50 00 50 1, 0 ≥1 500 0 00 25 50 50 00 75 00 ,0 ,0 ,5 ,0 (×106) (×103) Precision –7 –5 –2 , 0, 10 Variant length (bp) ≤– Fig. 4 | Pangenome graph evaluation. a,b, Precision and recall of autosomal (v.3.0) genomic context. c, Average SV precision, recall and frequency in the small variants (a) and SVs (b) in the pangenomes relative to consensus variant Dipcall-confident regions stratified by length in the MC (top) and PGGB sets. Small variants are compared to HiFi–DeepVariant calls. SVs are compared (bottom) graphs relative to consensus SV sets. The histogram bin size is 50 bp to the consensus of six reference-based SV callers (Methods). Comparisons are for SVs <1 kb and 500 bp for SVs ≥1 kb. restricted to the Dipcall-confident regions and then stratified by the GIAB singleton. The remaining 2.4% were missing from gnomAD. Simi- calculated the read depth of on-target and off-target edges based on lar results were obtained using the PGGB graph by repeating this the sample paths in the graph. On average, more than 94% of on-target exercise with small variant calls detected by pairwise alignment of edges were supported by at least 5 reads, and we observed 2 peaks in assemblies to GRCh38 using Dipcall38 and by calling small variants the read depth distribution of on-target edges (Extended Data Fig. 4, from the HiFi sequencing data using DeepVariant40. Given that 1KG middle): a minor peak corresponding to the edges in heterozygous samples are included in gnomAD, these missing variants are expected regions, and a major peak at twice the minor peak corresponding to the to be a mixture of false negatives in gnomAD and false positives in the edges in homozygous regions. By contrast, only 7% or fewer off-target pangenome. edges were supported by at least 5 reads (Extended Data Fig. 4, right). To further explore the quality of variant calls captured by assem- In addition to HiFi reads, we used ONT reads from 29 out of the 44 sam- bly and graph construction, we compared pangenome-decoded ples to perform the same analysis. Even though the data were lower in variants against GRCh38 to variant sets identified by conventional coverage, similar results were obtained (Supplementary Figs. 6 and 7). reference-based genotyping methods (Supplementary Fig. 5 and These data also show that the pangenome graphs performed better Methods). These reference-based call sets were generated from the at capturing genome variation than the above benchmarking results PacBio HiFi reads and haplotype-resolved assemblies using the follow- imply. For example, a mean of 89.3% of putative false-positive small vari- ing different discovery methods: DeepVariant40, PBSV41, Sniffles42 with ant calls were supported by ≥5 HiFi reads, and 75.3% by ≥10 reads (85.9% Iris43, SVIM44, SVIM-asm45, PAV5 and the Hall-lab pipeline (Methods). For and 73.8%, respectively, for SVs). This result suggests that most putative benchmarking small variants, we excluded regions that contained SVs errors are in fact real variants that were missed by the reference-based detected or implied by the alignment of the haploid assemblies of that callers used to create the truth set (Supplementary Fig. 8 and Sup- sample to GRCh38, as current benchmarking tools do not account for plementary Table 13). different representations of small variants inside or near SVs (Meth- To assess gene alignments in the pangenome, we used the Compara- ods). Comparing small variants (Fig. 4a) and SVs (Fig. 4b) from the tive Annotation Toolkit (CAT)46 to liftover GENCODE (v.38) annotations pangenomes to the reference-based sets, we observed a high level using the MC pangenome alignment onto the individual haplotype of concordance that varied, as expected, by the relative repeat con- assemblies. CAT lifted and annotated a median of 99.1% of 86,757 tent of the surrounding genome. Overall, variant calling performance protein-coding transcripts per assembly (Extended Data Fig. 5, Supple- was high in both the MC and PGGB graphs. For example, in relatively mentary Fig. 9, Supplementary Tables 14 and 15 and Methods), making it unique easy genomic regions constituting 75.42% of the autosomal comparable to the Ensembl-mapping-based pipeline (median of 99.4% genome, samples showed a mean of 99.64% recall and 99.64% preci- per assembly). This result supports the idea that the MC pangenome sion for small variants in the MC graph. Meanwhile, in high-confidence captures most transcript homologies. When comparing the CAT and regions (around 90% of autosomal genome), samples showed 97.91% Ensembl annotations per assembly, median Jaccard similarities of 0.99 recall and 96.66% precision (Fig. 4a). Performance was lower for SVs for both genes and transcripts were obtained (Methods). A median of than for small variants (Fig. 4b), as expected, but was still strong. Vari- 360 (0.4%) protein-coding transcripts per assembly mapped at differ- ant calling performance was lower in highly repetitive genome regions ent loci between the Ensembl and CAT annotations. (3.87% of autosomal genome; Fig. 4a,b), for which more work will be required to achieve high-quality variant maps. These values are likely Pangenomes represent complex loci to be significant underestimates of variant calling quality, consider- We next turned our attention to complex multiallelic SVs, which have ing known errors in the truth set owing to the inherent limitations of historically been difficult to map using reference-based methods. To reference-based variant callers (see below). Stratifying the insertion screen for complex SVs, we identified bubbles >10 kb from Minigraph and deletion SVs within the pangenome, we observed relatively high that exhibited at least five structural alleles among the assembled hap- levels of agreement with the reference-based methods regardless of lotypes (Methods). We found that 620 out of 76,506 total sites (0.81%) length (Fig. 4c). were complex, and 44 of these overlapped with medically relevant An independent measure of the quality of the pangenome graphs protein-coding genes47 (Supplementary Table 16). Some are well-known is the extent to which sample haplotype paths through the graph are complex SV loci, and all are known to be structurally variable based well supported by the raw sequencing data. When we calculated the on previous short-read SV mapping studies10,19,32. However, whereas number of supporting reads by aligning them to the MC graph using previous short-read SV calls at these loci are typically imprecise owing GraphAligner (Methods), more than 97% of HiFi reads were aligned to alignment issues and low-resolution read-depth analysis methods, to the MC graph after filtering (Extended Data Fig. 4, left). We further here we resolved their structure at single-base resolution. We selected 318 | Nature | Vol 617 | 11 May 2023 a RHD d HLA-Y HCG4B HLA-K HLA-U RHD HLA-A HCG9 HLA-H RHCE HLA-J HCP5B HLA-W b RDH;RHCE RHCE e HLA-A HLA-Y ins GRCh38 HG00733-2 GRCh38 HG00735-2 RHD-RHCE(10);inv;RHCE-RHD(10) RHD;RHD;RHCE-RHD(9)-RHCE HLA-H/HCG4B/HLA-K/HLA-U del HLA-W del HG03516-2 HG01106-1 HG00438-1 HG03453-1 c f Count Frequency Haplotype name Gene 43 0.48 RHD;RHCE RHD RHCE 22 0.24 RHD;RHCE-RHD(2)-RHCE TMEM50A 15 0.17 RHCE Gene Count Frequency Haplotype name HCP5B 3 0.03 RHD-RHCE(2-3)-RHD;RHCE HLA−H 57 0.63 HLA-A HLA−K 2* 0.02 RHD-RHCE(8)-RHD;RHCE 25 0.28 HLA-Y ins HLA−U HLA−A 1 0.01 RHCE-RHD(2)-RHCE HLA−W 7 0.08 HLA-H/HCG4B/HLA-K/HLA-U del HCG9 HLA−Y 1* 0.01 RHD-RHCE(2-9)-RHD;RHCE 1 0.01 HLA-W del HLA−J 1* 0.01 RHD;RHCE-RHD(9)-RHCE 1* 0.01 RHD-RHCE(10);inv;RHCE-RHD(10) 1* 0.01 RHD;RHD;RHCE-RHD(9)-RHCE Fig. 5 | Visualizing complex pangenome loci. a–c, Structural haplotypes of through the graph (b). The colour gradient and lines show the path of each RHD and RHCE from the MC graph. Locations of RHD and RHCE within the allele; red, start of a path; blue, end of a path. Frequency and linear structural graph (a). The colour gradient is based on the precise relative position of each visualization of all structural haplotypes called by the graph among 90 haploid gene; green, head of a gene; blue, end of a gene. The lines alongside the graph assemblies (c). Asterisks indicate newly discovered haplotypes. d–f, Structural are based on the approximate position of gene bodies, including exons and haplotypes of HLA-A from the PGGB graph, visualized using the same transcription start sites. Different structural haplotypes take different paths conventions as a–c. del, deletion; ins, insertion; inv, inversion. five clinically relevant complex SV loci for detailed structural analy- previously unreported. The long sequence (65 kb) inserted with HLA-Y sis: RHD–RHCE, HLA-A, CYP2D6–CYP2D7, C4 and LPA (Methods). For occurred at high frequency (28%) but has little homology to GRCh38. each locus and graph, we identified their locations within the graph We also compared the representation of these five loci in the MC and then annotated paths within this subgraph with known genes. and PGGB graphs (Supplementary Fig. 13). Each graph independently We traced the individual haplotypes through the subgraph to reveal recapitulated the same haplotype structures. In general, in the PGGB the structure of each assembly. In CYP2D6–CYP2D7 (Extended Data graph, many SV hotspots, including the centromeres, were transitively Fig. 6), C4 (Supplementary Fig. 10) and LPA (Supplementary Fig. 11), we collapsed into loops through a subgraph representing a single repeat recapitulated previously described haplotypes. For CYP2D6–CYP2D7, copy. This feature tends to reduce the size of variants found in repetitive our calls matched 96% of haplotypes of 76 assemblies called by Cyrius sequences. Assemblies that contained multiple copies of the homol- using Illumina short-read data48. Two discrepancies appeared to be ogous sequence traversed these nodes a corresponding number of caused by errors from Cyrius, and the third was a false duplication in times. By contrast, the MC graph maintained separate copies of these the HG01071-2 pangenome assembly revealed by Flagger. This com- homologous sequences. parison suggests that the pangenomes faithfully agree with existing knowledge of this complex locus. In RHD–RHCE (Fig. 5a–c), in addition to previously described haplotypes, we inferred the presence of five Applications of the pangenome new haplotypes, which included one duplication allele of RHD and Pangenome-based short variant discovery one inversion allele between RHD and RHCE that swaps the last exon of Our pangenome reference aims to broadly improve downstream analy- both genes. Around HLA-A (Fig. 5d–f and Supplementary Fig. 12), two sis workflows by removing mapping biases that are inherent in the use deletion alleles have been previously described—albeit with imprecise of a single linear reference genome such as GRCh38 or CHM13. As an breakpoints10—but an insertion allele carrying a HLA-Y pseudogene was initial test case, we studied whether mapping against our pangenomes Nature | Vol 617 | 11 May 2023 | 319 Article a c SV deletions (filtered set) SV insertions (filtered set) SV other (filtered set) 1.0 1.0 1.0 HPRC Giraffe 104 0.998 DeepTrio 103 F1: 0.9986 0.8 0.8 0.8 103 103 AF PanGenie AF PanGenie AF PanGenie DragenGRAPH DragenGRAPH Dragen DeepTrio 0.6 0.6 0.6 F1: 0.9966 F1: 0.9975 102 102 0.4 0.4 102 0.4 0.996 HPRC Giraffe DeepVariant F1: 0.9971 101 101 0.2 0.2 101 0.2 Recall 0 0 0 100 100 100 0.994 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 d AF assemblies AF assemblies AF assemblies BWAMEM BWAMEM DeepTrio DeepVariant F1: 0.9964 F1: 0.9951 25,000 HPRC PanGenie filtered Number of SV alleles per sample 0.992 HGSVC PanGenie lenient 1KG Illumina calls 20,000 15,000 0.997 0.998 0.999 1.000 Precision b 10,000 0.98 HPRC Giraffe 5,000 DeepTrio F1: 0.9836 HPRC Giraffe 0 DeepVariant F1: 0.9809 All Easy Low mappability Other difficult Tandem repeat 0.97 and SegDup DragenGRAPH DragenGRAPH Dragen DeepTrio e F1: 0.9733 F1: 0.9768 103 HPRC PanGenie filtered Recall HGSVC PanGenie lenient BWAMEM 1KG Illumina calls DeepTrio Number of SV alleles 0.96 F1: 0.9738 102 BWAMEM DeepVariant 101 F1: 0.9712 0.95 100 0.980 0.985 0.990 0.995 0 Precision 106 105 104 103 102 102 103 104 105 106 Deletions Insertions Fig. 6 | Performance gains for pangenome-aided analysis of short-read WGS (28,433 deletions, 84,755 insertions, 32,431 other alleles). d, Number of SVs data. a,b, Precision–recall curves showing the performance of different present (genotype 0/1 or 1/1) in each of the 3,202 1KG samples in the filtered combinations of linear reference and various mappers and variant callers HPRC genotypes (PanGenie) after merging similar alleles (n = 100,442 SVs), the evaluated against the GIAB (v.4.2.1) HG005 benchmark (a) and the challenging HGSVC lenient set (n = 52,659 SVs) and the 1KG Illumina calls (n = 172,968 SVs) medically relevant genes (CMRG; v.1.0) benchmark (b). Giraffe uses the MC in GIAB regions. In the box plots, lower and upper limits represent the first and pangenome graph, BWA-MEM uses GRCh38 and Dragen Graph uses GRCh38 third quartiles of the data, the white dots represent the median and the black with additional alternative haplotype sequences. c, Comparison of AFs lines mark minima and maxima of the data points. e, Length distribution of observed from the PanGenie genotypes for all 2,504 unrelated 1KG samples SV insertions and SV deletions contained in the filtered HPRC genotypes and the AFs observed across 44 of the HPRC assembly samples in the MC graph. (PanGenie), the HGSVC lenient set and the 1KG Illumina calls. Only variants with The PanGenie genotypes include all variants contained in the filtered set a common AF > 5% across the 3,202 samples were considered. could improve the accuracy of calling small variants from short reads. pangenome mapping, but regions with errors in GRCh38 and large We used Giraffe49 to align short reads from the GIAB benchmark sam- L1HS sequences benefitted the most from the pangenomic approach ples18 to the MC pangenome graph. For comparison, we aligned reads to (Supplementary Figs. 15 and 16). GRCh38 using BWA-MEM50 and to Dragen Graph51, which uses GRCh38 We next benchmarked variant calling using parent–child trios. Using augmented with alternative haplotypes at variant sites. We called SNPs DeepTrio52 resulted in better performance compared with DeepVariant and indels with DeepVariant40 and the Dragen variant caller51 (Methods). across all samples of the GIAB (Fig. 6a) and the challenging medi- Our pangenomic approach (Giraffe plus DeepVariant) outperformed cally relevant gene benchmarks (Fig. 6b and Supplementary Fig. 14). the other approaches for calling small variants (Fig. 6a), with gains Moreover, improvements appeared to be additive to those from the for both SNPs and indels (Supplementary Fig. 14 and Supplementary pangenome. For example, DeepTrio using Giraffe alignments gave Table 17). For example, it made 21,700 errors (false positives or false the highest calling accuracy, with the number of errors decreasing negatives) in the confident regions of the GIAB truth set using 30× reads from 21,700 (single sample calling) to 10,098 (trio calling) for HG005. from HG005. By contrast, 36,144 errors were made when DeepVariant used the reads aligned to GRCh38, and 26,852 errors when using the A pangenome variant resource Dragen pipeline. In challenging medically relevant genes47, the increase To create a community resource to aid the development of methods in performance was even larger for both SNPs (F1 score, defined as and the analyses of pangenome-based population genetics, we used the harmonic mean of precision and recall, of 0.985 for Giraffe plus Giraffe to align high-coverage short-read data from 3,202 samples of the DeepVariant compared with <0.976 for the other methods) and indels 1KG19 to our pangenome graph and DeepVariant to call small variants (F1 score of 0.961 for Giraffe plus DeepVariant compared with <0.958 (Methods). The Mendelian consistency computed across 100 trios from for the other methods) (Fig. 6b). Many regions benefitted from using those samples was comparable to the one computed across samples 320 | Nature | Vol 617 | 11 May 2023 from the GIAB truth set, which indicated that comparable call set qual- 9,227 for HGSVC and 6,099 for the 1KG calls (Supplementary Fig. 27); ity was obtained (Supplementary Figs. 17 and 18). The number of small that is, a gain of 104.0% for HPRC over 1KG and of 34.8% over HGSVC. variants called was consistently higher across different ancestries, with This result confirms that short-read-based SV discovery relative to a on average 64,000 more variants per sample compared with the 1KG linear reference genome misses a large proportion of SVs5,6,8. As antici- catalogue (Extended Data Fig. 7a). Given that our pangenome-based pated, the number of SVs per sample within ‘easy’ genomic regions was calls showed improved performance in challenging regions (Fig. 6b), consistent across all three call sets, particularly in low-mappability and this call set across the 1KG cohort now provides the genetics and genom- tandem repeat regions, and the use of our pangenome reference led ics communities with AF estimates for complex but medically relevant to substantial gains (Fig. 6d), including for common variants (Fig. 6e loci. For example, our approach was able to detect the gene conversion and Supplementary Fig. 28). Although the newly identified SVs were event covering the second exon RHCE, which was observed in about 25% harder to genotype because they are primarily located in repetitive of assembled haplotypes (Fig. 5c and Extended Data Fig. 8). Moreover, regions, genotype concordances were high and close to the ones for for KCNE1, we provide calls and frequencies in a 40 kb region, spanning known SVs (Supplementary Fig. 29). 3 exons, that could not be previously assessed owing to the presence of a false duplication in GRCh38 (Supplementary Fig. 19; see also an associ- Improved tandem repeat representation ated article53 for genome-wide analysis of interlocus gene conversion). VNTRs are particularly variable among individuals and are challenging to access with short reads. The gains in the number of genotyped SVs in SV genotyping VNTRs (Fig. 6d and Supplementary Fig. 28) prompted us to investigate The ability to represent polymorphic SVs is a key advantage of a whether our pangenome reference could also improve read mapping graph-based pangenome reference. To demonstrate the utility of the in VNTR regions. We first established orthology mapping between sequence-resolved SVs inherent to our pangenome, we used PanGenie54 haplotypes in our pangenome reference using danbing-tk55. The orthol- to genotype the bubbles in the MC graph. We decomposed bubbles ogy can be established for 94,452 out of the 98,021 VNTR loci (96.4%) into their constituent variant alleles (Supplementary Figs. 20 and 21) discovered by TRF56. When mapping simulated short reads to GRCh38 and found that 22,133,782 bubbles represented 20,194,117 SNP alleles, with BWA-MEM, the rate of unmapped reads was 6.6–8.5 times greater 6,848,115 indel alleles and 413,809 SV alleles (Supplementary Fig. 22 compared with mapping to the MC graph with Giraffe (Extended Data and Methods). Of these SV alleles that were non-reference (neither Fig. 10a, Supplementary Fig. 30 and Supplementary Table 19). The true GRCh38 nor T2T-CHM13), 17,720 were observed in biallelic contexts and negatives were on average 1.9% higher than the GRCh38 approach, and 396,089 at multiallelic loci with more than 1 non-reference allele, includ- the true positives were on average 0.087% higher. The graph approach ing extreme cases in which all 88 haplotypes showed distinct alleles also reduced false negatives by 2.1-fold. Read depth over a locus is cor- (Supplementary Fig. 22). To analyse the genotyping performance of related with the copy number of a duplication and we evaluated how PanGenie, we conducted a leave-one-out experiment in which we repeat- well length variants in VNTR regions can be estimated using either the edly removed one sample from the graph and re-genotyped it using the MC graph or GRCh38. The graph approach performed better for 80% of remaining haplotype paths in the graph and short-read data for the the loci (48,085 out of 60,386) and increased the median r2 from 0.58 left-out sample (Methods). In line with previous results5,54, we obtained to 0.70 (Supplementary Fig. 31). high genotype concordance across all variant types and genomic contexts (Extended Data Fig. 9). Furthermore, we used PanGenie Improved RNA sequencing mapping to genotype HG002 and evaluated genotypes based on SVs at chal- To evaluate the benefit of our pangenome reference on transcriptom- lenging medically relevant loci47. This analysis resulted in a precision ics, we simulated RNA sequencing (RNA-seq) reads and mapped them of 0.74 and an adjusted recall of 0.81 (Methods). to a pangenome and to a standard reference genome (Methods). The Next we genotyped the 3,202 samples from the 1KG19 (Methods). pangenome-based pipeline using vg mpmap57 achieved significantly We filtered the resulting SV genotypes using a machine-learning lower false mapping rates than a linear reference pipeline using either approach5,54 that assessed different statistics, including Mendelian vg mpmap or STAR58 (Extended Data Fig. 10b). Compared with the linear consistency and concordance, to assembly based calls. As a result, we reference pipelines, the pangenome pipeline also showed reduced produced a filtered, high-quality subset of SV genotypes containing allelic bias and increased mapped coverage on heterozygous variants, 28,434 deletion alleles, 84,752 insertion alleles and 26,439 other SV which could benefit studies of allele-specific expression (Supplemen- alleles (Supplementary Table 18 and Methods). Many of the alleles not tary Fig. 32). With real sequencing data, mapping rates were more dif- included in the filtered set stemmed from complex, multiallelic loci and ficult to interpret in the absence of a ground truth (Supplementary were enriched for rare alleles. As independent quality control measures Fig. 33). Instead, we focused on the correlation in exon coverage to for genotypes in the filtered set, we assessed the Hardy–Weinberg independent PacBio long-read Iso-Seq data. The analysis showed that equilibrium values (Supplementary Figs. 23–25) and compared AFs the correlation was highest when mapping to a spliced pangenome observed across the genotypes of all 2,504 unrelated samples to the graph derived from the MC graph (Supplementary Fig. 34). The pange- respective AFs of the 44 assembly samples (88 haplotypes) contained nome pipeline showed a modest increase in correlation over the linear in the graph. Pearson correlation values of 0.96, 0.93 and 0.90 for the reference pipelines (0.006–0.011). In addition, mapping the simulated deletion, insertion and other SV alleles, respectively, were observed reads to the MC graph led to improved gene expression estimates rela- (Fig. 6c), which indicated the high quality of the genotypes. To quan- tive to the linear GRCh38, regardless of whether alternative contigs tify our ability to detect additional SVs, we compared our filtered set were included in GRCh38 (Supplementary Fig. 35). of genotypes to the HGSVC PanGenie genotypes (v.2.0 ‘lenient’ set)5 and the Illumina-based 1KG SV call set19. We analysed the number of Improved chromatin immunoprecipitation and sequencing detected SV alleles in each sample (homozygous or heterozygous) and analysis stratified them by genome annotations from GIAB (Fig. 6d and Meth- We used the pangenome to re-analyse H3K4me1 and H3K27ac data from ods) as well as using our own more detailed annotations (Supplemen- chromatin immunoprecipitation and sequencing (ChIP-seq) and assay tary Fig. 26). Both of the PanGenie-based call sets detected more SVs for transposase-accessible chromatin with high-throughput sequenc- (HPRC, 18,483 SVs per sample; HGSVC, 12,997 SVs per sample) than the ing (ATAC-seq) of monocyte-derived macrophages from 30 individuals short-read-based 1KG call set (9,596 SVs per sample), with a particularly with African ancestry or European ancestry59. Overall, we observed a substantial advance for deletions <300 bp and insertions (Fig. 6e). The net increase in the number of peak calls, whereby, on average, 2–3% of respective average numbers of SVs per haplotype were 12,439 for HPRC, peaks were found only when using the MC pangenome (Extended Data Nature | Vol 617 | 11 May 2023 | 321 Article Fig. 10c). Moreover, the newly found peaks were replicated in more existing linear reference-based methods for variant discovery, particu- samples than expected by chance (Supplementary Fig. 36). In addition, larly within the non-repetitive fraction of the genome. Our study of we recovered epigenomic features that were specific to SV alleles not complex and medically relevant loci showed that the pangenomes faith- present in GRCh38 (termed non-reference). For example, across all fully recapitulated existing knowledge and will enable future efforts H3K4me1 datasets, we assigned 1,326 events to the non-reference SV to study the role of complex variation in human disease. Further work allele, 1,443 to the reference allele and 2,008 to both alleles within het- will be required to more comprehensively identify medically relevant erozygous SVs (Extended Data Fig. 10d), with some replicated multiple complex SVs and to ensure the accuracy of each allele represented in times across samples (Supplementary Fig. 37). Of these, there were 194 the pangenome. SVs with peaks that were observed only in African ancestry genomes, Where the pangenome graphs differ is principally in how they handle 150 that were observed only in European ancestry genomes and 216 CNV sequences. The PGGB method will frequently merge CNVs, whereas that were observed in both African and European ancestry genomes. the MC graphs represent CNV copies as independent subgraphs. Both As expected, rare alleles were enriched for ancestry-specific events approaches have merits, and which approach to favour will take further (Supplementary Fig. 37). experimentation and community input, and may vary by the specific application. The PGGB method retained all centromeric and satellite sequences, whereas the MC graph pruned much of this sequence. This Discussion made it practical with current methods to use the MC graphs for read We have publicly released 94 de novo haplotype assemblies from a alignment applications. However, pruning these sequences is not a diverse group of 47 individuals. This provides a large set of fully phased satisfactory solution. Longer term, more work is needed to determine human genome assemblies and outperforms earlier efforts on many how best to align and represent these large repeat arrays within pange- levels of assembly quality5,16,60. For example, compared with Ebert nomes, particularly as T2T assembly becomes commonplace and these et al.5, the average median base level accuracy is nearly an order of arrays therefore completed. Furthermore, although the PGGB graph magnitude higher, the N50 contiguity is nearly double and the struc- retained centromeric and satellite sequences, in principle, by ena- tural accuracy is higher33. These improvements are the result of recent bling analysis of previously inaccessible parts of the pangenome, our improvements in de novo assembly driven both by better sequencing initial population-genetic analysis of these regions (Methods) leaves technology and coordinated innovations in assembly algorithms20,21. open questions about assembly accuracy and alignment, especially in To validate assembly structural accuracy, we developed a new pipeline areas of the genome where mutation rates are thought to be an order of that maps low error, long reads to each diploid assembly to support magnitude greater66. This suggests that significant care must be taken the predicted haplotypes. This pipeline indicated that more than 99% when studying them, and new methods may need to be developed of each assembly, and greater than 90% of the assembled sequence to fully understand and characterize this component of the human representing highly repetitive arrays, was structurally correct. Some pangenome. challenges around loci that harbour copy number polymorphisms and/ A near-term application of pangenome references will be to improve or inversions remained33. Although the focus of this effort was to build reference-based sequence mapping workflows. In these workflows, the a reference resource, highly accurate haplotype-resolved assemblies pangenome can act as a drop-in replacement for existing references, enabled us to access previously inaccessible regions, highlighting new with the read mappings projected from the pangenome space back forms of genetic variation and providing new insights into mutational onto an existing linear reference for downstream processing. This is processes such as interlocus gene conversion53. how the Giraffe–DeepVariant workflow functions: DeepVariant, the Accompanying these assemblies are 94 sets of Ensembl gene anno- variant caller, never needs to consider the complexity of the pange- tations, representing a large collection of de novo assembled human nome, but the workflow benefits from a mapping step that accounts transcriptome annotations. Each transcriptome annotation is highly for sequences that are missing from the linear reference. Making the complete, particularly for protein-coding transcripts. These putative switch to using pangenome mapping is not significantly more com- transcriptome annotations enabled us to analyse sequence-resolved putationally expensive49 and resulted in an average 34% reduction CNVs. In detail, we assembled genic CNVs (mostly singletons) for 1,115 in false-positive and false-negative errors compared with using the different protein-coding genes, confirming earlier mapping-based standard reference methods (Supplementary Fig. 38). These benefits analyses that predicted that the majority of rare genic CNVs occur out- were also greatest at complex loci47. Pangenomes not only improve side known SDs31. These CNV genes accounted for 0.6–4.4 Mb of addi- variant calling but also improve transcript mapping accuracy57 and tional genic sequences per haplotype compared with GRCh38. These detection of ChIP-seq peaks67. contained genes known to have CNVs associated with human health, SVs have been mostly excluded from short-read studies because including amylase61 (four to ten copies), β-defensin62 (three to seven cop- methods to genotype them using a linear reference have limited ies, DEFB107A) and NOTCH2NLC–NOTCH2NLB63 (one additional copy). accuracy and sensitivity. Previous short-read, linear reference stud- The pangenomes presented here are both a set of individual haploid ies have discovered 7,500–9,500 SVs per sample19,68, whereas long-read genome assemblies and an alignment of these assemblies. The combi- sequencing efforts have routinely discovered around 25,000. Ebert nation can be efficiently described as a variation graph14,64. A new set of et al.5 showed that using PanGenie, a pangenomic approach, with 32 exchange formats for pangenomics, including extensions of the graphi- samples, a subset of these variants could be genotyped in short-read cal fragment format (GFA) that encode variation graphs, are emerging34. genomes (about 13,000 genotyped on average, ranging from 12,000 An associated article65 to this work demonstrated that the pangenomes to 15,000 per sample). Using the same PanGenie method, the HPRC presented here can be losslessly stored using a compressed, binary pangenome increases this to around 18,500 (ranging from 16,900 to representation of GFA in just 3–6 GB despite representing more than 24,900) per sample, enabling the genotyping of the substantial majority 282 billion bases of individual sequence, with strongly sublinear scaling of SVs discovered using long-reads per sample. The draft pangenome as new genomes are added. Creating pangenome graphs is an active therefore delivers better SV calling than previous approaches, extract- research topic, so we developed multiple pipelines, and details of these ing latent information from short-read samples that are already avail- methods are further explored in companion papers35,36. We found con- able. So, in the future, the pangenome will enable the inclusion of tens cordance between the different construction approaches used here, of thousands of additional SV alleles into genome-wide association whereby the MC and PGGB pangenomes contained nearly the same studies. Looking beyond short reads, in the future, the combination number of small variants and SVs of various types. Furthermore, these of the pangenome and low-cost long-read sequencing should prove to encoded pangenome variants showed high levels of agreement with be a potent combination for comprehensive SV genotyping. 322 | Nature | Vol 617 | 11 May 2023 These new pangenomic workflows could benefit individuals of differ- 1. 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature 526, 68–74 (2015). ent ancestries differently. For read mapping and small variant calling, 2. International Human Genome Sequencing Consortium. Initial sequencing and analysis of we observed a consistent improvement across individuals (Extended the human genome. Nature 409, 860–921 (2001). Data Fig. 7). Moreover, the pangenome might improve SV genotyping 3. Nurk, S. et al. The complete sequence of a human genome. Science 376, 44–53 (2022). 4. Aganezov, S. et al. A complete reference genome improves analysis of human genetic differently across individuals owing to the stronger divergence of the variation. Science 376, eabl3533 (2022). alleles from the reference. In the 1KG cohort, we observed that the 5. Ebert, P. et al. Haplotype-resolved diverse human genomes and integrated analysis of genotyped samples clustered by super-population labels (Extended structural variation. Science 372, eabf7117 (2021). 6. Chaisson, M. J. P. et al. Multi-platform discovery of haplotype-resolved structural variation Data Fig. 11), which would suggest different levels of detection bias in human genomes. Nat. Commun. 10, 1784 (2019). that are mitigated with the pangenome. However, we caution that the 7. Wenger, A. M. et al. Accurate circular consensus long-read sequencing improves variant composition of the samples underlying the pangenome relative to detection and assembly of a human genome. Nat. Biotechnol. 37, 1155–1162 (2019). 8. Zhao, X. et al. Expectations and blind spots for structural variation detection from long-read the composition of the set of samples genotyped could potentially assemblies and short-read genome sequencing technologies. Am. J. Hum. Genet. 108, influence these results; an analysis with more samples is warranted. 919–928 (2021). The openly accessible, diverse assemblies and pangenome graphs 9. Chiang, C. et al. The impact of structural variation on human gene expression. Nat. Genet. 49, 692–699 (2017). we present here form a draft of a pangenome reference. There are 10. Sudmant, P. H. et al. An integrated map of structural variation in 2,504 human genomes. many remaining challenges to growing and refining this reference. Nature 526, 75–81 (2015). For example, assembly reliability analysis revealed roughly an order 11. Computational Pan-Genomics Consortium. Computational pan-genomics: status, promises and challenges. Brief. Bioinform. 19, 118–135 (2018). of magnitude more erroneously assembled sequences in the HPRC 12. Paten, B., Novak, A. M., Eizenga, J. M. & Garrison, E. Genome graphs and the evolution of assemblies than in the T2T-CHM13 complete assembly. Similarly, in genome inference. Genome Res. 27, 665–676 (2017). a companion analysis, Strand-seq data from a subset of assemblies 13. Garrison, E. et al. Variation graph toolkit improves read mapping by representing genetic variation in the reference. Nat. Biotechnol. 36, 875–879 (2018). revealed 6–7 Mb of incorrectly oriented sequence per haplotype33, 14. Eizenga, J. M. et al. Pangenome graphs. Annu. Rev. Genomics Hum. Genet. 21, 139–162 which indicates that there is room to structurally improve the assem- (2020). blies. Furthermore, despite being predicted to have less than 1 base 15. Wang, T. et al. The Human Pangenome Project: a global resource to map genomic diversity. Nature 604, 437–446 (2022). error per around 200,000 assembled bases, base level sequencing 16. Shafin, K. et al. Nanopore sequencing and the Shasta toolkit enable efficient de novo errors are still an issue. For example, we identified more than a dozen assembly of eleven human genomes. Nat. Biotechnol. 38, 1044–1053 (2020). 17. Porubsky, D. et al. Fully phased human genome assembly without parental data using apparent frameshifts and nonsense mutations per genome annotation single-cell strand sequencing and long reads. Nat. Biotechnol. 39, 302–308 (2021). that are probably the result of sequencing errors. The cohort we present 18. Zook, J. M. et al. Extensive sequencing of seven human genomes to characterize benchmark is also relatively small notwithstanding the significant effort to generate reference materials. Sci. Data 3, 160025 (2016). 19. Byrska-Bishop, M. et al. High-coverage whole-genome sequencing of the expanded 1000 the underlying long-read sequencing resource. Our near-term goal is Genomes Project cohort including 602 trios. Cell 185, 3426–3440.e19 (2022). to expand the pangenome to a diverse cohort of 350 individuals (which 20. Cheng, H., Concepcion, G. T., Feng, X., Zhang, H. & Li, H. Haplotype-resolved de novo should capture most common variants), to push towards T2T genomes assembly using phased assembly graphs with hifiasm. Nat. Methods 18, 170–175 (2021). for this cohort (to properly represent the entire genome in almost all 21. Jarvis, E. D. et al. Semi-automated assembly of high-quality diploid human reference individuals) and to refine the pangenome alignment methods (so that genomes. Nature 611, 519–531 (2022). telomere-to-telomere alignment is possible, capturing more complex 22. Guarracino, A. et al. Recombination between heterologous human acrocentric chromosomes. Nature https://doi.org/10.1038/s41586-023-05976-y (2023). regions of the genome). This will give us a more comprehensive repre- 23. Garg, S. Computational methods for chromosome-scale haplotype reconstruction. sentation of all types of human variation. Genome Biol. 22, 101 (2021). We acknowledge that references generated from the 1KG samples 24. Garg, S. Pstools: a toolkit for fully phased sequences on chromosome level. GitHub https://github.com/shilpagarg/pstools (2020). alone are insufficient to capture the extent of sequence diversity in the 25. Rahman, A. & Pachter, L. CGAL: computing genome assembly likelihoods. Genome Biol. human population. To ensure that we are able to maximize our surveys 14, R8 (2013). of sample diversity while abiding by principles of community engage- 26. Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0, http://www.repeatmasker.org (2013–2015). ment and avoiding extractive practices14,15, we will broaden our efforts 27. Altemose, N. et al. Complete genomic and epigenetic maps of human centromeres. Science to recruit new participants to improve the representation of human 376, eabl4178 (2022). genetic diversity. A richer human reference map promises to improve 28. Alkan, C., Sajjadian, S. & Eichler, E. E. Limitations of next-generation genome sequence assembly. Nat. Methods 8, 61–65 (2011). our understanding of genomics and our ability to predict, diagnose and 29. Frankish, A. et al. GENCODE 2021. Nucleic Acids Res. 49, D916–D923 (2021). treat disease. A more diverse human reference map should also help 30. MacArthur, D. G. et al. A systematic survey of loss-of-function variants in human protein- ensure that the eventual applications of genomic research and precision coding genes. Science 335, 823–828 (2012). 31. Sudmant, P. H. et al. Diversity of human copy number variation and multicopy genes. medicine are effective for all populations. We recognize that the value of Science 330, 641–646 (2010). this project will partly be in the future establishment of new standards 32. Handsaker, R. E. et al. Large multiallelic copy number variations in humans. Nat. Genet. for how we capture variant diversity, the opportunity to disseminate 47, 296–303 (2015). 33. Porubsky, D. et al. Gaps and complex structurally variant loci in phased genome assemblies. science into diverse communities and continued efforts to engage with Genome Res. https://doi.org/10.1101/gr.277334.122 (2023). diverse voices in this ambitious goal to build a common global reference 34. Li, H., Feng, X. & Chu, C. The design and construction of reference pangenome graphs resource. The methods we are developing should prove valuable for with minigraph. Genome Biol. 21, 265 (2020). 35. Hickey, G. et al. Pangenome graph construction from genome alignment with Minigraph- other species. Indeed, other groups are pioneering such efforts69,70. In Cactus. Preprint at bioRxiv https://doi.org/10.1101/2022.10.06.511217 (2022). parallel with our efforts to obtain a more comprehensive collection of 36. Garrison, E. et al. Building pangenome graphs. Preprint at bioRxiv https://www.biorxiv. diverse and highly accurate human reference genomes, we anticipate org/content/10.1101/2023.04.05.535718v1 (2023). 37. Armstrong, J. et al. Progressive Cactus is a multiple-genome aligner for the thousand- further optimization and rapid improvement of the pangenome refer- genome era. Nature 587, 246–251 (2020). ence, enabling an increasingly broad set of applications and use cases 38. Li, H. et al. A synthetic-diploid benchmark for accurate variant-calling evaluation. Nat. for both the research and clinical communities. Methods 15, 595–597 (2018). 39. Tettelin, H., Riley, D., Cattuto, C. & Medini, D. Comparative genomics: the bacterial pan- genome. Curr. Opin. Microbiol. 11, 472–477 (2008). 40. Poplin, R. et al. A universal SNP and small-indel variant caller using deep neural networks. Online content Nat. Biotechnol. 36, 983–987 (2018). 41. Pacific Biosciences. PBSV: a suite of tools to call and analyze structural variants in diploid Any methods, additional references, Nature Portfolio reporting summa- genomes from PacBio SMRT reads. GitHub https://github.com/PacificBiosciences/pbsv ries, source data, extended data, supplementary information, acknowl- (2021). edgements, peer review information; details of author contributions 42. Sedlazeck, F. J. et al. Accurate detection of complex structural variations using single- molecule sequencing. Nat. Methods 15, 461–468 (2018). and competing interests; and statements of data and code availability 43. Kirsche, M. et al. Jasmine and Iris: population-scale structural variant comparison and are available at https://doi.org/10.1038/s41586-023-05896-x. analysis. Nat. Methods 20, 408–417 (2023). Nature | Vol 617 | 11 May 2023 | 323 Article 44. Heller, D. & Vingron, M. SVIM: structural variant identification using mapped long reads. 1 Department of Genetics, Yale University School of Medicine, New Haven, CT, USA. 2Center for Bioinformatics 35, 2907–2915 (2019). Genomic Health, Yale University School of Medicine, New Haven, CT, USA. 3Division of 45. Heller, D. & Vingron, M. SVIM-asm: structural variant detection from haploid and diploid Biology and Biomedical Sciences, Washington University School of Medicine, St. Louis, MO, genome assemblies. Bioinformatics 36, 5519–5521 (2020). USA. 4Genomics Institute, University of California, Santa Cruz, CA, USA. 5Institute for Medical 46. Fiddes, I. T. et al. Comparative Annotation Toolkit (CAT)—simultaneous clade and personal Biometry and Bioinformatics, Medical Faculty, Heinrich Heine University, Düsseldorf, genome annotation. Genome Res. 28, 1029–1038 (2018). Germany. 6Center for Digital Medicine, Heinrich Heine University, Düsseldorf, Germany. 47. Wagner, J. et al. Curated variation benchmarks for challenging medically relevant autosomal genes. Nat. Biotechnol. 40, 672–680 (2022). 7 Division of Oncology, Department of Internal Medicine, Washington University School of 48. Chen, X. et al. Cyrius: accurate CYP2D6 genotyping using whole-genome sequencing Medicine, St. Louis, MO, USA. 8Institute of Genetics and Biophysics, National Research data. Pharmacogenomics J. 21, 251–261 (2021). Council, Naples, Italy. 9Department of Data Sciences, Dana-Farber Cancer Institute, Boston, 49. Sirén, J. et al. Pangenomics enables genotyping of known structural variants in 5202 MA, USA. 10Department of Biomedical Informatics, Harvard Medical School, Boston, MA, USA. diverse genomes. Science 374, abg8871 (2021). 11 Department of Genetics, Genomics and Informatics, University of Tennessee Health Science 50. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Center, Memphis, TN, USA. 12McDonnell Genome Institute, Washington University School of Preprint at arXiv https://doi.org/10.48550/arXiv.1303.3997 (2013). Medicine, St. Louis, MO, USA. 13Department of Genetics, Washington University School of 51. Miller, N. A. et al. A 26-hour system of highly sensitive whole genome sequencing for Medicine, St. Louis, MO, USA. 14Novo Nordisk Foundation Center for Biosustainability, emergency management of genetic diseases. Genome Med. 7, 100 (2015). 52. Kolesnikov, A. et al. DeepTrio: variant calling in families using deep learning. Preprint at Technical University of Denmark, Copenhagen, Denmark. 15Quantitative Life Sciences, McGill bioRxiv https://doi.org/10.1101/2021.04.05.438434 (2021). University, Montréal, Québec, Canada. 16Genomics Research Centre, Human Technopole, 53. Vollger, M. R. et al. Increased mutation and gene conversion within human segmental Milan, Italy. 17Department of Genome Sciences, University of Washington School of Medicine, duplications. Nature https://doi.org/10.1038/s41586-023-05895-y (2023). Seattle, WA, USA. 18Quantitative Biology Center (QBiC), University of Tübingen, Tübingen, 54. Ebler, J. et al. Pangenome-based genome inference allows efficient and accurate Germany. 19Biomedical Data Science, Department of Computer Science, University of genotyping across a wide spectrum of variant classes. Nat. Genet. 54, 518–525 (2022). Tübingen, Tübingen, Germany. 20Tree of Life, Wellcome Sanger Institute, Hinxton, Cambridge, 55. Lu, T.-Y., Human Genome Structural Variation Consortium & Chaisson, M. J. P. Profiling UK. 21Northeastern University, Boston, MA, USA. 22Department of Quantitative and variable-number tandem repeat variation across populations using repeat-pangenome Computational Biology, University of Southern California, Los Angeles, CA, USA. 23European graphs. Nat. Commun. 12, 4250 (2021). 56. Benson, G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Res. 27, 573–580 (1999). Campus, Hinxton, Cambridge, UK. 24Coriell Institute for Medical Research, Camden, NJ, USA. 57. Sibbesen, J. A. et al. Haplotype-aware pantranscriptome analyses using spliced 25 Department of Computer Science, University of Pisa, Pisa, Italy. 26Center for Health Data pangenome graphs. Nat. Methods 20, 239–247 (2023). Science, University of Copenhagen, Copenhagen, Denmark. 27Division of Medical Genetics, 58. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013). University of Washington School of Medicine, Seattle, WA, USA. 28Google, Mountain View, CA, 59. Groza, C. et al. Genome graphs detect human polymorphisms in active epigenomic state USA. 29Barrett and O’Connor Washington Center, Arizona State University, Washington, DC, during influenza infection. Preprint at bioRxiv https://doi.org/10.1101/2021.09.29.462206 USA. 30Department of Ecology and Evolutionary Biology, University of California, Santa Cruz, (2022). CA, USA. 31Core Unit Bioinformatics, Medical Faculty, Heinrich Heine University, Düsseldorf, 60. Audano, P. A. et al. Characterizing the major structural variant alleles of the human genome. Cell 176, 663–675.e19 (2019). Germany. 32Vertebrate Genome Laboratory, The Rockefeller University, New York, NY, USA. 61. Falchi, M. et al. Low copy number of the salivary amylase gene predisposes to obesity. 33 National Institutes of Health (NIH)–National Human Genome Research Institute, Bethesda, Nat. Genet. 46, 492–497 (2014). MD, USA. 34Center for Computational and Genomic Medicine, The Children’s Hospital of 62. Mohajeri, K. et al. Interchromosomal core duplicons drive both evolutionary instability Philadelphia, Philadelphia, PA, USA. 35Institute for Society and Genetics, College of Letters and disease susceptibility of the chromosome 8p23.1 region. Genome Res. 26, 1453–1467 and Science, University of California, Los Angeles, CA, USA. 36Institute for Precision Health, (2016). David Geffen School of Medicine, University of California, Los Angeles, CA, USA. 37Division of 63. Fiddes, I. T. et al. Human-specific NOTCH2NL genes affect Notch signaling and cortical General Internal Medicine and Health Services Research, David Geffen School of Medicine, neurogenesis. Cell 173, 1356–1369.e22 (2018). University of California, Los Angeles, CA, USA. 38Department of Biomolecular Engineering, 64. Eizenga, J. M. et al. Efficient dynamic variation graphs. Bioinformatics 36, 5139–5144 (2021). 65. Sirén, J. & Paten, B. GBZ file format for pangenome graphs. Bioinformatics 38, 5012–5018 University of California, Santa Cruz, CA, USA. 39Dovetail Genomics, Scotts Valley, CA, USA. (2022). 40 Division of Oncology, Department of Medicine, Stanford University School of Medicine, 66. Logsdon, G. A. et al. The structure, function and evolution of a complete human Stanford, CA, USA. 41Institute for Genomic Health, Icahn School of Medicine at Mount Sinai, chromosome 8. Nature 593, 101–107 (2021). New York, NY, USA. 42Program in Bioethics and Institute for Human Genetics, University of 67. Groza, C., Kwan, T., Soranzo, N., Pastinen, T. & Bourque, G. Personalized and graph California, San Francisco, CA, USA. 43Genome Biology Unit, European Molecular Biology genomes reveal missing signal in epigenomic data. Genome Biol. 21, 124 (2020). Laboratory, Heidelberg, Germany. 44Genome Informatics Section, Computational and 68. Collins, R. L. et al. A structural variation reference for medical and population genetics. Statistical Genomics Branch, National Human Genome Research Institute, National Institutes Nature 581, 444–451 (2020). of Health, Bethesda, MD, USA. 45Computer Sciences Department, Barcelona Supercomputing 69. Tang, D. et al. Genome evolution and diversity of wild and cultivated potatoes. Nature 606, 535–541 (2022). Center, Barcelona, Spain. 46Departament d’Arquitectura de Computadors i Sistemes 70. Zhou, Y. et al. Graph pangenome captures missing heritability and empowers tomato Operatius, Universitat Autònoma de Barcelona, Barcelona, Spain. 47Material Measurement breeding. Nature 606, 527–534 (2022). Laboratory, National Institute of Standards and Technology, Gaithersburg, MD, USA. 48 Department of Public Health Sciences, University of California, Davis, CA, USA. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 49 Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD, USA. published maps and institutional affiliations. 50 Berlin Institute for Medical Systems Biology, Max Delbrück Center for Molecular Medicine in the Helmholtz Association, Berlin, Germany. 51National Center for Biotechnology Information, Open Access This article is licensed under a Creative Commons Attribution National Library of Medicine, National Institutes of Health, Bethesda, MD, USA. 52Al Jalila 4.0 International License, which permits use, sharing, adaptation, distribution Genomics Center of Excellence, Al Jalila Children’s Specialty Hospital, Dubai, UAE. 53Center and reproduction in any medium or format, as long as you give appropriate for Genomic Discovery, Mohammed Bin Rashid University of Medicine and Health Sciences, credit to the original author(s) and the source, provide a link to the Creative Commons licence, Dubai, UAE. 54Center for Computational Biology, Johns Hopkins University, Baltimore, MD, and indicate if changes were made. The images or other third party material in this article are USA. 55Department of Human Genetics, McGill University, Montréal, Québec, Canada. included in the article’s Creative Commons licence, unless indicated otherwise in a credit line 56 Canadian Center for Computational Genomics, McGill University, Montréal, Québec, to the material. If material is not included in the article’s Creative Commons licence and your Canada. 57Institute for the Advanced Study of Human Biology (WPI-ASHBi), Kyoto University, intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, Kyoto, Japan. 58Howard Hughes Medical Institute, Chevy Chase, MD, USA. 59Laboratory of visit http://creativecommons.org/licenses/by/4.0/. Neurogenetics of Language, The Rockefeller University, New York, NY, USA. 60These authors contributed equally: Wen-Wei Liao, Mobin Asri, Jana Ebler. ✉e-mail: egarris5@uthsc.edu; © The Author(s) 2023 tobias.marschall@hhu.de; ira.hall@yale.edu; hli@jimmy.harvard.edu; bpaten@ucsc.edu 324 | Nature | Vol 617 | 11 May 2023 Methods vials as follows: HiFi, 2 × 107 cells; Nanopore, 5 × 107 cells; Omni-C, 5 × 106 cells; Strand-seq, 1 × 107 cells; 10x Genomics, 4 × 106 cells; and Bionano, Sample selection 4 × 106 cells. Cells for Strand-seq were stored in 65% RPMI-1640, 30% FBS We identified parent–child trios from the 1KG in which the child cell and 5% DMSO and frozen as viable cultures. All other cells were washed line banked within the NHGRI Sample Repository for Human Genetic in PBS and flash-frozen as dry cell pellets. Cells used for ONT-UL produc- Research at the Coriell Institute for Medical Research was listed as tion were separately expanded from the original expansion culture having zero expansions and two or fewer passages, and rank-ordered lot to a bank of five vials of 5 × 106 cells. A single vial was subsequently representative individuals as follows. Loci with MAFs less than 0.05 expanded to a total culture size of 4 × 108 cells, which resulted in a total were removed. MAFs were measured in the full cohort (that is, 2,504 of eight passages. Cells were also reserved for g-banded karyotyping individuals, 26 subpopulations) regardless of each individual’s sub- and Illumina Omni2.5 microarray. population labelling. For each chromosome, principal component analysis (PCA) was performed for dimension reduction. This resulted Karyotyping and microarray. G-banded karyotype analysis was per- in a matrix with 2,200 features, which was then centred and scaled formed on 5 × 106 cells collected at passage five (for HiFi, Nanopore and using smartPCA normalization. The matrix was further reduced to Omni-C) and passage eight (for ONT-UL). For all cell lines, 20 metaphase 100 features through another round of PCA. cells were counted, and a minimum of 5 metaphase cells were analysed We defined the representative individuals of a subpopulation as and karyotyped. Chromosome analysis was performed at a resolution those who are similar to the other members in the group (which, in of 400 bands or greater. A pass/fail criterion was used before cell lines this scenario, is the subpopulation they belong to), as well as different proceeded to sequencing. Cell lines with normal karyotypes (46,XX or from individuals outside the group. Group is defined by previous 1KG 46,XY) or lines with benign polymorphisms that are frequently seen in population labels (for example, ‘Gambian in Western Division’). We apparently healthy individuals were classified as passes. Cell lines were did this as follows. For each sample, we first calculated the intragroup classified as failures if two or more cells harboured the same chromo- distance dintra, which is the average of L2-norms between the sample and somal abnormality. DNA used for microarray was isolated from frozen samples of the same subpopulation. The intergroup distance, dinter, was cell pellets (3 × 106 to 7 × 106 cells) using a Maxwell RSC Cultured Cells similarly defined as the average of L2-norms between the sample and DNA kit on a Maxwell RSC 48 instrument (Promega). DNA was geno- samples from all other subpopulations. The L2-norms were derived in typed at the Children’s Hospital of Philadelphia’s Center for Applied the feature space of the PCA. The score of this sample was then defined Genomics using an Infinium Omni2.5-8 v.1.3 BeadChip (Illumina) on as 10 × dintra + dinter/(n – 1), where n is the number of subpopulations. an iScan System instrument (Illumina). For each subpopulation, if fewer than three trios were available, all were selected. Otherwise, trios were sorted by ranking children with HiFi sequencing. PacBio HiFi sequencing was distributed between max(paternalrank, maternalrank), where paternalrank and maternalrank are two centres: Washington University in St. Louis and the University of the respective ranks of each parent’s score, selecting the three trios Washington. We describe the protocols used at each centre separately. with a maximum value. We ranked by parent scores because during the year 1 effort, the child samples did not have sequencing data and Washington University in St. Louis. High-molecular-weight DNA was therefore had to be represented by the parents. isolated from frozen cell pellets using a Qiagen MagAttract HMW DNA Ideally, we would have selected the same number of candidates from kit and sheared using a Diagenode Megaruptor I to 20 kb mode size. each subpopulation and have an equal number of candidates from both At all steps, DNA quantity was checked on a Qubit Fluorometer I with sexes. To correct for imbalances, we applied the following criteria for a dsDNA HS Assay kit (Thermo Fisher), and sizes were examined on a each subpopulation’s candidate set: (1) when the sex was unbalanced FEMTO Pulse (Agilent Technologies) using a Genomic DNA 165 kb kit. (that is, off by more than one sample), we tried to swap in the next-best SMRTbell libraries were prepared for sequencing according to the candidate of the less represented sex or did nothing if this was not protocol ‘Procedure & Checklist—Preparing HiFi SMRTbell Libraries possible; (2) if a subpopulation had fewer individuals than the desired using the SMRTbell Express Template Prep Kit 2.0’. After SMRTbell gen- sample selection size (that is, all candidates were selected), their unused eration, material was size-selected on a SageELF system (Sage Science) slots were distributed to other unsaturated subpopulations. The latter using the ‘0.75% 1-18 kb’ program (target 3,450 bp in well 12), and some choice is arbitrary but should have little impact on the overall results. combinations of fraction 3 (average size of 15–21 kb), fraction 2 (average The genetic information used in this study was derived from pub- size of 16–27 kb) and fraction 1 (average size of 20–31 kb) were selected licly available cell lines from the NHGRI Sample Repository for Human for sequencing, depending on the empirical size measurements and Genetic Research and the NIGMS Human Genetic Cell Repository at the available mass. The selected library fractions were bound with Sequenc- Coriell Institute for Medical Research. Therefore, this study is exempt ing Primer v.2 and Sequel II Polymerase v.2.0 and sequenced on Sequel II from human research approval as the proposed work involved the col- instruments (PacBio) on SMRT Cells 8M using Sequencing Plate v.2.0, lection or study of data or specimens that are already publicly available. diffusion loading, 2 h of pre-extension and 30 h of movie times. Sam- ples were sequenced to a minimum HiFi data amount of 108.5 Gbp Sequencing (35× estimated genome coverage) on four SMRT Cells. Cell line expansion and banking. Lymphoblastoid cell lines (LCLs) used for sequencing from the 1KG collection (Supplementary Table 1) University of Washington. High-molecular-weight DNA was isolated were obtained from the NHGRI Sample Repository for Human Genetic from frozen cell pellets using a modified Gentra Puregene method and Research at the Coriell Institute for Medical Research. HG002 (GM24385) sheared using gTUBE (Covaris) to 20 kb mode size. At all steps, DNA quan- and HG005 (GM24631) LCLs were obtained from the NIGMS Human tity was checked by fluorometry on a DS-11 FX instrument (DeNovix) with Genetic Cell Repository at the Coriell Institute for Medical Research. All a Qubit dsDNA HS Assay kit (Thermo Fisher), and sizes were examined on expansions for sequencing were derived from the original expansion a FEMTO Pulse (Agilent Technologies) using a Genomic DNA 165 kb kit. culture lot to ensure the lowest possible number of passages and to SMRTbell libraries were prepared for sequencing according to the pro- reduce overall culturing time. Cells used for HiFi, Nanopore, Omni-C, tocol ‘Procedure & Checklist—Preparing HiFi SMRTbell Libraries using Strand-seq, 10x Genomics and Bionano production and for g-banded the SMRTbell Express Template Prep Kit 2.0’. After SMRTbell generation, karyotyping and Illumina Omni2.5 microarray were expanded to a total material was size-selected on a SageELF system (Sage Science) using culture size of 4 × 108 cells, which resulted in a total of five passages after the ‘0.75% 1–18 kb’ program (target 3,400 bp in well 12), and fraction 2 cell line establishment. Cells were split into production-specific sized (average size of 17–20 kb) or fraction 1 (average size of 18–20 kb) was Article selected for sequencing, depending on the empirical size measurements and 3.2 μl 2% Triton X-100) and incubating at room temperature for and available mass. For some samples, the SageELF program ‘0.75% aga- at least 1 h. Tubes were gently tapped halfway through elution. DNA rose, 10 kb–40 kb’ (target 10,000 bp in well 10) was used, and fractions 6 was collected by transferring eluate with a P200 wide-bore pipette to and 7 were pooled together for sequencing (average size of 17–21 kb). The a new 1.5 ml microcentrifuge tube. Some liquid and DNA remained on selected library fractions were bound with Sequencing Primer v.2 and the Nanobind disk after pipetting. The tube containing the Nanobind Sequel II Polymerase v.2.0 and sequenced on Sequel II instruments disk was spun in a centrifuge at 10,000g for 5 s, and any additional liquid (PacBio) on SMRT Cells 8M using Sequencing Plate v.2.0, diffusion that came off the disk was transferred to the eluate tube. This process loading, 3–4 h of pre-extension and 30 h of movie times. Samples were was repeated if necessary until all DNA was removed. The samples sequenced to a minimum HiFi data amount of 96 Gbp (30× estimated were pipette mixed 5 times (approximately 10 s to aspirate and 10 s to genome coverage) on at least four SMRT Cells. dispense for each cycle) with a wide-bore P200 pipette to homogenize the sample. Samples were further allowed to rest at room temperature Comparisons of HiFi production methods. Although subtle differ- overnight to allow DNA to solubilize (disperse). ences in HiFi data production methods existed between the University of Washington and Washington University in St. Louis, the resulting Library preparation data were highly similar, with overlapping assembly statistics from DNA tagmentation and FRA. Circulomics EB+ (EB buffer with 0.02% most samples. These initial genomes were sequenced at a time when Triton X-100) was prepared, and 140.82 μl EB+ was aliquoted into a methods were being refined and optimized for HiFi sequencing, as it 1.5 ml Eppendorf DNA Lo-Bind tube. UHMW DNA (300 μl) from above was a relatively new process. The primary differences in protocols are was aliquoted into the same tube with a wide-bore P200 pipette. The part of the nucleic acid isolation, fragmentation and size selection, mixture was slowly pipetted up and down 3 times with a wide-bore P200 with the downstream sequencing-specific applications being more pipette set to 150 µl. In a separate 1.5 ml Eppendorf DNA Lo-Bind tube, consistent. Both teams were closely engaged with each other as well the following reagents were added in sequence: 144 μl of FRA dilution as with our company associates, including New England Biolabs (NEB), buffer, 9.18 μl of 1 M MgCl2 and 6 μl FRA. The tube was tapped to mix Qiagen, Diagenode and Sage Science, to provide optimal end products. and spun down using a microcentrifuge. The EB–Triton X-100–DNA mixture was added to the FRA dilution buffer–MgCl2–FRA mixture Nanopore ultra-long sequencing protocol. For the 18 additional sam- with a wide-bore P200 pipette. This mixture was then pipette mixed ples, we used the nanopore unsheared long-read sequencing protocol16. 15–20 times with a wide-bore P1000 pipette set to 600 µl. The mix- This generated about 60× coverage of unsheared sequencing from 3 ture appeared homogeneous when pipette mixing was finished. The PromethION flow cells and a N50 value of around 44 kb. For the 29 newly tube was then incubated for 15 min at room temperature. The mixture selected HPRC samples (Results), we used the protocol outlined below. was then pipette mixed 5 times with a wide-bore P1000 pipette set to 600 µl and incubated at room temperature for an additional 15 min. DNA extraction. Around 50 million cells in a pellet were resuspended The mixture was incubated at 30 °C for 1 min, followed by 80 °C for in 200 µl of PBS, and the resuspended cells were aliquoted (40 µl) into 1 min and then held at 4 °C. five 1.5 ml DNA Lo-bind Eppendorf tubes. The following procedure for DNA extraction was completed for each of the five aliquots. Each FRA clean-up. Clean-up used a Nanobind disk. A 5 mm Nanobind disk tube contained sufficient DNA for three libraries loaded onto one flow was added to the above-described reaction mixture followed by 300 µl cell. The following reagents were added in sequence to each tube with of Circulomics buffer NAF10. The tube was gently tapped 10–20 times pipette mixing (10 times up and down) using a P200 wide-bore pipette: to mix. The mixture was placed on a platform rocker at 20 r.p.m. for 40 μl of proteinase K, 40 μl of buffer CS and 40 μl of CLE3. The samples 2 min at room temperature. A DNA ‘cloud’ was visible on the Nanobind were then incubated at room temperature (18–25 °C) for 30 min. Next, disk. The tube was spun for 1–2 s using a benchtop microcentrifuge 40 μl of RNase A was added to each tube with pipette mixing (10 times) and placed on a magnetic rack. The binding solution was removed and with a P200 wide-bore pipette, and samples were incubated at room discarded. The Nanobind disk was washed by adding 350 µl ONT long temperature for 3 min. Two hundred microlitres of BL3 was mixed with fragment buffer (LFB) and gently tapped 5 times to mix. The tube was 200 μl PBS in a 1.5 ml Eppendorf tube. Four hundred microlitres of this spun for 1–2 s using a microcentrifuge and placed on a magnetic rack. BL3–PBS mixture was then added to each sample and the samples mixed The ONT LFB was removed and discarded. Care was taken to not pipette 10 times with a P1000 wide-bore pipette set to 600 μl. DNA attached to the Nanobind disk. This LFB wash was repeated. The Samples were incubated for 10 min at room temperature and then tube was then briefly spun (microcentrifuge) to move the Nanobind pipette mixed 5 times, then incubated at room temperature for 10 min disk to the bottom of the tube. DNA was eluted from the Nanobind and pipette mixed 5 times and then further incubated for 10 min at room disk by adding 125 µl of ONT EB to the tube. The tube was incubated temperature. A white precipitate may form after addition of BL3. This for 30 min at room temperature then gently tapped 5 times (mixing) is normal. A Nanobind disk was added to the cell lysate, then 600 μl of and incubated for an additional 30 min at room temperature. Fluid isopropanol was added. Mixing was performed by inversion of the tube was slowly aspirated 4 times over the Nanobind disk before removing 5 times. Tubes were further mixed on a tube rotator (9 r.p.m. at room the eluate from the tube. The eluate was transferred to a new 1.5 ml temperature for 10 min). The tubes were then placed on a magnetic Eppendorf DNA Lo-Bind tube using a wide-bore P200 pipette. The tube rack, and the Nanobind disk positioned closer to the top of the eluate was then pipette mixed 2 times with a wide-bore P200 pipette. tube to avoid inadvertent removal of the DNA bound to the Nanobind disk. The supernatant was discarded using a pipette and 700 μl of buffer Adapter attachment and rapid adaptor. Rapid adaptor (RAP) was CW1 was added to each tube. The tube in the magnetic rack was then added to the DNA preparation. To 120 µl of eluate (from above), 3 µl inverted 4 times for mixing. A second and third wash with 500 μl of of ONT RAP was added. The mixture was pipette mixed 8 times with buffer CW2 (inversion mix 4 times for each wash) was performed. After a wide-bore P200 pipette. The mixture was then incubated for 15 min the second CW2 wash, liquid was removed from the tube cap and the at room temperature and then again pipette mixed 8 times with a tubes spun on a mini-centrifuge for 2 s, and replaced on the magnetic wide-bore P200 pipette. rack. Residual liquid was removed from the bottom of the tube, taking care not to remove DNA associated with the Nanobind disk. Elution from RAP reaction clean-up with Nanobind. The final library clean-up step the Nanobind disk was accomplished by adding 160 μl Circulomics removes unligated adaptor. In brief, 120 µl Circulomics EB was added elution buffer (EB) plus 0.02% Triton X-100 (comprising 316.8 μl EB to 123 µl of the above-described RAP reaction mixture. The mixture was slowly pipette mixed 3 times with a wide-bore P1000 pipette set to running Hifiasm (https://dockstore.org/organizations/HumanPange- 240 µl. Each aspiration took about 10 s and each dispense took around nome/collections/Hifiasm). All assemblies were generated using this 10 s. A 5 mm Nanobind disk was added to the reaction mixture followed workflow, running on AnVIL71. Cleaning assemblies and fixing some by 120 µl Circulomics buffer NAF10. Mixing was accomplished by gentle structural issues were performed through a combination of automated tapping. The tube was incubated for 5 min at room temperature with- workflows and manual curation as described below. Manual curation out agitation or rotation. The tube was gently tapped 5 times (each was performed using Jupyter notebooks available at GitHub (https:// time for 2–3 times) during the 5 min of incubation. The tube was spun github.com/human-pangenomics/hpp_production_workflows/tree/ for 1–2 s using a microcentrifuge and placed on a magnetic rack. The master/assembly/y1-notebooks). binding solution was discarded. Next 350 µl of ONT LFB was added to the tube and mixed by gentle tapping 5 times. The tube was then spun Filtering adapter-contained reads and running Hifiasm. Before for 1–2 s using a microcentrifuge and placed on a magnetic rack. The producing the assemblies, we detected and removed the reads con- ONT LFB was removed and discarded. Next the Nanobind disk was taining PacBio adapters using a bash script from the HiFiAdapterFilt washed by adding 350 µl ONT LFB. The tube was gently tapped 5 times repository72 (commit 64d1c7b). This script first creates a database of to move LFB over the surface of the disk. The tube was then incubated the PacBio adapter sequences, as illustrated below: at room temperature for 5 min. The tube was then spun for 1–2 s using a microcentrifuge and placed on a magnetic rack. The ONT LFB was >gnl|uv|NGB00972.1:1-45 Pacific Biosciences Blunt Adapter removed and discarded. The tube was briefly spun using a microcen- ATCTCTCTCTTTTCCTCCTCCTCCGTTGTTGTTGTTGAGAGAGAT trifuge to move the Nanobind disk to the bottom of the tube. To elute >gnl|uv|NGB00973.1:1-35 Pacific Biosciences C2 Primer DNA from the Nanobind disk, 126 µl ONT EB was added to the tube. AAAAAAAAAAAAAAAAAATTAACGGAGGAGGAGGA The tube was incubated for 30 min at room temperature, then gently tapped 5–10 times and incubated for an additional 1–2 h at room tem- It then runs blastn with tuned parameters to detect adapter-containing perature. The eluate was then transferred to a new 1.5 ml Eppendorf reads as follows: DNA Lo-Bind tube using a wide-bore P200 pipette using the same tech- nique described above for passing the eluate over the Nanobind disk blastn -db ${DATABASE} -query ${HIFI_FASTA} -task blastn -reward before removing the eluate from the tube. The mixture was then pipette 1 -penalty -5 -gapopen 3 -gapextend 3 -dust no -soft_masking true mixed 2–3 times with a wide-bore P200 pipette. The library was stored -evalue 700 -searchsp 1750000000000 -outfmt overnight at 4 °C before sequencing to facilitate maximal dissolution of DNA. For 43 samples (out of 47), we removed less than 0.15% of the reads; all of the 29 HPRC-selected samples are among these 43 samples, Flow cell loading and sequencing. ONT sequencing buffer (SQB) indicating the low level of adapter contamination in the HiFi data (68 µl) was added to 82 µl of the eluate from above. The mixture was produced by the HPRC. HG005, which is one of the 18 additional sam- pipette mixed 4 times with a wide-bore P200 pipette set to 150 µl. ples, had the highest adaptor contamination percentage, at about 1%. Each aspiration of 150 µl took 10–20 s, and each dispense of 150 µl (Supplementary Fig. 39) took 10–20 s. Samples were then incubated at room temperature for The removed reads were then aligned to the T2T-CHM13 (v.2.0) refer- 10 min. Next the samples were pipette mixed 8 times with a wide-bore ence to ensure that there was no chromosomal or locus-specific bias P200 pipette set to 150 µl. Before loading the library, the flow cell was in the filtering process. Supplementary Fig. 40 shows a snapshot of primed with flush buffer/flush tether mixture per ONT directions. The the IGV browser73 illustrating the coverage of the adapter-containing library was then added to the flow cell. The mixture was viscous but reads along the genome. The locations of the reads were almost evenly loaded smoothly in about 1 min. Some samples took 2 min maximum distributed along the genome and, excluding centromeres, we barely to load. The sequencing run had a re-mux time set for every 6 h. Base found any region covered with more than two adapter-containing calling was performed using Guppy (v.4.0.11), with default parameters reads, even in HG005, which had the highest contamination and the high-accuracy PromethION model (dna_r9.4.1_450bps_hac_ percentage. prom.cfg). The trio-binning mode of Hifiasm needs haplotype-specific k-mers for trio phasing the assembly graph. To generate these k-mers, we Dovetail Omni-C. We prepared Omni-C libraries from each cell line used parental Illumina short reads for the 47 HPRC samples, which are using a Dovetail Omni-C kit (Dovetail Genomics) with the following publicly available from the 1KG dataset19. For each parental short-read modifications as. First, we aliquoted 1 million cells for fixation with sample, we used yak count (v.0.1)74 to generate the k-mer hash tables, formaldehyde and DSG. We digested chromatin with DNAse I until running it once for each of the paternal and maternal read sets: DNA fragments of a desired length were obtained. Per the protocol, we performed end repair on the chromatin, followed by ligation of a yak count -k31 -b37 -o pat.yak paternal.fq.gz biotinylated bridge oligonucleotide, followed by ligation of free chro- yak count -k31 -b37 -o mat.yak maternal.fq.gz matin ends. We reversed crosslinks and purified proximity ligated DNA. We converted the DNA into an Illumina sequencing library using a NEB The adapter-filtered HiFi reads along with the parental k-mer tables Ultra II library preparation kit (NEB) with a Y-adaptor. We enriched for were then given to Trio-Hifiasm (v.0.14) to produce haplotype-resolved ligation products using streptavidin bead capture on the final library. assembly graphs. Only the sample HG002 was re-assembled with Each capture reaction was then split into two replicates before the Trio-Hifiasm (v.0.14.1), which is explained in more detail in the next final PCR enrichment step to preserve complexity. All libraries were subsection. uniquely dual indexed and sequenced on an Illumina Novaseq Platform with read lengths of 2 × 150 bp. hifiasm -o ${SAMPLE_NAME} -t 48 -1 pat.yak -2 mat.yak hifi.fq.gz Phased assembly pipeline Hifiasm produces one graph per haplotype in GFA format. Each We describe the main automated and manual steps taken before, during haplotype-specific GFA file was then converted to FASTA format and after assembly. A combined set of workflow description language using Gfatools75. The assemblies produced by Trio-Hifiasm (v.0.14) are (WDL) formatted assembly workflows is available from Dockstore that released under v.2 after performing the three cleaning steps described captures each of the steps for filtering adapter-contained reads and at the end of this section. Article the bedtools maskfasta command (https://dockstore.org/workflows/ Manually fixing issues. We used paftools.js asmgene, from the mini- github.com/human-pangenomics/hpp_production_workflows/Mask map2 repository (https://github.com/lh3/minimap2/tree/master/ Assembly:master?tab=info). misc)76 to count the number of apparent gene duplications for each of the assemblies produced by Trio-Hifiasm (v.0.14). Asmgene does not bedtools maskfasta \ distinguish between true duplicates and errors. Looking at its results, -fi ${inputFastaFN} \ we were able to find the duplication trend and detect any outlier. This -bed ~{adapterBed} \ assessment acted as a proxy for detecting high-level duplication errors. -fo ~{outputFasta} We used the Ensembl v.99 cDNA sequences77 as the input gene set for running asmgene. In the second cleaning step, we used VecScreen to detect mitochon- drial contigs and contigs consisting of nonhuman sequences from other # Aligning genes to GRCh38 and each Hifiasm haploid assembly organisms, such as bacteria, viruses and fungi. These contigs were then minimap2 -cx splice:hq hs38.fa cdna.fa > hs38.paf dropped from the assemblies using a WDLized version of samtools minimap2 -cx splice:hq ${pat/mat}.fa cdna.fa > ${pat/mat}.paf faidx. It is worth noting that the contigs with nuclear mitochondrial DNAs within them were not dropped. # Detecting gene duplications paftools.js asmgene -a hs38.paf ${pat/mat}.paf samtools faidx \ $inputFastaFN \ Three samples were detected as outliers in terms of the number of `cat contigsToKeep.txt` | gzip \ gene duplications. To identify the cause of this issue, we aligned back the > ~{outputFasta} HiFi reads to those assemblies and checked the depth of coverages and mapping qualities. It showed that the samples HG01358, HG01123 and In the last cleaning step, we selected one contig as the best mito- HG002 contained false duplications in at least one haplotype of length chondrial contig per diploid assembly. To do this, selection we aligned around 55 Mb (in h1tg000058l contig), about 14 Mb (h1tg000013l) and the sequence of the mitochondrial DNA (with the RefSeq identifier around 70 Mb (h2tg000045l), respectively. In the assembly graphs of NC_012920.1) to each diploid assembly using minimap2 with the of HG01358 and HG01123, the duplicated HiFi reads that appeared parameters -cx asm5 --cs. Then we selected one contig with the high- multiple times were used as anchors to manually determine the exact est mapping score and the lowest number of mismatches as the best boundaries of the duplicated regions in the contigs. These two contigs mitochondrial contig (we selected one randomly if multiple best were then manually fixed by breaking the contigs at the duplication contigs existed). This contig was then rotated and flipped (if neces- start and end points and discarding the duplicated sequence from the sary) to match the start and orientation of NC_012920.1.fa and then assembly. In detail, for HG01123 for h1tg000013l, we discarded the added to the maternal assembly of the corresponding sample. Only interval [94439457, contig end]. For HG01358 for h1tg000058l, we the HG01071 sample did not produce any identifiable mitochondrial kept the interval [0, 95732608), renaming the contig to h1tg000058l_1, contig. we discarded the interval [95732608, 150395342) and kept the interval Masked and cleaned mitochondrial assemblies were then acces- [150395342, contig end], renaming it to h1tg000058l_2. To address sioned to GenBank, where they underwent another round of adapter the false duplication in HG002 we re-assembled it using a newer ver- masking and removal of contamination, which was mostly Epstein– sion of Trio-Hifiasm (v.0.14.1), which was reported not to have this Barr virus used to generate the LCLs. The final cleaned assemblies problem. in GenBank were downloaded, and the contig identifiers were We also evaluated the phasing accuracy of the assemblies by pre-pended with the sample name and haplotype integer (where 1 = using yak trioeval (see below). We detected a single large misjoin paternal and 2 = maternal). For example, a contig assigned the in a maternal contig of the HG02080 assembly. It contained an name JAGYVH010000025 in sample HG02257’s maternal assembly approximately 22-Mb-long paternal block in the middle of the con- was renamed to be HG02257#2#JAGYVH010000025. The renamed tig, which resulted in two switch errors at the edges of this block. assemblies were then released to our Amazon Simple Storage Service This block was manually discarded from the assembly and the con- (S3) and Google Cloud Platform (GCP) buckets. In the process of tig was broken into two smaller ones. In detail, in HG02080 for the downloading from GenBank, three of the assemblies (HG00733 h2tg000053l contig, we kept the interval [0, 41506503), renaming it paternal, HG02630 paternal and NA21309 maternal) had their down- to h2tg000053l_1, we discarded the interval [41506503, 63683095) loads prematurely stopped, which resulted in missing sequences. and kept the interval [63683095, contig end], renaming it to Notably, NA21309 is missing its mitochondrial contig. Details can h2tg000053l_2. be found on the HPRC’s year 1 assembly GitHub repository (https:// We finally searched for interchromosomal misjoins using the Mini- github.com/human-pangenomics/HPP_Year1_Assemblies). The graph pangenome (see below for construction details). An ‘interchro- assemblies held at the International Nucleotide Sequence Data- mosomal misjoin’ was defined by a chimeric Minigraph alignment (see base Collaboration are not truncated, but the truncated copies below) consisting of ≥1 Mb subalignments on different chromosomes. were retained in S3 and GCP as they were used in the construction of the pangenomes. Cleaning steps. To clean the raw assemblies, we performed three After submission to GenBank, the assemblies were aligned against additional steps: masking the remaining HiFi adapters, dropping the CHM13 using Winnowmap, and multiple contigs were found to be contigs that were contaminated in their entirety and removing any unmapped. These contigs were subjected to BLAST and found to be redundant mitochondrial contigs. almost exclusively Epstein–Barr virus sequence. GenBank confirmed In the first cleaning step, the sequence of the PacBio SMRTbell (personal communication) that these unmapped contigs should adapter was aligned to each assembly using minimap2 with the param- have been dropped as contamination, but since the genomes were eters -cxsr -f5000 -N2000 --secondary=yes. We extracted only the hits already in active use, they elected not to remove them at this time. with less than or equal to 2 mismatches and which were longer than A list of the contigs that should have been dropped can be found 42 nt. In addition, eukaryotic adapters in each assembly were identified on the year 1 assembly GitHub repository (https://github.com/ by VecScreen78. The combined minimap2 and VecScreen adaptor hits human-pangenomics/HPP_Year1_Assemblies/blob/main/genbank_ (when present) were hard-masked in the assemblies using a WDL of changes/y1_genbank_remaining_potential_contamination.txt). merqury.sh sample.meryl maternal.meryl paternal.meryl Assembly assessment pipeline .merqury Several steps in assembly assessment were managed through a StandardQC workflow written using WDL, run on AnVIL, and avail- GIAB-based assembly quality assessment. As a complementary able at Dockstore (https://dockstore.org/workflows/github.com/ and stratified assessment of assembly quality, we used the GIAB human-pangenomics/hpp_production_workflows/StandardQC). Indi- assembly benchmarking pipeline to compare assembly-based variant vidual tools within the workflow were run in Docker containers with calls to GIAB’s small variant benchmarks (v.4.2.1) for two GIAB samples specific tool versions installed for consistency and reproducibility. assembled in this work: HG002 and HG005. We evaluated the HG002 Details are available within the Dockstore-deposited workflow. The and HG005 HPRC assemblies aligned to GRCh38. Variants were called StandardQC workflow takes short-read data for parental and child from assemblies using Dipcall (v.0.3) (using mimimap2 (v.2.2.4))38. We samples, the two assembly haplotypes, and it produces an analysis over used -z200000,10000 parameter to improve alignment contiguity, as various quality metrics produced by the tools described below. For it was previously shown to improve variant recall in regions with dense each task, the workflow produced a small human-readable summary variation, such as the major histocompatibility complex81. Small vari- file, which is also easy to parse for summarizing steps, as well as the full ant evaluation was performed using hap.py (v.3.15)82, benchmarking output from the tool for manual inspection. Specific tool invocations against v.4.2.1 of high-confidence SNP, small indel and homozygous can be determined from the deposited workflow and are described in reference calls for the GIAB samples HG002 and HG005. Comparisons the subsequent sections. were performed with and without restriction to the associated dip- call region file (dip.bed) to assess recall within and outside assembled Measuring interchromosomal joins. Contigs were aligned to CHM13 regions. For better comparisons of complex variants, hap.py was run (v.2.0) with Minigraph (v.0.18) and processed with the following com- using vcfeval83. Variant calls were stratified using GIAB stratifications mand line: (v.3.0)84, stratifying true-positive, false-positive and false-negative vari- ant calls in challenging and targeted regions of the genome. minigraph -cxasm chm13v2.0.fa contigs.fa | paftools.js misjoin – Trio-based assembly phasing assessment. Assembly phasing was The ‘misjoin’ command reports an interchromosomal join if a contig assessed using Yak and described using two statistics: switch error has two ≥1 Mb alignments to two different chromosomes. and Hamming error rates. Switch error describes the number of times two adjacent phased variants incorrectly switch between maternal and Assembly contiguity assessment. Assembly contiguity was assessed paternal haplotypes. Hamming error rate relates to the total number for each haplotype using QUAST79. These statistics included total of misphased variants per assembled contig. Yak generates phasing sequence assembled, total assembled contigs and contig NG50 (assum- statistics separately for each haplotype using parental k-mers gathered ing a genome size of 3.1 Gb). All reference-based analyses were skipped. from Illumina short-read sequencing of the parents. QUAST was invoked with the following command: Yak generates k-mer databases for the sample and both parental haplotypes (as described above). We used Yak to generate phasing python /opt/quast/quast-5.0.2/quast-lg.py -t 16 -o .quast metrics with the following command: --large --est-ref-size 3100000000 --no-icarus yak trioeval -t 32 paternal.yak maternal.yak > ..yak_phasing.txt Assembly QV assessment. Assembly QV was determined using two separate k-mer-based tools. The first is Yak20. Yak’s QV estimation hap- Hi-C-based assembly phasing assessment. An alternative approach pens separately on each haplotype. The k-mer databases for Yak were for phasing evaluation is to use Hi-C reads that do not require trio infor- generated using the following command: mation. We computed the switch error rate for local phasing evaluation and the Hamming error rate for global phasing evaluation. We imple- yak count -t16 -b37 -o .yak <(cat ) <(cat ) maximum Hi-C read support to detect switch errors on heterozygous positions. In this procedure, we first identified heterozygous k-mers QV estimation using Yak was generated with the following command: (hets) from phased assemblies using 31-mers. Then we mapped Hi-C reads to the assemblies using these 31-mers. If there were >5 reads yak qv -t 32 -p -K 3.2g -l 100k .yak that supported a switch between consecutive hets in assemblies, we > ..yak.qv.txt considered a haplotype switch. For each het pair, we noted whether Hi-C reads supported or did not support the phase. We considered Assembly QV was also determined using Meryl and Merqury80. Meryl a switch error when a het site had a phase switched support relative generates k-mer databases and Merqury determines haplotype QV to that of the previous heterozygous site. The switch error rate is the jointly with both haplotypes. number of local switches divided by the number of heterozygous sites. The k-mer databases with Meryl were generated with the following We performed this operation for the entire contig over all contigs for commands. Databases were generated separately for each read file switch calculations. For the Hamming error calculations, we considered using meryl count and merge with meryl union-sum. Parental-specific Hamming distance on the entire contig level divided by the number of k-mers (hap-mers) were generated using merylu hapmer. heterozygous sites. This measure gives a global view of phasing errors and implicitly penalizes any long switches in contigs. meryl k=21 threads=64 memory=32 count output .meryl Assembly read-based evaluation of Flagger meryl union-sum output .meryl ments to the HPRC assemblies and running Flagger (v.0.1), a read-based bash hampers.sh maternal.meryl paternal.meryl sample.meryl pipeline for evaluating diploid/dual assemblies. All the WDL-based workflows for running these steps are deposited in the Dockstore col- QV estimation using Merqury was generated with the following lection (https://dockstore.org/organizations/HumanPangenome/ command: collections/Flagger-Secphase). Article filter_alt_reads -i ${INPUT_BAM} -o ${ALT_FILTERED_BAM} -f Preparing the HiFi alignments. We aligned back the HiFi reads of each ${ALT_BAM} -v ${SNPS_VCF} sample to its diploid assembly. The alignments were produced with winnowmap (v.2.03) using the following commands: Running the evaluation pipeline. To assess the read mappings result- # making the k-mer table with meryl ing from our diploid alignment process, we used the following five meryl count k=15 output merylDB asm.fa steps, which are combined into a pipeline that we refer to as Flagger. meryl print greater-than distinct=0.9998 merylDB > repetitive_k15.txt Flagger essentially fits a mixture model to successive coverage blocks of the read-to-diploid assembly alignment and then classifies each block # alignment with winnowmap to a category predicting the accuracy of the assembly at that location. winnowmap -W repetitive_k15.txt -ax map-pb -Y -L --eqx --cs -I8g <(cat pat_asm.fa mat_asm.fa) reads.fq.gz | samtools view -hb > read_ Step 1: calculating depth of coverage. First, after producing and alignment.bam cleaning the HiFi alignments, we calculated the depth of coverage for each assembly base by samtools depth -aa (the -aa option allows out- For all samples, we used the full HiFi read sets mentioned in Sup- putting the bases with zero coverage): plementary Table 1, except HG002, for which we downsampled the read set to 35×. samtools depth -aa -Q 0 read_alignment.bam > read_alignment.depth To exclude unreliable alignments, we removed all chimeric alignments and alignments shorter than 2 kb or with a gap-compressed mismatch The output of samtools depth was then converted into a more effi- ratios higher than 1%. As the assembly is diploid and the reads aligned cient format with the .cov suffix. This format is implemented specifi- to the homozygous regions are expected to have low mapping qualities, cally for Flagger and is more efficient, as the consecutive bases with the we did not filter alignments on the basis of their mapping qualities. In same coverage take only one line. We implemented a program called Supplementary Fig. 41, we plot the histograms of mapping qualities depth2cov for converting the output of samtools depth to the .cov and the distributions of alignment identities for one sample, HG00438, format. as an example. The statistics of three sets of alignments were plotted: the alignments to the diploid assembly and to each haploid assembly depth2cov -d read_alignment.depth -f asm.fa.fai -o read_alignment.cov (maternal and paternal) separately. It indicates that the reads have higher identities when the diploid assembly is used as reference but Step 2: fitting the mixture model. In the second step, the frequencies about 20% more reads have mapping qualities lower than 10. of coverages were calculated using cov2counts. The output file with Generally, in highly homozygous regions, the aligner may not be able the .count suffix is a two-column tab-delimited file: the first column to select the correct haplotype as the primary alignment because of shows coverages and the second column shows the frequencies of either read errors or misassemblies. To detect these cases, we searched those coverages. for secondary alignments for which the scores were almost as high as the primary alignment of the same read. For each such read, we made cov2counts -i read_alignment.cov -o read_alignment.counts a pseudo-multiple alignment of the read sequence and the assembly blocks captured by all secondary and primary alignments. Using this The python script fit_gmm.py takes a file .counts suffix, fits a Gauss- alignment, we searched for the read bases that were mismatched in at ian mixture model and finds the best parameters through expectation- least one alignment but not all alignments. We called such bases single maximization. This mixture model consists of four main components nucleotide markers. For each alignment, we calculated a consistency and each component represents a specific type of region: score by considering only the single nucleotide markers and taking the (1) Erroneous component, which is modelled by a Poisson distribution. summation of their base qualities with a negative sign. We then sorted To avoid overfitting, this mode only uses the coverages below 10 so the alignments (regardless of being primary or secondary) based on its mean is limited to be between 0 and 10. It represents the regions this score. If the best alignment was a secondary alignment, we assigned with very low read support. the primary tag to this alignment and removed the other alignments. (2) (Falsely) duplicated component, which is modelled by a Gaussian The percentage of the total reads with swapped alignments ranged distribution, the mean of which is constrained to be half of mean from 0.03% (HG03453) to 0.44% (HG005) across 47 HPRC samples. of the haploid component. It should mainly represent the falsely This result shows that only a small percentage of the reads needed to duplicated regions. be relocalized using this method. This step was performed through (3) Haploid component, which is modelled by a Gaussian distribution. the Secphase (v.0.1) workflow, which is available in the Dockstore col- It represents blocks with the coverages that we expect for the blocks lection (https://dockstore.org/organizations/HumanPangenome/ of an error-free assembly. collections/Flagger-Secphase). (4) Collapsed component, which is actually a set of components each By calling variants, it is possible to detect the regions that either need of which follows a Gaussian distribution, the mean of which is a mul- polishing (that is, are errors) or that have alignments from the wrong tiple of the mean of the haploid component. It represents regions haplotype because of mismappings. We used DeepVariant (v.1.3.0) that have additional copies present in the underlying genome that with the parameter --model_type=“PACBIO” to call variants on these have been ‘collapsed’ into a single copy. alignments. The variants were then filtered to include only the biallelic SNPs with a variant frequency higher than 0.3 and genotype quality It was noted that the model components may change for different higher than 10. regions owing to regional coverage differences and that the resulting systematic differences affect the accuracy of the partitioning process. bcftools view -Ov -f PASS -m2 -M2 -v snps -e ‘FORMAT/VAF < To make the coverage thresholds more sensitive to the local patterns, 0.3 || FORMAT/GQ < 10’ ${OUTPUT_VCF} > ${SNPS_VCF} the diploid assembly was split into windows of length (5–10 Mb) and Having the biallelic SNPs, we found the alignments with alternative a distinct model was fit for each window. Before fitting, we split the alleles and removed them from the bam file. For this step, we imple- whole-genome coverage file produced in the first step into multiple mented and used the program filter_alt_reads, running the following coverage files for each window. We implemented and ran split_cov_ command: by_window for splitting: split_cov_by_window -c read_alignment.cov -f asm.fa.fai -s 5000000 -p ${OUTPUT_PREFIX} Assessing T2T-CHM13 using Flagger. To estimate the false-positive rate of Flagger, we applied it to the T2T-CHM13 (v.1.1) reference. The This produced a list of coverage files, each of which ends with direct output of Flagger showed that about 12.77 Mb (around 0.41%) ${CONTIG_NAME}_${WINDOW_START}_${WINDOW_END}.cov of the T2T-CHM13 reference assembly was flagged as potentially unreliable. The HPRC assemblies were almost free of rDNA arrays, but We then repeated the above-described steps for each resulting there were modelled sequences for rDNA arrays in the T2T-CHM13 coverage file. (v.1.1) reference. These arrays were flagged as falsely duplicated in One important observation is that for short contigs, the coverage their entirety, which indicated that Flagger with HiFi reads may not distribution is generally too noisy to satisfactorily fit the mixture model. be able to correctly evaluate rDNA arrays. Therefore, to make a fair To address this issue, we performed the window-specific coverage comparison, we excluded rDNA arrays (about 9.92 Mb in total) from the analysis only for the contigs longer than 5 Mb. For the shorter contigs, reference evaluation, which decreased the number of unreliable bases we used the results of the whole-genome analysis. to 5.58 Mb (around 0.18%). We additionally identified about 2.76 Mb of a region beside chromosome 1–HSat2 that was mis-flagged as col- Step 3: extracting blocks of each component. Using the fitted lapsed. This mis-flagging was the impact of the systematic coverage model, we assigned each coverage value to one of the four components rise on the neighbouring HSat2 that altered the fitted mixture model. (erroneous, duplicated, haploid and collapsed). To do so for each cover- By manually fixing this mis-flagging, we had about 2.82 Mb (0.09%) of age value, we picked the component with the highest probability. For unreliable blocks in T2T-CHM13 (v.1.1). This number is about 9.3 times example, the coverage value 0 is frequently assigned to the errone- lower than the average for the HPRC assemblies. These unreliable ous component. In Supplementary Fig. 42, the coverage intervals are blocks are mainly a combination of ‘Unknown’ blocks, which could coloured based on their assigned component. not be properly assigned and the regions with HiFi-specific coverage drops. The results of this analysis are available in the HPRC S3 bucket Step 4: incorporating coverage biases in HSats. According to an (https://s3-us-west-2.amazonaws.com/human-pangenomics/index. article describing a complete human genome3, there are some satel- html?prefix=submissions/e9ad8022-1b30-11ec-ab04-0a13c5208311-- lite arrays (for example, HSat1, HSat2 and HSat3) for which the HiFi COVERAGE_ANALYSIS_Y1_GENBANK/FLAGGER/APR_08_2022/FINAL_ coverage drops or increases systematically owing to biases in sample HIFI_BASED/T2T-CHM13/). preparation and sequencing. Such platform-specific biases mislead the pipeline. As a result, the falsely duplicated component may contain Repeat masking a mixture of falsely duplicated and coverage-biased blocks. Similar Repeat masking on each assembly was iteratively performed using effects occur for the collapsed component. RepeatMasker (v.4.1.2-p1). The first step masked used the default To incorporate such coverage biases and to correct the results in human repeat library, and the second step used a repeat library aug- the corresponding regions, we first found the regions of each haploid mented by CHM13 satellite DNA sequences on the original assemblies assembly for which a coverage bias is expected. To find such regions, after hard masking the initial repeat masked DNA. The augmented we lifted over the CHM13 HSat1, HSat2 and HSat3 annotation to each repeat library (final_consensi_gap_nohsat_teucer.embl.txt) is avail- assembly by aligning the assembly contigs to the reference T2T-CHM13 able at Zenodo (https://doi.org/10.5281/zenodo.5537107), and a par- (v.1.1) and GRCh38 (chromosome Y) and projecting the HSat coordi- allelized repeat masking pipeline (RepeatMaskGenome.snakefile) is nates back to the assembly (using python script project_blocks.py). available at GitHub (https://github.com/chaissonlab/segdupanno- Then we ran fit_gmm.py to fit a mixture model for the blocks assigned tation). The union of the two steps generated the complete repeat to each HSat type and adjusted the parameter --coverage, the starting masking. point of the expectation-maximization process, based on the expected coverage in the corresponding HSat. For HSat1, HSat2 and HSat3 we set SD annotation --coverage to 0.75, 1.25 and 1.25 times the average sequencing coverage, SDs were annotated using sedef 85 after masking repeats in each assem- respectively. Finally, we decomposed each HSat based on the inferred bly. Repeats annotated with more than 20 copies corresponded to coverage thresholds and replaced the previous assigned component unannotated mobile elements and were excluded from the analysis. by the new one. The pipeline for annotating SDs is available at GitHub (https://github. com/ChaissonLab/SegDupAnnotation/releases/tag/vHPRC). Step 5: using high-quality alignments to correct spurious flags. In some cases, the duplicated component was mixed up with the haploid SD reliability one. It usually happens when the coverage in the haploid component The reliable and unreliable regions for all haplotype assemblies were drops systematically or the majority of a long contig is falsely dupli- aligned to T2T-CHM13 (v.2.0) and then subdivided into 5 kb windows cated. To address this issue, we used another indicator of a false duplica- and intersected with the SD annotations for T2T-CHM13. SD annotations tion, which is the accumulation of alignments with very low mapping were unavailable for chromosome Y on T2T-CHM13 (v.2.0) at the time quality (MAPQ). We produced another coverage file using only the of analysis. Furthermore, the chromosome Y added to T2T-CHM13 is alignments with MAPQ > 20. Whenever we found a region flagged as from the HPRC HG002 sample; therefore chromosome Y was excluded. duplicated with more than five high-quality alignments, we changed For each class of unreliable region (unknown, erroneous, duplicated, the flag to haploid. collapsed and haploid), we calculated the average number of base After the correction made in step 5, we merged blocks from each pairs overlapping SDs across the haplotype assemblies and annotated component closer than 1,000, and the overlap of any two components each 5 kb window with the most representative overlapping SD (the SD after merging was flagged as ‘unknown’ to show that this block could with the highest product of identity and length). Then using the most not be properly assigned. The BED files produced by Flagger are avail- representative SD, we calculated the average length and identity of SDs able in the HPRC S3 bucket (https://s3-us-west-2.amazonaws.com/ overlapping each class of unreliable region for all the 94 haplotypes and human-pangenomics/index.html?prefix=submissions/e9ad8022- compared the length of identity of SDs that overlapped the different 1b30-11ec-ab04-0a13c5208311--COVERAGE_ANALYSIS_Y1_GENBANK/ types of errors in the assembly. The code for this analysis is available FLAGGER/APR_08_2022/FINAL_HIFI_BASED/FLAGGER_HIFI_ASM_SIM- on GitHub (https://gist.github.com/mrvollger/3bdd2d34f312932c12 PLIFIED_BEDS/). 917a4379a55973). Article For each transcript that mapped to the target genome, the qual- Ensembl mapping pipeline for annotation ity of the mapping was assessed on the basis of aligning the original To create high-confidence annotations, a new Ensembl annotation pipe- reference sequence with the newly identified target sequence. Again, line was developed. The pipeline clusters and maps spatially proximal if the coverage or identity of the aligned sequence was <98% or <99%, genes in parallel (to help avoid issues with individually mapping near respectively, the reference transcript sequence was re-aligned to the identical paralogues) and attempts to resolve inconsistent mappings target region, this time using Exonerate87. Exonerate, although slower by both considering the synteny of the gene neighbourhood in rela- than minimap2, has the ability to handle very small exons and can incor- tion to the GRCh38 annotation and the identity and coverage of the porate CDS data to preserve the CDS (introducing pseudo-introns as underlying mappings. needed). The following command was used: A reference gene set was created from a subset of the GENCODE (v.38) genes29, which was mapped to the HPRC assemblies through a two-pass exonerate -options --model cdna2genome --forwardcoordinates alignment process. This excluded readthrough genes and genes on FALSE --softmasktarget TRUE --exhaustive FALSE --score 500 patches or haplotypes, and only included one copy of the genes on --saturatethreshold 100 --dnawordlen 15 --codonwordlen 15 the X/Y PAR region (only one copy, chromosome X, is modelled in the --dnahspthreshold 60 --bestn 1 --maxintron [max_intron_size] Ensembl representation of the PAR genes). -coverage_by_aligned 1 --querytype dna --targettype [target_type] First, to minimize the difficulty of mapping near identical paralogues, --query [query_file] --target [target_file] --annotation [annotation_file] > a jumping window of 100 kb in length was used to identify clusters of [output_file] genes to map in parallel (Supplementary Fig. 43). The initial window was positioned at the start of the most 5′ gene for each chromosome in the When more than one approach was used to model the transcript, the GRCh38 reference and extended 100 kb from the start of the gene. Any mapping with the highest combined identity and coverage was selected. genes fully or partially overlapping the window were then included in For genes not mapped through the initial regional anchors, a second the cluster. The next 3′ gene that did not overlap the previous window approach was used. The expected location of the gene was located using was then identified, and a new window was created and the process high-confidence genes mapped during the first phase. High-confidence repeated. This resulted in both clustered genes and non-clustered mappings were those for which there was a single mapped copy of genes (genes were considered not clustered when there was only one the gene, all the transcripts had mapping scores of 99% coverage and gene within the window). The regions to map were then identified on identity on average and the gene also had a similar gene neighbour- the basis of the start of the most 5′ gene and the end of the most 3′ gene hood to the neighbourhood in the reference (at least 80% of the of in each cluster (or simply the 5′ and 3′ end of the gene in the case of the same genes in common for the 100 closest neighbouring genes in non-clustered genes). the reference). After this step, the entire genome region underlying For each region defined in the previous step, anchor points were the missing gene, including a 5 kb flanking sequence, was mapped then selected to help map the region on the target genome. Two 10 kb against the target genome using minimap2: anchor points were created 5 kb from the 5′ and 3′ edge of the region, and a central 10 kb anchor was created around the midpoint of the minimap2 --cs --secondary=yes -x map-ont [genome_index] [gene_ region in the GRCh38 genome. The sequences of these anchors were genomic_sequence] > [alignment_file] them mapped against the target genome using minimap2 (ref. 86) with the following command: The resulting hits were then filtered on the basis of overlap with the expected region that the missing gene should lie in. If there was no minimap2 --cs --secondary=yes -N 10 -x map-ont [genome_index] expected region calculated (cases in which no pair of high-confidence [anchor_file] > [alignment_file] genes could be found to define the 5′ and 3′ boundaries of the expected location of the missing gene, for example, at the edge of a scaffold) or The resulting hits were examined to determine high-confidence no hit overlapping the expected region was found, the top reported hit regions in the target genome. High-confidence regions were ones in was used providing it passed an identity cut-off of 99%. The selected which all three anchors were on the same top-level sequence, in col- hit or hits were then extended on the basis of how much of the original inear order, with ≥99% sequence identity and ≥50% hit coverage, and reference gene they covered to ensure that minor local variants between with a similar distance between the anchors when compared with the the reference and target regions did not lead to the target region being reference genome. If no suitable candidate region was found with all truncated. Once extended, the remaining hits were then clustered on three anchors, pairs of mapped anchors were then assessed in a similar the basis of genomic overlap and merged into unique regions. The manner. missing genes were then attempted to be mapped to these regions The sequence-selected region or regions were then retrieved and using an identical process as described above for the initial mappings, aligned against the corresponding GRCh38 region using MAFFT. For involving MAFFT, minimap2 and Exonerate. each gene, the corresponding exons were retrieved and the coordinates To minimize the occurrences of mis-mapped paralogues, each gene were projected through the alignment of the two regions. Transcripts was checked for exon overlap in both the target and the reference. If were then reconstructed from the projected exons. For each transcript, the overlapping genes were not identical at a locus between the ref- the coverage and identity when aligned to the parent transcript from erence and the target, then a conflict was identified. For each gene GRCh38 were calculated. present, filtering was done to reduce or remove the conflict based on If the resulting transcript had either a coverage of <98% or an iden- a number of factors, including whether the genes were in the expected tity of <99%, the parent transcripts were aligned to the target region location, whether the genes were high-confidence mappings, the using minimap2 in splice-aware mode, with the high-quality setting average percent identity and coverage of the transcript for the genes for Iso-Seq/cDNA transcripts enabled. The maximum intron size was and the neighbourhood score. When it was not possible to resolve a set to 100 kb by default. For transcripts with reference introns larger conflict between two genes, both were kept. This concluded the primary than 100 kb, the maximum intron size was scaled and set as 1.5 times mapping process. the length of the longest intron (to allow some variability): After this process, potential recent duplications were identified. To search for recent duplications, the canonical transcript of each gene minimap2 --cs --secondary=no -G [max_intron_size] -ax splice:hq (the longest transcript in the case of noncoding genes, or the transcript -u b [expected_target_region] [transcript_sequences] > [sam_file] with the longest translation followed by the longest overall sequence for protein-coding genes) was selected and aligned across the genome file using bcftools view. The chromosome files for each sample were using minimap2 in a splice-aware manner: combined into one VCF using bcftools concat. minimap2 --cs --secondary=no -G [max_intron_size] -ax splice:hq Gene duplication analysis -u b [genome_index] [input_file] > [sam_file] Duplicated genes were detected as multi-mapped coding sequences using Liftoff88 supplemented by a complementary approach (gb-map) Mappings that had exon overlap with existing annotations from the with multi-mapped gene bodies. The combined set was formed by primary mapping process on the target genome were removed. For including all Liftoff gene duplications and duplicated genes detected new mappings that did not overlap existing annotations, the quality by gb-map. of the alignment was then assessed by aligning the mapped transcript sequence to the corresponding reference transcript to calculate the Liftoff. We ran Liftoff (v.1.6.3) to annotate extra gene copies in each of coverage and per cent identity of the mapping. Different coverage and the assemblies. Liftoff was run with the flag -sc = 0.90 to find additional per cent identity cut-offs were used for these mappings on the basis of copies of genes, with an identity threshold of at least 90%. An example the type of transcript mapped. Protein-coding and small noncoding command is below: transcripts used a coverage and identity cut-off of 95%, whereas long noncoding transcripts used a coverage and identity cut-off of 90%. liftoff -p 10 -sc 0.90 -copies -db -u Pseudogene transcripts had a lower coverage cut-off of 80%, but the -o -polish same identity cut-off of 90% as long noncoding transcripts. When looking for new paralogues, for cases in which multiple canoni- cal transcripts mapped to a locus, a single representative transcript was The additional copies of the genes were identified as such in the selected. This was based on the following hierarchy of gene biotype output gff3 with the field extra_copy_number (equal to anything groups: coding, long noncoding, pseudogene, small noncoding, and other than 0). For this analysis, we also only considered genes that miscellaneous or undefined. were multi-exon, protein-coding genes. The additional gene copies If there were multiple transcripts for the highest represented group, were further filtered to remove any genes outside the ‘reliable’ haploid the transcript with the longest sequence was selected as the repre- regions as determined by the Flagger pipeline. sentative. gb-map. The gene-body mapping pipeline identifies duplicated genes Gene annotation quality analysis by first aligning transcripts of protein-coding and pseudogenes Frameshifts. For the Ensembl and CAT gene annotation sets, we identi- (GENCODE v.38) to each assembly and then multi-mapping the geno­ fied the locations of frameshifting indels by iterating over the coding mic sequences of each corresponding gene. Alignments of at least sequence of each transcript and looking for any gaps in the alignment. 90% identity and 90% of the length of the original duplication were If the gap had a length that was not a multiple of 3, and its length was considered candidate duplicated genes. Candidates were removed <30 bp long (to remove probable introns from consideration), the gap if they overlapped previously mapped transcripts from other genes, was determined to be a frameshift and its location saved to a BED file. low-quality duplications and genes identified through CAT and Liftoff analysis. Nonsense mutations. We also analysed the number of nonsense muta- tions that would cause early stop codons in both the Ensembl and CAT Gene family analysis. To account for gene duplications in high-identity gene annotation sets. We identified nonsense mutations by iterating gene families, gene families were identified on the basis of sequence through each codon in the coding sequence of the predicted transcripts. alignments from gb-map. Genes that mapped reciprocally with 90% If there was an early stop codon before the canonical stop codon at the identity and 90% length were considered a gene family. A single gene end of the transcript, we saved the location in a BED file. was selected as the representative gene for the family, and any gene duplication in the family was counted towards that gene. Validation of mutations using Illumina. For both sets of mutations, we then lifted over the coordinates of the mutations to be on the GRCh38 Pangenome graph construction reference so that we could use existing variant call sets on GRCh38. Minigraph. Minigraph can rapidly perform assembly-to-graph map- We used halLiftover to lift over each set of coordinates, using the pings using a generalization of the minimap2 algorithm34. New SVs of GRCh38-based HAL file from the MC alignment. Then we used bedtools at least 50 bp detected in the mapping can then be added to the graph. intersect to intersect with the variant call file for each of the assemblies. To construct a pangenome graph, one chosen reference assembly, The following sample commands were used: GRCh38 in this case, was used as a starting graph, and the mapping and SV addition steps were repeated for each additional assembly, greedily. halLiftover GRCh38-f1g-90-mc-aug11.hal This iterative approach is analogous to partial order alignment (POA)89. GRCh38 Graphs constructed in this way describe the structural variation within the samples and provide a coordinate system across the reference and bedtools intersect -wo -a -b > will never align a portion of the reference assembly onto another por- tion of the reference assembly. In this way all reference positions have The VCF files used in this intersection were downloaded from the a unique location within the created pangenome. Minigraph (v.0.14) 1KG (https://urldefense.com/v3/__https://ftp.1000genomes.ebi.ac.uk/ was used with -xggs options. The input order was GRCh38, CHM13 then vol1/ftp/data_collections/1000G_2504_high_coverage/working/ the remainder in lexicographic order by sample name. 20201028_3202_raw_GT_with_annot/20201028_CCDG_14151_B01_GRM_ WGS_2020-08-05_chr$i.recalibrated_variants.vcf.gz__;!!NLFGqXoFfo MC. Graphs constructed by Minigraph only contain structural vari- 8MMQ!r6nD4EtteJ7k2BauOrREfgIrlxEI2Upx35sNHiqyI8Did-a6UUUz ation (≥50 bp) by default. The aim of the MC pipeline is to refine the zxGVQkwYkb-bE_rlHQN2Jw2cBdlw7te_-Q$). output of Minigraph to include smaller variants, down to the SNP Where $i was replaced with each chromosome number. From there, level. Doing so allows the graph to comprehensively represent most each chromosome VCF was split so that each sample was in its own variation, as well as to embed the input haplotypes within it as paths, Article which is important for some applications49. To remove noisy alignments (4) Cactus base alignment: Cactus is a tool that uses a graph-based from the MC pangenome, long (≥100 kb) non-reference sequences approach to combine sets of pairwise alignments obtained from identified as being satellite, unassignable to a reference chromosome lastz into a multiple genome alignment37. When aligning different or which appear unaligned to the remainder of the assemblies were species, it uses a phylogenetic tree to progressively decompose the removed from the graph. This resulted in a pangenome with significantly alignment into a subproblem for each ancestral node in the tree. We reduced complexity that nevertheless maintained all sequences of adapted it to also accept chromosome-scale sequence-to-Minigraph the starting reference assembly and the large majority of those in the mappings as produced above, and improved its runtime on align- additional haplotypes. The MC pipeline is composed of the five steps ments of many sequences by replacing its base aligner with abPOA91. described below and in more detail in ref. 35. The script and commands The core algorithm described in ref. 37 remains unchanged, whereby to reproduce this process can be found at GitHub (https://github.com/ the pairwise alignments were used to induce a sequence graph, ComparativeGenomicsToolkit/cactus/blob/81903cb82ae80da342515 then filtered using the Cactus alignment filtering algorithm, and 109cdee5a85b2fde625/doc/pangenome.md#hprc-version-10-graphs). components of unaligned sequence were then processed by the A newer, simpler version of the pipeline that no longer requires satellite base alignment and refinement algorithm. The resulting graph masking can be found at GitHub (https://github.com/ComparativeGe- was used to infer an ancestral sequence (not explicitly used in this nomicsToolkit/cactus/blob/5fed950471f04e9892bb90531e8f63be91 work) and then exported to a hierarchical alignment (HAL) file92. We 1857e1/doc/pangenome.md#hprc-graph). implemented a converter, hal2vg93 that converts the HAL alignment Paths from the reference, GRCh38, are acyclic in the MC graph. Paths into a sequence graph in VG format. from any other haplotypes can contain cycles (as a result of different (5) Post-processing and whole-genome indexing: the following query segments mapping to the same target), but they are relatively rare. post-processing steps were performed on each chromosome (1) Satellite masking: Minigraph is unable to map through highly graph. First, unaligned sequences >10 kb in length, including repetitive sequences such as centromeres and telomeres and, as sequences not aligned to Minigraph, were removed to filter out these regions are also enriched for misassemblies (see the section any under-alignment artefacts that might later be mistaken for ‘Assembly assessment’ of Results), we decided to exclude them insertions. Next, GFAffix94 was used to normalize the graphs by from the MC graphs used in this work. dna-brnn is a tool that uses merging together redundant node prefixes and suffixes. Nodes a recurrent neural network to quickly identify alpha satellites as were flipped as necessary to ensure that reference paths always well as human satellites 1 and 2 (ref. 90). We ran it with its default visit their forward orientations. The chromosomes were combined parameters on all input sequences and cut out any identified regions into a whole-genome graph, indexed and exported to VCF, all using ≥100 kb, except on the reference. The three satellite families that vg. Patched versions of both the GRCh38-based and CHM13-based dna-brnn detects account for the majority of satellite sequences, graphs were created when it was discovered that short contigs but not all. As such, gaps ≥100 kb in Minigraph mappings were also split-mapping to distant locations had induced large deletions. removed. They were detected by mapping each assembly, after The deletions were removed using vg clip -D 10000000 (and the having removed the dna-brnn regions, to the Minigraph (using the pipeline has since been corrected to no longer produce them). procedure described below). Overall, an average of 188.6 Mb of Allele-filtered graphs, used for short-read mapping, were pro- sequence from each (non-reference) assembly was excluded from duced (from the patched graphs) by removing all nodes traversed the graph. by fewer than 9 haplotype paths (minimum AF = 10%) using vg clip (2) Assembly-to-graph mapping: Minigraph generalizes the fast seeding -d 9 -m 10000. The chromosome HAL files were also combined and chaining algorithms of minimap2, but it does not currently pro- into a whole-genome HAL file using halMergeChroms, and clipped duce exact alignments in cigar strings or otherwise. For this work, sequences added back (to facilitate running CAT) using halUnclip. an option, --write-mz, was added to report chains of minimizers, which in this case are 15 bp exact matches, and all assemblies were PGGB. The PGGB uses a symmetric, all-to-all comparison of genomes to mapped to the Minigraph graph using it. The resulting minimizers generate and refine a pangenome. We applied it to build a pangenome were then converted into PAF files with cigars representing exact graph from all genome assemblies and references (both GRCh38 and pairwise alignments between the query contigs and Minigraph node CHM13). The resulting PGGB graph represents all alignment relation- sequences, and all mappings with MAPQ < 5 were excluded, as were ships between input genomes in a single graph. The PGGB graph is a overlapping query regions >10 kb. lossless model of the input assemblies that represents all equivalently. (3) Chromosome decomposition: the Minigraph graphs do not contain This arrangement enables all of our pangenome assemblies to be used inter-chromosomal rearrangements, but the mappings performed as reference systems, a property that we used to explore the scope of in the previous steps can imply them. That is, a contig can partially pangenome variation in a total way. Owing to ambiguous placement map to multiple chromosomes. In most cases, these mappings of variation in all-to-all pairwise alignments, many SV hotspots, includ- involve similarity across different acrocentric short arms. To avoid ing the centromeres, are transitively collapsed into loops through a introducing misleading interchromosomal events, and because it is subgraph representing a single repeat copy, a feature that tends to necessary to run the subsequent steps individually on chromosomes reduce the size of variants found in repetitive sequences. In contrast owing to memory requirements, the mappings were divided by ref- to MC, PGGB does not filter rapidly evolving satellite sequences or the erence chromosome. This was done by splitting the Minigraph into regions that do not reliably align. This increases its size and complex- connected components and using the RN tags to determine their ity relative to the MC graph and adds a significant amount of single- corresponding chromosome names. The PAF mappings were used to ton sequences relative to the Minigraph and MC graphs. However, determine the coverage of each query contig with each chromosome this property enables annotations and coordinates of all assemblies component. This coverage was used to assign each query contig to a in the pangenome to be related to the graph structure and utilized single chromosome by choosing the chromosome with the highest in subsequent downstream analyses. We applied the PGGB model to coverage. Contigs with insufficient coverage to any chromosome investigate the full pangenome and integrate annotations established (<90% for contigs with lengths in [1,10 kb); <80% for [10 kb,100 kb), de novo on the diverse assemblies into a single model for analyses of <75% for [100 kb,1 Mb) and <70% for ≥mb.) were considered ambigu- pangenome diversity and of complex structurally variable loci (MHC ous and not included in the graph. In the GRCh38-based graph, all and 8p inversion). unplaced and random contigs were grouped together into the same PGGB generates a pangenome graph in three phases. (1) Alignment: component. in the first phase, the wfmash aligner95 is used to generate all-vs-all alignments of input sequences. This method, wfmash, applies the map- (v.0.1.0)100 with options -l 0 -r 10000000. Next, variant sites were ping algorithm of MashMap2 to find homologies at a specified length classified into small variant (<50 bp) and SV (≥50 bp) sites. The SV sites and per cent identity. It then derives base-level alignments using a were then annotated as described in the methods section of article that high-order version of the WFA algorithm (wflign), which first aligns describes Minigraph34. In brief, the longest allele sequence of each SV sequences in segments of 256 bp, then patches up the base-level align- site was extracted and stored in the FASTA format. The interspersed ment with local application of WFA. wfmash was designed and devel- repeats, low-complexity regions, exact tandem repeats, centromeric oped specifically for the problem of building all-to-all alignments for satellites and gaps in the longest allele sequences were then identified large pangenomes. (2) Graph induction: the input FASTA sequences using RepeatMasker (v.4.1.2-p1) with the NCBI/RMBLAST (v.2.10.0) and PAF-format alignments produced by wfmash are converted to a search engine and Dfam (v.3.3) database, SDUST (v.0.1)101, ETRF (com- graph (in GFA format) using seqwish96. This losslessly transforms the mit fc059d5)102, dna-brnn (v.0.1)90 and seqtk gap (v.1.3)103, respectively. input alignments and sequences into a graph. (3) Graph normaliza- SDs were identified if the total node length in a site was ≥1,000 bp and tion: we applied a normalization algorithm—smoothxg97—to simplify ≥20% of bases of these nodes were annotated as SD in the reference complex motifs that occur in STRs and other repetitive sequences, or in individual assembly (‘SD annotation’ subsection). To find hits as well as to mitigate underalignment. The graph is first sorted using to the GRCh38 reference genome, minimap2 (v.2.24) with options a path-guided stochastic gradient descent method98 that organizes -cxasm20 -r2k --cs was used to align the longest allele sequences to the graph in one-dimension to optimize path distances and graph the reference genome. Based on the identified features, SV sites were distances. This sort provides a way to partition the graph into smaller classified into various repeat classes using mgutils.js anno (https:// pieces over which we applied a multiple sequence alignment algorithm github.com/lh3/minigraph/blob/master/misc/mgutils.js) with minor (abPOA)91. These pieces were laced back into a final graph. We iter- modifications to enable it to work with the files derived from the MC and ated this process twice using different target POA lengths to remove PGGB graphs. boundary effects caused at the borders of the MSA problems. Finally, we applied GFAffix94 to remove redundant furcations from the topology of Pangenome size and growth. We use the heaps tool of the odgi pange- the graph. nome analysis toolkit98 to estimate how the euchromatic autosomal To build the HPRC PGGB graph, we used both the CHM13 and GRCh38 pangenome grows with each additional genome assembly added. Here references as a target and mapped all contigs against these with wfmash, we approximated euchromatic regions by non-satellite DNA, which was requiring a full length mapping at 90% total identity, collecting all identified by dna-brnn in the construction of the MC graph (see the contigs that mapped to a given chromosome. Contigs that did not map ‘MC’ subsection). Although the MC non-reference haplotypes of the under this arrangement were then partitioned using a split-mapping MC graph do not contain satellite DNA, the PGGB graph does. Conse- approach, requiring 90% identity over 50 kb to seed the mappings, and quently, we subset the PGGB graph to segments contained in the MC putting the contig into the chromosome bin for which it had the best graph. We additionally excluded reference haplotypes (GRCh38 and split mapping. We thus initially partition the data into 25 chromosome CHM13) from the analysis. We then sampled permutations of the 88 sets: one for each autosome, one for each sex chromosome, and finally non-reference (neither GRCh38 nor T2T-CHM13) haplotypes. In each the mitochondria. permutation, we calculated the size of the pangenome after adding We then applied PGGB (v.0.2.0+531f85f) to each partition to build a the first 1, 2, …, N haplotypes in both graphs. This produced a collec- chromosome-specific graph. Run in parallel over 6 PowerEdge R6515 tion of saturation curves from which we derived a median saturation AMD EPYC 7402P 24-core nodes with 384 GB of RAM, this process curve onto which we fitted a power law function known as Heaps’ Law. requires 22.49 system days, or around 3.7 days wallclock. To develop The exponent of this function is generally understood to represent a robust process to build the HPRC graph, the PGGB team iterated the degree of openness—or diversity—of a pangenome39. Summing the build 88 times. The final chromosome graphs were compacted up, we called odgi heaps -i -S -n200 to generate pange- into a single ID space using vg ids -j, then for each reference (GRCh38 nome saturation curves for 200 permutations. Next to calculating a and CHM13) a combined VCF file was generated from the graph with non-permuted cumulative base count, we also counted the number vg deconstruct (v.1.36.0/commit 375cad7). of common (≥5% of all non-reference haplotypes) and core (≥95% of A handful of key parameters defined the shape of the resulting graph. all non-reference haplotypes) bases in the pangenome graphs. To this First, in wfmash, we required >100 kb mappings at 98% identity. We end, we used a tool called panacus104 and supplied a list of the samples mapped each HPRC assembly contig and reference chromosome (both in which they are grouped according to their assigned superpopula- GRCh38 and CHM13) to all the other 89 input haplotypes. To reduce tion (pangenome-growth -m -t bp ). We complexity, and false-positive SNPs resulting from misaligned regions, repeated the count, this time including only segments of depth ≥2, that we applied a minimum match length filter (in seqwish) of 311 bp. This is, contained at least twice in any haplotype sequence. meant that the graph that we induced was relatively ‘underaligned’ locally, and only through normalization in smoothxg did we compress Decomposing pangenome graphs based on allele traversals. Pange- the bubble structures that are produced. For smoothxg, our first iter- nome graphs were decomposed topologically into a set of nested sub- ation attempts to generate 13,033 bp-long POA problems, whereas graphs, termed snarls, that each correspond to one or a collection of the second is 13,177 bp. These lengths provided a balanced trade-off genetic variants. These snarls were then converted to VCF format using between run time and variant detection accuracy. vg deconstruct99. Large (>100 kb) deletions in MC and PGGB graphs In addition to a graph (in GFA), PGGB generates visualizations of were removed using vcfbub (v.0.1.0)100 with options -r 100000. To ease the graph in one and two dimensions, which show both the topology the comparison of variants with other call sets for each individual, the (two dimensional) and path-to-graph relationship (one dimensional). multi-sample VCF files were converted to per-sample VCF files using A code-level description of the build process is provided at GitHub bcftools view -a -I -s , and the multiallelic sites were (https://github.com/pangenome/HPRCyear1v2genbank). split into biallelic records using bcftools norm -m -any. Owing to the limitations of snarl decomposition, snarls may contain multiple vari- Pangenome graph assessment ants that cannot be further decomposed into nested snarls using vg Annotating variant sites in pangenome graphs. Variant sites in Mini- deconstruct. If snarls of this kind are compared with truth calls, the graph and in MC and PGGB graphs were discovered using gfatools evaluation will not be accurate. We solved this problem by compar- bubble (v.0.5)75 and vg deconstruct99, respectively. Large (>10 Mb) ing reference and alternate allele traversals for each snarl to infer the spurious deletions in MC and PGGB graphs were removed using vcfbub minimalist representation of variants (Supplementary Fig. 44). Article For PAV (v.0.9.1)5, assemblies were aligned to the GRCh38 human Annotating small variants in pangenome graphs with AFs from reference genome with no alternatives using minimap2 (v.2.21)86 with gnomAD. Variant sites in MC and PGGB graphs were discovered using options -x asm20 -m 10000 -z 10000,50 -r 50000 --end-bonus=100 vg deconstruct99. The resulting VCF files were then decomposed on the --secondary=no -a --eqx -Y -O 5,56 -E 4,1 -B 5. These alignments are basis of allele traversals (‘Decomposing pangenome graphs based on then trimmed to reduce the redundancy of records and to increase the allele traversals’ subsection). The multi-nucleotide polymorphisms contiguity of alignments. SVs, indels and SNPs were called by using cigar and complex indels were further decomposed into SNPs and simple string parsing of the trimmed alignments. Inversion calling in PAV uses indels using vcfdecompose --break-mnps --break-indels from RTG a new k-mer density assessment to resolve inner and outer breakpoints tools (v.3.12.1)83, so that they could be annotated with gnomAD later. For of flanking repeats, which does not rely on alignment breaks to identify comparison, variants called from PacBio HiFi reads using DeepVariant inversion sites. This is designed to overcome limitations in alignment and from haplotype-resolved assemblies using Dipcall were also used. methodologies and to expand inversion calls, which result in duplica- For each discovery method, small variants (<50 bp) were extracted and tions and deletions of sequence on the boundaries. normalized using bcftools norm -c s -f -m -any. Next, all per-sample VCF files were combined into one github.com/hall-lab/competitive-alignment/blob/master/call_assem- VCF file using bcftools concat -a -D after dropping individual genotype bly_variants.wdl) (commit 830260a). In brief, the maternal and paternal information using bcftools view -G. To annotate small variants with assemblies were aligned to the GRCh38 human reference genome using AFs from gnomAD105, the gnomAD (v.3.1.2) per-chromosome VCF files minimap2 (v.2.1)86 with options -ax asm5 -L --cs. Large indels (>50 bp) were downloaded and concatenated into one VCF file using bcftools were detected using the ‘call_small_variants’ task, based on paftools concat. The VCF file was then compressed into a file in the gnotate (v.2.17-r949-dirty). For large SV, breakpoints were mapped based on format using make-gnotate from slivar (v.0.2.7)106 with options --field split alignments of assembly contigs to the reference genome and AC:gnomad_ac --field AN:gnomad_an --field AF:gnomad_af --field classified as SVs using a series of custom python scripts in the ‘call_sv’ nhomalt:gnomad_nhomalt. The small variants were annotated with task. The breakpoint-mapped SVs were then filtered on the basis of the gnomAD using slivar expr --gnotate . coverage of the reference genome by the assembly contigs (calculated using bedtools genomecov (v.2.28.0)). For each haplotype assem- Variant benchmarking bly, a BED file of ‘excluded regions’ was defined comprising genomic Calling variants from PacBio HiFi reads. The PacBio HiFi reads were regions covered by more than one distinct contig or with more than aligned to the GRCh38 human reference genome with no alternatives 3× coverage by a single contig. Breakpoint-mapped SVs where either using Winnowmap2 (v.2.03)107 with -x map-pb -a -Y -L --eqx --cs. The breakpoint or >50% of the outer span intersected an excluded region MD tags required by Sniffles were calculated using samtools calmd. were filtered. The resulting BAM files were sorted and indexed using SAMtools. For small variants, the two-pass mode of DeepVariant (v.1.1.0)107 with Merging SV call sets. To integrate per-sample VCF files generated by WhatsHap (v.1.1)108 was used to call SNPs and indels from the PacBio three HiFi-based and three assembly-based SV callers, svtools109 was HiFi read alignments. The resulting VCF files were used as truth sets used. For each individual, VCF files from the six callers were jointly for small variant benchmarking. sorted and then merged using svtools lsort and lmerge, first using Three discovery methods were used to call SVs from the PacBio a strict criterion (svtools lmerge -f 20), followed by a more lenient HiFi read alignments. For PBSV (v.2.6.2)41, SV signatures were identi- second merge (svtools lmerge -f 100 -w carrier_wt). The autosomal SV fied using pbsv discover with --tandem-repeats to improve the calling performance in repetitive regions. SVs SV call set for comparison. were then detected using pbsv call with --ccs --preserve-non-acgt -t DEL,INS,INV,DUP,BND -m 40 from the signatures. For SVIM (v.2.0.0)44, Defining confident regions for variant benchmarking. For SVs, con- SVs were called using svim alignment with --read_names --zmws fident regions were generated using Dipcall. Although useful for small --interspersed_duplications_as_insertions --cluster_max_distance variants, current benchmarking tools such as hap.py/vcfeval cannot 0.5 --minimum_depth 4 --min_sv_size 40. In contrast to PBSV and properly compare different representations of small variants in and Sniffles, SVIM outputs all calls no matter their quality. To determine around SVs. Therefore, for each sample, the confident regions from the threshold used for filtering low-quality calls, a precision–recall Dipcall were further processed as follows: curve was generated across various quality scores by comparing with (1) Exclude any SD, self-chain, tandem repeat longer than 10 kb or the GIAB (v.0.6) Tier 1 SV benchmark set for HG002 (Supplementary satellite DNA if there are any breaks in the Dipcall BED file in the Fig. 45). Consequently, SVIM calls with a quality score lower than ten repeat region +15 kb flanking sequence on each side. The rationale were excluded. For Sniffles (v.1.0.12b)42, SVs were discovered with is that breaks in the Dipcall BED file are generally caused by missing -s 4 -l 40 -n -1 --cluster --ccs_reads. Unlike PBSV and SVIM, Sniffles sequence or errors in the assembly or reference or by large SVs or does not generate consensus sequences of insertions from aggregating CNVs for which we do not have tools to benchmark small variants multiple supporting reads. Therefore, Iris (v.1.0.4)43 was used to refine in these regions. the breakpoints and insertion sequences with --hifi --also_deletions (2) Exclude 15 kb around all breaks in the Dipcall BED file for the same --rerunracon --keep_long_variants. All resulting VCF files were sorted reason as noted above. and indexed using BCFtools. (3) Exclude 15 kb around all gaps in GRCh38 because alignments are unreliable. Calling SVs from haplotype-resolved assemblies. Three discovery (4) Exclude variants >49 bp in the Dipcall VCF file and any tandem methods were used to call SVs from the haplotype-resolved assemblies repeats overlapping these SVs +50 bp on each side. generated using Trio-Hifiasm. For SVIM-asm (v.1.0.2)45, assemblies were aligned to the GRCh38 Benchmarking variants. Variant sites in MC and PGGB graphs were human reference genome with no alternatives using minimap2 (v.2.21)86 discovered using vg deconstruct99. Variant sites with alleles larger with -x asm5 -a --eqx --cs and then sorted and indexed using SAM- than 100 kb in MC and PGGB graphs were then removed using vcfbub tools. SVs were called using svim-asm diploid with --query_names (v.0.1.0)100 with options -l 0 -a 100000. The resulting VCF files were --interspersed_duplications_as_insertions --min_sv_size 40. The result- further processed using vcfwave from vcflib110 with option -I 1000. In ing VCF files were sorted and indexed using BCFtools. brief, vcfwave realigned alternate alleles against the reference allele for each variant site using the bidirectional wavefront alignment (BiWFA) appropriate haplotypes. The parameters used were default parameters, algorithm111 to decompose complex alleles into primitive ones. The except as shown below. An example CAT command run is: multi-sample VCF files were then converted to per-sample VCF files luigi --module cat RunCat --hal=CHM13-f1g-90-mc-aug11.hal --ref- using bcftools view -a -I -s and the multiallelic sites genome=GRCh38 --workers=8 --config=cat-hprc.gencode38.auto- were splitted into biallelic records using bcftools norm -m -any. Next somes.config --work-dir work-hprc-gencode38-chm13 --out-dir out- the autosomal small variants (<50 bp) from a given pangenome graph hprc-gencode38-chm13 --local-scheduler --assembly-hub --maxCores (query set) were compared with the HiFi-DeepVariant call set (truth 8 --binary-mode local > cat.hprc.gencode38.autosomes.chm13.log set) using vcfeval from RTG tools (v.3.12.1)83 with options -m annotate --all-records --ref-overlap --no-roc. Note that the multi-nucleotide Comparisons were made between the resulting CAT annotations and polymorphisms and complex indels were reduced to SNPs and sim- those from the Ensembl pipeline by looking at the parent GENCODE ple indels using vcfdecompose --break-mnps --break-indels from identifiers for each gene and transcript in the sets. Numbers of shared RTG tools (v.3.12.1)83 The comparison was performed independently and unique identifiers between the sets were tabulated. Because the for each individual. Recall and precision were calculated within the two annotation sets used slightly different versions of GENCODE (v.38), refined Dipcall confident regions (‘Defining confident regions for only those identifiers attempted to be mapped by both pipelines were variant benchmarking’ subsection) and then stratified using the considered. Additionally, features were considered to be at the same GIAB (v.3.0) genomic context. To evaluate the SV (≥50 bp) calling per- locus if their genomic intervals overlapped. formance, the autosomal SVs from a given pangenome graph (query set) were compared to the consensus SV call set (truth set) for each Identifying medically relevant sites individual using truvari bench (v.3.2.0)112 with options --multimatch SV sites in Minigraph were discovered using gfatools bubble. To -r 1000 -C 1000 -O 0.0 -p 0.0 -P 0.3 -s 50 -S 15 --sizemax 100000 obtain the number of observed alleles per site, per-sample alleles --includebed . Recall and precision were called using minigraph -cxasm --call (Supplementary Table 21). were then stratified using the GIAB (v.3.0) genomic context and by SV sites with alleles larger than 10 kb and at least five observed alleles variant length. were selected as complex SV sites. The complex SV sites were further filtered on the basis of whether they overlapped with medically rel- Alignment of long reads to pangenomes evant protein-coding genes47 using bedtools intersect. To understand PacBio HiFi reads. PacBio HiFi reads from 44 HPRC samples (excluding whether the medically relevant complex SV sites are known in previous the held out samples) were aligned to the MC graph using GraphAligner studies, the coverage of SVs from the 1KG call sets10,19 was computed (v.1.0.13)113 with option -x vg and stored in the GAF format34. For each using bedtools coverage -F 0.1 (Supplementary Table 16). All complex read that aligned to multiple places in the graph, the alignment with SVs were examined using Bandage115 and visually compared to previous the highest score was retained. To remove low-quality alignments, a short-read SV call sets10,19,68,116 using IGV. read with <80% of read length aligned to the graph was discarded. After filtering the read-to-graph alignments, the read depth of each edge was Analysis of five complex loci calculated using vg pack (v.1.33.0)114 with options -Q -1 -D. Note that the Visualization of graph structures of five loci. We extracted sub- resulting GAF files did not contain a mapping quality (encoded as 255 graphs and paths for five loci in the MC and PGGB graphs using gfabase for missing) for each alignment, therefore the option -Q -1 was given to (v.0.6.0)117 and odgi (v.0.6.2)98 with the following example commands: vg pack to ensure that these alignments were used during read-depth gfabase sub GRCh38-f1g-90-mc-aug11.gfab GRCh38. calculation. Next, the edges of each sample were classified into either chr1:25240000-25460000 --range --connected --view --cutpoints on-target or off-target depending on whether they were on the sample 1 --guess-ranges -o RH_locus.walk.gfa paths (encoded as W-lines in MC GFA files) or not. odgi extract -i chr1.pan.smooth.og -o chr1.pan.RH_locus.og -b chr1.RH_locus.bed -E -P Oxford Nanopore reads. ONT reads obtained from 29 HPRC samples (samples labelled HPRC in Supplementary Table 1) were aligned against We then visualized the graph structures of the subgraphs using Band- the MC graph. The alignments were produced using GraphAligner age (v.0.8.1)115. (v.1.0.13) with parameter settings -x vg --multimap-score-fraction 1 --multiseed-DP 1. The number of reads in these datasets range between Alignment of genes to graphs. We aligned Ensembl (release 106)77 1 million and 5.4 million and have an average read length of 28.4 kb. On GRCh38 version gene sequences to the MC graph and PGGB graph using average, 99.68% of the reads received hits from one or more locations in GraphAligner (v.1.0.13)113 with parameter settings -x vg --try-all-seeds the graph. For each read, we determined its best hit based on alignment --multimap-score-fraction 0.1 to identify the gene positions within the score and discarded all its lower-scoring alignments in subsequent graphs. To show locations of genes on Bandage plots, we applied colour analyses. The alignment identities of these best hits peaked above gradients from green to blue to the nodes of each gene. Lines alongside 95%, with an average ratio of alignment-length-to-read-length (ALRL) the Bandage plots showing approximate gene positions, exons and of 0.880 (s.d. = 0.302) and average MAPQ value of 59.35. The alignment transcription start sites based on Ensembl Canonical transcripts were set was further quality-pruned by discarding alignments that either had drawn by hand. an ALRL lower than 0.8 or a MAPQ value lower than 50. The surviving alignments had an overall average ALRL of 0.968 (s.d. = 0.047) and Structural haplotypes identification. Sequences of each assembly are effectuate an overall genome coverage between 10.5-fold and 43-fold represented by paths in a GFA file. We identified SVs in each assembly by across the 29 samples (Supplementary Fig. 46). tracing these paths through different ‘big’ bubbles (>5,000 bp) in either the MC graph or PGGB graph within those gene regions. We selected the Annotating genes within pangenomes 5,000 bp bubble size based on manual inspection of Bandage plots. An We ran CAT46 to annotate each of the genomes within a pangenome example command to identify big bubbles at a RH locus is as follows: graph. CAT projects a reference annotation, in this case GENCODE (v.38), to each of the haplotypes using the underlying alignments within bcftools filter hprc-v1.0-mc-grch38.vcf.gz -r chr1:25240000- the graph. CAT (commit eb2fc87) was run on both of the GRCh38-based 25460000 | grep LV=0 | awk '{OFS=“\t”; print $1,$2,$3,$4,$5}' | tr “,” and the CHM13-based MC graphs . For each graph, the autosomes were “\t” | awk '{OFS=“\t”; for (i=1;i<=NF;i++) {len=length($i); if (len>5000) first run all together, and then the sex chromosomes were run on the {print $1,$2,$3; next}}}' Article To identify gene conversion events (as gene conversions are not For DeepVariant, we trained on the HG002, HG004, HG005, HG006 shown as bubbles in the graphs), we identified nodes that were different and HG007 samples, with HG003 held out. The trained DeepVariant between a gene and its homologous gene (for example, RHD and RHCE) model is available at https://s3-us-west-2.amazonaws.com/human- based on the GraphAligner alignments described above. We refer to pangenomics/index.html?prefix=publications/PANGENOME_2022/ these as paralogous sequence variants. A gene conversion event was DeepVariant/models/DEEPVARIANT_MC_Y1. detected if a path of a gene goes through more than four paralogous For DeepTrio, we trained two sets of models: one on HG002, HG003, sequence variants of its homologous gene in a row. HG004, HG005, HG006 and HG007, with HG001 held out; and one on HG001, HG005, HG006 and HG007, with the HG002, HG003 and Visualization of linear gene structures. We counted the number of HG004 trio held out. Each DeepTrio model set included parental assemblies for each structural haplotype and computed their frequency. and child models. The two trained child deeptrio models are avail- We visualized linear haplotype structures (for example, in Fig. 5c) using able at https://s3-us-west-2.amazonaws.com/human-pangenomics/ gggenes (v.0.4.1)118 based on the structural haplotypes determined for index.html?prefix=publications/PANGENOME_2022/DeepTrio/ each assembly from the pangenome graphs. The length of intervals models/deeptrio/child and https://s3-us-west-2.amazonaws.com/ between genes is fixed (except for TMEM50A and RHCE, because those human-pangenomics/index.html?prefix=publications/PANGENOME_ two genes are immediately next to each other). Lengths of genes are 2022/DeepTrio/models/deeptrio-no-HG002-HG003-HG004/child, shown as proportional to gene lengths in GRCh38. respectively. The two trained parental DeepTrio models are available at https://s3-us-west-2.amazonaws.com/human-pangenomics/ Pangenome point genotyping index.html?prefix=publications/PANGENOME_2022/DeepTrio/ Alignment of reads to the pangenome. The short reads were first models/deeptrio/parent and https://s3-us-west-2.amazonaws.com/ split into chunks to parallelize the read mapping to the ‘allele-filtered human-pangenomics/index.html?prefix=publications/PANGENOME_ graph’ pangenome, as defined in the ‘MC’ subsection. This pangenome 2022/DeepTrio/models/deeptrio-no-HG002-HG003-HG004/parent, is included within the dataset accompanying this paper and can be iden- respectively. tified as ‘clip.d9.m1000.D10M.m1000’. Mapping was performed with Giraffe49 from vg release v.1.37.0. For trio-based runs, the trio-sample Variant calling with DeepVariant. DeepVariant (v.1.3) was evaluated sets of short reads were mapped to the pangenome using Giraffe from on HG003, using the model we trained with HG003 held out (see ‘Model vg release v.1.38.0. Note that the core vg algorithms for Giraffe map- training’). We used the --keep_legacy_allele_counter_behavior flag ping and surjection (conversion from graph space to linear space) (introduced to support this analysis) and a minimum mapping quality are the same in both vg v.1.37.0 and v.1.38.0. The output alignments, of 1 in the make_examples step, before calling the variants with call_ surjected to GRCh38 in BAM format as explained below, are available variants. Both VCFs and gVCFs were produced. The WDL workflow at https://s3-us-west-2.amazonaws.com/human-pangenomics/index. used for single sample mapping and variant calling was deposited into html?prefix=publications/PANGENOME_2022/DeepTrio/samples in the Zenodo (https://doi.org/10.5281/zenodo.6655968). bam directory of each sample’s directory, and are organized by aligner. Variant calling on GRCh38 with BWA-MEM and DeepVariant. Small Surjection to GRCh38 and indel realignment. To perform variant variants were also called using a more traditional pipeline. We aligned calling, GAM alignments were surjected onto the chromosomal paths reads with BWA-MEM50 to GRCh38 with decoys but no ALTs. DeepVariant from GRCh38 (chromosomes 1–22, X and Y) using vg surject and the then called small variants from the aligned reads. The same version and --prune-low-cplx option to prune short and low-complexity anchors parameters were used for DeepVariant. Only the model was changed during realignment. The BAM files were sorted and split by chromo- to the default DeepVariant model. some using SAMtools (v.1.3.1)119. The reads were realigned, first using bamleftalign from FreeBayes (v.1.2.0)120, and then with ABRA (v.2.23)121 Variant calling with DeepTrio. Small variants were also called using on target regions that were identified using RealignerTargetCreator DeepTrio (v.1.3). For HG001, we used the DeepTrio models we trained from GATK (v.3.8.1)122 and expanded by 160 nucleotides with bedtools with HG001 held out (see ‘Model training’). For the HG002, HG003 and slop (v.2.21.0)123. HG004 trio and HG005, HG006 and HG007 trio, we used the models trained with the HG002, HG003 and HG004 trio held out; the HG005, Model training. To perform variant calling with DeepVariant and Deep- HG006 and HG007 trio (except for chromosome 20) was still included Trio, we trained machine-learning models specific to our graph refer- in the training set. We used the --keep_legacy_allele_counter_behavior ence and Giraffe alignment pipeline based on our alignments. For all and a minimum mapping quality of 1 in the make_examples step models, chromosome 20 was entirely held out from all input samples before calling the variants with call_variants. Both VCFs and gVCFs were to provide a control. produced and are available at https://s3-us-west-2.amazonaws.com/ Training was performed on Google’s internal cluster, using unre- human-pangenomics/index.html?prefix=publications/PANGENOME_ leased Google tensor processing unit (TPU) accelerators, from a cold 2022/DeepTrio/samples in the vcf directory of each sample’s directory, start (that is, without using a pre-trained model as input). We believe and are organized by mapping and calling condition. The WDL workflow that nothing about the way in which we executed the training is essential used for trio-based mapping and variant calling was deposited into to the results obtained. Cold start training is estimated to be feasible Zenodo (https://doi.org/10.5281/zenodo.6655962). outside the Google environment; therefore the claims we present here are falsifiable, but it is not expected to be cost-effective. Research- Variant calling on GRCh38 with BWA-MEM, Dragen Graph and ers looking to independently replicate our training should consider DeepTrio. For DeepTrio, small variants were also called using a more doing warm start training from a base model trained on other data, traditional pipeline and a graph-based implementation of Illumina’s using commercially available graphics processing unit (GPU) accelera- Dragen platform (v.3.7.5). The conditions evaluated were each a tors. An example procedure can be found in the DeepVariant training combination of a mapper and a reference. The Giraffe-HPRC condi- tutorial at GitHub (https://github.com/google/deepvariant/blob/r1.3/ tion used Giraffe (v.1.38.0)49 to align reads to the HPRC reference. The docs/deepvariant-training-case-study.md). We predict that this more BWA-MEM condition used BWA-MEM (v.0.7.17-r1188)50 to align reads accessible method would produce equivalent results. to the hs38d1 human reference genome with decoys but no ALTs. For both DeepVariant and DeepTrio, the true variant calls being The Dragen-DeepTrioCall condition used Illumina’s Dragen platform trained against came from the GIAB benchmark (v.4.2.1). (v.3.7.5)51 against their default graph, which was constructed using the same GRCh38 reference with decoys but no ALTs, and population a top-level bubble with a reference length exceeding 100 kb; that is, contigs, SNPs and liftover sequences from datasets internal to their top-level bubbles longer than that are replaced by their child nodes in platform. DeepTrio then called small variants from the aligned reads. the snarl tree. The VCF also contained the haplotypes for all 44 assembly The same version and parameters were used for DeepTrio (v.1.3). Only samples, representing paths in the pangenome graph. We additionally the default model was used for these conditions. We also applied the removed all records for which more than 20% of all 88 haplotypes car- native Dragen caller and joint genotyper to the Dragen-Graph-based ried a missing allele (“.”). This resulted in a set of 22,133,782 bubbles. In a alignments for comparison purposes, referred to as Dragen-DragenCall next step, we used PanGenie (v.1.0.0)54 to genotype these bubbles across and Dragen-DragenJointCall, respectively. Dragen-DragenCall imple- all 3,202 samples from the 1KG based on high-coverage Illumina reads19. ments a single-sample based method and is what is the default use-case for processing Dragen-Graph-mapped data. Dragen-DragenJointCall Decomposition of variants. We genotyped all top-level bubbles across uses a pedigree-backed implementation that informs which variants all 1KG samples. Whereas biallelic bubbles can be easily classified repre- are likely to be de novo and which are erroneous given the genotype senting SNPs, indels or SVs, this becomes more difficult for multiallelic information of the parents. To make a fairer comparison with Dragen, we bubbles contained in the VCF. In particular, larger multiallelic bub- tested these configurations to assess what implementation of Dragen bles can contain a high number of nested variant alleles overlapping variant calling produced the best results given the available trio data. across haplotypes, represented as a single bubble in the graph. This is especially problematic when comparing the genotypes computed for Evaluation using the GIAB benchmark. The small variant calls were the entire bubble to external call sets, as coordinates of the bubble do evaluated using HG001–HG007 with the GIAB benchmark (v.4.2.1)124, not necessarily represent the exact coordinates of individual variant on HG002 in challenging medically relevant autosomal genes47, and alleles carried by a sample in such a region (Supplementary Fig. 20). on HG002 using a preliminary draft assembly-based benchmark. For To tackle this problem, we implemented a decomposition approach the draft assembly-based benchmark, we used Dipcall38 to align a scaf- that aimed to detect all variant alleles nested inside multiallelic top-level folded, high-coverage Trio-Hifiasm assembly21,38 to GRCh38 and call records. The idea was to detect variants from the node traversals of the variants, and then we excluded structurally variant regions from the reference and alternative alleles of all top-level bubbles. Given the dip.bed file as described above for the benchmarking of small variants node traversals of a reference and alternative path through a bubble, from the pangenome graph. The comparison between the call sets and our approach was to match each reference node to its leftmost occur- truth set was made with RTG’s vcfeval83 and Illumina’s hap.py tool82 on rence in the alternative traversal, resulting in an alignment of the node confident regions of the benchmark. We used high-coverage read sets traversals (Supplementary Fig. 21a). Nested alleles could then be deter- of the GIAB HG001, HG002 and HG005 trio child samples and evalu- mined based on indels and mismatches in this alignment. As the node ated performance within the held-out chromosome 20 for the GIAB traversals of the alternative alleles can visit the same node more than (v.4.2.1) truth set, or the entire genome for the reduced truth set of the once (which is not the case for the reference alleles of the MC graph), challenging medically relevant autosomal genes. The evaluation was this approach is not guaranteed to reconstruct the optimal sequence also stratified using the set of regions provided by the GIAB at GitHub alignment underlying the nodes in these repeated regions. (https://github.com/genome-in-a-bottle/genome-stratifications)125. As an output, the decomposition process generated two VCF files. The first one is a multiallelic VCF that contains exactly the same variant Variant calls across samples from the 1KG. We applied our small records as the input VCF, just that annotations for all alternative alleles variant calling pipeline to the high-coverage read sets for the 3,202 of a record were added to the identifier tag in the INFO field. For each samples of the 1KG19. The output alignments, in the GAM format, and alternative allele, the identifier tag contains identifiers encoding all the VCFs were saved in public buckets at https://console.cloud.google. nested variants it is composed of, separated by a colon. The second VCF com/storage/browser/brain-genomics-public/research/cohort/1KGP/ is biallelic and contains a separate record for each nested variant identi- vg/graph_to_grch38. We selected 100 trios among those samples to fier defining reference and alternative allele of the respective variant further evaluate the quality of the calls. We tested all variants that have (Supplementary Fig. 21b). Both VCFs are different representations of at least one alternative allele in a trio for Mendelian consistency. In ad- the same genomic variation, that is, before and after decomposition. We dition, for each variant, we only considered trios for which the child’s applied this decomposition method to the MC-based VCF file, used the genotype was different from the genotype of at least one of the parents multiallelic output VCF as input for PanGenie to genotype bubbles, and to minimize bias created by systematic calls (for example, all homozy- used the biallelic VCF as well as the identifiers to translate PanGenie’s gous or all heterozygous). We looked at the fraction of variants–trios genotypes for bubbles to genotypes for all individual nested variant that failed Mendelian consistency in the entire genome and in sites alleles. All downstream analyses of the genotypes are based on this that do not overlap simple repeats as defined by the ‘simpleRepeat’ biallelic representation (that is, after decomposition). track downloaded from the UCSC Genome Browser. The results were Although the majority of short bubbles (<10 bp) are biallelic, par- compared with Mendelian consistency of calls provided by the 1KG that ticularly large bubbles (>1,000 bp) tend to be multiallelic. Sometimes used GATK HaplotypeCaller on the reads aligned to GRCh38. We also each of the 88 non-reference (neither GRCh38 nor T2T-CHM13) hap- repeated this analysis on the two trios of the GIAB (v.4.2.1) benchmark lotypes contained in the graph covered a different path through such (HG002–HG007) and across the different methods of our evaluation a bubble (Supplementary Fig. 22a), leading to a VCF record with 88 described above (BWA-MEM and DragenGraph mappers; DeepVariant, alternative alleles listed. We determined the number of variant alleles DeepTrio and Dragen variant callers). located inside biallelic and multiallelic bubbles in the pangenome after decomposition. As expected, the majority of SV alleles was located SV genotyping with PanGenie inside of the more complex, multiallelic regions of the pangenome VCF preprocessing. We used a VCF file created on the basis of snarl (Supplementary Fig. 22b). traversal of the MC graph as a basis for genotyping. The records con- tained in this VCF represent bubbles in the underlying pangenome Genotyping evaluation based on assembly samples. We conducted graph and their nested variants, derived from the snarl tree. Each variant a leave-one-out experiment to evaluate PanGenie’s genotyping perfor- was marked according to their level in this tree. Variants annotated by mance for the call set samples. For this purpose, we repeatedly removed ‘LV=0’ correspond to the top-level bubbles. We used vcfbub (v.0.1.0)100 one of the panel samples from the MC VCF and genotyped it using only with parameters -l 0 and -r 100000 to filter the VCF. This removed all the remaining samples as an input panel for PanGenie. We later used non-top-level bubbles from the VCF unless they were nested inside the genotypes of the left-out sample as ground truth for evaluation. Article We repeated this experiment for five of the call set samples (HG00438, defined a filtered set of variants that included the positive set and all HG00733, HG02717, NA20129 and HG03453) using 1KG high-coverage variant alleles with a score of ≥ –0.5. Illumina reads19. PanGenie is a re-genotyping method. Therefore, like We show the number of variant alleles contained in the unfiltered set, every other re-typer, it can only genotype variants contained in the the positive set and the filtered set in Supplementary Table 18. As our input panel VCF, that is, it is not able to detect variants unique to the focus is on SVs and as 65% of all SNPs and indels are already contained genotyped sample. For this reason, we removed all variant alleles (after in the positive set, we applied our machine-learning approach only decomposition) unique to the left-out sample contained in the truth to SVs. We found that 50%, 33% and 26% of all deletion, insertion and set for evaluation. To evaluate the genotype performance, we used other alleles, respectively, were contained in the final, filtered set of the weighted genotype concordance54. Extended Data Fig. 9 shows variants. Note that these numbers take all distinct SV alleles contained the results stratified by different regions. Extended Data Fig. 9a shows in the call sets into account. Especially for insertions and other SVs, concordances in biallelic and multiallelic regions of the MC VCF. The many of these alleles are highly similar, with sometimes only a single biallelic regions include only bubbles with two branches. The multi- base pair differing. Therefore, it is probable that many of these actu- allelic regions include all bubbles in which haplotypes cover more ally represent the same events. Our genotyping and filtering approach than two different paths. Extended Data Fig. 9b shows the same results helped to remove such redundant alleles. stratified by genomic regions defined by GIAB that we obtained from To evaluate the quality of the PanGenie genotypes, we compared the following genotypes: easy (https://ftp-trace.ncbi.nlm.nih.gov/ the AFs observed for the SV alleles across all 2,504 unrelated 1KG sam- ReferenceSamples/giab/release/genome-stratifications/v3.0/GRCh38/ ples to their AFs observed across the 44 assembly samples in the MC union/GRCh38_notinalldifficultregions.bed.gz); low-mappability call set. Supplementary Figs. 23–25 show the results for SV deletions, (https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/ insertions and other SV alleles. We observed that the AFs between genome-stratifications/v3.0/GRCh38/union/GRCh38_alllowmapand- both sets matched well, resulting in Pearson correlations of 0.93, segdupregions.bed.gz); repeats (https://ftp-trace.ncbi.nlm.nih.gov/ 0.87 and 0.81 for deletions, insertions and other alleles, respectively, ReferenceSamples/giab/release/genome-stratifications/v3.0/GRCh38/ contained in the unfiltered set. For the filtered set, we observed cor- LowComplexity/GRCh38_AllTandemRepeats_gt100bp_slop5.bed.gz); relations of 0.96, 0.93 and 0.90, respectively. We also analysed the and other-difficult (https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSam- heterozygosity of the PanGenie genotypes across all 2,504 unrelated ples/giab/release/genome-stratifications/v3.0/GRCh38/OtherDifficult/ 1KG samples and observed a relationship close to what is expected GRCh38_allOtherDifficultregions.bed.gz). by Hardy–Weinberg equilibrium (Supplementary Figs. 23–25, lower Here and in the following, we considered results for SNPs, indels panels). (1–49 bp), SV deletions, SV insertions and other SV alleles, defined as follows: SV deletions include all alleles for which length(REF) ≥50 bp Number of SVs per sample. We compared our filtered set of vari- and length(ALT) = 1 bp; SV insertions include all alleles for which ant alleles to the HGSVC PanGenie genotypes (v.2.0 lenient set)5 and length(REF) = 1 bp and length(ALT) ≥ 50 bp; all other alleles with a Illumina-based SV genotypes19. A direct comparison of the three call length ≥50 bp are included in ‘others’. sets is difficult. The HGSVC and HPRC call sets are based on variant calls Overall, weighted genotype concordances were high for all vari- produced from haplotype-resolved assemblies of 32 and 44 samples, ant types. In particular, variant alleles in biallelic regions of the graph respectively5. For each call set, variants were re-genotyped across all were easily genotypable. Alleles inside multiallelic bubbles were more 3,202 1KG samples. Note that the call set samples for HPRC and HGSVC difficult to genotype correctly as PanGenie needs to decide between are disjoint. As re-genotyping cannot discover new variants, both call several possible alternative paths, whereas there are only two such sets will miss variants carried by 3,202 samples that were not seen in paths for biallelic regions (Extended Data Fig. 9a). Furthermore, geno- the assembly samples. By contrast, the 1KG call set contains short-read typing accuracy depended on the genomic context (Extended Data based variant calls produced for each of the 3,202 1KG samples. Another Fig. 9b). Regions with low mappability, repetitive regions and other difference between the HGSVC and HPRC call sets is that in the HGSVC difficult regions were harder to genotype than regions classified as call set, highly similar alleles are merged into a single record to correct ‘easy’ by GIAB. for representation differences across different samples or haplotypes. The HPRC call set, however, keeps all these alleles separately even if Creating a high-quality subset. We generated genotypes for all 3,202 there is only a single base pair difference between them. To make the 1KG samples with PanGenie and defined a high-quality subset of SV call sets better comparable, we merged clusters of highly similar alleles alleles that we could reliably genotype. For this purpose, we applied in the HPRC filtered set before comparisons with other call sets. This a machine-learning approach similar to what we have previously pre- was done with Truvari (v.3.1.0)112 using the command: truvari collapse sented5,54. We defined positive and negative subsets of variants based -r 500 -p 0.95 -P 0.95 -s 50 -S 100000. on the following filters: ac0_fail, a variant allele was genotyped with an To be able to properly compare the call sets despite their differences, AF of 0.0 across all samples; mendel_fail, the mendelian consistency we counted the number of SV alleles present in each sample (genotype across trios is <80% for a variant allele. Here, we use a strict defini- 0/1 or 1/1) in each call set and plotted the corresponding distributions tion of Mendelian consistency, which excludes all trios with only 0/0, stratified by genome annotations from GIAB (same as above, Fig. 6d). only 0/1 and only 1/1 genotypes; gq_fail, <50 high-quality genotypes We also generated the same plot including only common SV alleles with were reported for this variant allele; self_fail, genotyping accuracy an AF > 5% across all 3,202 samples (Supplementary Fig. 28). Both plots of a variant allele across the panel samples is <90%; nonref_fail, not showed that both assembly based call sets (HPRC and HGSVC) were a single non-0/0 genotype was genotyped correctly across all panel able to access more SVs across the genome than the short-read-based samples. 1KG call set, especially deletions <300 bp and insertions (Fig. 6e). This The positive set included all variant alleles that passed all five filters. result confirms that SV callers based on short reads alone miss a large The negative set contained all variant alleles that passed the ac0_fail proportion of SVs located in regions inaccessible by short-read align- filter but failed at least three of the other filters. We trained a support ments, which has been previously reposed by several studies5,8. In the vector regression approach based on these two sets that used multiple ‘easy’ regions, the number of SVs per sample was consistent across all features, including AFs, Mendelian consistencies or the number of three call sets. For the other regions, however, results indicated that the alternative alleles transmitted from parents to children. We applied HPRC-filtered genotypes provided access to more variant alleles than this method to all remaining variant alleles genotyped with an AF > 0, the HGSVC lenient set, especially insertions and variants in regions of resulting in a score between –1 (bad) and 1 (good) for each. We finally low mappability and tandem repeats (Fig. 6d,e). To evaluate the new SVs in our filtered HPRC call set, we revisited non-overlapping loci. Using the raw VNTR coordinates on GRCh38 the leave-one-out experiment we had previously performed on the as input, VNTR regions across 96 haplotypes (including GRCh38 and unfiltered set of variants (see above). We restricted the evaluation to CHM13) were annotated using the build module in danbing-tk (v.1.3)55 the following subset of variants: (1) those that are in our filtered set (dist_scan=700, dist_merge=1, TRwindow=100000, MBE_th1=0.3, but not in the 1KG Illumina calls (novel); (2) those in our filtered set and MBE_th2=0.6). in the 1KG Illumina call set (known); and (3) all variants in our filtered Whole-genome paired-end error-free short reads were simulated at set. To find matches between our set and the Illumina calls, we used a around 30× for each genome, or equivalently about 15× for each hap- criterion based on reciprocal overlap of at least 50%. Results are shown lotype. A read pair was generated for every 20 bp with fragment size of in Supplementary Fig. 29. We generated two versions of this figure: the 500 bp and read length of 150 bp. Paired-end read mapping to the MC first one (top) excludes variants that are unique to the left-out sample graph was done using Giraffe (v.1.39.0)49 using the command vg giraffe and therefore not typable by any re-genotyping method, and the second -x $pref.xg -g $pref.gg -H $pref.gbwt -m $pref.min -d $pref.dist -p -f one includes these variants (bottom). In general, genotype concord- <(zcat $h1 $h2) -i -t 16, whereas mapping to GRCh38 was done using ances of all lenient variants (brown, dark purple) were slightly higher BWA-MEM (v.0.7.17-r1188)50 using the command bwa mem -t 16 -Y -K compared with the concordances we observed for the unfiltered set 100000000 -p $ref <(zcat $h1 $h2). For a fair comparison, GRCh38 (Supplementary Fig. 29). Furthermore, concordances of the known plus decoy minus ALT/HLA contigs were used as reference to match variants were highest. This is expected, as these variants tended to the paths included in the MC graph. be in regions easier to access by short reads. Concordances for novel variants were slightly worse. This was also expected, as these variants Evaluating read mapping accuracy at VNTR regions. To evaluate tended to be located in more complex genomic regions that are gener- the performance of read mapping using the MC graph plus Giraffe, the ally harder to access. However, even for these variants, concordances VNTR information from danbing-tk were used to annotate each node were still high, which indicated that the PanGenie genotypes for these in the graph by traversing each haplotype path. Every node that covers variants are of high quality. a VNTR region has a tuple that denotes the intersected interval; any aligned reads overlapping with the interval were considered mapped Evaluation based on medically relevant SVs. In addition to all 3,202 to the VNTR. Similarly, a read simulated from an interval overlapping 1KG samples, we genotyped sample HG002 based on Illumina reads with a VNTR was considered derived from the VNTR. To evaluate the per- from ref. 18. We used the GIAB CMRG benchmark containing medically formance of GRCh38 plus BWA-MEM, the mapped region by each read relevant SVs47 downloaded from https://ftp-trace.ncbi.nlm.nih.gov/Ref- was obtained using the bamtobed submodule in BEDTools (v.2.30.0)123. erenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/ The VNTR annotations on GRCh38 were used to determine whether a CMRG_v1.00/GRCh38/StructuralVariant/ for evaluation. Similar to read was mapped to a VNTR. the 1KG samples, we used the MC-based VCF (see above) containing For each read, we tracked its source and mapped VNTR and VNTRs, variant bubbles and haplotypes of 44 assembly samples as an input and used this information to compute accuracy. Only VNTRs pre- panel for PanGenie. We extracted all variant alleles with a length ≥50 bp sent in danbing-tk’s annotations were tracked; otherwise they were from our genotyped VCF (biallelic version, after decomposition). We labelled untracked, the same as non-VNTR regions. A true positive converted the ground truth VCF into a biallelic representation using denotes mapping from a VNTR to its original VNTR. An exogenous bcftools norm -m -any and kept all alleles with length ≥50 bp. We false positive denotes mapping from untracked regions to a VNTR. An used Truvari (v.3.1.0)112 with parameters --multimatch --includebed endogenous false positive denotes mapping from a VNTR to another -r 2000 --no-ref a -C 2000 --passonly VNTR. A true negative denotes mapping from untracked regions to to compare our genotype predictions to the medically relevant SVs. untracked regions. A false negative denotes mapping from a VNTR Results are shown in Supplementary Table 22. As PanGenie is a re-typing to untracked regions (Supplementary Fig. 30). Any alignments in method, it can only genotype variants provided in the input and there- the JSON output of Giraffe that did not contain the mapping field fore cannot detect novel alleles. As HG002 is not among the panel were considered unmapped. The two ends of a read pair that did not samples, the input VCF misses variants unique to this HG002 sample. map to the same chromosome by BWA-MEM were also considered Thus, these unique variants cannot be genotyped by PanGenie and unmapped. will be counted as false negatives during evaluation. Therefore, we computed an adjusted version of the recall that excluded SV alleles Estimating VNTR length variants from read depths. The WGS samples unique to HG002 (that is, alleles not in the graph) from the truth set for 35 genomes (HG00438, HG00621, HG00673, HG00733, HG00735, for evaluation. To identify which SV alleles were unique, we compared HG01071, HG01106, HG01109, HG01175, HG01243, HG01258, HG01361, each of the 44 panel samples to the ground truth VCF using Truvari to HG01891, HG01928, HG01952, HG01978, HG02055, HG02080, HG02145, identify the false negatives for each sample. Then we computed the HG02148, HG02257, HG02572, HG02622, HG02630, HG02717, HG02723, intersection of false-negative calls across all samples. The resulting HG02818, HG02886, HG03098, HG03453, HG03486, HG03492, set then contained all variant alleles unique to the HG002 ground truth HG03579, NA18906 and NA19240) were mapped to the MC graph using set. We found 15 such unique SV alleles among the GIAB CMRG variants. Giraffe as described in the ‘Pangenome point genotyping’ subsection. We removed these alleles from the ground truth set and recomputed Using the VNTR annotations described in the previous section, the precision–recall statistics for our genotypes. Adjusted precision–recall number of reads mapped to each VNTR region in the MC graph was values are shown in Supplementary Table 22. calculated as a proxy for VNTR length. VNTRs with invariant length All genotyping results produced by PanGenie are available at Zenodo across the 35 genomes were removed from analysis, leaving a total of (https://doi.org/10.5281/zenodo.6797328). 60,861 loci. The ground truth for a VNTR in a genome was computed from the number of bases spanned by the VNTR, averaged from the Read mapping at VNTR regions two haplotypes. Simulating and mapping VNTR reads. Raw VNTR coordinates As a baseline control, the read depth of each VNTR region for the 35 on GRCh38 (chromosomes 1–22 and sex chromosomes only) were WGS samples produced by mapping reads to GRCh38 was also com- generated using TRF (v.4.09)56 with command trf hg38.fa 2 7 7 80 puted with mosdepth (v.0.3.1)126 using the command mosdepth -t 4 -b 10 50 500 -f -d. Only repeats with a period size between 6,000 $VNTR_bed -x -f $hg38 $pref $cram. To be able to compare with the and 10,000 bp, total length >100 bp and not overlapping with cen- graph-based approach, VNTRs with missing annotation on GRCh38 tromeric regions were selected, leaving a total number of 98,021 were further removed, leaving a total of 60,386 VNTRs. Article We also used the results of the mapping experiment to quantify allelic RNA-seq mapping evaluation bias. We used vg deconstruct to identify sites of variation in the MC We augmented the allele-filtered graph (see the ‘Pangenome point graph of haplotypes from HG002, with deletions greater than 10 kb genotyping’ subsection) with edges for splice junctions to create a removed to avoid spurious variants. Variants overlapping exons were spliced pangenome graph using the rna subcommand in the vg toolkit selected and normalized using BCFtools (v.1.16)130. Next we counted (v.1.38.0) with a maximum node length set to 32 (vg rna -k 32)57. The the number of mapped reads from each of the two haplotypes that transcript annotations that were used to define the splice junctions overlapped each heterozygous exonic variant57. Specifically, the read consisted of the CAT transcript annotations on each assembly together count for each allele was calculated as the average count across the two with splice junctions from the GENCODE (v.38) annotation29. Tran- breakpoints of an allele. This was done to treat different variant types scripts from the GENCODE (v.38) annotation were further added as and lengths equally. We then tested all variants with a read coverage paths to the spliced pangenome graph. For comparison, we also cre- of at least 20 for allelic bias using a two-sided binomial test at α = 0.01. ated a spliced reference constructed from the reference sequence, All tests that reject the null hypothesis are false positives, as the reads once again using the GENCODE (v.38) transcript annotation. For both were simulated without allelic bias. The results were split into different graphs, we created the indices needed for mapping using the vg toolkit classes of variants and plotted against the number of sites that achieved (v.1.38.0) with default parameters, except when pruning, for which a coverage of at least 20, which is a rough indicator of mapping sensitiv- edges on embedded paths were restored (vg prune -r). Furthermore, ity. Indels larger than 50 bp were excluded. for the spliced HPRC pangenome graph, it was necessary to use stricter We also compared gene expression inference using the mapped pruning parameters (vg prune -r -k 64 -M 64). For the spliced refer- reads. For the vg mpmap (v.1.43.0) graph mappings, we used RPVG ence, we also created the index needed by the RNA-seq mapper STAR (commit 1d91a9e)57 in transcript inference mode to quantify expression. using default parameters. Two different transcript annotations were used as input to RPVG. The We simulated RNA-seq reads with a pipeline that was designed to pre- MC pantranscriptome was created from the CAT transcript annota- serve complex genome variation in the simulated data. The transcript tions on each assembly and the GENCODE (v.38) transcriptome. The sequences used for the simulation were derived from the GENCODE mpmap-RPVG pipeline was compared with Salmon (v.1.9.0)131 and RSEM (v.38) transcript annotations projected onto assembled haplotypes (v.1.3.3)132. These methods were provided the GENCODE (v.38) tran- from HG002. Specifically, we used MC to create an alignment between scriptome as input, both with and without transcripts on the GRCh38 GRCh38 and the two HPRC HG002 assembly haplotypes, which were alternative contigs. Any genes unique to the alternative contigs were held out of the main pangenome graph for benchmarking. We then filtered. For Salmon, the GRCh38 reference was used as a decoy and used CAT to lift the transcript annotations over to these haplotypes. duplicate transcripts were kept. Bowtie2 (v.2.4.5)133 was used as a map- We constructed a spliced personal genome graph using the vg rna sub- per for RSEM. Gene expression values were calculated by summing the command, and then we simulated reads using vg sim (commit 2cea1e2) corresponding transcript level values for each gene. GffRead (v.0.12.7) using an Illumina NovaSeq cDNA read set (SRR18109271) to fit model was used to create a table of gene names and transcript identifiers from parameters. This essentially amounts to simulating directly from the the transcript annotation134. The accuracy was measured in two ways: (1) projected transcript sequences. The transcripts were simulated with the Spearman correlation between the simulated and inferred expres- uniform expression, split evenly between the two haplotypes, keeping sion values, and (2) the mean absolute relative difference between the the reads from each haplotype separate. This expression profile is not simulated and inferred expression values. biologically realistic, but it avoids the difficulty of choosing a particular The scripts that were used for graph construction, read simulation, expression profile as representative for all tissues and stages of develop- mapping and evaluation are available at GitHub (https://github.com/ ment. Moreover, existing estimated profiles would be biased towards jonassibbesen/hprc-rnaseq-analyses-scripts). the tools that were used to estimate them. We simulated 5,000,000 paired-end 150 bp RNA-seq reads. ChIP-seq analysis Both simulated and real Illumina RNA-seq reads were mapped to the We aligned H3K4me1, H3K27ac and ATAC-seq obtained from monocyte- graphs using vg mpmap (commit c0c4816) with default parameters. In derived macrophages from 30 individuals59 using vg map13 to the addition, the reads were mapped to the spliced reference using STAR GRCh38 reference genome graph and to the HPRC genome graph. Then (v.2.7.10a) with default parameters58. For the real data, we used previ- we called peaks using Graph Peak Caller (v.1.2.3)135 on both sets of align- ously published data NA12878 RNA-seq data (SRR1153470)127 and the ments for each of the 30 H3K4me1, H3K27ac and ATAC-seq samples. ENCODE project (ENCSR000AED, replicate 1)128,129. To identify HPRC-only peaks, we projected HPRC coordinates to the We used the same approach as previously described57 to evaluate the GRCh38 path using Graph Peak Caller and compared intervals using alignments. In brief, for the simulated data, the graph alignments were BEDTools123. We named HPRC peaks that overlapped GRCh38 peaks as compared with the truth alignments by estimating their overlap on the common peaks and those that do not as HPRC-only. We calculated the reference genome paths. The graph alignments were projected to the expected frequency (as inverse cumulative distributions) of common reference paths using vg surject -S. A uniquely aligned read (one with and HPRC-only peaks among the 30 samples by resampling the peaks primary MAPQ ≥ 30) was considered correct if it overlapped 90% of the of each sample from the peaks of all the samples and re-counted the truth alignment. A multi-mapped read was considered correct if any of number of overlaps. We repeated this simulation 100 times and plot- the multi-mappings was correct under the same criterion. For the real ted the average curves. We determined heterozygous variants in our data, the average read coverage of each exon on the reference path samples by aligning WGS datasets for each sample to the HPRC graph calculated from the projected graph alignments were compared with using vg map and genotyping the variants using vg call -a. We narrowed the corresponding coverages estimated from long-read alignments. For the list of heterozygous SVs above 50 bp in each sample with the aim the long-read data, we used PacBio Iso-Seq alignments from the ENCODE of looking for allelic-specific peaks. For each epigenomic sample, we project (ENCSR706ANY, all replicates), which come from the same cell obtained allelic-specific read counts within peaks that lie on the pre- line as the Illumina data. The long-read alignments were used to define viously identified loci by running vg call -a on the epigenomic HPRC the exons, and only primary long-read alignments with a mapping qual- alignments, which outputs the numbers of reads on each path in a ity of at least 30 were used. The alignments for the four Iso-Seq replicates bubble (DP and AD fields in the VCF output). We then assigned peaks (ENCFF247TLH, ENCFF431IOE, ENCFF520MMC and ENCFF626GWM) to the SV or reference allele, or both alleles with a two-tailed binomial were combined and filtered using SAMtools (v.1.15)130. BEDTools test parameterized on the sum of reads on both alleles and P = 0.05. Any (v.2.30.0) was used to convert the alignments to exons coordinates123. peak with reads on one allele, but not on the other was assigned to the allele with the reads. Read counts were proportionally adjusted for the the differences between cluster count distributions in metacentric difference in length between the reference and SV alleles. versus acrocentric chromosomes across the entire chromosome, the Processed data, scripts and code for the above steps are available at q-arm and the p-arm. There were non-significant differences between Zenodo (https://doi.org/10.5281/zenodo.6564396). the distributions between acrocentric and metacentric chromosomes at a chromosome scale and in the q-arms, but a significant difference Population genetic analyses (Wilcoxon P = 0.013) in the case of acrocentric p-arms (Supplementary Although the size of the population sample represented in our pange- Fig. 49). nome is small, it provides access to previously under-ascertained This analysis indicates that significant challenges remain for the use regions of the genome. We sought to understand the potential utility of these new regions in population genetic studies. Patterns observed of these regions for future population genetic studies using regional in PCA projections of the pangenome across all chromosomes suggests PCA based on variants called compared to the CHM13 and GRCh38 a distinct process of variation sharing between populations within the references. For these analyses, we considered both the PGGB (whole short arms of the acrocentrics. In effect, we observed a more homog- pangenome, combined) and MC (reference-based, distinct CHM13 enous population in these regions when using the CHM13 assembly as and GRCh38) graphs. For both graph models, the CHM13 VCFs provide a reference. This reference contains real sequences in these regions, access to regions that were not previously observed by studies based whereas GRCh38 contains gaps, which render analysis impossible. The on GRCh38, for which short-read-based studies may have difficulty apparent population homogenization could be driven by error. We reliably aligning and calling variants. In combination, these two graphs mitigated this issue by utilizing only SNPs found in Flagger-confident provide cross-validation of implied population genetic patterns in regions, but this does not guard against potential sources of alignment these new regions, which we explore here. error that are likely to be amplified by the repetitive sequences in these To understand chromosome-specific patterns of variation, we loci. It is also possible that the chromosome-specific partitioning pro- applied PCA to each autosomal chromosome independently to the cess applied by both graph models is failing to correctly partition con- VCFs from PGGB (PGGB-CHM13, PGGB-GRCh38). tigs on these short arms. The known homology between the short arms To ensure that observed patterns were not derived from higher rates bolsters the possibility of ongoing sequence information exchange of assembly error in the repetitive regions of acrocentric p-arms, we between non-homologous chromosomes3, which would be consistent used our Flagger-confident region annotations to prune the PGGB with the patterns we observed. In summary, this analysis shows that graph (using odgi inject to inject the confident regions as subpaths when using CHM13 as a reference, the behaviour of sequences on the and then odgi prune to remove the full original paths that were includ- short arms of the acrocentrics in the PGGB graph is not similar to that ing unreliable regions) to only confident regions of assemblies. We of other sequences in the pangenome. then reapplied vg deconstruct to this graph to obtain a new set of SNPs (the code for the PGGB graphs pruning and variant calling on Reporting summary the pruned graphs can be found at the following link: https://github. Further information on research design is available in the Nature Port- com/pangenome/HPRCyear1v2genbank/blob/main/workflows/confi- folio Reporting Summary linked to this article. dent_variants.md). Genome-wide, we found that pruning reduced the number of called SNPs by only 1.188% (previous N = 23,272,652, pruned N = 22,996,113). The total reduction in the acrocentrics was higher, Data availability with 6.29% fewer SNPs (previous N = 3,735,605, pruned N = 3,676,746), Sequencing data, assemblies and pangenomes produced by the HPRC which indicated the difficulty in assembling these regions. We note are available at AnVIL (https://anvilproject.org) in the AnVIL_HPRC that the PCA sample distributions remained almost identical (data not workspace. Data are also available as part of the AWS Open Data Pro- shown), which indicated that the patterns observed in the full graph gram (https://registry.opendata.aws) in the human-pangenomics S3 are maintained despite assembly issues. In these filtered PGGB-CHM13 bucket (https://s3-us-west-2.amazonaws.com/human-pangenomics/ and PGGB-GRCh38 VCFs, we considered all biallelic SNPs relative to the index.html). In addition, data have been uploaded to the International chosen reference, regardless of variant nesting level (data not shown; Nucleotide Sequence Database Collaboration for long-term storage filtering for only SNPs LV = 0 or LV > 0 produced nearly identical results). and availability. Supporting information about the data (including A qualitative evaluation suggested no significant differences in PCA index files with S3 and GCP file locations) can be found in our GitHub patterns across the metacentric chromosomes (Supplementary Fig. 47). repositories (see below). Sequencing data for 29 selected HPRC samples However, in the p-arms of the acrocentrics (chromosomes 13, 14, 15, 21 from the 1KG cohort (Results) are uploaded to BioProject PRJNA701308 and 22), which are accessible in the PGGB-CHM13 VCF, we observed a (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA701308). Sequencing reduction in population differentiation and a higher rate of variance data created by the HPRC for samples in the cohort of 18 additional sam- explained in the lowest principal component. ples are uploaded to BioProject PRJNA731524 (https://www.ncbi.nlm. To investigate this quantitatively, we measured the number of nih.gov/bioproject/PRJNA731524). Both sets of assemblies are grouped clusters implied by the PCA for the PGGB-CHM13 VCFs using k-means by an Umbrella BioProject PRJNA730822. Data used in this paper have clustering to automatically determine the optimal number of clusters additional information available at AnVIL and at GitHub (https://github. for each PCA (gap_stat clustering in the fviz_nbclust function of the com/human-pangenomics/HPP_Year1_Data_Freeze_v1.0). Assem- factoextra R package) (analysis code available at: https://github.com/ blies along with assembly annotations, such as RepeatMasker and SilviaBuonaiuto/hprcPopGenAnalysis). Applying this approach to three Ensembl gene annotations, can be viewed in an assembly hub in the PCAs per chromosome VCF, we obtained optimal cluster counts for the UCSC Genome Browser (http://hprc-browser.ucsc.edu). The data and p-arm, q-arm and entire chromosome. In metacentric chromosomes, annotations can also be accessed through the Ensembl Rapid Release we usually observed optimal numbers of clusters approximately cor- Genome Browser (https://rapid.ensembl.org) and a dedicated Ensembl responding to the number of expected world population groupings in project page for centralized access to HPRC data (https://projects. the input genomes (3–5, as in Supplementary Fig. 48). However, in the ensembl.org/hprc), which includes links to download the data locally. p-arms of the acrocentrics, we observed much fewer, in general only File locations for the assemblies (and select annotation files) that are one cluster, indicative of reduced population differentiation compared stored in S3 and GCP can also be found at AnVIL or GitHub (https:// with other parts of the acrocentric chromosomes. This pattern was only github.com/human-pangenomics/HPP_Year1_Assemblies). Pange- apparent in the PGGB graph based on CHM13. To evaluate the differ- nomes were uploaded to the European Nucleotide Archive as analysis ence quantitatively, we applied a Wilcoxon rank-sum test to compare objects and are organized under Umbrella BioProject PRJNA850430 Article (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA850430). File loca- 102. Li, H. ETRF: exact tandem repeat finder. GitHub https://github.com/lh3/etrf (2019). 103. Li, H. Seqtk: a toolkit for processing sequences in FASTA/Q formats. GitHub https:// tions for pangenomes, indices and variant calls derived from the github.com/lh3/seqtk (2018). pangenomes can be found at AnVIL or GitHub (https://github.com/ 104. Doerr, D. Panacus: calculate growth statistics for pangenome graphs. GitHub https:// human-pangenomics/hpp_pangenome_resources). Variant calls pro- github.com/marschall-lab/panacus (2022). 105. Karczewski, K. J. et al. The mutational constraint spectrum quantified from variation in duced by the analysis performed in this paper are available at Amazon 141,456 humans. Nature 581, 434–443 (2020). S3 (https://s3-us-west-2.amazonaws.com/human-pangenomics/index. 106. Pedersen, B. S. et al. Effective variant filtering and expected candidate variant yield in html?prefix=publications/PANGENOME_2022) or where indicated in studies of rare human disease. NPJ Genom. Med. 6, 60 (2021). 107. Jain, C., Rhie, A., Hansen, N. F., Koren, S. & Phillippy, A. M. Long-read mapping to repetitive the relevant sections below. PanGenie genotypes produced for the 1KG reference sequences using Winnowmap2. Nat. Methods 19, 705–710 (2022). samples based on the MC graph are available at Zenodo (https://doi. 108. Martin, M. et al. WhatsHap: fast and accurate read-based phasing. Preprint at bioRxiv org/10.5281/zenodo.6797328). https://doi.org/10.1101/085050 (2016). 109. Larson, D. E. et al. svtools: population-scale analysis of structural variation. Bioinformatics 35, 4782–4787 (2019). 110. Garrison, E., Kronenberg, Z. N., Dawson, E. T., Pedersen, B. S. & Prins, P. A spectrum of free Code availability software tools for processing the VCF variant call format: vcflib, bio-vcf, cyvcf2, hts-nim and slivar. PLoS Comput. Biol. 18, e1009123 (2022). The list of all tools, including versions and/or code commits, used for 111. Marco-Sola, S. et al. Optimal gap-affine alignment in O(s) space. Bioinformatics 39, this study are available in Supplementary Table 20. btad074 (2023). 112. English, A. C., Menon, V. K., Gibbs, R. A., Metcalf, G. A. & Sedlazeck, F. J. Truvari: refined structural variant comparison preserves allelic diversity. Genome Biol. 23, 271 (2022). 71. Schatz, M. C. et al. Inverting the model of genomics data sharing with the NHGRI Genomic 113. Rautiainen, M. & Marschall, T. GraphAligner: rapid and versatile sequence-to-graph Data Science Analysis, Visualization, and Informatics Lab-space. Cell Genom. 2, 100085 alignment. Genome Biol. 21, 253 (2020). (2022). 114. Hickey, G. et al. Genotyping structural variants in pangenome graphs using the vg toolkit. 72. Sim, S. HiFiAdapterFilt: remove CCS reads with remnant PacBio adapter sequences and Genome Biol. 21, 35 (2020). convert outputs to a compressed .fastq (.fastq.gz). GitHub https://github.com/sheinasim/ 115. Wick, R. R., Schultz, M. B., Zobel, J. & Holt, K. E. Bandage: interactive visualization of de HiFiAdapterFilt (2021). novo genome assemblies. Bioinformatics 31, 3350–3352 (2015). 73. Robinson, J. T., Thorvaldsdóttir, H., Wenger, A. M., Zehir, A. & Mesirov, J. P. Variant review 116. Abel, H. J. et al. Mapping and characterization of structural variation in 17,795 human with the Integrative Genomics Viewer. Cancer Res. 77, e31–e34 (2017). genomes. Nature 583, 83–89 (2020). 74. Li, H. Yak: yet another k-mer analyzer. GitHub https://github.com/lh3/yak (2020). 117. Lin, M. F. Gfabase: GFA insert into GenomicSQLite. GitHub https://github.com/mlin/ 75. Li, H. Gfatools: tools for manipulating sequence graphs in the GFA and rGFA formats. GitHub gfabase (2021). https://github.com/lh3/gfatools (2021). 118. Wilkins, D. Gggenes: a ggplot2 extension for drawing gene arrow maps. GitHub https:// 76. Li, H. New strategies to improve minimap2 alignment accuracy. Bioinformatics 37, github.com/wilkox/gggenes (2022). 4572–4574 (2021). 119. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 77. Cunningham, F. et al. Ensembl 2022. Nucleic Acids Res. 50, D988–D995 (2022). 2078–2079 (2009). 78. National Center for Biotechnology Information. VecScreen: screen a sequence for vector 120. Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. contamination. Natl Cent. Biotech. Inf. https://www.ncbi.nlm.nih.gov/tools/vecscreen Preprint at arXiv https://doi.org/10.48550/arXiv.1207.3907 (2012). (2016). 121. Mose, L. E., Wilkerson, M. D., Hayes, D. N., Perou, C. M. & Parker, J. S. ABRA: improved 79. Mikheenko, A., Prjibelski, A., Saveliev, V., Antipov, D. & Gurevich, A. Versatile genome coding indel detection via assembly-based realignment. Bioinformatics 30, 2813–2815 assembly evaluation with QUAST-LG. Bioinformatics 34, i142–i150 (2018). (2014). 80. Rhie, A., Walenz, B. P., Koren, S. & Phillippy, A. M. Merqury: reference-free quality, 122. Poplin, R. et al. Scaling accurate genetic variant discovery to tens of thousands of samples. completeness, and phasing assessment for genome assemblies. Genome Biol. 21, 245 Preprint at bioRxiv https://doi.org/10.1101/201178 (2017). (2020). 123. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic 81. Chin, C.-S. et al. A diploid assembly-based benchmark for variants in the major features. Bioinformatics 26, 841–842 (2010). histocompatibility complex. Nat. Commun. 11, 4794 (2020). 124. Wagner, J. et al. Benchmarking challenging small variants with linked and long reads. Cell 82. Krusche, P. et al. Best practices for benchmarking germline small-variant calls in human Genom. 2, 100128 (2022). genomes. Nat. Biotechnol. 37, 555–560 (2019). 125. Olson, N. D. et al. PrecisionFDA Truth Challenge V2: calling variants from short and long 83. Cleary, J. G. et al. Comparing Variant Call Files for performance benchmarking of next- reads in difficult-to-map regions. Cell Genom. 2, 100129 (2022). generation sequencing variant calling pipelines. Preprint at bioRxiv https://doi.org/10.1101/ 126. Pedersen, B. S. & Quinlan, A. R. Mosdepth: quick coverage calculation for genomes and 023754 (2015). exomes. Bioinformatics 34, 867–868 (2018). 84. Zook, J. Genome in a Bottle—v3.0 genome stratifications. Natl Inst. Stand. Technol. 127. Tilgner, H., Grubert, F., Sharon, D. & Snyder, M. P. Defining a personal, allele-specific, and https://doi.org/10.18434/mds2-2499 (2021). single-molecule long-read transcriptome. Proc. Natl Acad. Sci. USA 111, 9869–9874 85. Numanagic, I. et al. Fast characterization of segmental duplications in genome assemblies. (2014). Bioinformatics 34, i706–i714 (2018). 128. Davis, C. A. et al. The Encyclopedia of DNA Elements (ENCODE): data portal update. 86. Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, Nucleic Acids Res. 46, D794–D801 (2018). 3094–3100 (2018). 129. ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human 87. Slater, G. S. C. & Birney, E. Automated generation of heuristics for biological sequence genome. Nature 489, 57–74 (2012). comparison. BMC Bioinformatics 6, 31 (2005). 130. Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10, giab008 88. Shumate, A. & Salzberg, S. L. Liftoff: accurate mapping of gene annotations. (2021). Bioinformatics 37, 1639–1643 (2020). 131. Patro, R., Duggal, G., Love, M. I., Irizarry, R. A. & Kingsford, C. Salmon provides fast and 89. Lee, C., Grasso, C. & Sharlow, M. F. Multiple sequence alignment using partial order graphs. bias-aware quantification of transcript expression. Nat. Methods 14, 417–419 (2017). Bioinformatics 18, 452–464 (2002). 132. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-seq data with or 90. Li, H. Identifying centromeric satellites with dna-brnn. Bioinformatics 35, 4408–4410 without a reference genome. BMC Bioinformatics 12, 323 (2011). (2019). 133. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 91. Gao, Y. et al. abPOA: an SIMD-based C library for fast partial order alignment using adaptive 9, 357–359 (2012). band. Bioinformatics 37, 2209–2211 (2021). 134. Pertea, G. & Pertea, M. GFF Utilities: GffRead and GffCompare. F1000Res. 9, 304 (2020). 92. Hickey, G., Paten, B., Earl, D., Zerbino, D. & Haussler, D. HAL: a hierarchical format for 135. Grytten, I. et al. Graph Peak Caller: calling ChIP-seq peaks on graph-based reference storing and analyzing multiple genome alignments. Bioinformatics 29, 1341–1342 (2013). genomes. PLoS Comput. Biol. 15, e1006731 (2019). 93. Hickey, G. Hal2vg: convert HAL to vg-compatible sequence graph. GitHub https://github. com/ComparativeGenomicsToolkit/hal2vg (2021). 94. Doerr, D. & Marijon P. GFAffix: identify walk-preserving shared affixes in variation graphs Acknowledgements We would like to acknowledge S. Bidwell and other members of the and collapses them into a non-redundant graph structure. GitHub https://github.com/ GenBank staff at the National Center for Biotechnology Information (NCBI; NLM/NIH) for their marschall-lab/gfaffix (2022). work to release the assemblies into GenBank. Certain commercial equipment, instruments or 95. Guarracino, A., Mwaniki, N., Marco-Sola, S. & Garrison, E. wfmash: a pangenome-scale materials are identified to adequately specify experimental conditions or reported results. aligner. Zenodo https://doi.org/10.5281/zenodo.6949373 (2021). Such identification does not imply recommendation or endorsement by the National Institute 96. Garrison, E. & Guarracino, A. Unbiased pangenome graphs. Bioinformatics 39, btac743 of Standards and Technology nor does it imply that the equipment, instruments or materials (2023). identified are necessarily the best available for the purpose. Computational infrastructure 97. Garrison, E. et al. smoothxg: normalization of variation graphs with local partial order and support were provided by the Centre for Information and Media Technology at Heinrich realignment. Zenodo https://doi.org/10.5281/zenodo.7239231 (2022). Heine University Düsseldorf. This work was funded in part by the National Human Genome 98. Guarracino, A., Heumos, S., Nahnsen, S., Prins, P. & Garrison, E. ODGI: understanding Research Institute of the National Institutes of Health under award numbers U41HG010972, pangenome graphs. Bioinformatics 38, 3319–3326 (2022). 1U01HG010973, U41HG007234, 1R01HG011274, R01HG010485, U24HG010262 U01HG010963, 99. Paten, B. et al. Superbubbles, ultrabubbles, and cacti. J. Comput. Biol. 25, 649–663 (2018). U24HG007497 and R01HG011649. This work was funded in part by the National Institutes of 100. Garrison, E. vcfbub: popping bubbles in vg deconstruct VCFs. Zenodo https://doi.org/ Health under award numbers U01HG010961, OT2OD033761, U24HG011853, R01-HG006677, 10.5281/zenodo.7239225 (2022). R35-GM130151, R01HG002385, R01HG010169, U01HG01973, 5U01HG010971, R01GM123489, 101. Li, H. SDUST: symmetric DUST for finding low-complexity regions in DNA sequences. U24HG009081 and 1ZIAHG200398. This work was funded in part by the Intramural Research GitHub https://github.com/lh3/sdust (2019). Program of the National Human Genome Research Institute, National Institutes of Health. The work of F.T.-N. and V.A.S. was supported by the National Center for Biotechnology W.T.H., P.P., E.E.E., T.M., H.P.J. and H.M. Pangenome graph creation: D.D., G.H., A.G., S.H., Information of the National Library of Medicine (NLM), National Institutes of Health. This work M.N.M., F.V., E.G., Y.G. and S.M.-S. Data coordination and management: M.H., J.K.L., R.S.F., was funded in part by the USDA National Institute of Food and Agriculture, grant number 2018- W.T.H., M.J., C.T., A.M.P., E.D.J., K.H.M., T.W., L.L.A.-F., S.C., M.D., S.F. and R.E.G. Transcriptome 67015-28199, and the National Science Foundation (NSF), grant IOS-1744309, and NSF PPoSS and annotation: M.H., J.M.E., F.J.M., M.R.V., M.J.P.C., K.B., M.D., A.F., C.G.G., L.H., T.H. and F.F.T. award number 2118709 (E.G. and P.P.). This work was funded in part by the Natural Sciences Pangenome applications (small variants): G.H., J. Monlong, C.M., A.M.N., P.P., J.M.Z., G. Baid, and Engineering Research Council of Canada (NSERC). G.Bourque is supported by a Canada A.B., A.C., P.-C.C., D.E.C., H.P.J., A.K., M.N., K.S. and J.W. Pangenome visualization and complex Research Chair Tier 1 award, a FRQ-S, Distinguished Research Scholar award and by the World loci analysis: S.L., J.C., C.F. and A.G. Population genetic analysis: S.B., A.G. and V.C. Sample Premier International Research Center Initiative (WPI), MEXT, Japan. J.S. was supported by the selection: X.F., K.M.M., A.M.P., E.E.K., E.E.E., K.H.M., S.F. and J.O.K. Sequencing: R.S.F., M.J., Carlsberg Foundation. This work was funded in part by intramural funding at the National M.W.M., K.M.M., H.E.O., A.M.P., E.E.E., E.D.J., K.H.M., C.A.B., O.F., R.E.G., K. Hoekzema, J.O.K., Institute of Standards and Technology. E.E.E., D.H. and E.D.J. are investigators of the Howard J.K., A.P.L., J. Mountcastle, S.S. and A.D.S. Pangenome applications (ChIP-seq analysis): C.G. Hughes Medical Institute. This work was funded in part by an Oxford Nanopore Research grant and G. Bourque. Principal investigator and laboratory organizer within the HPRC: M.J., M.W.M., (SC20130149) awarded to M. Akeson, University of California Santa Cruz. This work was funded M.J.P.C., P.F., E.E.K., A.M.P., E.E.E., D.H., E.D.J., K.H.M., T.W., T.M., B.P., R.E.G. and V.A.S. in part by Wellcome Trust award numbers WT104947/Z/14/Z, WT222155/Z/20/Z and WT108749/ Pangenome applications (VNTR analysis): T.-Y.L. and M.J.P.C. Pangenome applications Z/15/Z. This work was funded in part by a Juan de la Cierva fellowship grant (IJC2020-045916-I) (RNA-seq analysis): J.A.S. and J.M.E. Development of algorithms and software: J.S., H.P.J., H.L., funded by MCIN/AEI/ 10.13039/501100011033 and by the European Union NextGenerationEU/ G.H., A.M.N., P.P., A.G., S.H., D.D., M.N.M., F.V., E.G., Y.G., S.M.-S., H.L., J.M.E., J.A.S., B.P., T.M. PRTR. This work was funded in part by the Novo Nordisk Foundation (NNF21OC0069089). S.H. and X.H.C. Ethical, legal and social implications: R.M.C.-D., N.A.G., B.A.K., A.M. and A.B.P. acknowledges funding from the Central Innovation Programme (ZIM) for SMEs of the Federal Programme organization: A.L.F., B.I.S., M.W.S. and H.J.S. Ministry for Economic Affairs and Energy of Germany. This work was supported by the BMBF- funded de.NBI Cloud within the German Network for Bioinformatics Infrastructure (de.NBI) Competing interests E.E.E. is a scientific advisory board (SAB) member of Variant Bio. P.F is a (031A532B, 031A533A, 031A533B, 031A534A, 031A535A, 031A537A, 031A537B, 031A537C, member of the SABs of Fabric Genomics and Eagle Genomics. E.E.K. is a member of the SAB 031A537D and 031A538A). This work was funded in part by the German Federal Ministry of of Encompass Biosciences, Foresite Labs and Galateo Bio and has received personal fees Education and Research (BMBF) (031L0184A) and the European Commission, Innovative from Regeneron Pharmaceuticals, 23&Me and Illumina. A.B., A.C., P.-C.C., D.E.C., G.Baid, A.K., training network (ITN) (956229). W.-W.L. was supported in part by the Government Scholarship M.N. and K.S. are employees of Google and own Alphabet stock as part of the standard to Study Abroad (GSSA) from the Ministry of Education of Taiwan. compensation package. Author contributions Pangenome empirical analysis and pangenome quality control: W.-W.L., D.D., M.H., G.H., C.M., J. Monlong, H.J.A., J.M.Z., E.E.E., T.M., I.M.H., P.M. and J.W. Paper writing: Additional information W.-W.L., M.A., J.E., D.D., M.H., G.H., S.L., J. Monlong, R.S.F., S.G., T.-Y.L., M.W.M., A.M.N., H.E.O., Supplementary information The online version contains supplementary material available at T.P., J.A.S., M.R.V., G. Bourque, K.H.M., E.G., T.M., I.M.H., B.P., R.E.G. and L.H. Paper editing: https://doi.org/10.1038/s41586-023-05896-x. W.-W.L., D.D., M.H., G.H., X.H.C., H.C., A.G., A.M.N., P.P., A.M.P., E.E.E., E.D.J., K.H.M., E.G., E.E.K., Correspondence and requests for materials should be addressed to Erik Garrison, Tobias T.M., I.M.H., H.L., B.P., O.E.C., P.E., G.F., A.N.A.T. and A.V.Z. Assembly creation: M.A., J.K.L., H.C., Marschall, Ira M. Hall, Heng Li or Benedict Paten. A.M.P., H.L., D. Puiu, A.A.R. and A.V.Z. Assembly quality control and assembly reliability Peer review information Nature thanks Melissa Gymrek, Martin Mascher and the other, analysis: M.A., J.K.L., H.C., J.C., S.G., K. Howe, T.P., D. Porubsky, C.T., M.R.V., A.M.P., J.M.Z., E.E.E., anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer K.H.M., H.L., R.E.G., S.K., J. McDaniel, S.N., N.D.O., D. Puiu, M.R., A.A.R., A.R., V.A.S., K.S., F.T.-N., reports are available. J.W., B.W., J.M. D.W. and A.B.P. Pangenome applications (structural variants): J.E., G.H., H.J.A., Reprints and permissions information is available at http://www.nature.com/reprints. Article Extended Data Fig. 1 | Characterizing uncovered reference bases using satellite 3 comprising ~50% of these bases, mainly due to their highly repetitive peri/centromeric annotation and evaluating the completeness of different composition and also higher frequency compared to other satellites. Other satellite families. We characterized the regions not covered by the assembly centromeric satellites, centromeric transition regions, and rDNA arrays alignments to the T2T-CHM13 (v.2.0) reference and also investigated the accounted for another ~40% of the uncovered bases on average. Panels c and d completeness of the peri/centromeric satellites across all HPRC assemblies. display the average lengths of uncovered regions located within each satellite We characterized these regions using the peri/centromeric annotation family. Panels e and f show what percentage of each satellite family was covered available for the T2T-CHM13 (v.2.0) reference. We made separate bar plots for by at least one assembly alignment. The most complete centromeric regions male and female samples to exclude chromosome X for the paternal assemblies (~90% coverage) are divergent/monomeric alpha satellites, gamma satellites of male samples and exclude chromosome Y for all other assemblies. Panels and centromeric transition regions. The rDNA arrays have been covered by ~8% a and b indicate that on average ~90% of the uncovered bases are located in on average, which made them the least completely assembled repeat arrays. peri/centromeric regions with the active/inactive alpha satellites and human Extended Data Fig. 2 | Segmental duplication reliability. a, Average number of Mbp per haplotype of correctly or incorrectly assembled SDs lifted from T2T-CHM13 (v.2.0). b, The features of the most identical and longest overlapping SDs for each type of assembly error calculated in 5 kbp windows. Article Extended Data Fig. 3 | The differences in pangenome graph construction pangenome graph constructed by Minigraph. PGGB used a symmetric, methods for Minigraph, MC, and PGGB. a, Two haplotypes (H1 and H2) vary in all-by-all alignment of haplotypes to build a pangenome graph whose structure copy number of a chromosomal segment S. The S1, S2, and S3 segments are is not affected by the order of inputs (unlike Minigraph and MC). The critical highly similar with only a SNP or a small indel. b, Pangenome graph structures difference in graph construction is that, due to ambiguous pairwise for Minigraph, MC, and PGGB. Minigraph used H1 as an initial backbone and relationships of paralogs, PGGB tends to collapse copy-number polymorphic then augmented with SVs (≥50 bp) from H2, such that the SNP in S2 is not loci like segmental duplications and VNTRs into a single copy through which represented in the pangenome graph. MC added small variants (<50 bp) to the haplotypes loop, while Minigraph and MC do not. Extended Data Fig. 4 | HiFi read depth of on- and off-target edges in the MC Right: read depth distribution of off-target edges. Samples are sorted by graph. Left: fraction of reads aligned to the pangenome graph after filtering sequencing coverage (Supplementary Table 1). low-quality alignments. Middle: read depth distribution of on-target edges. Article Extended Data Fig. 5 | Gene mapping in the pangenome graphs. The first graph based on GRCh38, and CAT run on the PGGB graph. The second three three show the percentage of protein-coding genes from GENCODE (v.38) able show the percentage of noncoding genes from GENCODE (v.38) able to be to be mapped in the gene annotation sets from Ensembl, CAT run on the MC mapped on the same annotation sets. Extended Data Fig. 6 | Structural haplotypes of CYP2D6 and CYP2D7 from paths through the graph. The colour gradient and lines show the path of each the MC graph. a, Locations of CYP2D6 and CYP2D7 within the graph. The colour allele; red, start of a path; blue, end of a path. c, Frequency and linear structural gradient is based on the precise relative position of each gene; green, head of a visualization of all structural haplotypes called by the graph among 90 haploid gene; blue, end of a gene. b, Different structural haplotypes take different assemblies. Article Extended Data Fig. 7 | Performance comparison of pangenome-based coloured by their super-population label from the 1KG. b, The proportion of variant calling and read mapping across populations. a, Number of variants mapped reads that align perfectly (y-axis) is shown for a subset of samples from with at least one alternate allele (i.e. excluding homozygous for the reference the 1KG, ordered by the number of variants called (x-axis). Two mapping allele) for each in the 1KG samples. The number of variants in the 1KG callset approaches are compared: mapping short reads to GRCh38 with BWA (green); (x-axis) are compared to the variants found when aligning reads to the HPRC mapping to the HPRC pangenome with Giraffe (orange). The samples were pangenome and calling variants with DeepVariant (y-axis). Points (samples) are selected to span the x-axis. Extended Data Fig. 8 | Improved genotyping in the challenging medically- The HPRC-Giraffe-DeepVariant calls show higher frequencies. In particular, relevant gene RHCE. a, Gene annotation of part of the RHCE gene. b, Genotyping the gene-converted alleles, at about 25.406-25.410 Mbp, are observed at ~25% performance in this region for three approaches (horizontal panels). The top frequency, similar to estimates from the HPRC haplotypes (Fig. 5a–c). panel, using the HPRC pangenome, shows the best performance with most d,e, A pangenomic view of the gene-converted region showing 1 of 4 haplotypes variants being true positives (TP, blue points) based on the CMRG (v.1.0) truth in the HPRC pangenome supporting the non-reference alleles. The inclusion of set while more other methods have a higher number of false negatives (FN, this haplotype in the HPRC pangenome enables short sequencing reads, here red points). c, Allele frequency across 2,504 unrelated individuals of the 1KG. from HG002, to map along this gene-converted haplotype. Article Extended Data Fig. 9 | Leave-one-out experiment. A leave-one-out complexity: biallelic regions of the MC graph include only bubbles with two experiment was conducted by repeatedly removing one of the assembly- branches, and multiallelic regions include all bubbles with > 2 different samples from the panel VCF and genotyping it based on the remaining samples. alternative paths defined by the 88 haplotypes. b, results of the same Plots show the resulting weighted genotype concordances for different variant experiment stratified by different genomic regions defined by the GIAB. allele classes. a, weighted genotype concordances are stratified by graph Extended Data Fig. 10 | Additional applications supported by the and STAR was used to map to the linear reference. c, Proportion of all ChIP-seq pangenome reference. a, Performance of read alignment in VNTR regions peaks that are called only in the MC graph. Each data point represents samples using the MC graph versus GRCh38. All statistics are expressed relative to that were assayed for H3K4me1, H3K27ac histone marks or chromatin the total number of reads simulated from each genome. b, Performance of accessibility using ATAC-seq. d, H3K4me1 peaks that overlap an SV for which RNA-seq read alignment. Mapping rate and false discovery rate are stratified the sample is heterozygous. The reads within the peak are partitioned between by mapping quality producing the curves shown. The MC graph is compared the SV or reference allele. The red boundary represents regions where a to a graph derived from the 1KG variant calls and to GRCh38. Each reference is binomial test assigns a peak to the SV allele, both alleles, or the reference allele. augmented with splice junctions. vg mpmap was used to map to the graphs, Article Extended Data Fig. 11 | Number of SVs per sample in the HPRC PanGenie superpopulation, while the right plot shows the same results including african filtered set as well as the 1KG Illumina calls for all 3,202 1KG samples. samples and including the assembly samples present in the graph (marked by a Samples are coloured by superpopulation. The left plot excludes the african black circle).