Long read concordance

Many key pharmacogenomics genes have homology to other parts of the reference genome. This can confound short read alignments and reduce variant calling accuracy. Long reads derived from these regions cover more unique sequences on one or both sides of the homology, leading to more confident mapping, alignment, and variant calling.

For selected All of Us samples that have both srWGS and long-read whole genome sequencing (lrWGS) data available, single-sample VCFs containing SNP and indel calls for each were input to Stargazer. Calls were compared for each gene for each sample between the two technologies and results were aggregated per gene. Low concordance is indicative of potentially decreased accuracy with short reads. Additional notes are in the gene-specific details.

The 1005 samples represent the intersection of the All of Us CDRv7 dataset and the CDRv8 srWGS data. DPYD has fewer samples than other genes because of pipeline issues.

Table 3. Long read concordance results

Per-allele cosine similarity

As an additional check, population allele frequencies were compared between our results and the published results from CPIC available through PharmGKB (retrieved June 2024). We used a cosine similarity metric based on vectors of per-population allele frequencies for computed versus published. These tables are available in the PDF at the bottom of this article. Cosine similarity can be underestimated in cases where the published results are derived from a comparatively small cohort size. There are also discrepancies in how the populations are defined for All of Us versus CPIC, especially for African genetic ancestry and Americas samples. CPIC allele frequencies are often compiled from multiple studies and assign participants according to genetic ancestry categories as described on the PharmGKB site.

Note that when the allele frequency across all the intersecting populations is zero or not available for one dataset or the other then the cosine similarity is zero.

Data description

Star allele calls are provided as per-gene TSVs. The fields that are within each TSV vary depending on the gene and are described in Table 4. All TSVs contain the fields person_id, gene, genotype, and phenotype. The TSVs can be concatenated easily, but are provided separately for memory usage considerations. We recommend researchers start by looking at the phenotype field, which contains a pre-computed predicted outcome of the genotype. The genotype field is quite variable depending on the gene.

The diplotype is used to predict the phenotype using CPIC guidelines. Drug transporters and non-drug metabolizing enzymes may not have easily inferred phenotypes. For some genes, especially in cases where the relevant CPIC drug guidelines include additional genes, phenotype information may be given as "variant present" or "variant absent". Example phenotypes are given in gene-specific details.

For some pharmacogenes, a net activity score quantifies how gene activity compares to two reference alleles. They may vary based on the role of the enzyme that they are associated with. The activity score quantifies the activity of the enzyme. For example, activity scores are only appropriate for certain drug-metabolizing enzymes, such as those in the cytochrome P450 family where activity scores are defined by CPIC. There can be a possible phenotype and possible genotype field for some genes. There are two additional columns unique to CYP2D6 output, which are "copy_number" and "cyrius_filter".

Table 4. Pharmacogenomics possible fields file description

Column Explanations:

• Field name -- The name of the field. In tsv files, this will appear on the first row of the file.
• Type -- Data type.
• Key? -- Whether this field makes up a unique key for the row. Note that all key fields together make a unique key for the row.
• Notes -- Any other relevant information.

High concordance genes

CYP2C_CLUSTER

CYP2C_CLUSTER refers to the non-coding SNP rs12777823 upstream of the CYP2C cluster of genes. Only 11 samples were available to compare with the GeT-RM results, but given the perfect long read concordance for those samples and good coverage over this region, we include the 2C cluster in our high concordance results.

Gene name: 2C_CLUSTER. Genotype values for CYP2C_CLUSTER are formatted as “rs12777823 reference (G)/rs12777823 reference (G)”. Phenotypes are denoted as "Variant Present" or "Variant Absent".

ABCG2

GeT-RM validation is not available for this gene. However, ABCG2 only contains a single SNP. Due to high long read call concordance and lack of phasing challenges, this gene is included in the high concordance release.

Gene name: ABCG2. Genotype values for ABCG2 are formatted as “rs2231142 reference (G)/rs2231142 reference (G)”. Phenotypes are denoted as "Normal Function", "Decreased Function", or "Poor Function".

CACNA1S

GeT-RM validation is not available for this gene. Non-reference alleles in CACNA1S are rare with the most common allele reported at 0.01% allele frequency in the East Asian genomic ancestry population in gnomAD. c.3257G>A variant comparison allele frequencies are based on gnomAD for better correspondence.

Gene name: CACNA1S. Genotype values for CACNA1S are formatted as “Reference/Reference”. Phenotypes are denoted as "Malignant Hyperthermia Susceptibility" or "Uncertain Susceptibility".

CFTR

Phenotypes are given as Ivacaftor-responsive or not, in accordance with CPIC guidelines. CFTR alleles can also be categorized by class according to Table 1 of the CPIC guidelines in order to relate variants to treatment strategies. Some CFTR allele frequency cosine similarities are high though AFs appear low because of rounding in the PDF output of the eur.cpic column.

Gene name: CFTR. Genotype values for CFTR are formatted as “Reference/Reference”. Phenotypes are denoted as "normal_function", "decreased_function", or "unknown_function". An activity score is included for CFTR.

CYP2C9

The only discrepancy in the GeT-RM validation is that Stargazer calls a *61 allele that was not included in any of the GeT-RM assays. In sample HG01190 that allele is instead called as *2 in GeT-RM. Concordance is still high enough to merit this gene being included in our high concordance list.

Six CYP2C9 alleles have no sensitivity in our callset due to filtered SNPs in the joint-called VDS. Of these, *7, *59, *66, and *79 are classified as unknown function and *37 and *46 are classified as decreased function, but very rare. CPIC has no observations of these alleles and gnomAD shows missing decreased function alleles to be very rare:

• *37 is AF=0.0002229 in EAS pop
• *46 is AF=0.00008913 in EAS, AF=0.00001098 in SAS

Gene name: CYP2C9. Genotype values for CYP2C9 are formatted as “*1/*2”. Phenotypes are denoted as "Normal Metabolizer", "Intermediate Metabolizer", "Poor Metabolizer”, or “Indeterminate”. An activity score is included for CYP2C9.

CYP2D6

The Cyrius publication reports 99% diplotype accuracy on 144 GeT-RM samples. However, that validation is not reproduced here.

Gene name: CYP2D6. Genotype values for CYP2D6 are formatted as “*2/*4”. An example of a phenotype value is "CYP2D6 Normal Metabolizer". An activity score is included for CYP2C9. The columns "copy_number" and "cyrius_filter" are unique to CYP2D6 output. Filters are applied by the Cyrius tool used to produce CYP2D6 diplotype calls. Filters include: PASS, Not_assigned_to_haplotypes, More_than_one_possible_genotype, and LowQ_high_CN.

G6PD

Note that G6PD is located on chrX and diplotypes are modified based on the participant's X chromosome count derived in the genomic metrics auxiliary data (see the Controlled CDR Directory). If multiple alleles are reported for a sample genotyped as XY, then that sample's diplotype and phenotype are given as indeterminate.

G6PD concordance with GeT-RM and lrWGS each exceed 90%, so this gene is included in the high concordance list.

Gene name: G6PD. Genotype values for G6PD are formatted as “A”. Phenotypes are given as "Normal", "Variable", "Deficient", or "Indeterminate". As of this publication, phenotype terms for G6PD have not been standardized by CPIC.

NUDT15

The Pratt 2022 paper with NUDT15 and TPMT results only genotyped NUDT15 in 18 samples. Only 4 of those samples overlapped with NYGC 1000 Genomes samples. However, concordance with long read calls is 100% and allele frequency cosine similarity is very high for the three populations and nine alleles that overlap CPIC. This gene is included in the high concordance list.

Gene name: NUDT15. Genotype values for NUDT15 are formatted as “*1/*1”. Phenotypes are denoted as "Normal Metabolizer", "Intermediate Metabolizer", "Poor Metabolizer”, or “Indeterminate”.

RYR1

The called genotypes have no sensitivity to c.14512C>G, which does show a phenotype of increased activity. However this allele appears to be very rare. CPIC does not have population allele frequency estimates for this allele, but this variant is absent from gnomAD 4.1.

RYR1 has excellent long read concordance and good allele frequency similarity for many alleles.

Gene name: RYR1. Genotype values for RYR1 are formatted as “Reference/Reference”. An example of a phenotype value for RYR1 is "Uncertain Susceptibility".

TPMT

The called genotypes have no sensitivity to *28. This allele is not in the GeT-RM validations, is classified as unknown function, absent in CPIC, and very low AF in gnomAD. The AF in gnomAD in the EUR genetic ancestry group is 0.000002542.

TPMT has excellent GeT-RM concordance and long read concordance and good allele frequency similarity for many alleles.

Gene name: TPMT. Genotype values for TPMT are formatted as “*1/*1”. Phenotypes are denoted as "Normal Metabolizer", "Intermediate Metabolizer", "Poor Metabolizer”, or “Indeterminate”.

VKORC1

GeT-RM validation assessed VKORC1 as three separate assays for the three known VKORC1 SNPs. Stargazer run on GVCF input (as used in GeT-RM validations) fails to phase the rs72547529/V66M and rs61742245/D36Y SNPs. GeT-RM validations were compared only for the rs9923231/*S1 allele calls. The All of Us participants are only genotyped for the rs9923231
allele. With respect to the rs9923231, validation concordance is high and population allele frequency is in very good agreement so we include this in the list of high concordance genes.

Gene name: VKORC1. Genotype values for VKORC1 are formatted as “rs9923231 reference (C)/rs9923231 reference (C)". Phenotypes are given as "Variant Absent", "Variant Present", or "Indeterminate".

Low concordance genes

CYP2B6

CYP2B6 concordance was evaluated in comparison with updated GeT-RM calls from the Aldy publication (Hari et al. 2023). GeT-RM calls were updated according to the logic detailed in the supplemental notes. For CYP2B6, several diplotype calls were updated to alleles that were called by a consensus of evaluated tools where the called alleles included all the variants in the GeT-RM allele. A few other CYP2B6 diplotypes were updated because none of the tools evaluated confirmed the GeT-RM allele.

CYP2B6 long read concordance is poor. This may be a result of the high degree of self-similarity of the reference genome for parts of this gene, as evidenced by the "self-chain" track in the UCSC Genome Browser. The most common error mode seen in the lrWGS versus srWGS comparison is that lrWGS prefers *6 but srWGS calls *9.
A diagram from the CAVE viewer comparing CYP2B6 alleles at PharmVar shows some of the distinctions.

Population allele frequency comparisons are poor for *19, *20, *21, and *28, which are only present in CPIC in the sub-saharan African population. The All of Us genetic ancestry groups do not distinguish sub-saharan African from African ancestry.
Gene name: CYP2B6. Genotype values for CYP2B6 are formatted as “*1/*6”. An example of a phenotype value for CYP2B6 is “normal_metabolizer”. CYP2B6 also includes an activity_score.

CYP2C19

CYP2C19 concordance was evaluated in comparison with updated GeT-RM calls from the Aldy publication, similar to CYP2B6. CYP2C19 GeT-RM concordance is low due to data format issues with the 1000 Genomes validation data that don't apply to the All of Us pipeline, but lrWGS concordance is ~90%. This gene is included in the low concordance list.

Note that *38 is the "reference" allele with no variants compared to the GRCh38 reference, not *1, which has a SNP at chr10:94842866:A>G.

Gene name: CYP2C19. Genotype values for CYP2C19 are formatted as “*1/*38”. Phenotypes are denoted as "Normal Metabolizer", "Intermediate Metabolizer", "Poor Metabolizer”, or “Indeterminate”.

CYP3A5

This gene is included in the low concordance list.

Gene name: CYP3A5. Genotype values for CYP3A5 are formatted as “*3/*3”. Phenotypes are denoted as "Normal Metabolizer", "Intermediate Metabolizer", "Poor Metabolizer”, or “Indeterminate”. CYP3A5 includes a possible_genotype field.

CYP4F2

Variants on the *2 and *3 alleles are much farther apart than the insert size for short read sequencing, so phasing is challenging.

According to PharmVar "Although CYP4F2 is part of a CPIC guideline, function has not been assigned by CPIC."

Gene name: CYP4F2. Genotype values for CYP4F2 are formatted as “*1/*1”. An example of a phenotype value is "CYP4F2 Variant Absent".

DPYD

The DPYD gene is too large to phase variants across its entirety. As such, the nomenclature typically doesn't follow the star allele convention and variants are not mutually exclusive.

DPYD concordance was evaluated in comparison with update GeT-RM calls from the Aldy publication, similar to CYP2B6. Cosine similarity metrics are based on the first two called alleles, which may result in underestimating the similarity. DPYD calls have no sensitivity to c.46C>G, c.1403C>A, or c.2582A>G due to filtered SNPs, however these are all classified as normal function alleles.

Activity scores are given to quantify DPYD function. In scenarios where a sample carries >2 alleles for DPYD, CPIC guideline recommends using the 2 alleles with the lowest activity to assign activity score for the gene. Therefore a genotype can contain >2 alleles. Qualitative phenotypes are given based on activity scores.

Gene name: DPYD. Genotype values for DPYD are formatted as “c.496A>G/Reference”. An example of a phenotype value is "DPYD Normal Metabolizer”. DPYD includes an activity score formatted as “DPYD Activity Score 2.0”.

SLCO1B1

SLCO1B1 is a large gene containing variants that can be tens of thousands of bases apart. As such, the phasing of variants can be inconclusive. The PCR- and array-based methods used to generate the GeT-RM SLCO1B1 data do not convey any sample-specific phasing information. More information about the GeT-RM assays is available at the CDC website.

The GeT-RM validation performed here allowed for matches if the sample contained the correct variants, but the phasing differed from the GeT-RM call. Specifically:

• *1/*15 is cis version of *5/*37
• *1/*14 is cis version of *4/*37
• *21 became suballele of *20 and then *1/*20 is cis version of *19/*37

Phasing differences do not affect the predicted phenotype in these cases. Additionally, Stargazer reported *1A allele is updated to *1 and *1B becomes *37 to match the CPIC definitions.

Even with all these improvements, the GeT-RM concordance is about 55%. However, concordance with lrWGS calls is nearly 90%, so we believe most of the GeT-RM issues to be attributable to the older assay data.

Gene name: SLCO1B1. Genotype values for SLCO1B1 are formatted as “*1/*15”. Phenotypes are denoted as "Normal Function", "Decreased Function", or "Poor Function". SLCO1B1 includes a possible_genotype field and a possible_phenotype field.

UGT1A1

This gene has not been curated by PharmVar. The SNP defining the *60 allele is far upstream of the UGT1A1 transcript start site and is not included in Stargazer or PharmCAT. *60 GeT-RM calls were converted to reference for the purposes of the evaluation. *80 allele calls present in Stargazer were only produced by one GeT-RM assay and not present in any GeT-RM consensus calls. *80 calls were converted to reference for the purposes of validation. These changes helped UGT1A1 calls achieve a GeT-RM concordance of 98%. However, long read concordance is less than 90%. This may be due to low long-read sensitivity caused by coverage dropout, especially given the TATA repeats in this gene. Due to the low lrWGS concordance, UGT1A1 is included in the low concordance release.

Gene name: UGT1A1. Genotype values for UGT1A1 are formatted as “*1/*1”. Phenotypes are denoted as "Normal Metabolizer", "Intermediate Metabolizer", "Poor Metabolizer”, or “Indeterminate”. UGT1A1 includes a possible_genotype field and a possible_phenotype field.