An efficient genotyper and star-allele caller for pharmacogenomics

Abstract

High-throughput sequencing provides sufficient means for determining genotypes of clinically important pharmacogenes that can be used to tailor medical decisions to individual patients. However, pharmacogene genotyping, also known as star-allele calling, is a challenging problem that requires accurate copy number calling, structural variation identification, variant calling, and phasing within each pharmacogene copy present in the sample. Here we introduce Aldy 4, a fast and efficient tool for genotyping pharmacogenes that uses combinatorial optimization for accurate star-allele calling across different sequencing technologies. Aldy 4 adds support for long reads and uses a novel phasing model and improved copy number and variant calling models. We compare Aldy 4 against the current state-of-the-art star-allele callers on a large and diverse set of samples and genes sequenced by various sequencing technologies, such as whole-genome and targeted Illumina sequencing, barcoded 10x Genomics, and Pacific Biosciences (PacBio) HiFi. We show that Aldy 4 is the most accurate star-allele caller with near-perfect accuracy in all evaluated contexts, and hope that Aldy remains an invaluable tool in the clinical toolbox even with the advent of long-read sequencing technologies.

Figure 1.

An example of an incorrect long read alignment to the reference genome and its correction. If a donor genome (above) contains two copies of CYP2D6 pharmacogene, any long read (gray rectangle) that spans both copies will get aligned to the reference genome (below) that contains only a single CYP2D6 copy. However, this read will get its second half (containing CYP2D6 sequence) incorrectly aligned to the CYP2D7 pseudogene owing to the high sequence similarity between these genes. The final result is the overabundance of coverage in the pseudogene region compared with the CYP2D6 region (an Integrative Genomics Viewer [IGV; Robinson et al. 2011] coverage plot is shown above the reference genome).

Figure 2.

A sample decomposition of aggregate coverage into individual structural configurations. (A) An example database of CYP2D6 structural configurations containing three such configurations (v_CYP2D6, v_CYP2D7, and v_CYP2D6*13). Regions on top of the configurations that were defined (i.e., e₁, i₁, etc.) are shaded with lighter color. In this example, v_CYP2D6 corresponds to the g₁. (B) Sample decomposition of the aggregate coverage vector cn, observed after aligning the reads originating from the donor genome (above) to the reference genome (below). As can be seen, cn can be expressed as the sum of four structural configuration vectors from the database.