Locityper

Software

Clone Locityper to access additional scripts:

git clone https://github.com/tprodanov/locityper

Install additional dependencies (not necessarily via conda):

conda install agc seqtk brotli

Data preparation

Download test data from Zenodo (opens in a new tab) (107 Mb) and decompress it:

xz -d data.tar.xz

Please, also download AGC (opens in a new tab)-compressed HPRC2 assemblies here (3.1 Gb) (opens in a new tab).

Download (936 Mb) (opens in a new tab) or copy CHM13 reference genome into the working directory:

wget https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0.fa.gz
gzip -d chm13v2.0.fa.gz
mv chm13v2.0.fa CHM13.fa
samtools faidx CHM13.fa
jellyfish count --canonical --lower-count 2 --out-counter-len 2 --mer-len 25 \
    --threads 8 --size 3G --output CHM13.jf CHM13.fa

Workflow

Constructing reference panels

While it is possible to extract target haplotypes from a Minigraph-cactus VCF file (see instructions), we recommend extracting haplotypes directly from the assemblies since some variation may be missing from the VCF file. Please run

locityper/extra/extract-targets.sh \
    -i CHM13.fa -i HPRC_r2_assemblies_0.6.1.agc \
    -n data/assembly_aliases.txt \
    -c data/targets.CHM13.bed -r CHM13.fa \
    -o haplotypes

The script will map extract target sequences (-c <BED file>) from the corresponding reference genome (-r <FASTA>); and map them to each of the assemblies (-i). In this case, it's the CHM13 reference genome and the assemblies from the HPRC2 AGC file. To speed things up, we only specified GRCh38 and two other assemblies in the -n data/assembly_aliases.txt, setting others to skip. On a server, you can run multiple instances of this command at the same time, as they will automatically process different assemblies. In total, this process will take around 2-3 minutes per assembly. You can find more information about the output files here.

Constructing target database

Now we need to construct Locityper target database on top of the extracted haplotypes.

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

We also provided local haplotypes for all HPRC2 and HGSVC3 assemblies in data/panels. Let's rerun database construction:

locityper target -d db -r CHM13.fa -j CHM13.jf -L data/panels/targets.bed --force

Next step is recommended, but not required: we augment the database with pairwise haplotype alignments:

locityper augment -d db

If necessary, you can prune the database to remove very similar haplotypes up to a divergence threshold (-t) or to retain some number of haplotypes (-n):

locityper prune -d db -o db-pruned -t 0.0002

Preprocessing and genotyping

First, we need to preprocess WGS data (you can find more information here):

mkdir preproc
locityper preproc -i data/HG00733.*.fastq.gz -o preproc/HG00733 -r CHM13.fa -j CHM13.jf

Afterwards, we genotype the data (more information here):

mkdir gt
locityper genotype -i data/HG00733.*.fastq.gz -p preproc/HG00733 -d db -o gt/HG00733 --leave-out HG00733

In this example we added --leave-out HG00733 because this sample is already present in the Locityper database.

Reformatting results

Locityper has written results into gt/HG00733/loci/{locus}/res.json.gz. We can summarize them with

locityper/extra/into_csv.py -i gt/* -o summary.csv

Results may differ due to different versions and non-deterministic nature of Locityper. In our case, results look like this (explanations can be found here):

sample   locus   genotype             quality  total_reads  unexpl_reads  weight_dist  warnings
HG00733  NDUFA6  NA20762.1,HG01784.2  14.5     3231         0             0.03476      *
HG00733  NDUFB1  HG00099.2,HG00731.2  29.8     3103         0             0.00209      *

Finally, since this sample is taken from the pangenome, we can evaluate the genotyping accuracy with

locityper/extra/eval_accuracy.py -i summary.csv -d db -o eval.csv

which in our case generated

sample   locus   genotype             quality  total_reads  unexpl_reads  weight_dist  warnings  query_type  edit  size   div          qv
HG00733  NDUFA6  NA20762.1,HG01784.2  14.5     3231         0             0.03476      *         hap         0     25226  0.000000000  inf
HG00733  NDUFA6  NA20762.1,HG01784.2  14.5     3231         0             0.03476      *         hap         1     25252  0.000039601  44.022957807
HG00733  NDUFA6  NA20762.1,HG01784.2  14.5     3231         0             0.03476      *         gt          1     50478  0.000019811  47.031021393
HG00733  NDUFB1  HG00099.2,HG00731.2  29.8     3103         0             0.00209      *         hap         0     25589  0.000000000  inf
HG00733  NDUFB1  HG00099.2,HG00731.2  29.8     3103         0             0.00209      *         hap         0     25589  0.000000000  inf
HG00733  NDUFB1  HG00099.2,HG00731.2  29.8     3103         0             0.00209      *         gt          0     51178  0.000000000  inf

The output shows that the sample was genotyped almost perfectly: one of the haplotypes in NDUFA6 had edit distance of 1 to the true haplotype; all other haplotypes were predicted perfectly.

Installation Preparing targets