GoKawiilPangenome reference
The pangenome reference used in this study is a sequence graph augmented with phased small variants. It is stored in two data files: a sequence graph in reference rGFA format16 containing the structural variants and a supplementary phased VCF33 file containing small variants. The rGFA sequence graph is constructed first from haplotype assemblies, then the small variants are extracted from the assemblies and added to the VCF file30.
The sequence graph G(V,E) has vertices v ϵ V and edges e ϵ E. A segment is a sub-sequence s of one or more input haplotype assemblies. Each segment is made of two vertices, having respectively the forward and reverse directions. Traversing the forward vertex spells out s while traversing the reverse vertex spells out its reverse complement. An edge e connects two segments if their respective sub-sequences are adjacent in one or more input assemblies. The edges are bidirected such that each end of an edge is either pointing forward or reverse. The direction determines the traversal directions of the connected segments. Thus, there are four forms of edges: forward–forward, forward–reverse, reverse–forward and reverse–reverse.
The rGFA graph is incrementally constructed from a set of input haplotype assemblies (Supplementary Methods). The first input assembly serves as the reference assembly, which is the basis of the coordinate system. When another haplotype assembly is inserted in the graph, existing segments might be split up and additional vertices with new sequences introduced. The sequences spelled out by paths in the graph from a previous insertion can always be spelled out in the updated graph. Since these sequences do not change, they are denoted stable. To keep the sequences stable across insertions, we give each edge and vertex a rank corresponding to the input haplotype assembly that it is originally from. Traversal of a path of rank r spells out a sub-sequence in the r-th haplotype assembly sequence.
Weaver
Weaver maps paired-end short reads in FASTQ format to a pangenome reference, such as HPRC-ICE. The reference can either be linear in FASTA format or a graph in rGFA format with an optional phased VCF file containing small variants. It outputs alignments in Sequence Alignment/Map (SAM) format53, typically compressed as a BAM or CRAM file. Before mapping, Weaver preprocesses the pangenome reference once in an indexing step.
Indexing
The Weaver index has two components. The first component is a seed index that stores the minimizers found in the pangenome reference alongside their locations. A minimizer is a sub-sequence of length k (23 by default), called a k-mer, with the minimum hash value in a sliding window containing w overlapping k-mers (11 by default). We use the minimap2 (ref. 54) hash function, which returns the same hash value for a k-mer and its reverse complement. As a result, the same set of minimizers is extracted regardless of the orientation of sequences in the pangenome reference. The minimizers are extracted from all haplotype sequences represented in the graph and in the VCF. They are then added to a seed index along with their graph locations (Fig. 3a). The minimizers may therefore contain bases from alternative alleles in the VCF file. The second component is the ICU index, which stores all pairs of vertices that see each other. We define a vertex u as seeing vertex v if there exists a path from u to v in the graph, which is d bp (d = 1,500 by default) or less.
Mapping
During Weaver read mapping, minimizers are extracted from the read sequence in the same way as during indexing. The minimizers are then anchored to the graph using the seed index (Fig. 3b). Weaver links together seed index hits to form chains (Fig. 3c) when the seeds anchor onto vertices of same rank and are interspersed by the same distances in the read and the graph. If both reads in a pair have long chains that see each other, other small chains that do not see chains on the other read are discarded (Supplementary Methods 2). The chains are then extended along the stable sequence by traversing paths of the same rank in the graph. Then we perform a pairwise alignment between the stable sequence and the read sequence (Fig. 3d). An alignment score S is calculated as
... continue reading