In order to genotype a dataset, please run

locityper genotype -i/-a/-I input -d db -p bg/SAMPLE -o gt/SAMPLE

You can limit genotyping to a subset of loci using --subset-loci argument. Additionally, one can specify genotype priors using --priors.

Please find descriptions of the input files here.

Genotyping in a weighted mode

Since v0.17 Locityper allows users to specify weights for subregions of locus haplotypes. This can be used to upweight exons and downweight introns and intergenic sequence; or to generally upweight specific subregions, known to have high functional importance.

To do so, please run

locityper genotype ... --reg-weights paths.txt

with a file with two columns <locus> <path to BED file>. Corresponding BED files should have columns <haplotype> <start> <end> <weight>, where the weight ranges between 0 and 1. Note that all haplotypes must be fully covered by the BED file, otherwise Locityper will raise an error.

For a single locus, you can use the following bash tricks:

locityper genotype ... --reg-weights <(echo "LOCUS PATH")

Custom BED file for recruiting reads

By default, when genotyping an indexed BAM/CRAM file, Locityper recruits reads based on the BED file in db/loci/<LOCUS>/ref.bed. However, since v1.2.0 you can specify a custom BED file for read recruitment (--recr-bed), which would cover, for example, other repeat copies of the locus, or specify coordinates for a different reference genome. When corresponding BED files are located in db/loci/*/name.bed you can use --recr-bed @@name.bed and use --recr-bed @name.bed when they are located at the root of the database (db/name.bed). Additionally, Locityper always examines unmapped reads and optionally examines reads mapped to short contigs (--recr-alt-len). Please set --recr-alt-len 0 if you already included all relevant contigs in the BED file.

Solver parameters

To speed up genotyping, Locityper uses multiple increasingly accurate solving steps. First, during a filtering step, locus is genotyped based on read alignments, without accounting for read depth. This is a very fast, but inaccurate procedure. Afterwards, by default, 5000 genotypes are processed using Greedy Solver (1 attempt per genotype). Finally, Simulated Annealing is run on 20 top genotypes with 20 optimization attempts per genotype. This setup can be reproduced using arguments -S greedy:i=5k,a=1 -S anneal:i=20,a=20. Each -S argument specifies the solver name: greedy; anneal; highs; or gurobi (last two require the additional installation steps). Then i=N specifies the number of top genotypes, processed during this step and a=N specifies the number of optimization attempts.

Solvers can have additional arguments, which can specified through comma with key=value format:

• Stochastic greedy:
  • x0/start: possible values best or random [default: best];
  • s/sample: sample size during each greedy iteration [default: 10];
  • p/plato: how many iterations without improvement are allowed [default: 100].
• Simulated annealing:
  • n/steps: how many steps of simulated annealing (allow switch to a worse likelihood) [default: 20k];
  • p/plato: after or during annealing, how many steps without improvement are allowed [default: 10k],
  • P/prob: initially, significantly bad read reassignments happen at this probability [default: 0.5],
• HiGHs:
  • m/mode: possible values choose/simplex/ipm (see HiGHs manual) [default: simplex].

Output files

Suppose, Locityper output directories are organized as gt/SAMPLE. Then, you can find results in gt/SAMPLE/loci/LOCUS/res.json.gz. These files can be summarized using

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

If sample name is not the end of the path, use . before it (e.g. -i gt/./*/v1). Output CSV file will contain similar information to the JSON file, but without the secondary genotype predictions. In the example workflow output file was

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      *

Here, columns are:

• genotype: predicted genotype.
• quality: quality of a predicted genotype (Phred-probability of error). If there are many very similar genotypes, quality may be low. We advise users to instead focus on the number of unexplained reads and weighted distance (see below).
• total_reads: total number of reads, used to identify locus genotype,
• unexpl_reads: number of reads that map well to some haplotypes, but not to the predicted haplotype,
• weight_dist: sum distance from the primary prediction to other genotypes. Distances are weighted by the predicted probabilities. Since v1.7.0 if target database contained PAF alignments between haplotypes, distances will be calculated as edit distance, more accurate than previously used minimizer distances. In the JSON files this will be marked with dist_type: "edit".
• warnings: optional field that specifies heuristic warnings, raised by Locityper after genotyping.
• options (only in JSON): several top genotype predictions, along with their probabilities, likelihood mean and standard deviation.

Filtering Locityper predictions

At the moment, we do not have a universal filtering recommendation. In general, examine weight_dist and unexpl_reads values for the samples for your locus of interest, if some of them are too high compared to others, this may indicate incorrect predictions.

Converting to VCF

Locityper predictions can be converted into a VCF file based on an existing pangenome VCF. For that, please run

/path/to/locityper/extra/into_vcf.py -i gts.csv -d db \
    -v pangenome.vcf -g GRCh38 -o out-dir

Converting to FASTA

In addition to a VCF file, you can generate predicted haplotypes in a FASTA file. To do so, please run

/path/to/locityper/extra/into_fasta.py -i gts.csv -d db -o out-dir

Each entry in the output files will be described in the following way:

>SAMPLE HAPLOTYPE qual=XX ...