# Impute Genotypes using Sequences

This file is tab-delimited, note the column names:

model   usebx   bxlimit   k       s       nGen
diploid   TRUE    50000    10      5       50
diploid   TRUE    50000   15      10      100
pseudoHaploid   TRUE    50000   10      1       50

This is the table view of the tab-delimited file, shown here for clarity.

See the section below for detailed information on each parameter. This table serves as an overview of the parameters.

STITCH uses one of three "methods" reflecting different statistical and biological models:

• diploid: the best general method with the best statistical properties
  • run time is proportional to the square of k and so may be slow for large, diverse populations
• pseudoHaploid: uses statistical approximations that makes it less accurate than diploid
  • run time is proportional to k and may be suitable for large, diverse populations
• diploid-inbred: assumes all samples are completely inbred and as such uses an underlying haplotype based imputation model
  • run time is proportional to k

Each model assumes the samples are diploid and all methods output diploid genotypes and probabilities.

The useBX parameter is given as a true/false. Simulations suggest including linked-read information isn't helpful in species with short haploblocks (it might makes things worse). So, it's worth trying both options if you aren't sure about the length of haplotype blocks in your species.

The bxlimit parameter is an integer that informs STITCH when alignments with the same barcode on the same contig should be considered as originating from different molecules. This is a common consideration for linked-read analyses and 50kb (50000) is often a reasonable default. A lower value is considered more strict (fewer reads per moleucle) and a higher value is considered more generous (more reads per molecule). You can/should change this value if you have evidence that 50kb isn't appropriate. See haplotag data for a more thorough explanation.

The k parameter is the number of ancestral haplotypes in the model. Larger K allows for more accurate imputation for large samples and coverages, but takes longer and accuracy may suffer with lower coverage. There's value in in trying a few values of k and assess performance using either external validation, or the distribution of quality scores (e.g. mean / median INFO score). The best k gives you the best performance (accuracy, correlation or quality score distribution) within computational constraints, while also ensuring k is not too large given your sequencing coverage (e.g. try to ensure that each ancestral haplotype gets at least a certain average _X of coverage, like 10X, given your number of samples and average depth).

The s parameter controls the number of sets of ancestral haplotypes used and which final results are averaged over. This may be useful for wild or large populations, like humans. The s value should affect RAM and run time in a near-linearly.

The nGen parameter controls recombination rate between the sequenced samples and the ancestral haplotypes. It's probably fine to set it to \frac {4 \times Ne} {k} given some estimate of effective population size {Ne} . If you think your population can be reasonably approximated as having been founded some number of generations ago or reduced to 2 \times k that many generations ago, use that generation time estimate. STITCH should be fairly robust to misspecifications of this parameter.

STITCH is a genotype imputation software developed for use in the R programming language. It has quite a few model parameters that can be tweaked, but HARPY only focuses on a small handful that have the largest impact on the quality of the results. Imputation is performed on a per-contig (or chromosome) level, so Harpy automatically iterates over the contigs present in the input variant call file. Using the magic of Snakemake, Harpy will automatically iterate over these model parameters.

graph LR
    subgraph Inputs
        v[VCF file]:::clean---gen[genome]:::clean
        gen---bam[BAM alignments]:::clean
    end
    B([split contigs]):::clean-->C([keep biallelic SNPs]):::clean
    Inputs-->B & C & G
    C-->D([convert to STITCH format]):::clean
    D-->E([STITCH imputation]):::clean
    E-->F([merge output]):::clean
    G([create file list]):::clean-->E
    style Inputs fill:#f0f0f0,stroke:#e8e8e8,stroke-width:2px
    classDef clean fill:#f5f6f9,stroke:#b7c9ef,stroke-width:2px

The default output directory is Impute with the folder structure below. contig1 and contig2 are generic contig names from an imaginary genome.fasta for demonstration purposes. The directory model/params is a generic name to reflect the corresponding parameter row in the stitch parameter file, which would have explicit names in real use (e.g. modelpseudoHaploid/useBXTrue/k10_s1_nGen50/). The resulting folder also includes a workflow directory (not shown) with workflow-relevant runtime files and information.

Impute/
└── model/params
    ├── variants.imputed.bcf
    ├── variants.imputed.bcf.csi
    ├── contigs
    │    ├── contig1.vcf.gz
    │    ├── contig1.vcf.gz.tbi
    │    ├── contig2.vcf.gz
    │    └── contig2.vcf.gz.tbi
    ├── logs
    └── reports
        ├── data
        ├── contig1.stitch.html
        ├── contig2.stitch.html
        └── variants.imputed.html

While you are expected to run STITCH using your own set of configurable parameters as described in the section below, Harpy also runs STITCH with a few fixed parameters:

STITCH(
    ..., 
    bxTagUpperLimit      = 50000,
    niterations          = 40,
    switchModelIteration = 39,
    splitReadIterations  = NA
)

These are the summary reports Harpy generates for this workflow. You may right-click the images and open them in a new tab if you wish to see the examples in better detail.

STITCH Reports

Imputation Metrics

Aggregates the various outputs of a STITCH run into a single report along with bcftools stats.

Reports how effective STITCH was at genotype imputation.