Inverted genomic regions between reference genome builds in humans impact imputation accuracy and decrease the power of association testing

GWAS datasets and statistical analysis

We used four prostate cancer GWAS datasets to examine the effect of imputation accuracy affected by SNV sites that fall within the BBIS regions. These four datasets were ELLIPSE OncoArray (African ancestry, 4,231 cases/3,953 controls on Illumina Consortium-OncoArray, abbreviated as “ONCO-AAPC”), AAPC GWAS (4,822 cases/4,642 controls on Illumina Human1M-Duo BeadChip, abbreviated as “AAPC1M”), California and Uganda Prostate Cancer Study (1,590 cases/1,048 controls on Illumina H3Africa consortium array, abbreviated as “AAPC-H3”), and ELLIPSE OncoArray (Hispanic ethnicity, 1,192 cases/1,052 controls on Illumina Consortium-OncoArray, abbreviated as “ONCO-LAPC”). Details of the study description, data processing and quality control filtering, and models for association testing can be found in a previous publication.¹³ Briefly, association with prostate cancer risk was estimated in each dataset using logistic regression adjusting for study, age, and top 10 principal components. Per-allele ORs and standard errors from three AAPC association results were meta-analyzed using fixed-effect inverse-variance weighting. The analyses of these datasets for association with disease phenotype has been approved by the institutional review boards of University of Southern California.

Imputation via TOPMed imputation server

Imputation-ready GWAS datasets in either GRCh37 or GRCh38 genome build were submitted to the TOPMed imputation server (https://imputation.biodatacatalyst.nhlbi.nih.gov/) for quality control (QC) and imputation. The imputation pipeline used was michigan-imputationserver-1.5.7, the imputation software was minimac4-1.0.2, and the phasing software was eagle-2.4. Imputation results were then downloaded from the server to assess the imputation quality for common variants in the BBIS region. Imputation quality for each variant was assessed using the metric Rsq, the estimated value of the squared correlation between imputed genotypes and the true, unobserved genotypes,¹⁴ as recorded in the info files output by the imputation server.

Triple-liftOver script

Triple-liftOver is a PERL script that takes an input file in PLINK bim format and converts the genomic coordinate for three consecutive bases at each chromosomal position between genome builds using UCSC’s tool liftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver). It outputs the new positions in the destination build with a category column indicating whether a variant is found in a BBIS region (Figure 1). Variants that cannot be lifted over (in the unmapped output of liftOver) or are no longer on the same chromosome are excluded from the output file. For a focal SNV to qualify as an inverted site, the heuristic requires either the succeeding base in the old build to become the preceding base in the new build or the preceding base in the old build to become the succeeding base in the new build. The triple-liftOver program does not use any allele information, nor does it verify these alleles against the reference alleles in the source and destination genome reference builds. This procedure finds the inverted sites between source and target genome reference builds based on the relative position change between the focal SNV and its two neighboring bases.

Merging SNVs inverted between genome builds into contiguous stretches

To better characterize the extent of these stretches of inverted sequence between genome builds, we merged consecutive inverted SNV sites identified by triple-liftOver on each chromosome that are within 250 kb of each other. The genomic coordinates for each stretch span the first to the last SNV. A stretch defined by a single SNV flanked by two SNVs that are not found to be in regions inverted between genome builds would be assigned a length of 1 bp. Note that the span of these stretches calculated in this manner is likely to be underestimated due to potential errors near the boundaries and would differ across array platforms due to differences in SNV content.

Trait-associated SNVs

The GRCh37 build files from GWAS catalog (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/gwasCatalog.txt.gz) and the Global Biobank Engine (https://biobankengine.stanford.edu/downloads) were used to identify the inverted sites for trait-association SNVs. For GWAS catalog file, we downloaded 331 M variant-trait pairs and examined 162,001 unique single-nucleotide sites across chromosome 1 (chr1) to chr22, X and Y. For the Global Biobank Engine, there were genome-wide summary statistics for ∼3,800 traits, and we retained all SNVs with p < 1 × 10⁻⁸ for any trait. In total, there are 2,912 traits with at least 1 SNV associated with p < 1 × 10⁻⁸, for a total of ∼2.1 M SNV-trait pairs and 298,556 unique chromosomal locations.

An approach to identify SNVs located in the BBIS regions

As of freeze 8 (r2) of the TOPMed imputation server (last accessed on August 31st, 2022), the server allows internal conversion of genome build so users can submit GWAS datasets for imputation without explicitly converting the coordinate from GRCh37 to GRCh38. We observed that this practice results in a number of directly genotyped SNVs in the input dataset being labeled “typed only” after the server converts the genome build, suggesting that these SNVs are not found in the TOPMed imputation reference panel despite being common in the population. This also created instances of imputed variants immediately adjacent to the “typed only” variant in the output with complementary alleles (Table 1).

If the genome build is updated using UCSC’s tool liftOver prior to submission to the imputation server, we found that liftOver maps these SNVs to a coordinate that is off by 1 bp compared with the mapped coordinate from the imputation server. Therefore, because of the 1-bp shift in the implemented build conversion in TOPMed server, these SNVs could not be matched to their counterparts in the reference panel and would need to be imputed back by the server.

Genome-wide examinations of this occurrence using several GWAS datasets indicated that this is not a general phenomenon across the genome. The vast majority (>99.8%) of the SNVs found on the GWAS arrays that we examined would be correctly mapped from GRCh37 to GRCh38 by the imputation server. However, we noticed when the conversion does fail, these SNVs are not randomly scattered along the genome but tend to cluster into regions (Figure S1). Since the imputed versions of these SNVs tend to have complementary alleles (e.g., rs28549707 from Table 1), we suspected that in these regions, orientations are reversed between GRCh37 and GRCh38. Specifically, the designated reference alleles for these SNVs in GRCh37 reside on the forward (plus) strands, but in GRCh38, they are on the reverse (minus) strands. Previous genome-wide alignments of GRCh38 to GRCh37 reference revealed 11 Mb (0.37% of total length) of inverted sequences. Though we could not locate a complete record of all inversions, the 10 longest inverted regions previously reported in Schneider et al.¹ coincided with our mapping of SNVs erroneously converted by the imputation server. We thus termed these regions as BBIS regions.

In principle, one approach to rescue SNVs in BBIS region is to perform in-house conversion of input data from GRCh37 to GRCh38 prior to imputation. Since liftOver does not utilize the allele information, the converted dataset in GRCh38 can be submitted to TOPMed imputation server for QC checks, followed by inspection and correction of SNVs with strand flips, as needed. Although this approach corrects non-palindromic SNVs where potential strand flips are easy to detect, it does not systematically detect and correct palindromic SNVs.

Having recognized this specific issue, we devised a heuristic that we termed triple-liftOver that could identify SNVs that fall within BBIS regions. In essence, we convert the genome build using liftOver not just for the bp of the SNV of interest (i) but also for the bp before (i − 1) and after the SNV (i + 1). If the sequence is inverted around this SNV location (j) in the new build, the succeeding base will become the preceding base, and the preceding base will become the succeeding base (Figure 1). To account for the rare event that the 3-bp sequence is no longer contiguous in the new genome reference build, we lessen the constraint to flag a site as inverted if either of the neighboring bases show the orientation change.

Validation of triple-liftOver

We validated the triple-liftOver approach using three GWAS datasets for prostate cancer¹³ that were genotyped on three different Illumina platforms: Human1M-Duo (AAPC1M), Consortium-OncoArray (ONCO-AAPC), and H3Africa (AAPC-H3) (material and methods). Using all non-palindromic SNVs on each array platform, we aimed to identify the proportion of server-identified strand flips that would be detected using the triple-liftOver approach. Across the three Illumina arrays, our approach identified all non-palindromic SNVs (733 for Human1M-Duo, 410 for Consortium-OncoArray, and 1,325 for H3Africa) flagged by the imputation server as strand flips. Given the 100% sensitivity in identifying non-palindromic SNVs residing in the BBIS regions, we used our approach to further identify 33, 53, and 59 palindromic SNVs for the Human1M-Duo, Consortium-OncoArray, and H3Africa arrays, respectively, that would escape detection by the imputation server (Table S1). The detected palindromic and non-palindromic SNVs in BBIS regions cluster into 36, 25, and 51 stretches spanned by consecutive markers on each of Human1M-Duo, Consortium-OncoArray, and H3Africa, respectively, together covering approximately 2.7–5.4 Mb in length, consisting of 501–1,578 markers found on an array (Table S2; material and methods).

We further validated the detected BBIS-region SNVs by checking the reference allele in the human reference GRCh37 and GRCh38. If a SNV identified by triple-liftOver is not actually found in the BBIS region, we would expect either of the two alleles for the SNV in GRCh37 to match the reference allele in GRCh38. In practice, we observed that for all 1,927 unique SNV sites found in BBIS regions across the three arrays, either the reference allele or the alternative allele of the SNV in GRCh37 is complementary to the reference allele in GRCh38. We note that the Illumina array designs are biased toward non-palindromic SNVs to avoid strand confusions; hence, we detected relatively fewer palindromic SNVs compared with the non-palindromic SNVs that were localized in BBIS regions. We expect our approach to identify a greater number of palindromic SNVs that could impact downstream analyses when the data are genotyped on other array platforms without this bias.

Triple-liftOver takes a variant-centric approach: an input list of SNV coordinates is used to identify a subset exhibiting behaviors consistent with the SNV falling within the BBIS regions. To evaluate how many SNVs across the genome may fall into BBIS regions, we interrogated each biallelic SNV from IMPUTE2 1000 Genomes phase 3 legend file (https://mathgen.stats.ox.ac.uk/impute/1000GP_Phase3.html) with triple-liftOver. Out of 80,978,435 SNVs with unique chromosome locations in GRCh37 that can be lifted over to the same chromosome in GRCh38, our program identified 208,930 as inverted sites (0.258%) (Figure S1; Table S3 and S4). Of these 208,930 inverted sites, we found that for 99.99% of them (208,903 out of 208,930), one of the two alleles of the SNV in GRCh37 is complementary to the reference allele in GRCh38. For the remaining 27 SNVs (0.01%), neither the reference allele nor the alternative allele of the SNV is complementary to the reference allele in GRCh38, but sequence alignment search using BLAT¹⁵ of nearby sequences shows that the sequences are indeed inverted, suggesting an annotation error for the SNV allele or a sequencing error in one of the reference genome builds.

Impact of undetected SNVs in BBIS regions on genotype imputation

We next examined the impact on imputation quality due to these SNVs falling within BBIS regions. BBIS regions are sparse, and the proportion of SNVs falling in these regions are low (∼0.26%–0.37% of the human genome¹); therefore, uncorrected SNVs that are inverted between genome builds do not result in an appreciable impact on the overall imputation quality (Figure S2). Instead, their impact on imputation, particularly for the common variants, is confined more locally. To better visualize the local impact on imputation accuracy, we grouped nearby (not necessarily contiguous) inverted SNVs that were detected by triple-liftOver into 18, 15, and 27 regions across the genome for each of the three arrays we examined, spanning ∼2.2–3.8 Mb and comprising 475–1,375 SNVs. We then extended each region by 500 kb to systematically examine the imputation quality locally within the BBIS region and the surrounding non-inverted regions.

We evaluated three possible approaches to prepare the input genetic data for imputation: (1) submitting the array data in GRCh37 and relying on the TOPMed imputation server to convert the coordinates into GRCh38, (2) manually converting the coordinates into GRCh38 using liftOver followed by manual correction of detectable strand issues for non-palindromic SNVs, or (3) using triple-liftOver to systematically identify SNVs falling in BBIS regions and flip the strand for all detected SNVs, including both palindromic and non-palindromic SNVs. We examined Rsq output by the imputation server to compare imputation accuracy across the three approaches.

Using the region around 46 Mb on chr10 from OncoArray (ONCO-AAPC) as an illustrative example, we describe in more detail the impact of three approaches to dealing with SNVs in BBIS regions when imputing a dataset in GRCh37. In approach 1, all genotyped SNVs (n = 69) falling within this chr10 BBIS region were dropped and then re-imputed with relatively poor Rsq (Figure 2A), as the TOPMed server genome build conversion in these regions would result in a 1-bp shift, as described earlier (Table 1). By correcting for strand flips at non-palindromic SNVs (n = 58), approach 2 improved imputation quality for directly genotyped SNVs as they became recognized by the imputation server. However, the imputation quality for other SNVs and the overall imputation quality in the region remained poor (Figure 2B). The first approach yielded a mean Rsq 0.74 compared with 0.67 in the second approach (p < 2.2 × 10⁻¹⁶ by Wilcoxon rank-sum test), which failed to correct palindromic SNVs (n = 12; Figure 2D). The lower Rsq overall in approach 2 suggests that the uncorrected palindromic SNVs disrupt local haplotype structures due to their genotypes being incompatible for these SNVs. In approach 3, after fixing the strand issue for these 12 SNVs as informed by triple-liftOver, the overall imputation quality is increased to mean Rsq of 0.93 across the 3,903 common variants within 1.1 Mb surrounding a BBIS region (Figure 2C).

Improvement in imputation accuracy was generally observed for each of the regions we examined across the three array platforms, though the magnitude of the improvement depends on the number of palindromic and non-palindromic SNVs falling within an inverted region and how well the imputation algorithm can impute these SNVs without their actual genotypes. We also observed qualitatively similar patterns using the OncoArray dataset from Latinos (ONCO-LAPC), suggesting that the adverse effect caused by these inverted SNVs is not a phenomenon unique to the African and African American populations we examined (Figure S3–S6).

Impact of undetected SNVs in BBIS regions on downstream association testing

Despite the low prevalence of SNVs found in BBIS regions, their impact on phenotypic associations is appreciable. The BBIS region on chr10 includes the MSMB gene (chr10:46033313–46046269 on GRCh38), which contains a consistently replicated susceptibility signal, rs10993994, for prostate cancer¹³ (Figure 2). We tested a palindromic SNV within this locus, rs10763546, for its association with prostate cancer in 10,643 cases and 9,643 controls. This variant was directly genotyped on one of the three Illumina arrays (ONCO-AAPC) but had to be imputed on the other two array platforms. When submitted for imputation in hg19 (approach 1), rs10763546 was imputed poorly in all three array platforms (due to the 1-bp error in genome build conversion; Rsq ∼0.57–0.66 for rs10763546). As a result, we found only little evidence of association (odds ratio [OR] = 1.11, p_meta = 0.0011 after meta-analysis across ONCO-AAPC, AAPC1M, and AAPC-H3; Table S5), which would not be declared genome-wide significant in a standard GWAS. If we converted the genome build manually using liftOver prior to submission for imputation (approach 2), the palindromic SNV rs10763546 would not be detected as a strand flip, resulting in relatively poor imputation, even for a directly genotyped variant on the OncoArray (Rsq = 0.80). Association testing in ONCO-AAPC resulted in a null association for rs10763546 (OR = 1.04, p = 0.42; Table S5) and only a marginally improved meta-analysis association signal (OR = 1.10, p_meta = 0.00020; Table S5). Finally, when we applied triple-liftOver to identify and correct all SNVs requiring strand flips in this region (approach 3), the imputation accuracy for this SNV and the local region is improved (Rsq = 1 for rs10763546), resulting in a nearly genome-wide significant meta-analysis result of OR = 1.14 and p = 2.86 × 10⁻⁷. We repeated the same analysis for this SNV in ONCO-LAPC (1,192 cases and 1,052 controls). The version with fixing all strand issues (approach 3) yields an OR = 1.41 and p = 6.52 × 10⁻⁸ compared with OR = 1.19 and p = 0.045 from the server liftOver version (approach 1; Table S5). Taken together, our proof-of-principle analysis demonstrated that the ability to detect or further fine-map this associated locus would be hampered simply due to a change in reference genome build.

We also examined all of the trait-associated SNVs found in the GWAS catalog¹⁶ and the Global Biobank Engine,¹⁷ which is primarily based on UK Biobank data (material and methods). For the GWAS catalog, we examined 162,000 unique SNVs, among which 151 SNVs are found in the BBIS regions (∼0.09%, in a total of 277 variant-trait pairs; Table S6). In the Global Biobank Engine, we obtained 2.1 M SNV variant-trait pairs that would collapse into 300,000 unique SNVs, and 409 SNVs would be in the BBIS regions (∼0.14%, in a total of 3,385 variant-trait pairs; Table S7). These estimates are broadly consistent with the estimated proportion of SNVs falling into the BBIS regions (∼0.26%–0.37% of the human genome¹). Functional annotation of SNVs in BBIS regions revealed 1,961 non-synonymous variants and 678 substitutions in 54 genes that were predicted to be deleterious by Sift and PolyPhen. Among non-coding SNVs, 187 variants had CADD scores greater than 20, corresponding to the top 1% most deleterious substitutions in the genome (Table S8; Figure S7). Taken together, we demonstrate that BBIS regions harbor variants of significance to multiple phenotypes and that ignoring allelic errors within these regions could lead to missed opportunities in identifying a genetic association.