An exploratory phenome wide association study linking asthma and liver disease genetic variants to electronic health records from the Estonian Biobank

The Estonian Biobank, governed by the Institute of Genomics at the University of Tartu (Biobank), has stored genetic material/DNA and continuously collected data since 2002 on a total of 52,274 individuals representing ~5% of the Estonian adult population and is increasing. To explore the utility of data available in the Biobank, we conducted a phenome-wide association study (PheWAS) in two areas of interest to healthcare researchers; asthma and liver disease. We used 11 asthma and 13 liver disease-associated single nucleotide polymorphisms (SNPs), identified from published genome-wide association studies, to test our ability to detect established associations. We confirmed 2 asthma and 5 liver disease associated variants at nominal significance and directionally consistent with published results. We found 2 associations that were opposite to what was published before (rs4374383:AA increases risk of NASH/NAFLD, rs11597086 increases ALT level). Three SNP-diagnosis pairs passed the phenome-wide significance threshold: rs9273349 and E06 (thyroiditis, p = 5.50x10-8); rs9273349 and E10 (type-1 diabetes, p = 2.60x10-7); and rs2281135 and K76 (non-alcoholic liver diseases, including NAFLD, p = 4.10x10-7). We have validated our approach and confirmed the quality of the data for these conditions. Importantly, we demonstrate that the extensive amount of genetic and medical information from the Estonian Biobank can be successfully utilized for scientific research.


Introduction
Genetic data are an important resource for scientific research and potential drug target identification [1] and genome-wide association studies (GWAS) have identified many disease-associated genetic variants [2]. Complementary to GWAS are phenome-wide association studies (PheWAS) which use a genotype-to-phenotype approach, testing for associations between specific genetic variants over a wide spectrum of phenotypes [3]. Combining genetic data with phenotypes defined and validated in electronic health records (EHRs) permits associations between genetic variants and disease outcomes, including diagnoses and procedures not commonly found in GWAS studies. To date, despite a relative abundance of EHR and genetic data becoming available, few large scale PheWAS studies have linked these data [4]. Understanding the full range of associations along with understanding the functional mechanisms of causal genetic variants will have important implications for the design of novel therapies across indications to reduce morbidity and mortality.
The Estonian Biobank, governed by the Institute of Genomics at the University of Tartu (Biobank) and launched in 2007, aimed to create a biobank with biological samples and genetic material with linkage to EHRs to investigate the genetic, environmental and behavioral background of common diseases in the Estonian population [5]. To explore the utility of data in the Biobank, we conducted a PheWAS in two areas of interest to healthcare researchersasthma and liver disease-to determine whether established and novel disease associations together with early disease indicators could be detected.
Asthma is a chronic inflammatory disease affecting the airways and lung characterized by recurring symptoms including breathlessness and wheezing caused by inhalation of environmental particulates such as allergens. Asthma requires both acute and long-term management for alleviation of symptoms which vary in severity between individuals. It is estimated that 235 million people worldwide suffer from asthma [30]. A number of asthma associated SNPs have been identified by GWAS and include variants near TSLP, BIRC3, IL18R1,  HLA-DQB1, IL33, GSDMB, GSDM1, IL2RB, SLC22A5, IL13, and RORA gene regions [31][32][33].

Objectives
In an effort to assess the utility of the genetic and phenotypic data available in the Estonian Biobank, our first two objectives were to show our ability to detect previously reported GWAS associations of (1) asthma-associated SNPs with biomarkers (lab measures) and clinical diagnoses of asthma, and (2) liver disease-associated SNPs and clinical diagnoses of NAFLD/ AstraZeneca provided support in the form of salaries for author GJ, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. STACC and Quretec provided support in the form of salaries for authors SR and JV, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. GlaxoSmithKline provided support in the form of salaries for authors NG, DW and MA, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. Pfizer provided support in the form of salaries for author AKL, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. European Regional Development Funds for the Center of Excellence of Estonian ICT research EXCITE and Estonian Research Council grant IUT34-4: These funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The funders provided support in the form of salaries for authors JV, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the 'author contributions' section.
NASH within the Estonian Biobank. Both objectives are first, to test the Biobank suitability for validation/replication studies, and second, to confirm no systematic errors in the Biobank data before moving on to objective three, to conduct a PheWAS to test the association of the selected 24 SNPs with all other ICD-10 clinical diagnoses. Where significant associations between SNPs and ICD-10 diagnostic codes are found, our fourth objective was to determine associations of corresponding SNPs with lab/biomarker measures (quantitative traits) prior to the date of diagnosis.

Database
The Estonian Biobank has had continuous data collection since 2002 on a total of 52,274 voluntary adult individuals, of which 50,000 are currently active (alive and data continuously updated) in the database, representing approximately 5% of the Estonian adult population. The Biobank has the ability to access Estonian primary care and secondary/inpatient care EHRs via linkage to a central e-health database where EHR information is sent and uploaded by healthcare professionals. Additionally, participant information is updated through linkage to other databases/registries including the Estonian population registry, cancer registry, and death registry. This linkage ensures prospective and retrospective capture of patient and disease characteristics in the form of ICD-10 codes [5], demographics, lifestyle information, laboratory measurements, biomarker measurements, diagnoses and drug utilization (prescriptions and fill data). Additional information is captured by questionnaires, including lifestyle data e.g. smoking and alcohol consumption. Questionnaire information is usually captured during patient enrolment; however, where new questionnaires are introduced, existing Biobank participants are invited to complete these to increase completeness of data. Number of participants in the Biobank with genetic data available is expected to increase over 150K by the end of 2019.

Case ascertainment
To be eligible for the study, patients were required to have genetic and EHR data linkage for at least one year (365 days) after recruitment. Patients were excluded if age, gender or genetic information was missing or if patients had no EHR linkage. ICD-10 diagnostic codes were used to determine clinical diagnoses.
Where primary SNPs of interest were not available on this array, the Broad Institute SNP Annotation and Proxy Search (SNAP) tool was used to identify proxy SNPs (SNPs which can represent the primary SNP of interest) with a minimum linkage disequilibrium r 2 �0.6 [39]. After these steps, 11 asthma (1 proxy) and 13 liver disease-associated (8 proxies) SNPs of interest were included in our analyses ( Table 1). For 4 of these SNPs, associations with one or more lab measurements (6 in total) were reported.

PheWAS codes
ICD-10 diagnostic codes were mapped to "PheWAS codes" using methodology from Neuraz et al. [40]. This involved converting individual ICD-10 codes into a higher order of grouped codes for a single disease, e.g. grouping codes A00.0, A00.9, A00.1, A00 for cholera into a three-character code A00 ("PheWAS code") and using broader A00-A09 range (intestinal infectious diseases) as an exclusion criterion for the control group (S1 Table). If a significant association was observed, we repeated the analysis using the individual ICD-10 code to define the case group while keeping the control group criteria unchanged, e.g A00.0 as the case group and A00-A09 as exclusion criteria.

Other covariates, laboratory measures and biomarkers
To describe the patient population, we extracted data on a range of covariates including birth year, gender, ethnicity, smoking status and body mass index (BMI). Laboratory and biomarker measures included sodium, potassium, blood urea, blood glucose, Creatinine, c-reactive protein, haemoglobin, platelet count, B-type natriuretic peptide / brain natriuretic protein, troponin I, troponin T, white blood cell count, neutrophil count, eosinophil count, serum uric acid, fibrinogen, A-fetoprotein, gamma-glutamyl transferase, alanine transaminase (ALT), alkaline phosphatase, albumin, aspartate transaminase, bilirubin, total cholesterol, high-density lipoprotein cholesterol and low density lipoprotein cholesterol.

Data extraction and analysis
Biobank participant information, e.g. lifestyle information, was extracted directly from the Estonian Biobank questionnaires completed during recruitment. Incident diagnoses were extracted from local copy of central e-health database, and laboratory measures were extracted using both e-health database and local copies of two main Estonian hospitals' databases. In total, there were more than 1.6 million diagnoses recorded of 12,845 different ICD-10 codes, which were grouped into 1,888 higher level ICD-10 (PheWAS) codes used in this analysis. There were over 640,000 lab/biomarker measurements available; however, for 42% of the participants no measurements were recorded. Summary statistics for each laboratory measurement can be found in S1 Text. SNP data was called from the Illumina Infinium Global Screening chip array. After variant calling, the data was filtered using PLINK software [41] sample-wise: call rate>95%, no sex mismatches between phenotype and genotype data, heterozygosity within 3 standard deviations of average heterozygosity over all samples to eliminate possible inbreeding and DNA contamination; and marker-wise: HWE p-value>1x10 -6 , callrate>95%, and Illumina Geno-meStudio GenTrain score >0.6, Cluster Separation Score >0.4. SNPs of interest were then extracted from the data. Analysis of this study was conducted in R version 3.4.3 using the Phe-WAS R package [42] customized to allow use of ICD-10 codes. The high-performance computing center at the University of Tartu was used for analysis. Calculations for statistical power can be found in S2 Text.
Logistic regression models were used to evaluate association of SNP/proxy SNP variation with ICD-10 diagnoses. Odds ratios (ORs) and 95% confidence intervals (CI) were estimated. Linear regression models were used to evaluate associations between SNPs and laboratory measurements. Most of the laboratory measurements (15 out of 26) were log-transformed due to positive skewness assessed graphically (see For objective 1, cases for asthma were defined as ever having been diagnosed with ICD-10 codes J45-J46, and controls were defined from the remaining individuals as never (prospectively and retrospectively in the patients' medical records) having been diagnosed with ICD-10 codes J40-J47. Similarly, for objective 2, cases for NASH/NAFLD were defined as having ICD-10 codes K75.8 or K76.0 and controls as individuals never having been diagnosed with ICD-10 codes K70-K77. Full list of J40-J47 and K70-K77 diagnoses with their counts is given in S4 Text. Altogether, 21 previously reported SNP-disease and 6 SNP-lab measurement (covering 4 SNPs out of 24) association validation tests were performed in the afore-mentioned objectives. For objective 3 (PheWAS), each PheWAS code was tested individually, and only PheWAS codes with at least 20 cases and 20 controls were included in the analyses. To explore what drives the detected significant associations, an exact ICD-10 code instead of PheWAS code was used in the subsequent analysis.
A p-value significance threshold 0.05 was used for objectives 1 and 2 (confirmation of liver and asthma SNP associations). For objective 3, a p-value of 2.0x10 -6 (�0.05/1,000/24 where 1,000 is the effective number of ICD-10 diagnoses tested and 24 is the number of SNPs tested) was used to reduce the likelihood of chance associations and identify the most prominent differences. To search for evidence of systematic bias, a QQ-plot of PheWAS p-values was used to evaluate whether the observed distribution was different from what would be expected under the null hypothesis.

Protection of human subjects
Research at Estonian Biobank is regulated by Human Gene Research Act and all participants have signed a broad informed consent. IRB approval for current study was granted by Research Ethics Committee of University of Tartu, approval nr 236/T-23.

Results
After applying the exclusion criteria, 26,766 (51.2%) Caucasian individuals remained from the original 52,274 ( Table 2). The main reason for the drop in the number of samples is missing genotype data (not genotyped with the given genotyping array). In this study, most participants were female (71.8%), 41.6% of individuals had a body mass index (BMI) of between 18.5-25 (normal weight), 59.2% were never smokers and 75.9% were of Estonian nationality ( Table 3).

Preselected SNPs and risk of asthma and fatty liver disease
Two asthma-associated SNPs (rs11071559 (RORA) and rs1837253 (TSLP)) passed the significance threshold of <0.05 ( Table 4) and were directionally consistent with previous studies. While Moffat et al. showed that rs1837253 significantly reduces severe asthma in one of their datasets, it did not replicate in the other and the association direction with asthma was the opposite (increase) in the full dataset [31]. We also observe the increased effect of allele C of rs1837253 in our data. From the lab measurement tests, we confirmed that each additive effect allele T of rs2846848 (BIRC3) had a significant effect on reducing neutrophil levels ( Table 4). All other associations were non-significant.
Five fatty liver disease-associated SNPs (rs780094 (GCKR), rs2281135 (PNPLA3), rs8418 (PNPLA3), rs58542926 (TM6SF2), rs2980875 (TRIB1)) passed the significance threshold of <0.05 ( Table 5) and were directionally consistent with previous studies. Additionally, we found that the AA genotype of rs4374383 (MERTK) significantly increases the chance of developing NASH/NAFLD, which is the opposite direction of effect to what was shown by Patin et al. [29]. For rs9992651 (HSD17B13) and rs11597086 (ERLIN1) we also assessed the association to ALT levels. In our analysis, both SNPs are significantly associated with increased ALT levels, but rs11597086 was in the opposite direction of what has been shown by Yuan et al. [21]. We could only replicate some previously reported associations with asthma or fatty liver disease-associated SNPs, but this is likely due to having only a small number of cases/measurements for each of these conditions (Tables 4 and 5) and limitations of EHRs and ICD-10 codes which might not always result in clear signals for all diseases (see Discussion). However, while many of the associations we attempted to confirm were not strong enough in our analysis to pass the significance threshold, almost all of the effect directions were concordant between our study and the original studies, and the QQ-plot of p-values from our confirmatory analysis showed clear enrichment of small p-values (Fig 1). We can consequently expect the Estonian Biobank data to be suitable for the PheWAS of these SNPs, though we might lack sufficient power to detect modest association signals for specific diseases.

PheWAS association between pre-selected SNPs and ICD-10 codes
We conducted a PheWAS to test the association of 24 investigated asthma and liver disease SNPs with all ICD-10 clinical diagnoses. Three SNP-diagnosis pairs passed the PheWAS significance threshold: rs9273349 (HLA-DQB1) and E06 (thyroditis); rs9273349 (HLA-DQB1) and E10 (type-1 diabetes) ( Table 6, Fig 2); and rs2281135 (PNPLA3) and K76 (non-alcoholic liver diseases, including NAFLD) ( Table 6). The QQ-plot of all the association p-values (Fig B in S2  Text) does not indicate any systematic biases in our study (e.g. due to population stratification) as some inflation is expected due to analyzing known disease-associated SNPs.
In order to explore what exact diagnoses drive the detected associations, we repeated the analysis with the individual ICD-10 diagnosis codes instead of PheWAS codes ( Table 7). We found that rs9273349 is associated with autoimmune thyroiditis (E06.3) and insulin-dependent diabetes mellitus without complications (E10.9), with ophthalmic complications (E10.3) and with multiple complications (E10.7). rs2281135 is associated with fatty (change of) liver, not elsewhere classified, including NAFLD (K76.0).

Associations between significant PheWAS results and laboratory measurements
Using the rs9273349-E06, rs9273349-E10 and rs2281135-K76 pairs, the lab/biomarker measures (quantitative traits) closest in time prior to diagnosis were identified but no significant effects were observed under a significance threshold 9.6x10 -4 (�0.05/2/26 where 2 is the number of distinct SNPs and 26 is the number of lab/biomarker measures) due to very small sample sizes after discarding all measurements on patients without the underlying diagnosis. Using all the measurements regardless of whether a patient had been diagnosed with a disease or not, only rs9273349 (HLA-DQB1) passed the PheWAS significance threshold for decreasing average (over all patient's measurements) levels of cholesterol per addition of an effect allele (p = 1.1x10 -6 ) (Fig 2).

Discussion
This is the first PheWAS study using Estonian Biobank data linked to EHRs. Large scale Phe-WAS are scarce [44] with other PheWAS utilising mostly small cohorts [42,[45][46][47][48][49][50][51][52]. Recently large scale PheWAS have become possible in the UK Biobank, but other biobanks will still be required for further validation or replication. To assess the Estonian Biobank's suitability for such a study or association validation/replication purposes, we focused on asthma and liver disease associated SNPs only and as a first task, investigated whether the effects previously reported in the literature (GWAS) could also be detected in our data. That is also to confirm that the data have no systematic errors and are suitable for PheWAS analysis. We replicated previous GWAS results, reporting a significant association of TSLP and RORA gene variants with asthma [31,33] (Table 4) and GCKR, PNPLA3, TRIB1 and TM6SF2 gene variants with the risk of developing liver diseases, notably NAFLD/NASH [8,[11][12][13][14][15]20,[24][25][26]29,53] (Table 5).
We observed decreasing neutrophil levels per addition of an effect allele in the BIRC3 gene variant, consistent with previously reported results. In addition, variants in HSD17B13 and ERLIN1 influence alanine transaminase (ALT) levels. The HSD17B13 association is consistent with a recent study reporting that a loss-of-function variant associated with decreased levels of ALT and aspartate aminotransferase (AST) and reduced the risk of liver disease and progression from NAFLD to NASH [22]. Our observed association between rs11597086 (ERLIN1) and ALT level is in the opposite direction of what has been previously reported [21]. The limited replication of our results with previously published GWAS results in asthma and liver disease is likely due to smaller case sizes in our data than in the original studies, leading to low power, but with continuous enrolment and data collection, statistical power will improve in time or could be enhanced by meta-analysis with other studies.
From the PheWAS, we identified the association of asthma-associated genetic variant rs9273349 with type 1 diabetes, and autoimmune thyroiditis. rs9273349 is a variant of the major histocompatibility complex, class II, DQ beta 1 (HLA-DQB1) gene which has an important role in the immune system. HLA-DQB1 is anchored to the cell surface membrane and functions to present extracellular proteins into the cell [54]. A number of studies have identified the association of HLA-DQB1 with thyroiditis (notably autoimmune Hashimoto's thyroiditis) and type 1 diabetes [55][56][57][58][59]. Recently, a study by Verma et al. detected an association between acquired hypothyroidism (usually caused by thyroiditis) and rs17843604, a SNP in the same HLA-DQA1/B1 region [44]. Thyroiditis is a complex immune disorder of unknown aetiology where an infiltration of T and B lymphocytes occurs as a reaction to thyroid antigens. These B and T lymphocytes then produce thyroid autoantibodies resulting in clinical hypo-or  hyper-thyroidism [55]. Type 1 diabetes is also a T lymphocyte driven disease which results in the destruction of insulin producing pancreatic islet cells. As suggested in [31,60,61], the amino acid variation of HLA-DQB1 variants could cause differential and incorrect binding of peptides, altering cellular sensitization resulting in cellular hypo-or hyperactivity, disrupting homeostatic immune function. Additionally, Mosaad et al. observed microalbuminuria in type 1 diabetic patients, suggesting that alterations in HLA-DQB1 expression effect homeostatic regulation [59]. This study has a number of strengths which include; a large sample size of ICD-10 diagnoses and quantitative measures linked to 26,766 genotyped individuals, relatively long follow-up time, integration of questionnaire data and linkage of EHR data with laboratory and other databases. However, there are a number of limitations to consider. The large proportion of females may introduce bias and reduce generalizability to the general population. However, this is a general problem of all voluntary biobank cohorts as women tend to enroll more actively than men [5]. Furthermore, EHR-linked data are not collected for research purposes, and without strict standards for data collection and format, the quality of the data may vary broadly (missing data, non-coherent format, errors on data insertion etc.). Even when an association between a SNP and an ICD-10 diagnostic code or biomarker/lab value is clearly identified, this does not imply a causal variant. ICD-10 coding may be a limitation itself as in some  An exploratory PheWAS linking asthma and liver disease genetic variants to EHR from the Estonian Biobank cases it can be difficult to ascertain the exact condition. The limitation of using a biomarker/lab measure closest to diagnosis is that it may not be the most sensitive or specific predictor of disease. Additionally, these measures do not take into account potential confounders e.g. response to medication, age and time-period effects. Furthermore, small sample sizes for specific conditions limit statistical power. For example, many diseases are rare and due to the fine granularity of ICD-10 codes, using the exact codes results in very small sample sizes for case groups. With the small number of cases, the power to detect an association is limited. Therefore, it is helpful to conduct a two-step PheWAS by first using higher level diagnosis codes to screen for any association signals, and subsequently evaluating the more detailed codes to understand where specifically the association is coming from. In our study, after detecting the association between rs9273349 and E06 (thyroiditis), an exact diagnosis analysis revealed that the association was driven by E06.3 (autoimmune thyroiditis). Similarly, the association between rs2281135 and K76 "Other diseases of liver, non-alcoholic fatty liver disease (NAFLD)" is actually driven by more specific code K76.0 "Fatty (change of) liver, not elsewhere classified".
In conclusion, this is the first PheWAS conducted using the Estonian Biobank demonstrating the extensive amount of genetic and medical information which can be successfully utilized for scientific research. The Estonian Biobank has 50K participants, all have signed informed consent, allowing regular data update from all health databases throughout their lives. Notably, its value is increasing over time-not only because of the continuous EHR data addition and vast amount of genetic data available, but also the participation count of Estonian Biobank (with genetic data available) is expected to rise to 150K by the end of 2019. We showed that Biobank data can be effectively used as a validation/replication database. We replicated 9 GWAS associations for asthma and liver disease-associated SNPs and found the opposite effect directions for 2 associations (rs4374383 increases the risk of NAFLD, rs11597086 increases ALT level). Furthermore, this PheWAS exploring the association of 11 asthma-associated and 13 liver disease-associated SNPs with other diseases and biomarkers is one of few studies to use ICD-10 diagnostic codes to link genetic data with EHRs. Although we did not detect any novel associations in our study, we were able to confirm 3 phenome-wide significant associations based on our data-rs9273349 and thyroiditis, rs9273349 and type-1 diabetes, rs2281135 and non-alcoholic liver diseases, including NAFLD. Considering also the continuous addition of the new data, this highlights the usability of Estonian Biobank for using it effectively as a validation database and conducting extended PheWAS studies with much larger sets of SNPs in the future.
Supporting information S1 Table.