First, index reference genome using samtools faidx (opens in a new tab) and Jellyfish (opens in a new tab) (e.g. using k = 25):

samtools faidx reference.fa
jellyfish count --canonical --lower-count 2 --out-counter-len 2 --mer-len 25 \
    --threads 8 --size 3G --output counts.jf reference.fa

Database with target loci can be constructed using locityper target (synonym: locityper add).

locityper target -d db -r reference.fa -j counts.jf -L loci.bed

where loci.bed is a five-column BED file: chrom start end name path/to/haplotypes.fa. Haplotypes can be extracted using the script from the following section. Alternatively, Locityper can extract haplotypes from a pangenome VCF file (see here), in that case the fifth column is not needed. Nevertheless, note that the pangenome VCF may not contain all the relevant variation.

In addition to -L you can provide regions with -l name chrom:start-end [haplotypes.fa]. It is even possible to use multiple -l/-L arguments together:

locityper target -d db -r reference.fa -j counts.jf \
    -l locus1 chr:start-end haps.fa -L loci2.bed -L loci3.bed

For a leave-out evaluation you can discard some haplotypes using --leave-out argument, but this is also possible to do during genotyping without modifying the database.

If you already have local haplotypes, you can provide them directly using -l/-L arguments. If you don't know exact region coordinates, you can align locus alleles to the reference genome with

minimap2 -cx asm20 genome.fa alleles.fa | \
    awk '{ printf("%s:%d-%d\n",$6,$8+1,$9) }' | \
    sort | uniq -c | sort -k1,1nr | head -n20

This will produce a list of possible reference coordinates ranked by the number of alleles. You can then select the most common coordinates / largest region at your discretion.

You can freely add more loci to the existing database using the same commands.

Extracting local haplotypes from assemblies

First, we need to obtain target haplotypes from the pangenome assemblies. To do that, you can use the following command:

locityper/extra/extract-targets.sh \
    -i CHM13.fa -i GRCh38.fa -i HPRC2.agc -i folder_with_HGSVC3 \
    -n assembly_aliases.txt \
    -c targets.CHM13.bed -r CHM13.fa \
    -o haplotypes

The script will first extract target sequences from the reference genome (-r) based on the coordinates from -c. Then, it will map these sequences to all assemblies, specified in (-i). Here, we provided two reference genomes (CHM13 and GRCh38), AGC (opens in a new tab)-compressed HPRC2 assemblies, and a folder with HGSVC3 assemblies. Assemblies can be renamed according to the two-column alias file (-n); some assemblies can be skipped when the second column is "skip".

This script can be executed multiple times simultaneously, as multiple instances will automatically select different assemblies to process. Each instance will consume approximately 20-30 Gb RAM and will process each assembly in 2-4 minutes.

Now, output haplotypes directory will contain the following files:

File ref_targets.fa contains target sequences from the reference genome (CHM13 in this case). copy_num.csv.gz contains approximate copy numbers of each target among the assemblies. by_assembly directory will contain assembly-specific files: paf.gz with target-to-assembly mappings and bed.gz with merged mappings. Combined haplotypes will be located in panels/{locus}.fa.gz.

After removing problematic loci from haplotypes/targets.bed, you can supply the file to Locityper via

locityper target -d db -r reference.fa -j counts.jf -L haplotypes/targets.bed

Extracting haplotypes from a VCF

As a shortcut to extracting haplotypes from the pangenomes, you can use a pangenome VCF file. For a Minigraph-Cactus (opens in a new tab) VCF file, we first need to transform it using vcfbub (opens in a new tab):

vcfbub -l 0 -i raw.vcf.gz | bgzip > vcfbub.vcf.gz
tabix -p vcf vcfbub.vcf.gz
locityper target -d db -r reference.fa -j counts.jf -L targets.bed -v vcfbub.vcf.gz

Reference haplotype name is either inferred, or needs to be provided with -g NAME.

By default, Locityper forbids overlapping variants, unless --ignore-overl is used. Nevertheless, ignoring overlapping variants leads to lost information, as only the first of the overlapping variants will be used.

Note that there are two problems associated with VCF files:

• Some variation may be absent, then all haplotypes will be more similar to the reference genome than they truly are.
• Pangenome VCF file may contain extremely large variants (see next section).

Locus extension

When constructing locus database from a pangenome VCF file, Locityper tries to extend the locus such that the new boundaries do not overlap any input variants. This process may fail if there are very long pangenomic bubbles overlapping the locus. In such cases, we recommend increasing allowed extension size (-e) or manually providing locus sequences / manually specifying bigger input region.

Since v1.1.0, Locityper allows multiple boundary extensions (default 20k 50k 200k). When one extension iteration fails, a bigger threshold is checked. Setting -e 0 will force Locityper to never extend loci, however, it will fail to add any regions overlapping boundary variants.

After locus extension it may happen that several target loci start overlapping. To check that, you can run

/path/to/locityper/extra/check_overlaps.py db [--move]

This command will list all overlapping target regions. With --move argument, the script will move all redundant loci (completely covered by other target locus) to db/redundant.

Augmentation and pruning

Locityper target database can be augmented by pairwise alignments between haplotypes using

locityper augment -d db

Afterwards, these pairwise alignments will be used to reconstruct missing read alignments and improve genotyping accuracy; as well to evaluate distance between predicted genotypes to mark potentially incorrect predictions.

In addition, database can be pruned such that clusters of similar haplotypes are discarded:

locityper prune -d db -o db-pruned

Pruning can be performed using divergence threshold (-t 0.0002 by default), or to retain at most N haplotypes (-n). Pruning is performed using hierarchical clustering. Therefore, pruned set may contain pairs of haplotypes with smaller divergence than the threshold, as the retained haplotypes belong to different clusters. For this reason, it is not recommended to run pruning on top of an already pruned database.

Pruned database will lead to faster genotyping but smaller accuracy since there will be fewer haplotypes that can be predicted by Locityper.