Diversity and longitudinal records: Genetic architecture of disease associations and polygenic risk in the Taiwanese Han population

HiGenome features

HiGenome is a pan-Taiwanese genomic database established by CMUH. This database was launched in 2018 with the primary objective of determining the associations between genes and diseases prevalent among the Taiwanese Han population and calculating the PRSs for disease incidence prediction in asymptomatic individuals, with the ultimate goal of achieving precision in health care. In this study, to complete the genomic database, we used the TPMv1 SNP array, which is specifically designed for the Taiwanese Han population. Whole-genome sequencing data pertaining to Taiwanese Han individuals, obtained from the TWB, were used to develop a platform for imputation analysis (18). Various platforms were established for pharmacogenomic analysis (19, 21, 22), human leukocyte antigen (HLA) typing (23), kinship verification (24), ancestral analysis, and PRS modeling (25) (Fig. 1).

CMUH is an extensive academic medical center located in Taichung, Taiwan. It collaborates with a network of hospitals across the northern, central, and southern regions of the country. In this study, patients were enrolled from highly populated towns and districts across Taiwan from 2018 to 2021 (Fig. 2A). Because both genotypic and phenotypic (clinical) data are essential for identifying disease-gene associations in GWASs, genotypic data were obtained using the TPMv1 SNP array, supplemented by whole-genome sequencing data obtained from the Taiwanese population, and enhanced using imputation algorithms. This approach expanded our dataset from ~680,000 SNPs to nearly 14 million reference points, consistent with the reference sequence of the Taiwanese Han population. Clinical data were collected from patient EMRs, which were matched with relevant PheCodes. After the participants were categorized into case and control groups on the basis of 1085 phenotypes, preliminary stratification was performed by age and sex, with each disease assigned to either a base or target dataset. After genetic and phenotypic data were combined, the base dataset was used to conduct GWASs and PheWASs for all phenotypes, and the target dataset was used to calculate the PRSs of 238 diseases (fig. S1).

Clinical characteristics of the HiGenome cohort

The ages of our participants ranged from 0 to 111 years, with a male-to-female ratio of 45.3:54.7. In table S1, the mean age of male (47.89 ± 21.72) participants was slightly higher than that of female (46.37 ± 21.07) participants. Retrospective analysis of patient EMRs revealed that ~85.9% of the participants were followed up for more than 1 year, 65.3% were followed up for more than 5 years, 46.3% were followed up for more than 10 years, and 27.9% were followed up for more than 15 years (Fig. 2B, left). In 2003, the total number of diagnostic instances reached 800,000; this number increased to ~7 million by 2021 (annual average: 3 million; Fig. 2B, right). Notably, 43% of the participants received treatment in the hospital’s internal medicine departments, including the pediatric (14%) and surgical (10%) departments (Fig. 2C). Diagnostic analysis revealed that the patients primarily sought treatment for neoplasms and for diseases affecting their circulatory, endocrine, metabolic, genitourinary, or digestive system (Fig. 2D). Figure 2E depicts the age distribution of all 1085 traits. In this figure, the median age of the disease group is plotted on the x axis, and the median age of the control group is plotted on the y axis. As observed, few traits were distributed along the reference line, indicating a consistent age distribution for these traits. Notably, most of the traits were associated with a higher median age in the disease group than in the control group, confirming that the incidence of most diseases increased with age and time. Figure 2F depicts the gender ratio of the traits. At the macro level, traits exclusive to males (coordinates at 1-1) or females (coordinates at 0-0) were observed (Fig. 2F, left). For most traits, the gender ratio of male in the control group consistently ranged between 0.49 and 0.42, reflecting our cohort’s overall gender distribution. In the case group, traits in the endocrine or metabolic categories were less prevalent among male participants, and certain neoplasm traits were more prevalent among male participants (Fig. 2F, right). For more details, please refer to table S1. Overall, these observations prompted us to integrate age and gender adjustment into our subsequent GWAS analyses. Notably, our data included long-term follow-up information. Patients with chronic hepatitis B virus (HBV) infection had a higher incidence of liver cancer compared with those without HBV infection (fig. S2A). In addition, long-term diabetes was associated with an increased risk of diabetic retinopathy. Moreover, female participants were more susceptible than were male participants to diabetic retinopathy, with a highly significant P value (6 × 10⁻⁶⁸; fig. S2B).

Ancestral distribution in the HiGenome cohort

To analyze ancestry, we conducted a principal components analysis (PCA), which resulted in the extraction of principal components 1 and 2. Data from the 1000 Genomes Project were used as a reference. After aligning the HiGenome cohort with the East Asian (EAS) cohort from the 1000 Genomes Project (Fig. 3A), we found that a subset of our participants exhibited values exceeding the three quartiles of principal components 1 and 2. These participants were characterized as a potentially non-EAS subset. To ensure homogeneity, these individuals were excluded from further analyses. As shown in Fig. 3B, the quality-controlled subset of the HiGenome cohort was consistent with the reference cohort in terms of principal component 1 and 2 distribution. Further exploration of the distribution of EAS individuals revealed a subset of participants exhibiting substantial deviations from the typical EAS distribution (Fig. 3C). However, after excluding this cohort on the basis of our established criteria, we found that the distribution of most participants was consistent with that of EAS individuals, thereby confirming the robustness of our quality control process (Fig. 3D). Predominant ancestral lineages were mapped to Southern Han Chinese individuals, followed by (in descending order of frequency) Han Chinese individuals from Beijing, Chinese Dai individuals from Xishuangbanna, Kinh individuals from Ho Chi Minh, Vietnam, Japanese individuals from Tokyo, Japan, and a small subset of individuals resembling the residents of Utah with northern or western European ancestry (Fig. 3E). Although most of the participants exhibited >50% single-source ancestry (Fig. 3E, top), some did not. These participants were of mixed EAS descent, and they were retained in subsequent analyses (Fig. 3E, bottom). These findings prompted us to include PCA adjustment in our subsequent GWAS analyses.

Analysis of HLA distribution and pharmacogenomic associations

To further explore the distribution of HLA and its subsequent integration into pharmacogenomic analysis, we established a predictive model for HLA. We used the results of HLA typing and chip genotyping from the TWB as input for model training. This model is capable of identifying various HLA types, including HLA-A with 35 types (1824 SNPs), HLA-B with 75 types (2007 SNPs), HLA-C with 34 types (2151 SNPs), HLA-DPA1 with 6 types (816 SNPs), HLA-DPB1 with 33 types (1713 SNPs), HLA-DQA1 with 17 types (2190 SNPs), HLA-DQB1 with 19 types (1778 SNPs), and HLA-DRB1 with 44 types (2277 SNPs). Across these eight subtypes, an average out-of-bag accuracy of 96.86% was detected, ranging from a minimum of 92.14% (HLA-B) to a maximum of 99.69% (HLA-DPA1), with a high accuracy observed in HMC class II (fig. S3A and table S2). Haplotype analysis of HLA revealed that HLA-A*11:01 (33.16%), HLA-B*40:01 (26.50%), and HLA-C*07:02 (20.87%) were the most frequent combinations, and diplotype analysis revealed that HLA-A*11:01/24:02 (12.94%), HLA-B*40:01/46:01 (8.47%), and HLA-C*01:02/07:02 (8.68%) were the most frequent combinations (fig. S3, B and C, and table S3). Our model also revealed varying association distances between each HLA subtype. For instance, HLA-A*01:01 and HLA-A*68:01 exhibited the most distant relationship, with HLA-B*15:13 and HLA-B*27:05 being the farthest apart. Additional associations are detailed in fig. S4. Given our previous ancestry analysis results (Fig. 3), we observed significant differences and a higher proportion of certain HLA subtypes in Southern Han Chinese individuals than in Han Chinese individuals from Beijing, including HLA-A*11:01(P < 7 × 10−8) and HLA-B*40:01(P < 7 × 10−70) (fig. S5 and table S4).

Analysis of drug metabolism genes revealed that the intermediate metabolizers for CYP2C19 (49.72%) and CYP3A5 (43.10%) were the most common genes, with the least common genes being the ultrarapid metabolizers for CYP2C19 and MT-RNR1:m.1494C>T, both at a rate of 0.003% (fig. S6 and table S5). In many genetic databases, results pertaining to drug metabolism genes are rarely presented, primarily because questionnaire-based clinical data do not cover the extensive use of various drugs. Overall, our genetic database is capable of tracking changes in drug dosages throughout the course of treatment with drugs such as warfarin (21) and aminoglycosides (22), making it particularly valuable.

Atlas of GWASs conducted using PheCodes available in HiGenome

To examine the genetic bases of common diseases, we analyzed our patients’ diagnostic data. The case group comprised patients with diseases confirmed by three or more diagnostic instances conforming to the PheCode definition, whereas the control group comprised patients with at least a single diagnosis not conforming to the PheCode definition. After the extracted phenotypes were stratified (with a distribution of 8:2) into base and target datasets, the base dataset was used for GWASs on disease associations, whereas the target dataset was used for PRS calculation. For the case group, a patient threshold of more than 150 was exceeded; therefore, we examined the disease-gene associations for all 1085 traits (table S6), depicted in a circular Manhattan plot (Fig. 4, outer ring). In this plot, the 10 most significant genes for each trait are highlighted. Among these associations, the most significant pertained to traits associated both with the musculoskeletal, hematopoietic, circulatory, endocrine, or metabolic systems and with mental disorders (P < 1 × 10⁻⁷⁰). Of the 187 traits with a P value of <1 × 10⁻¹⁰, the most prevalent disease classifications were related to the circulatory system, neoplasms, and the endocrine/metabolic. The predominant disease groups for which the number of patients exceeded 10,000 (across 13 traits) were found to be those affecting the circulatory, endocrine, metabolic, digestive, or genitourinary systems (Fig. 4, middle ring). Traits with a significant gene count of 10 or more were predominantly related to the endocrine, metabolic, circulatory, or integumentary (dermatology) systems (Fig. 4, inner ring). Table S2 presents detailed information on each trait. Each of 101 traits was associated with four or more significant genes. Diseases for which the sample was larger were associated with a larger number of significant genes. Given these findings, we stratified significant gene loci (P < 5 × 10⁻⁸) into three categories on the basis of the number of genes: strong association (more than eight loci), moderate association (four to eight loci), and weak association (fewer than four loci). Predominant associations were observed for traits related to neoplasms and the endocrine, metabolic, circulatory, integumentary, or genitourinary systems (fig. S7). However, weak associations were observed for most traits included in our database.

GWASs of diseases prevalent in the HiGenome cohort

The HiGenome database contains data on the following four traits: type 2 diabetes (T2D, 250-2), chronic renal failure or chronic kidney disease (CKD, 585-3), gout (274-1), and alcoholic liver damage (ALD, 317-11; Fig. 5). These traits are prevalent in Taiwan. For T2D, 57 significant gene loci were identified (P < 5 × 10⁻⁸; Fig. 5A, top). The most significant variant was rs2237897 (KCNQ1, P = 2.9 × 10⁻⁹³), with an MAF of 0.34. As indicated by the region plot, a strong association was observed between rs2237897 and its adjacent variants in the EAS population (r² > 0.6; Fig. 5A, middle). This variant was primarily associated with diseases affecting the endocrine or metabolic systems, such as diabetes mellitus and hyperlipidemia (Fig. 5A, bottom). For CKD, nine significant gene loci were identified (P < 5 × 10⁻⁸; Fig. 5B, top). The most significant variant was rs56094641 (FTO, P = 9.3 × 10⁻¹²), with an MAF of 0.13. As indicated by the region plot, a strong association was observed between rs56094641 and its adjacent variants in the EAS population (r² > 0.8; Fig. 5B, middle). This variant was primarily associated with diseases affecting the circulatory, endocrine, metabolic, or genitourinary systems, such as hypertension, diabetes mellitus, or CKD, respectively (Fig. 5B, bottom). For gout, 11 significant gene loci were identified (P < 5 × 10⁻⁸; Fig. 5C, top). The most significant variant was rs4148155 (ABCG2, P = 9.7 × 10⁻¹⁸⁷), with an MAF of 0.32. As indicated by the region plot, a strong association was observed between rs4148155 and its adjacent variants in the EAS population (r² > 0.8; Fig. 5C, middle). This variant was primarily associated both with diseases affecting the endocrine, metabolic, or genitourinary systems and with various symptoms, such as gout, abnormal blood chemistry, CKD, and calculus (Fig. 5C, bottom). For ALD, four significant gene loci were identified (P < 5 × 10⁻⁸; Fig. 5D, top). The most significant variant was rs3782886 (BRAP, P = 1.2 × 10⁻⁴³), with an MAF of 0.32. As indicated by the region plot, a strong association was observed between rs3782886 and its adjacent variants in the EAS population (r² > 0.8; Fig. 5D, middle). This variant was primarily associated with mental disorders and diseases affecting the endocrine, metabolic, or circulatory systems, such as ALD, hypertension, and gout (Fig. 5D, bottom). Although only four diseases are presented here, we conducted GWASs for all 1085 traits by using the PheCode definition. The results are available on HiGenome for interested researchers.

Comparative genetic analysis of disease variability in the Taiwanese Han and European populations

To explore the differences between Taiwanese Han and European populations, we conducted a meta-analysis of two databases to examine the weighted effects [odds ratio (OR)] on four diseases. Data from the PheWeb database of the UKBB were used. With the same PheCode definitions used, T2D (250-2) and gout (274-1) were found to exhibit similar outcomes, unlike CKD (585-3) and ALD (317-11; fig. S8A). According to our meta-analysis, the significance of the genetic associations for T2D (250-2) and gout (274-1) increased. However, the outcomes of CKD (585-3) and ALD (317-11) were similar to those of CMUH, and they were not rendered significant by the UKBB data (fig. S8B). As shown in table S7, no SNPs were significant in our CKD data, in the UKBB, or after our meta-analysis. Although certain SNPs were significant across all three datasets for ALD, their number was smaller than that associated with gout and T2D. For common diseases such as T2D, the data of CMUH exhibited the same significant associations with the CDKAL1 and FTO genes as in the UKBB. However, the RSPO3 and AUTS2 genes exhibited significant associations only in the CMUH database. For more details, please refer to table S8 (A to D). Together, our findings confirm the distinct genetic SNP profiles associated with diseases in our population and in European populations, particularly in those with CKD and ALD.

PRSs for diseases prevalent in the HiGenome cohort

To calculate PRSs, we refined our initial group of 1085 traits on the basis of the following criteria. More than 1000 patients had to have a specific trait and at least one significant gene locus (P < 5 × 10⁻⁸). This screening procedure yielded 238 traits, whose PRSs were calculated. Area under the curve (AUC) values were calculated for these PRSs. Even when a prediction was exclusively made using a PRS, most traits did not achieve an AUC of >0.6. This threshold was exceeded by only 15 traits, including ankylosing spondylitis, type 1 diabetes, and psoriasis. Notably, when the regression models were adjusted for confounders such as age, sex, and PCA results, the number of traits with an AUC of >0.8 reached 51 (table S9). We identified 31 traits for which neither the PRS alone nor its combination with clinical features achieved an AUC value of >0.6 (Fig. 6, left). By contrast, we identified nine traits for which the PRS alone achieved an AUC value of >0.6, with its combination with clinical features achieving an AUC value of >0.7. These traits were predominantly related to neoplasms and the endocrine, metabolic, circulatory, and musculoskeletal systems (Fig. 6, right).

A total of 19,486 variants were selected for the T2D PRS model. Differences were observed between the case and control groups in terms of PRS distribution. Specifically, the median PRS was significantly higher in the case group than in the control group (Fig. 7A, top left). A forest plot revealed that the ORs for PRS, sex, and age were 1.31 [95% confidence interval (CI), 1.28 to 1.33], 0.65 (95% CI, 0.62 to 0.67), and 1.06 (95% CI, 1.06 to 1.06), respectively, all with significant between-group differences (P < 0.001). No significant between-group difference was observed in any principal component (principal components 1 to 4; Fig. 7A, bottom). Cross-validation revealed median AUC values of 0.57 [interquartile range (IQR) = 0.565 to 0.575], 0.789 (IQR = 0.786 to 0.793), and 0.793 (IQR = 0.793 to 0.8) for PRS, clinical features, and their combination, respectively (Fig. 7A, top right). A total of 23,191 variants were selected for the CKD PRS model. The median PRS was significantly higher in the case group than in the control group (Fig. 7B, top left). The ORs for PRS, sex, and age were 1.17 (95% CI, 1.15 to 1.19), 0.61 (95% CI, 0.59 to 0.64), and 1.07 (95% CI, 1.06 to 1.07), respectively (all P < 0.001). No significant between-group difference was observed in any principal component (Fig. 7B, bottom). The AUC values were 0.537 (IQR = 0.532 to 0.542), 0.812 (IQR = 0.809 to 0.816), and 0.814 (IQR = 0.811 to 0.818) for PRS, clinical features, and their combination, respectively (Fig. 7B, top right). A total of 32 variants were selected for the gout PRS model. Significant differences were observed in the distribution patterns of the PRS between the case and control groups. The median PRS was significantly higher in the case group than in the control group (Fig. 7C, top left). The ORs for PRS, sex, and age were 1.38 (95% CI, 1.35 to 1.4), 0.26 (95% CI, 0.25 to 0.27), and 1.04 (95% CI, 1.04 to 1.04), respectively (Fig. 7C, bottom). The AUC values were 0.599 (IQR = 0.594 to 0.604), 0.771 (IQR = 0.767 to 0.775), and 0.783 (IQR = 0.78 to 0.787) for PRS, clinical features, and their combination, respectively (Fig. 7C, top right). A total of 23 variants were selected for the ALD PRS model. Differences were observed between the case and control groups in terms of PRS distribution. Specifically, the median PRS was significantly higher in the case group than in the control group (Fig. 7D, top left). The ORs for PRS, sex, and age were 1.14 (95% CI, 1.12 to 1.16), 0.26 (95% CI, 0.25 to 0.27), and 1.04 (95% CI, 1.04 to 1.04), respectively (Fig. 7D, bottom). The AUC values were 0.539 (IQR = 0.534 to 0.543), 0.718 (IQR = 0.714 to 0.722), and 0.722 (IQR = 0.718 to 0.727) for PRS, clinical features, and their combination, respectively (Fig. 7D, top right). For the four diseases, the AUC values for the PRS models were ~0.6. Among all clinical features, only age and sex had significant effects; however, no contributions from principal components were observed. Notably, the AUC value for clinical features reached 0.8, and the AUC value for the combination of PRS and clinical features occasionally exceeded this threshold, indicating the strong predictive ability of the combined model. For more information on the PRS results, please refer to HiGenome. Together, our findings underscore the benefits of incorporating clinical features with PRSs to increase the accuracy of predicting diseases prevalent among the Taiwanese Han population. This integrated approach holds promise for precision medicine.

Databases

The individuals in this study were divided into two cohorts. One cohort comprised individuals who had participated in a precision medicine project at CMUH. The primary objective of this project, which was launched in 2018, was to explore genetic predisposition to common diseases within the Taiwanese Han population and establish a refined system for predicting and preventing these diseases. The project focused on patients treated at CMUH and was approved by the Institutional Review Board of CMUH. The second cohort comprised patients whose data were extracted from EMRs covering the period from 2003 to 2021. This dataset included patients’ anamnesis and laboratory data, such as DNA SNP microarray findings, which are essential for investigating drug-induced side effects. At the Department of Laboratory Medicine of CMUH, which is accredited by the American College of Pathologists, the TPMv1 array is used to identify single-nucleotide variants, such as HLA types, associated with drug-induced side effects. This array is also used to facilitate GWASs of common diseases.

As of the time of this study, a total of 413,210 patients were enrolled. After excluding twin pairs, first-degree relatives, and individuals who were not of EAS ancestry, the final study cohort comprised 323,397 participants. Recruitment is still ongoing.

To ensure patient confidentiality, personal medical details were encrypted, and patient data were used exclusively for research purposes. In our case-control study, EMR data obtained from CMUH were used as the foundational dataset. This dataset included patient demographics, laboratory results, medical procedures, and diagnostic codes as outlined by the International Classification of Diseases, Ninth Revision, Clinical Modification (ICD-9-CM) and International Classification of Diseases, Tenth Revision, Clinical Modification (ICD-10-CM). CMUH archives disease data by using ICD-9-CM and ICD-10-CM codes, with the ICD-9-CM codes automatically converted into their corresponding ICD-10-CM codes. In this study, EMR data corresponding to the period from 2003 to 2021 were analyzed. Medical diagnoses were established in accordance with the PheCode criteria, which were applied on at least three distinct occasions. To establish a control group, we applied the PheCode criteria, including individuals who did not have PheCode-defined diseases (20). This study was approved by the Institutional Review Board of CMUH (approval nos. CMUH110-REC3-005 and CMUH111-REC1-176).

Genotyping

Blood samples were collected from all participants. Genomic DNA was extracted from 200-μl blood samples by using a MagCore Genomic DNA Whole Blood Kit (RBC Bioscience, New Taipei City, Taiwan) in accordance with the manufacturer’s instructions, followed by elution at a final volume of 60 μl. An Affymetrix Axiom genotyping platform was used to obtain genetic information from the samples; specifically, we used a custom Axiom TPMv1 SNP array (Thermo Fisher Scientific, Santa Clara, CA, USA), which includes 714,457 SNPs across the entire human genome. Genotype analysis was conducted using PLINK 1.9. After SNPs with a call rate of <0.95 (–geno 0.05) were excluded, samples and SNPs with missing rates (–geno 0.02 for SNPs and –mind 0.02 for samples) were removed. Subsequently, monomorphic SNPs with a count of <10 (–mac 10) and multiallelic SNPs were eliminated. Variants with Hardy-Weinberg equilibrium (HWE) P values of <1 × 10⁻⁶ and MAF values of <1 × 10⁻⁴ were also excluded. We incorporated the following analysis criteria into our study methodology: heterozygous outliers exceeding a standard deviation value of 5, PCA outliers exceeding an IQR of 3, and mismatches between genotypic and actual sex. We also used the KING-robust kinship estimator (PLINK 2.0) (53) to remove duplicate samples and first-degree relatives from our cohort and ensure that the genetic data were free from inflationary effects. SHAPEIT4 was used to phase the TPM arrays. Imputation was performed using Beagle 5.2, which is more effective and accurate than other imputation tools. The imputed data were filtered on the basis of the following criteria: an R² alternate allele dosage of <0.3 and a genotype posterior probability of <0.9 (18). Ultimately, a total of 14,181,206 applicable variants were obtained.

Genome-wide association studies

To identify disease-associated genetic variants, summary statistics were derived from the dataset of Taiwanese Han individuals. These summary statistics were derived using PLINK 2.0 (54). After relevant data were extracted from patient EMRs, the participants were divided into a case group and a control group on the basis of their PheCode classification. For details on participant distribution, please refer to table S1. Logistic regression was used to determine the correlations between various traits. Regression models were adjusted for several confounders, such as age, sex, and PCA results. To minimize the influence of pronounced linkage disequilibrium, which can lead to overestimation, we examined the most significant variant within each genomic region. A stringent P value (<5 × 10⁻⁸) was adopted to identify significant associations between the case and control groups. The results were visualized using Manhattan, region, and quantile-quantile plots. A regional visualization highlighting the variants of interest was constructed using the PheWeb platform (55).

PCA and ancestry analysis

After integrating our dataset with that of the 1000 Genomes Project, we performed SNP pruning. For this purpose, we used PLINK, adopted a window of 200 SNPs and an r² threshold of >0.1, and advanced in steps of 100 SNPs. This process enabled us to retain 66,728 variants that were common between our dataset and that of the 1000 Genomes Project. Ancestry analysis was conducted using ADMIXTURE in a supervised manner (56).

HLA imputation

To impute HLA alleles at a four-digit resolution, HLA genotype imputation was conducted using HiBAG (57) with attribute bagging analysis. After combining the results of the TWB2 chip from the TWB (approval no. TWBR11208-04) with targeted sequencing data for the HLA region, we used a cohort of 861 individuals to train our model. During the training phase, quality control measures were applied to the TWB2 chip, and 41,857 SNPs located in the HLA gene region on chromosome 6 were extracted. These SNPs were subsequently input into the HiBAG software to generate an HLA training model. After quality control, the SNP results of the TPMv1 chip were also input into the model to obtain corresponding HLA predictions. Default settings were used for analysis. HLA genes including HLA-A, HLA-B, HLA-C, HLA-DRB1, HLA-DQA1, HLA-DQB1, HLA-DPA1, and HLA-DPB1 were imputed, with an imputation posterior probability of >0.9 regarded as reliable (23).

Pharmacogene allele construction

After imputation, variants with an INFO R² of ≥0.3, an MAF of >0.0001, and an HWE P value of >1 × 10⁻⁷ in variant call format were included for analysis. With the default settings for PGxPOP used (58), pharmacogenomic phenotypes were generated from the imputed variant call format file, with GRCh38 as the reference genome. These tools are typically used to extract gene variants and infer star allele haplotypes for each individual in accordance with Clinical Pharmacogenetics Implementation Consortium (CPIC) allele definitions. A total of 10 genes, namely, CYP2B6, CYP2C19, CYP2C9, CYP3A5, CYP4F2, DPYD, NUDT15, SLCO1B1, TPMT, and VKORC1, were included for analysis. In cases where multiple pharmacogenomic phenotypes were inferred by PGxPOP for a single individual, the phenotype was recorded as “NA.” To enhance the results of drug metabolism for HLA genes, we used the results of HiBAG analysis.

Meta-analysis using UKBB data

A meta-analysis was conducted using data from TOPMed-imputed PheWeb (59) (https://pheweb.org/UKB-TOPMed/), with corresponding PheCode matching. This analysis was conducted using METAL software (60) (https://genome.sph.umich.edu/wiki/METAL_Documentation; version: generic-metal-2011-03-25).

PRS calculation and model construction

To calculate patient PRSs, data from the CMUH cohort was randomly divided into base and target datasets. The base dataset was used to investigate the associations between study variables and diseases delineated by PheCodes; this investigation was conducted using PLINK 1.9. After variants with MAF values of >0.01 were filtered, PRSs were compiled from the target dataset by using PRSice2, with the dataset of the 1000 Genomes Project (Phase 3), which is specific to the EAS population, selected as the reference standard. PRSs were calculated using z-score normalization. These scores, along with their clinical features, were used to construct logistic regression models, which were subsequently adjusted for confounding factors such as age, sex, and PCA results.

Phenome-wide association studies

A total of 58,257,251 ICD-9-CM or ICD-10-CM diagnostic codes were combined into 1791 PheCodes (61). Because of the limited variation in data and the insufficient number of participants in certain categories, the final categorization was narrowed down to 1085 PheCodes for subsequent analyses. Logistic regression was used to analyze the associations between variants and each PheCode by using the “PheWAS” package in R software (R Foundation for Statistical Computing, Vienna, Austria) (62). All summary statistics from PheWASs are publicly available on the HiGenome website, which provides interactive Manhattan plots, region plots, and PheWAS visualizations from the PheWeb platform (55).

Statistical analysis

Baseline continuous and categorical variables were analyzed using Wilcoxon rank sum test and Pearson’s chi-square test, respectively. Between-group comparisons were conducted using one-way analysis of variance (ANOVA) followed by Tukey’s post hoc test. A two-sided P value of <0.05 was considered statistically significant. The false discovery rate was adjusted using the Benjamini-Hochberg method. All statistical analyses were conducted using IBM SPSS Statistics version 22 (IBM, Armonk, NY, USA) and R software version 4.1.0.

Ethics statement

The study protocol was approved by the Institutional Review Board of China Medical University Hospital (CMUH111-REC1-176 and CMUH110-REC3-005). Deidentified genetic and clinical data were collected after informed consent was obtained from the patients.

This PDF file includes:

Other Supplementary Material for this manuscript includes the following: