Cyrius: accurate CYP2D6 genotyping using whole-genome sequencing data

Abstract

Responsible for the metabolism of ~21% of clinically used drugs, CYP2D6 is a critical component of personalized medicine initiatives. Genotyping CYP2D6 is challenging due to sequence similarity with its pseudogene paralog CYP2D7 and a high number and variety of common structural variants (SVs). Here we describe a novel bioinformatics method, Cyrius, that accurately genotypes CYP2D6 using whole-genome sequencing (WGS) data. We show that Cyrius has superior performance (96.5% concordance with truth genotypes) compared to existing methods (84-86.8%). After implementing the improvements identified from the comparison against the truth data, Cyrius's accuracy has since been improved to 99.3%. Using Cyrius, we built a haplotype frequency database from 2504 ethnically diverse samples and estimate that SV-containing star alleles are more frequent than previously reported. Cyrius will be an important tool to incorporate pharmacogenomics in WGS-based precision medicine initiatives.

Fig. 1. WGS data quality in CYP2D6/CYP2D7 region.

Mean mapping quality (red line) averaged across 2504 1kGP samples plotted for each position in the CYP2D6/CYP2D7 region (GRCh38). A median filter is applied in a 200 bp window. The nine exons of CYP2D6/CYP2D7 are shown as orange (CYP2D6) and green (CYP2D7) boxes. Two 2.8 kb repeat regions downstream of CYP2D6 (REP6, chr22:42123192-42125972) and CYP2D7 (REP7, chr22:42135344-42138124) are near-identical and essentially unalignable. The purple dashed line box denotes the unique spacer region (chr22:42138124-42139676) between CYP2D7 and REP7. Two major homology regions within the genes are shaded in pink and highlight areas of low mapping accuracy.

Fig. 2. Cyrius workflow, using NA12878 (*3/*68 + *4) as an example.

A CN(CYP2D6 + CYP2D7) is derived by counting and modeling all reads that align to either CYP2D6 or CYP2D7. The histogram shows the distribution of normalized CYP2D6 + CYP2D7 depth in 2504 1kGP samples, showing peaks at CN2, 3, 4, 5, 6, and 7. The red vertical line represents the value for NA12878, corresponding to CN5 that indicates an additional copy (could be CYP2D6 or hybrid). B SVs are called by examining the CNs of CYP2D6/CYP2D7 differentiating bases. Exons are denoted by yellow boxes. Blue dots denote raw CYP2D6 CNs, calculated as CN(CYP2D6 + CYP2D7) multiplied by the ratio of CYP2D6 supporting reads out of CYP2D6 and CYP2D7 supporting reads. The red diamond denotes the CN of genes that are CYP2D6-derived at the 3′ end (can be complete CYP2D6 or CYP2D7–CYP2D6 hybrid), calculated as CN(CYP2D6 + CYP2D7) minus CN(spacer). The CYP2D6 CN is called at each CYP2D6/CYP2D7 differentiating site and a change in CYP2D6 CN within the gene indicates the presence of a hybrid. In NA12878, the CYP2D6 CN changes from 2 to 3 between Exon 2 and Exon 1, indicating a CYP2D6-CYP2D7 hybrid (*68). C Supporting read counts of the star-allele defining protein-changing small variants are used to call the CN of each variant. The y axis shows the read counts for all queried small variant positions. Six variants are called in NA12878, one of which, g.100C>T, is called as two copies (one copy belongs to *4 and the other belongs to *68). Finally, star alleles are called based on detected SVs and small variants.

Fig. 3. Depth patterns in samples with different types of SVs.

Depth plots as described in Fig. 2B. CYP2D6 CN is called at each CYP2D6/CYP2D7 differentiating site and a change in CYP2D6 CN within the gene indicates the presence of a hybrid. The depth profiles for different SV patterns are shown in NA19239 (no SV), HG02465 (deletion, *5), HG01624 (duplication), HG01161 (CYP2D7–CYP2D6 hybrid, *13), NA24631 (CYP2D6-CYP2D7 hybrid, *36), NA12878 (CYP2D6-CYP2D7 hybrid, *68), HG00290 (tandem arrangement *13 + *2), and NA19982 (two different SVs, *13 and *68, one on each haplotype). The hybrids in NA24631 and NA12878 are confirmed with PacBio reads in Fig. 4.

Fig. 4. Structural variants validated by PacBio HiFi reads.

PacBio reads supporting CYP2D6-CYP2D7 hybrid *36 and *68, confirming SVs called in NA24631 and NA12878 (third row, Fig. 3). PacBio reads were realigned against modified sequence contigs representing the hybrids and plotted using sv-viz2 [42]. The black vertical lines mark the boundaries of the duplicated sequences, represented by the gray region. The red and blue regions represent flanking sequences.

Fig. 5. CYP2D6 allele frequencies across five ethnic populations.

A Ten most common haplotypes with altered CYP2D6 function. Those with increased function are labeled in red, those with no function in black and those with decreased function in blue. B Comparison between 1kGP and PharmGKB frequencies. Each dot represents a haplotype with a frequency ≥ 0.5% in either 1kGP or PharmGKB. SV-related haplotypes are marked in red, including the two haplotypes with the largest deviation (*36 + *10 in East Asians and *68 + *4 in Europeans). Other haplotypes with deviated values are annotated in blue. A diagonal line is drawn for each panel. Correlation coefficients are listed for each population.