Leveraging deep phenotyping from health check-up cohort with 10,000 Korean individuals for phenome-wide association study of 136 traits

The expanding use of the phenome-wide association study (PheWAS) faces challenges in the context of using International Classification of Diseases billing codes for phenotype definition, imbalanced study population ethnicity, and constrained application of the results in research. We performed a PheWAS utilizing 136 deep phenotypes corroborated by comprehensive health check-ups in a Korean population, along with trans-ethnic comparisons through using the UK Biobank and Biobank Japan Project. Meta-analysis with Korean and Japanese population was done. The PheWAS associated 65 phenotypes with 14,101 significant variants (P < 4.92 × 10–10). Network analysis, visualization of cross-phenotype mapping, and causal inference mapping with Mendelian randomization were conducted. Among phenotype pairs from the genotype-driven cross-phenotype associations, we evaluated penetrance in correlation analysis using a clinical database. We focused on the application of PheWAS in order to make it robust and to aid the derivation of biological meaning post-PheWAS. This comprehensive analysis of PheWAS results based on a health check-up database will provide researchers and clinicians with a panoramic overview of the networks among multiple phenotypes and genetic variants, laying groundwork for the practical application of precision medicine.

Statistical and computational analyses. Phenome-wide association study. We used PLATO 11 to run logistic regression analysis on 65 categorical outcomes and linear regression analysis on 71 continuous outcomes, incorporating 6,860,342 genetic variants in an additive model. We included age, sex, and the first three principal components to adjust for any potential confounding bias due to these variables. To identify significant results, we implemented multiple test correction through LD-aware Bonferroni correction. The conventional Bonferroni test assumes that the association tests for all SNPs are independent and thus divides the alpha by the total number of tests. For our study, instead of correcting p-values with the total number of SNPs, we use LD pruning to identify independent SNPs 12 . The threshold we used for association between SNPs was r 2 = 0.3, which is provided by Sobota et al. for the East Asian population 13 . We established genome-wide significance at P < 4.92 × 10 -10 .
Further exploratory analyses were performed using the associated 260,923 loci with a less stringent P < 1 × 10 -4 . Though we used the LD pruning method for Bonferroni correction, the p-value was still stringent. Thus, in addition to analyzing associations with a stringent p-value cutoff, this exploratory threshold allowed us to further expand the boundaries of research by involving a much wider PheWAS landscape 12 .
To perform systematic analysis of the PheWAS results, we leveraged cross-phenotype associations, in which one locus is associated with multiple phenotypes 14 . Such associations include polygenic inheritance, where a phenotype is influenced by more than one gene 15 (Fig. S4A); and pleiotropy, where a locus or a gene affects more than one phenotype 16 (Fig. S4B). To further explore and understand polygenicity and pleiotropy, we constructed two networks: a bipartite phenotype network, connecting phenotypes that shared at least one locus 14 (Fig. S4C) and a bipartite gene network, connecting genes that shared at least one phenotype 14 (Fig. S4D). In these connections  (PheWAS) was performed for 136 phenotypes adjusting for age, sex, and PC1-PC3. (C) We leveraged cross-phenotype associations to perform systematic analysis of the PheWAS results, which were polygenicity, pleiotropy, a bipartite gene network, and a bipartite phenotype network. The details are described in Fig. 4. (D) To ensure robustness of the PheWAS results, we further dissected the results to suggest applicable interpretations, the heritability for each phenotype; Correlation between phenotype heritability and the effect of the loci on genes and protein sequences associated with phenotypes. (E) Using cross-phenotype association information, we constructed phenotype-phenotype and phenotype-genotype networks. (F) We visualized the comparison of obesity indices (body mass index, waist circumference, visceral adipose tissue, and total adipose tissue amount). (G) We constructed cross-phenotype mappings, which have a core phenotype (Pheno-1 in the figure) and branches of connected phenotypes that share loci. These were partitioned by color according to the biological system involved. (H) We estimated causal inferences in the phenotype pairs from cross-phenotype associations using Mendelian randomization and constructed a causal inference map. (I) We performed transethnic and trans-nationality analysis among Korean, European, and Japanese populations. (J) We compared phenotype-phenotype pairs generated from SNP-based cross phenotype-association in the Biobank analysis with those generated from correlation analysis in the EHR-based H-PEACE cohort. We evaluated the overlap or exclusiveness of pairs for each phenotype by phenotype degree. www.nature.com/scientificreports/ Gene annotations. We mapped the variants to genes using Ensembl Variant Effect Predictor (VEP) 18 annotations (RefSeq). By default, VEP annotates variants in 5000 bp upstream and downstream. So, the variants in 5000 bp regions were mapped to the nearest genes.
Functional annotations (p value < 1 × 10 -4 ). We mapped genetic associations using VEP 18 in order to annotate the functional relevance of significant loci. Using the VEP annotation, we classified the biological consequences of loci in coding regions (stop-gained variant, slice acceptor variant, splice donor variant, and missense variant) and in non-coding regions. We also annotated UKBB and BBJ variants with VEP to conduct trans-ethnic and trans-national comparisons as described in a later section. www.nature.com/scientificreports/ Estimated heritability. To determine the contributions of genetic variants to the risk of certain phenotypes, we estimated the heritability of each phenotype. We estimated heritability using LD Score regression with LDSC (version 1.0.1) 19 on summary statistics from the PheWAS for all phenotypes. For this analysis, we used the East Asian LD Scores from 1000 Genomes as reference LD Score, which served as the independent variable in the LD Score regression (ref-ld-chr) and regression weights (w-ld-chr). General instructions and the East Asian LD Scores from 1000 Genomes are provided here: https:// github. com/ bulik/ ldsc.

Comparison in different populations.
To compare results across diverse populations, we performed a transethnic comparison utilizing PheWAS results from a European population and a trans-national comparison utilizing results from a Japanese population. For European population, data from the UK Biobank (UKBB) 9 was used; for the Japanese population, data from the Biobank Japan Project (BBJ) 7 was used. We downloaded the summary statistics and estimated heritability results of the phenotypes of these results from the following URLs: http:// www. neale lab. is/ uk-bioba nk/ and http:// jenger. riken. jp/ en/ result. We tabulated lists of the phenotypes in the UKBB and BBJ and searched for those that were most similar to phenotypes in our database. The manually curated overlapping phenotypes among GENIE, UKBB, and BBJ are given in Table S3.
Mendelian randomization. To better understand the causal inferences in cross phenotype mapping, we performed Mendelian randomization (MR) analysis on the phenotype pairs connected in the bipartite phenotype network. To avoid potential bias due to sample overlap between exposure and outcome, we split our dataset into two equal sets by random assignment of samples. PheWAS was conducted on each dataset separately to generate the summary statistics that were used as input to MR. Additionally, significant SNPs (P < 1 × 10 -4 ) from the initial PheWAS with all samples were used as instrument variables (IV). Furthermore, all IVs that were significant in outcome (P < 0.01) were removed, as IVs should not be directly associated with outcome. The SNPs were LD-clumped using very strict cutoff of clump kb = 10,000 and r2 = 0.001. We calculated p-values using the inverse-variance weighted (IVW) method from the MendelianRandomization package in R 20 . We adjusted for multiple testing using FDR correction. We also performed sensitivity analysis using MR-egger and the medianbased method.

Meta-analysis of PheWAS.
We performed meta-analysis using our PheWAS results and the BBJ results for all phenotypes that were available in both datasets. The BBJ summary statistics came from different studies, requiring harmonization of the files. Phenotype matches between GENIE and BBJ are listed in Table S3. Some of the phenotypes from GENIE matched to multiple phenotypes in BBJ; in such cases, we carried out meta-analysis separately for each BBJ phenotype. The meta-analysis was implemented using METAL 21 . The overall scheme of our study is shown in Fig. 1.

Results
After QC, the study population of the GENIE cohort included 9742 participants, comprising 5696 males and 4046 females, with average age 50.7 + / − 10.0 years. The characteristics of the study population are given in Table S4. See the Supplementary Appendix for detailed description of the results.
Phenome-wide association analysis. From the PheWAS on 136 phenotypes, we found significant associations for 65 phenotypes and 14,101 SNPs (P < = 4.92 × 10 -10 ). The counts of significant loci and genes associated with each phenotype are given in Table S5 and most significant variants are shown in Tables S6, S7 and Fig. S5. Approximately 1% of variants were in coding regions and 98.885% were in non-coding regions (Fig. S6,  Tables S8, S9). We systematically compared the significant associations of loci and their genes with phenotypes (P < 1 × 10 -4 ) to results from the BBJ and UKBB to determine if our results were replicated in other populations and also to look for novel findings (Fig. S7, Tables S10, S11). In the comparison between Korean and UK populations, fewer overlapping loci were identified, with the highest overlap ratio being 9.15% in fatty liver disease; 42 phenotypes did not have any overlap (Fig. 2, Fig. S8).
Population comparisons were further investigated for body mass index (BMI) in particular. 34 genes were unique in our populations relative to both Japanese and European populations (Table S12, Fig. S9). Of those unique genes, 23 have previously reported associations with obesity or body weight; the corresponding literature review and references are given in Table S13. The other 11 genes have not been previously reported as associated with obesity in humans, and could be candidate novel genes for BMI or obesity. The details of the genes are described in the Supplementary Appendix.

Systematic analysis of the PheWAS results.
To perform a systematic analysis of the PheWAS results, we leveraged cross-phenotype associations, where one locus is significantly associated with multiple phenotypes. For this analysis, significant loci were filtered by a less-stringent threshold, P < 1 × 10 -4 . The schematic structure for this analysis is shown in Fig. S4. Possible polygenicity (Fig. S4A, Table S14); possible pleiotropy ( Fig. S4B, Table S15); bipartite phenotype network (Fig. S4C, Table S16); and a bipartite gene network (Fig. S4D) were drawn from PheWAS results.
The bipartite phenotype network comprised 23,580 loci (2902 genes) with 135 phenotypes. There were 1926 distinct pairs of phenotypes. We calculated the degree properties of core phenotypes in this network (Table S17), where core phenotypes were those nodes connected to several phenotypes by shared variants (Fig. S4C) www.nature.com/scientificreports/ highest possible polygenicity was observed for mean corpuscular hemoglobin concentration (MCHC), with 782 genes. The bipartite gene network comprised 14,907 genes, which were connected through sharing associations with the same phenotypes. Table S18 give the gene degree and phenotype degree values for this network. The three genes with the highest phenotype degrees were; CUB and Sushi Multiple Domains 1 Protein (CSMD1), RNA-binding Fox-1 Homolog 1 (RBFOX1), and Protein Tyrosine Phosphatase Receptor Type D (PTPRD); this could be due to possible pleiotropy.
We compared the bipartite phenotype networks of the GENIE (Korea), BBJ (Japanese), and UKBB (European) cohorts. Fig. S10 visualizes the phenotype-phenotype pairs observed in each population; 288 pairs were simultaneously observed in all three populations (Table S19). Notably, these included the pairing of red blood cell count (RBC) and brain vascular atherosclerosis. There are reports of RBC having relation to coronary artery disease 22 and stroke mortality 23 , but not directly to brain vascular atherosclerosis . Secondary analysis of the PheWAS results. Heritability analysis. Heritability was calculated for each of the 136 phenotypes by regression of LD scores (Table S20). The top heritability values were obtained for compression fracture, spondylolisthesis, and height. In terms of biological categories and body systems, the highest heritability value was obtained for the musculoskeletal system (Table S21).
We further compared the heritability in our population with that in the Japanese and European populations (Table S20). Comparisons to each of the Japanese and UK populations are shown in Fig. S11, while the three-way comparison among Korean, Japanese, and UK populations is shown in Fig. S12. The prothrombin time (PT) and activated partial thromboplastin time (aPTT), which are biomarkers of coagulation function, showed similar Network analysis. Using cross-phenotype association information, we constructed phenotype-phenotype and phenotype-genotype networks. First, a network representation of gene-phenotype associations related to metabolic syndrome was constructed (Fig. 3A). 132 genes associated with metabolic syndrome and 128 phenotypes sharing 102 genes with metabolic syndrome were used to construct the network. In the metabolic syndrome sub-network, five genes had high degrees of connection and could be considered hub genes: PTPRD, DCC Netrin 1 Receptor (DCC), Proprotein Convertase Subtilisin/kexin Type 6 (PCSK6), Unc-13 Homolog C (UNC13C), and Contactin 4 (CNTN4). The phenotypes in this network comprised: of cardiovascular diseases, of metabolic diseases, used as markers for obesity, and other various disease. The phenotypes in this network comprised: of cardiovascular diseases, of metabolic diseases, used as markers for obesity, and other various disease. The phenotype nodes included triglyceride (TG), HDL cholesterol (HDL), hypertension, diabetes, and waist circumference (WC). These results give a genetic rationale for the definition of metabolic syndrome in the PheWAS perspective. Network analysis A network representation of gene-phenotype associations related to metabolic syndrome was constructed from 102 genes associated with metabolic syndrome and 128 phenotypes sharing those genes. Each edge is a phenotype-gene association, with genes for significant loci (P < 10 -4 ) being annotated by VEP. Node size is proportional to degree, which is the number of connections. Pink nodes correspond to phenotypes and green nodes to genes. (B) Relationships among obesity indices We visualized the comparison among the obesity indices such as body mass index (BMI), waist circumference (WC), visceral adipose tissue (VAT) and total adipose tissue (TAT) amount by drawing a the venn-diagram for cross phenotype association of phenotypes or genes. (C) Cross-phenotype mapping Crossphenotype mappings were generated based on the bipartite phenotype network, in turn constructed from the connections among phenotypes sharing at least one locus. Coffee consumption, which is one of the lifestyle phenotypes, had 31 phenotype degrees in the bipartite phenotype network. (D) Causal inference mapping We estimated causal inferences in phenotype pairs based on cross-phenotype associations using Mendelian randomization (MR), and constructed a causal inference map. The direction of the arrow is the causality result from MR (Blue arrows, skeletal muscle mass as outcome; Red arrows, skeletal muscle mass as exposure; Green arrows, bidirectional). Pairs observed in the bipartite phenotype network but insignificant in MR have straight black lines without arrows. www.nature.com/scientificreports/ We also constructed a phenotype-phenotype network using 1,926 phenotype pairs based on shared loci (P < 1 × 10 -4 ). Fig. S13 shows the phenotype-phenotype network for the whole dataset, and an interactive visualization tool of the phenotype-phenotype network is available (https:// hdpm. biome dinfo lab. com/ ddn/ genie/).

Relationships among obesity indices.
Obesity is a disease entity, which the interest in and research into, has been growing 24,25 . However, definitions of pathological obesity make inconsistent use of variable traits such as body mass index (BMI), waist circumference (WC), total adipose tissue area (TAT), and visceral adipose tissue area (VAT). The defining parameter for obesity also varies between researchers and with respect to the target disease. We visualized 26 the overlap or exclusiveness among BMI, WC, TAT, and VAT based on the bipartite phenotype network and pleiotropy potential of genes. As shown in Fig. 3B, connections were observed as quadrant intersections among BMI, WC, TAT, and VAT for seven phenotypes: CA19-9, GOT, GPT, body fat mass, body fat percent, weight, and metabolic syndrome. There were 15 phenotypes connected exclusively with VAT and WC, of which, most were crucial intermediate phenotypes that link obesity with diseases. Accordingly, it can be postulated that when defining obesity, VAT or WC would better represent the characteristics of pathogenic obesity. The two genes that are exclusively overlapped between VAT and WC (Fig. S14) could be candidate genes for explaining the pathogenic role of obesity (Table S22).
Cross-phenotype mapping. Cross-phenotype mappings were generated based on the bipartite phenotype network, in which the connected phenotypes shared at least one locus.
First, we constructed a cross-phenotype mapping focused on tumor markers. Table S23 shows the respective connected phenotypes we obtained for tumor markers. Fig. S15 shows the cross-phenotype mapping for CEA, which could be considered during oncological practice in order to take into consideration all the possible effects of phenotypes other than colorectal cancer progression itself.
Second, we constructed a cross-phenotype mapping focused on lifestyle factors. In this study, we visualized the cross-phenotype mapping for the coffee consumption. Coffee consumption had 27 phenotypes connected through sharing of significant loci (Fig. 3C). The results of these cross-phenotype mappings could provide the genetic background to explain interactions between environmental factors and disease, and might further provide basic knowledge necessary to conduct gene-environment interaction analysis.
Mendelian randomization analysis. We estimated the causal inferences in phenotype pairs based on crossphenotype associations using Mendelian randomization (MR) ( Table S24). As shown in Fig. 3D, we drew a causal inference mapping centered on skeletal muscle mass. The Mendelian randomization analysis yielded nine significant phenotypes, of which one was causal for skeletal muscle mass, two phenotypes were outcomes from skeletal muscle mass, and six had bidirectional relationships with skeletal muscle mass. This analysis revealed that skeletal muscle mass was a significant causal factor for metabolic syndrome and alcohol consumption.
We also performed Mendelian randomization with a focus on lifestyle factors that were causal exposures in cross-phenotype associations, such as alcohol consumption, coffee consumption, exercise amount, and smoking history (Table S25). Alcohol consumption was a significant causal exposure for ten phenotypes, coffee consumption for three phenotypes, exercise amount for six phenotypes, and smoking history for two phenotypes. Coffee consumption was also a significant causal exposure for three anthropometric measurements: body fat mass, visceral adipose tissue area, and waist circumference.
Comparison of the phenotype-phenotype pairs between PheWAS-driven versus EHR-driven. "Penetrance" in genetics is the proportion of those individuals carrying a certain genetic variant who also exhibit the associated phenotype, while "expressivity" measures the proportion of individuals that are carriers of a certain variant and show the associated phenotype to a certain extent 27 . As an indirect method to investigate the penetrance or expressivity of the significant loci identified in our study, we repeated bipartite phenotype network construction using an electronic health records (EHR)-driven method in H-PEACE cohort. Among the phenotypes used in PheWAS analysis, 76 phenotypes were also recorded for this cohort. PheWAS-driven pairs (1164 pairs) were selected based on shared SNPs with association P < 1 × 10 -4 , and EHR-driven pairs (1938 pairs) were selected based on correlation analysis with multi-test corrected P < 0.05. We compared these phenotype-phenotype pairs (Table S26) and evaluated the overlap or exclusiveness of the pairs for each phenotype. Of the 1164 pairs identified in the PheWAS-driven approach, 834 (71.65%) also manifested significance in the EHR-driven analysis. As shown in Fig. 4 and Table S27, high ratios of overlap were identified for skeletal muscle mass (95%) and alkaline phosphatase (93.48%), and low ratios for thyroid cancer (0%) and alpha fetoprotein (8%). When viewed in terms of biological category, the highest average % replication was obtained for anthropometric measurement (86.43%).

Meta-analysis of PheWAS from Korean and Japanese populations.
We performed a PheWAS meta-analysis by incorporating our data with the BBJ data (Japanese population). The results are given in Table S28, Figs. S16 and S17. All 51 phenotypes used in the meta-analysis had an increased number of significant variants in the Korean population, while 37 phenotypes had variants uniquely significant in the meta-analysis. Furthermore, height, diabetes and body mass index had more than 100 variants that were uniquely identified as significant in the meta-analysis.

Discussion
With the advancements in healthcare research that are being driven by big data, increasing efforts are being made to carry out data-wide association studies. PheWAS is one of the tools in that paradigm. However, previous studies faced major challenges in terms of deep phenotyping due to generally using ICD codes, which have limited clarity in their definitions; making the results robust by expanding its application; and the characteristics of population genetics, being highly affected by race and ethnicity. Here, we carried out PheWAS in a Korean population using comprehensive health check-up data linked with genotype data, and furthermore aimed to derive the biological meaning by performing secondary analysis of the PheWAS results. We also compared the results of PheWAS studies conducted in different populations to evaluate trans-ethnic differences. Finally, our bipartite phenotype network analysis of phenotypes using shared genetic association revealed hidden patterns between phenotypes. The deep phenotypes we used in our studies were corroborated during comprehensive health check-up by various confirmatory methods such as laboratory tests, endoscopy, CT scans, MRI, interview questionnaires, and so on. For each participant, all tests were done in the same institute and on the same day. This process of generating deep phenotypes makes for data quality that is well controlled and consistent when compared to results from phenotypes based on ICD codes, which can be discrepant with actual clinical diagnoses due to biases in billing pattern 28 . As we were able to use the raw data produced by the test, our analysis included a lot of endophenotypes. Endophenotype (intermediate phenotype) is a quantitative biological trait 29 that is reported to reliably reflect the function of the categorical biological system 29,30 and has reasonable heritability 31 . As such, an endophenotype could be more closely related to the genetic basis and cause of a clinical trait than would be a broad clinical phenotype such as an ICD code 32 .
We compared our PheWAS results with studies done in European (UK Biobank) and Japanese (Biobank Project Japan) populations and found several novel loci, replicated loci, replicated phenotype-phenotype pairs. We furthermore compared estimated heritability among the populations. Significant variants in the Korean population were partly replicated in both European and Japanese populations, though the replication rate was higher in the Japanese population. We also identified SNP-phenotype associations that were unique to the Korean population when compared to not only the European but also the Japanese population. Noticeably, in the comparison of significant variants associated with body mass index (BMI), the Korean population had novel unique variants (Fig. S8) associated with TERF2IP, ATRNL1, and BANF1. The results from these trans-ethnic and trans-nationality comparisons seemingly emphasize the importance of considering genetic differences among ethnicities, and also race. Koreans are generally included in the East Asian population; however, study of the human Y-chromosome 33 suggests that compared to other populations from Asia, the Korean population has characteristics of a distinct, mostly endogamous ethnic group, and living in a confined peninsula area has preserved these monogenic nationality traits. In a study comparing genetic structure and divergence among Han Chinese, Japanese, and Korean populations those three East Asian populations were shown to have distinct genetic makeup and could be distinguished based on their genetic characteristics 34 . In the meta-analysis of our population and the Japanese population, 72.5% of phenotypes had variants that were uniquely significant in the metaanalysis. Our study shows that the common and exclusive genetic associations of phenotypes should be taken into consideration when performing a population-based clinical study. Furthermore, meta-analysis of PheWAS studies in populations of the same ethnicity but different nationalities can discover uniquely significant variants.
In the comparison of the estimated heritability among different populations, the heritability in the Korean population of biomarkers for coagulation function, such as PT and aPTT, showed similar trends with that of the Japanese population, but manifested relatively high heritability when compared to the UK population. This Figure 4. Comparison of phenotype-phenotype pairs between PheWAS driven and EHR-driven analysis. There were 76 phenotypes also recorded in the EHR-driven cohort (H-PEACE cohort). PheWAS-driven pairs (1164) were based on shared SNPs with association P < 1 × 10 -4 , and EHR-driven pairs (1938) on correlation analysis with multi-test corrected P < 0.05 (Table S26). Skeletal muscle mass (95%) and alkaline phosphatase (93.48%) had high ratios of overlap, while thyroid cancer (0%) and alpha fetoprotein (8%) had low ratios. In terms of biological categories, the average replication % was highest for anthropometric measurement (86.43%) Of the 1164 pairs from the PheWAS-driven approach, 834 (71.65%) also manifested significance in the EHR-driven analysis. www.nature.com/scientificreports/ indicates that the contribution of genetic variants to variation in coagulation traits is affected by ethnical differences. Evaluating heritability difference by ethnicity will be important supportive information in the development of drugs as an aspect of precision medicine. We also leveraged the cross-phenotype association results to provide a panoramic overview of the network connections among multiple phenotypes and genetic variants. Specifically, we generated a phenotype-genotype network focused on metabolic syndrome (Fig. 3A). Metabolic syndrome is a cluster of metabolic abnormalities that are known to be associated with visceral adipose obesity 35 . A large number of epidemiological studies have been conducted on metabolic syndrome because it is a crucial target for healthcare, imposing an increased risk of developing conditions such as cardiovascular disease 35 , malignant disease 36 , depression 37 , and metabolic disease 35 . Early diagnosis is important to prevent the negative consequences of metabolic and this may be done by modifying the lifestyle and risk factors. The network we constructed provided a rationale for defining metabolic syndrome by phenotypes of TG, HDL, hypertension, diabetes, and WC, and for using the characteristics of metabolic syndrome to collectively integrate heterogeneous and complex disease status. The network included phenotypes of cardiovascular disease (coronary calcium score, cardiac ischemia, brain atherosclerosis, malignant disease (thyroid cancer, gastric cancer), and depression and metabolic disease (fatty liver, uric acid), which are known to be complications of metabolic syndrome. Other phenotypes in the network related to obesity, specifically visceral obesity indicator and visceral fat amount; obesity is a well-known cause of metabolic syndrome 35 . Furthermore, lifestyle factor phenotypes such as alcohol consumption, smoking habit, and exercise amount were also part of the network. These suggest modifiable targets for preventing the complications of metabolic syndrome. Finally, the network suggested hub genes associated with metabolic syndrome. Similar network analysis of PheWAS results might provide genotype-based evidence of connections among phenotypes or variants, which to date have been assumed from epidemiological research, and can also provide novel insights into connections that have not been previously reported or recognized.
We additionally used the bipartite phenotype network to perform cross-phenotype mapping. Table S23 shows the cross-phenotype mapping constructed for tumor markers. Tumor markers are highly used in clinical practice for tasks such as oncological screening and monitoring recurrence after treatment. The marker carcinoembryonic antigen (CEA) is recommended by the National Comprehensive Cancer Network (NCCN) guidelines for colon cancer and American Society of Clinical Oncology (ASCO) to test a diagnosis of colon cancer as a baseline for monitoring and then to regularly monitor for recurrence or metastasis of the colon cancer 38,39 . Testing for the marker PSA is recommend by the American Cancer Society (ACS) for men aged > 50 years, after an informed decision-making process 40 . Regular testing for another marker, serum alpha-fetoprotein (AFP), is recommended by the NCCN guideline in the follow-up of hepatocellular carcinoma 41 .
However, while testing for tumor markers is essential in the surveillance of malignant disease, their usage faces problems in the form of low sensitivity and specificity and the potential that they could be affected by factors other than the cancer itself. Thus, providing a cross-phenotype mapping for tumor markers could support an oncologist in interpreting the results of each tumor marker test. For instance, hemoglobin was included in our CEA cross-phenotype mapping. Thus, if a colorectal cancer patient has severe anemia, we should be cautious about interpreting a change in CEA; the anemia could attenuate or exaggerate its reflection of the patient's cancer status 42 . There are several reports that have used not only one tumor marker but a combination of tumor markers to monitor malignancies [43][44][45] . In Table S23, each tumor marker has pairs with multiple other tumor markers, which provide supporting evidence for combining tumor markers as a means to improve their utility in malignancy surveillance.
We also built cross-phenotype mappings for environmental factors. Figure 3C shows the cross-phenotype mapping for coffee consumption in particular. Similar visualization of the correlations between environmental factors and other phenotypes could provide insight into which disease should be considered for the investigation of the benefits or hazards of given environmental factors, and what also connections could provide a candidate model for gene x environment interactions.
In our study, we applied Mendelian randomization analysis to cross-phenotype networks in order to generate corresponding causal inference networks. To the best of our knowledge, this is the first approach to utilize MR in network-based analysis. MR enables the estimation of causal inference by evaluating the relationship between genetic susceptibility to the causal factor and the outcome in question 22 . As shown in Fig. 3D, we specifically drew a causal inference map for skeletal muscle mass. We visualized this map because skeletal muscle mass is regarded as an endocrine and paracrine organ, and is also suggested as a marker in diseases such as metabolic syndrome, diabetes, and more 46 . The analysis revealed skeletal muscle mass as having significant causal inference for metabolic syndrome. In the network, Skeletal muscle mass had six bidirectional associations. Bidirectional association means the "A" phenotype could cause "B", and at the same time "B" phenotype could cause "A", whether is in forward or reverse way 47,48 . Skeletal muscle mass had a bidirectional association with pulmonary function (FEV1, FVC). There are several epidemiological studies for this association [49][50][51] . In one of the studies, individuals with reduced skeletal muscle mass amount have caused a decrease in FVC and FEV1, because of weakened ability to inflate and deflate their lungs 49 . In another study, patients with chronic obstructive pulmonary disease (COPD) are a risk factor for skeletal muscle atrophy by complex combination of various pathophysiological alteration leading to suboptimal muscle work 50 . Though the effect of sarcopenia on pulmonary function is mainly emphasized in clinical practice, muscle recovery measures in poor pulmonary function patients should also be well understood. By the information provided by the bidirectional association network, it will raise alerts for researchers to focus on the reversed direction of causality, which is not well reported, by referring to our results. Thus, by performing MR, we can suggest which phenotype could be causal or an outcome in relation with a trait and also begin to elucidate the mechanism or pathophysiology for a disease of interest.
There are a couple of limitations for the Mendelian randomization (MR) analysis results. First, since there was high dimensional degree of significant association pairs (1767 significant pairs by threshold FDR < 0.05), we were www.nature.com/scientificreports/ not able to provide externally replicated analysis results. For increased confidence, all the significant results in the networks warrants further examination and replication in an external cohort. However, it was difficult to find a cohort with various deep phenotypes, especially those uniquely measured in comprehensive health check-up. For instance, in Fig. 3D, the core phenotype of the network is skeletal muscle amount. Though there are genetic studies regarding sarcopenia52, we couldn't find any cohort that simultaneously had skeletal muscle mass summary statistics as well as pulmonary function or bone mineral density test summary statistics results. But, as our study utilized a comprehensive health check-up cohort, we were able to evaluate the association between various phenotypes from apparently unrelated body systems. Second, there were several MR associations, difficult to be explained by currently reported epidemiological studies. We were able to provide epidemiologic evidence for some of the associations in Fig. 3D, such as skeletal muscle amount with metabolic syndrome46; pulmonary function49-51; liver function53; and bone mineral density54. But there were several MR associations, which are not biologically explainable. For instance, causal associations between PSA and gastric cancer; right and left intraocular pressure; exercise and hepatitis C virus carrier; skeletal muscle amount and alcohol consumption. There should be several reasons for these findings. First, in the analysis, end-phenotypes, endo-phenotypes and environmental phenotypes are altogether incorporated in the MR analysis. Gene by environmental interaction was not considered in the analysis. Second, our cohort is from the health check-up cohort, which could contain samples with positive morbidity relatively few. Third, we did not have an external validation population but analyzed in one sample population. This could have led to spurious correlations between phenotypes that are unrelated to genetics. These limitations should be further analyzed by performing gene by environment interaction; future analysis in a larger set of study population; and in an external validation study population. Though the MR analysis results are not validated in other cohorts and several significant associations were not biologically reasonable, it could raise a necessity to validate certain associations by other researchers to focus on certain pairs of phenotypes and organize a cohort for that purpose. In our study, we provide an openly accessible web-based phenotype-phenotype network for the whole dataset (Fig. S13), which is an interactive visualization tool of the phenotype-phenotype network (https:// hdpm. biome dinfo lab. com/ ddn/ genie). This tool will allow other researchers to easily access our results and pick up where the unreported associations or limitations are and build up future research.
Our study has several advantages. First, to the best of our knowledge, this is the first PheWAS study performed in the Korean population. As described above, several loci in this population differ from the Japanese population. We were also able to carry out trans-nationality analysis for the PheWAS. Second, we defined phenotypes directly using results from health check-ups and questionnaire responses from personal participants. This makes the resolution, clarity, and reliability of this study's results better than those of a PheWAS based on ICD codes or personal self-reports. These billing codes or personal memory can bring an underlying bias into the data registry. So, the phenotypes are objectively and precisely defined by these check-ups. Third, since all tests were performed in the same institute, under the same conditions, and by using the same machines, protocols, and chemicals, the produced data is consistent and its quality is highly controlled. Fourth, we performed secondary analysis of the PheWAS results in ways to derive the biological meaning, so that the results could be highly applicable and utilized more practically. We constructed a phenotype-phenotype network using all the phenotypes in our study (Fig. S13). Similarly constructing a phenotype-phenotype network based on comprehensive, deep phenotypes could provide clinicians and researchers with a detailed landscape of the interconnections between phenotypes and enable better understanding of their underpinnings. Furthermore, the phenotype-phenotype network not only includes disease status but also contains information on genes, environment, and lifestyle. Precision medicine pursues prevention and treatment strategies that take individual variability 1 , such as in genes, environment, and lifestyle, into account 2 . Accordingly, the networks generated by PheWAS would provide fundamental information for realizing precision medicine. Fifth, we provide summary statistics, which are significant. This will help other researchers to explore the phenotypes for making headway in further study.
Our study has several limitations. First, we did not have a set of Korean replication population because it was not possible to find such datasets with the variety of deep phenotypes incorporated in our study. However, we instead introduced the UKBB and BBJ as replication sets, and consequently identified multiple replicated loci. We also replicated the phenotype-phenotype pairs using a larger EHR-driven database of Korean samples to investigate whether the genetic connection was reflected at the actual phenotype level. Second, the study population was collected from those who had regular health check-ups, and therefore samples with positive morbidity were relatively few. Accordingly, the significance of the loci was low for some phenotypes. We tried to overcome this lack of statistical power by performing a meta-analysis with the UKBB and BBJ summary statistics, in which we were able to pick up additional significant loci. In a future study, we will incorporate diverse disease cohorts from the Korean population to increase the study power. Third, phenotype-phenotype networks were constructed from a single sample, because it was not possible to find an external set of population. This might have led to spurious correlations between phenotypes that are unrelated to genetics. Further, for most of our analysis we used suggestive p-value cutoff of 10 -4 , and there were 260,922 variants that passed the threshold across the phenotypes. However, if FDR < 0.05 cutoff is considered, the number is lower with 114,677 passing the threshold (Table S5), which could have led to inclusion of more false positive associations. Forth, in network analysis, it was based on a permissive p-value threshold, which can be associated with false-positive associations.
In conclusion, our study highlights the capacity for understanding the biological insights post-PheWAS by comprehensively exploiting the results. With the information generated by PheWAS, we attempted to provide a landscape that integrated an individual's genetic, lifestyle, and environmental factors along with health status. We provided several samples of actionable applications such as constructing a gene-phenotype association network related to metabolic syndrome; constructing cross-phenotype mappings; and visualizing causal inference mappings. Through analysis in the context of differences in ethnicity and nationality, our study shows that some phenotypes are common or exclusive in their genetic associations, and this should be taken into consideration www.nature.com/scientificreports/ when performing a population-based clinical study. The paradigm of PheWAS suggested in our study will eventually be the cornerstone for applying the core concepts of precision medicine to research and healthcare practice.

Data availability
Complete summary statistics of the GENIE cohort are not publicly available due to restrictions (institutional policy to protect the privacy of research participants), but are available from the corresponding author on reasonable request. However, all other data are contained in the article and its supplementary information or are available upon reasonable request. The summary statistics from UK biobank and Biobank Japan are available at from the following URLs: http:// www. neale lab. is/ uk-bioba nk/ and http:// jenger. riken. jp/ en/ result. The codes used in analysis for this paper is available at https:// github. com/ dokyo onkim lab.