1. Computational and Systems Biology
  2. Epidemiology and Global Health
Download icon

Linking glycemic dysregulation in diabetes to symptoms, comorbidities, and genetics through EHR data mining

  1. Isa Kristina Kirk
  2. Christian Simon
  3. Karina Banasik
  4. Peter Christoffer Holm
  5. Amalie Dahl Haue
  6. Peter Bjødstrup Jensen
  7. Lars Juhl Jensen
  8. Cristina Leal Rodríguez
  9. Mette Krogh Pedersen
  10. Robert Eriksson
  11. Henrik Ullits Andersen
  12. Thomas Almdal
  13. Jette Bork-Jensen
  14. Niels Grarup
  15. Knut Borch-Johnsen
  16. Oluf Pedersen
  17. Flemming Pociot
  18. Torben Hansen
  19. Regine Bergholdt
  20. Peter Rossing  Is a corresponding author
  21. Søren Brunak  Is a corresponding author
  1. University of Copenhagen, Denmark
  2. Odense University Hospital, Denmark
  3. Steno Diabetes Center Copenhagen, Denmark
  4. Rigshospitalet, Denmark
  5. Holbæk Hospital, Denmark
  6. Herlev-Gentofte Hospital, Denmark
  7. Technical University of Denmark, Denmark
Research Article
  • Cited 0
  • Views 373
  • Annotations
Cite this article as: eLife 2019;8:e44941 doi: 10.7554/eLife.44941

Abstract

Diabetes is a diverse and complex disease, with considerable variation in phenotypic manifestation and severity. This variation hampers the study of etiological differences and reduces the statistical power of analyses of associations to genetics, treatment outcomes, and complications. We address these issues through deep, fine-grained phenotypic stratification of a diabetes cohort. Text mining the electronic health records of 14,017 patients, we matched two controlled vocabularies (ICD-10 and a custom vocabulary developed at the clinical center Steno Diabetes Center Copenhagen) to clinical narratives spanning a 19 year period. The two matched vocabularies comprise over 20,000 medical terms describing symptoms, other diagnoses, and lifestyle factors. The cohort is genetically homogeneous (Caucasian diabetes patients from Denmark) so the resulting stratification is not driven by ethnic differences, but rather by inherently dissimilar progression patterns and lifestyle related risk factors. Using unsupervised Markov clustering, we defined 71 clusters of at least 50 individuals within the diabetes spectrum. The clusters display both distinct and shared longitudinal glycemic dysregulation patterns, temporal co-occurrences of comorbidities, and associations to single nucleotide polymorphisms in or near genes relevant for diabetes comorbidities.

Introduction

Electronic Health Records (EHRs) contain patient characteristics from different data layers including text narratives, assigned diagnosis codes, biochemical values, and prescription data. These data types display a high degree of complementarity, providing an excellent basis for deep phenotyping and patient stratification. Recent studies have shown how structured data derived from EHRs can be used to assess phenotypic variability of different disease areas (Li et al., 2015; Dahlem et al., 2015; Doshi-Velez et al., 2014; Kho et al., 2011; Kho et al., 2012). While the use of structured EHR data in many instances resembles traditional registry- or biobank-based research, the inclusion of unstructured data such as clinical narratives allows for the definition of even more fine-grained phenotypes, which could lead to novel subgroup stratifications (Li et al., 2015; Roque et al., 2011; Miotto et al., 2016).

A vast amount of information on symptoms, lifestyle, complications, and comorbidities is available from clinical narratives in unstructured EHR data. Text mining applying natural language processing (NLP) algorithms is one strategy, but simpler approaches have also been shown to be valuable in the context of clinical text, for reviews see Jensen et al. (2012) and Denny (2012). These methods work across language barriers and have been successfully implemented in for example adverse drug reaction detection (Warrer et al., 2012), subgrouping of chronic obstructive pulmonary disease (Fu et al., 2015), cancer subgrouping (Chen et al., 2015), and classification of epileptic children (Pereira et al., 2013). Such studies show the possibilities of using and integrating different parts of EHRs for matching phenotypically similar subgroups to biomarker data, which is key to improved treatment and characterizing etiological differences.

Several large initiatives have been established for utilizing EHRs, including the Electronic Medical Records and Genomics (eMERGE) consortium of DNA biorepositories that links genetic data with electronic medical records (McCarty et al., 2011; Gottesman et al., 2013), and EMR-driven nonnegative restricted Boltzmann machines (eNRBM) which use unsupervised learning for analyzing EHRs (Tran et al., 2015). Furthermore, other studies have used general approaches for finding direct and inverse comorbidities (Doshi-Velez et al., 2014; Roque et al., 2011; Gligorijevic et al., 2016).

Diabetes Mellitus (DM) is a difficult disease to stratify (American Diabetes Association, 2017). DM covers etiologically different metabolic disorders that exhibit the same phenotype, hyperglycemia, due to either insufficient insulin production relative to insulin demand or insulin resistance. Although DM is classified into different major subtypes, it has been hypothesized to represent a disease continuum rather than strict distinct disease subtypes (American Diabetes Association, 2017; Flannick et al., 2016). One recent data-driven study used five subgroups of adult-onset diabetes and clustered six parameters from the structured data of the EHR (Ahlqvist et al., 2018). DM is a complex disorder associated with several comorbidities and organ complications. These can be classified as macrovascular complications that is cardiovascular disease, and microvascular complications resulting in eye, kidney, and nerve damage. Cardiovascular complications alone are responsible for 50–80% of all-cause mortality in diabetes patients (Laakso, 2001). The severity of complications is affected by glycemic dysregulation, that is increased or fluctuating blood glucose levels (Stratton et al., 2000; UK Prospective Diabetes Study Group, 1998a; UK Prospective Diabetes Study Group, 1998b; Nathan et al., 1993), and successful reduction and prevention of diabetic complications have been observed when the glycemic dysregulation is reduced or removed (Stratton et al., 2000; UK Prospective Diabetes Study Group, 1998a). Therefore, risk factors for glycemic dysregulation are crucial to diabetes progression (Ahlqvist et al., 2015). Known risk factors for complications include age, diabetes duration, polypharmacy, comorbidities (Juarez et al., 2012), increased levels of circulating triglyceride and LDL-cholesterol, and lower levels of HDL-cholesterol (Saudek et al., 2006; Giannini et al., 2011; Bitzur et al., 2009). Finding new risk factors that can help classify poorly regulated versus well-regulated diabetes, such as other biochemical variables or genetic variants, could improve treatment and reduce diabetic complications.

In this study, we utilized the unstructured data of EHRs and performed a deep phenotypic characterization of a Danish diabetes cohort of 14,017 individuals, aged 18 to 101 at the end of the study, using vocabularies comprising both diagnosis codes and ‘exposome’ related terms. We used text-mined and assigned diagnosis codes to stratify the cohort and described it using both physiological and genetic variation data. The unstructured EHR data enabled us to classify patients based on their level of glycemic dysregulation and to identify potential biochemical and genetic markers associated with dysglycemia.

Results

Text mining the EHR corpus

The general aim of the text mining effort was to obtain a richer phenotypic characterization of each patient. Initially, each patient had in 4.9 assigned codes on average. Applying text mining with two vocabularies (ICD-10 and SDC-custom) resulted in a 4-fold increase to 18.6 codes per patient. Moreover, the distribution of codes across ICD-10 chapters changed considerably when adding the text-mined codes, with chapters I, VII, XVIII and XIX showing the largest increases (6, 15, 25 and 22-fold increase, respectively) (Figure 1). This illustrates the difference between the assigned diagnosis codes from the structured data and the much more symptom-rich codes detected by text mining.

Figure 1 with 6 supplements see all
Comparison of distributions of ICD-10 diagnosis codes with and without text mining.

(A) Percentage of diagnosis codes belonging to the different ICD-10 chapters and the relative increase in diagnosis codes from the different chapters when combining the text-mined and assigned codes. (B) Age distributions of text-mined and assigned ICD-10 diagnosis codes from the SDCC corpus divided into the 21 ICD-10 chapters.

Comorbidity clustering based on text-mined and assigned diagnosis codes

For each patient, the assigned and text-mined ICD-10 codes were combined to create a patient-specific diagnosis-vector where the primary diabetes type (E10 or E11) was not included. Contrary to cancer for example, where the ICD-10 diagnoses are quite reliable and highly detailed, the primary codes in a multi-organ disease like diabetes are used in a fuzzier way, as the knowledge on robust diabetes subtypes and their characteristics in the context of comorbidities is quite limited. We do therefore not want the clustering to be driven by the broad, less etiology-relevant primary codes from the endocrinology chapter, but rather by more objectively observed symptoms, other diseases and lifestyle features. Following code-abundance normalization and BM-25 correction the vectors were clustered using MCL producing 172 clusters (mean = 65 patients, min = 11, max = 979, median = 40), in which 11,208 patients (80.47%) were included Figure 2A. The remaining 2720 patients (19.53%) were in clusters with ten or less patients and were therefore omitted from subsequent analyses.

Figure 2 with 5 supplements see all
Phenotypic clusters found in the SDCC cohort.

The clustering was created with diagnosis vectors of 13,928 patients (with text in the record) comprising both text-mined and assigned ICD-10 codes. A total of 172 clusters were created, where 11,208 patients (80.47%) were captured in the clustering (clusters with five or less patients were discarded for statistical reasons). (A) Each node represents a patient within the corpus colored by the association to one of the 172 unique clusters. (B) The 71 clusters with at least 50 patients colored with the same palette as in (A).

Even though codes for the primary diabetes type were not part of the diagnosis vectors, specific clusters were significantly enriched for T1D patients (cluster 1: N = 506, adj. p-value=9.3e-51 and cluster 9: N = 101 adj. p-value=1.2e-10). Other clusters had significantly more T2D patients than expected (cluster 3: N = 233, adj. p-value=9.1e-10, cluster 5: N = 170, adj. p-value=3.8e-13 and cluster 6, N = 158 adj. p-value=8.4e-17). In addition, we observed a cluster significantly enriched with the ICD-10 term E13: other diabetes (cluster 25, N = 93, adj. p-value=1.8e-142), which includes diabetes due to genetic defects, post-pancreatectomy diabetes and post-procedural diabetes. Several other clusters had a mix of T1D and T2D patients according to the assigned codes. Further characteristics of the laboratory data and prescription data as well as the clusters regarding sex, age, observational time, years with diabetes etc. can be found in Supplementary files 13 and in Figure 1—figure supplement 1, Figure 1—figure supplement 2, Figure 1—figure supplement 3, Figure 1—figure supplement 4, Figure 2—figure supplement 1, Figure 2—figure supplement 2, Figure 2—figure supplement 3, Figure 2—figure supplement 4. The robustness of the clustering was found to be high (see description in Materials and methods and Figure 2—figure supplement 5). To maintain power in subsequent analyses we focused on clusters with at least 50 patients (71 clusters comprising 8652 patients, Figure 2B).

Enriched comorbidity and symptom patterns in diabetes patient clusters

The 71 clusters (Figure 2B) were grouped by hierarchical clustering, using distances obtained from cluster specific symptoms from the ICD-10 chapter XVIII (level 1). Six main groups and an outlier (cluster 70) were found containing 5, 8, 21, 11, 7 and 18 of the original clusters, respectively. The symptom groups are illustrated by the branch colors in Figure 3. The nodes represent the 71 clusters each depicted as a pie chart displaying the comorbidities and symptoms that are significantly enriched (adj. p-value≤0.05), see Supplementary file 4 for details on the enrichment and p-values.

Hierarchical clustering based on enriched comorbid ICD-10 diagnoses.

The comorbidities present in a minimum of 10 patients and significantly enriched (adj. p-value<=0.05) in each cluster are shown in the pie charts. The number of significant codes ranges from 1 to 10. Each color corresponds to an ICD-10 code chapter as listed in the legend of Figure 1. Six main groups and an outlier (cluster 70) resulted, and the colors of the dendrogram branches indicate to which hierarchical groups the clusters belong. The size of the pie charts represents the average diabetes duration (years with diabetes) divided into six bins. The 21 clusters where at least 50% of the patients have three or more HbA1c severity parameters are marked with a red line surrounding the pie chart.

The 71 clusters were defined based on the associated comorbidities, excluding DM without complications, and from the pie charts we observed that distinct diagnoses do indeed characterize the clusters. For example, ICD-10 code N40: Benign prostatic hyperplasia for cluster 56, L40: Psoriasis for cluster 16, F20: Schizophrenia for cluster 47, K29: Functional intestinal disorders for cluster 17, and Z94: Transplanted organ and tissue status for cluster 42. Using Fisher’s exact test, we found that: Symptoms related to skin and subcutaneous tissue (adj. p-value<0.001) characterized symptom group five and Symptoms related to digestive system and abdomen; cognition, perception, emotional state and behavior; and general symptoms and signs (adj. p-value<0.001 for all) characterized symptom group 3. These results correspond well to the enriched codes observed in Figure 3, as was the case for the other enriched codes across the 71 clusters within the six symptom groups.

Genomic characterization by SNP association of phenotypically determined clusters

We evaluated the 71 clusters in the six symptom groups, plus the outlier cluster, for SNPs that could characterize the different groups (details on the genetic data can be found in the ‘Genomic characterization’ section under Materials and methods). The five highest association signals (independent) for each group are shown in Supplementary file 5. Only results from analyses with more than 15 cases and a well-calibrated QQ-plot (visual inspection and a lambda inflation factor >0.96) are reported. Accordingly, clusters 1–5, 7–9, 12, 15–18, 21–23, 26, 31, 35, 39, 45, 46, and 66, as well as all aggregated symptom groups, met the criteria. The median coverage of the symptom clusters was 31% [range: 10–67]. SNPs characterizing the symptom groups were found in several instances and association signals to disease-associated genes were also found for several of the clusters (Figure 3). Most frequently found, unsurprisingly, were genes associated by GWAS to diabetes or diabetes-related cardio-metabolic traits (cluster 3: MYO3B, cluster 4: DAPK1, cluster 5: LPIN2, cluster 7: SAMD4A and FHIT, cluster 8: ERG and PLCB1, cluster 12: MYT1L, cluster 15: UBE2WP1, cluster 16: ADARB2, CDKAL1, and CLIP1, cluster 17: C8orf37-AS1, cluster 21: FHOD3 and MCF2L, cluster 24: MTCL1, cluster 26: NTM, cluster 31: PCDH15, CDH4, and DCTD, cluster 31: KLF12, cluster 39: FHOD3, cluster 45: IGF1R, BCAS3, and TENM4, cluster 46: NRXN3). Cluster eight is characterized by cardiovascular complications, and three of the top ranking genes for this cluster have been associated with LDL peak particle diameter (THBS4; Rudkowska et al., 2015), abdominal aortic aneurysm (ERG; Jones et al., 2017), pulse pressure (ERG; Warren et al., 2017), and diastolic blood pressure (PLCB1; Warren et al., 2017). Cluster 21 is enriched for the ICD-10 diagnosis foot ulcer (L97), and MCF2L, one of the top ranking genes for cluster, has been associated with both end-stage coagulation (Williams et al., 2013) and prothrombin time (Tang et al., 2012). In total, of the top five association signals that were mapped to genes (n = 103) we found five (CDKAL1, DCDC2C, KLF12, LPIN2, TLE1) to be related with diabetes.

Comorbidity pairs and patterns within symptom related clusters

We detected codes occurring significantly more or less together within and across the symptom groups (Fischer’s test with Bonferroni adjusted p-values<=0.01) defining distinct comorbidity pairs. If the comorbidity pairs covered more than 100 unique codes (symptom groups 4 and 7) we extracted only the most significant pairs until these pairs consisted of 100 unique codes.

Figure 4A illustrates the comorbidity correlations for the six main symptom groups where each pairwise interaction has a comorbidity score (see Material and methods). To characterize whether a diagnosis occurred significantly more before or after another, we made this analysis in a temporal manner. Figure 4B illustrates the comparison of the first diagnosis (row) to the second diagnosis (column). We found that especially the diagnoses related to diabetes (E13, O24), diabetes with complications (shortened to E10 and E11), obesity (E66), diseases of the pancreas (K86), poly- and proteinuria (R35 and R80), and to some extent hypertension and ischemic heart disease (I10, I20, I21, I25) are observed before other diagnoses (blue indicates that the row diagnosis is observed prior to the column diagnosis more than expected, and red indicates the opposite). Focusing on the different symptom groups, we detected which comorbidity pairs were unique in the different groups, and Figure 4C displays these unique comorbidity interactions.

Comorbidity patterns within the six symptom groups.

(A) Comorbidity correlations between the combined symptom groups. (B) Asymmetric comorbidity matrix for observing row diagnosis codes before column diagnoses. First, we calculated Bonferroni corrected p-values for diagnosis pair directionality, second, we extracted the top 100 unique diagnosis codes pairs with lowest adjusted p-values and lastly, we calculated a comorbidity score (CS) by using the log2 of observing the pair more or less than expected. The heat-map colors reflect the CS quantification. (C) Comorbidity pairs unique for each of the symptom groups. All interactions are observed significantly more (blue) or less (red) than expected (adj. p-value<=0.01). Arrows indicate that the diagnoses are observed in the particular order (Fischer’s exact test with Bonferroni correction p-value<=0.01). Node size indicates in how many symptom groups the diagnosis code is observed in, ranging from one group (the diagnosis is unique for the group, largest nodes) to six groups (all groups have the code, smallest nodes).

In symptom group two we found that L84: corns and callosities is observed significantly more together within patients with T1D than T2D (CS = 1.24, adj. p-value=4.06e-15 and CS = −1.58, adj. p-value=1.25e-03, respectively). Temporal analysis of diagnosis occurrence showed that T1D is observed before L84 (Figure 4B, mean time difference = 8.3 years, adj. p-value=1.01e-39). Corns or callosities are unproblematic in healthy people, but in diabetes patients they can cause skin defects that increase the risk for additional complications, for example foot ulcers which can lead to amputations (Apelqvist et al., 2000; Hunt, 2011).

Although not observed significantly together within any clusters the temporal analysis showed that the time between T2D and elevated blood glucose levels (R73) is significantly shorter in symptom group two than in groups 4, 5 and 6 (mean time = 2.2 days; adj. p-value=6.45e-04, 3.29e-06 and 2.73e-06, respectively).

In symptom group 5, five of the eleven clusters are enriched with ICD-10 codes from chapter XIII: Diseases of the musculoskeletal system and connective tissue, especially dorsopathies, spondylopathies and soft tissue disorders. Further, these diagnoses are observed exclusively in this group and show unique disease co-occurrence patterns, for example M48-M54 (other spondylopathies and dorsalgia, CS = 1.01, adj. p-value=1.8e-04) and M43-M47 (deforming dorsopathies and spondylosis, CS = 1.54, adj. p-value=1.91e-06). One of the top ranked genetic associations for this cluster (rs76548985, p-value=1.43e-06) is LINC00351, associated with sporadic amyotrophic lateral sclerosis (Xie et al., 2014). It is worth noting that clusters 8, 22, 33, 35, 45 within symptom group five are all enriched for drugs from ATC chapter A10B: blood glucose lowering drugs, excluding insulin (adj. p-value<0.05), and all but cluster eight are associated with glycemic dysregulation.

Within symptom group 7, we observed two diagnosis pairs less than expected: E11-E13 (CS = −1.46, adj. p-value=1.57e-04), and K86-E78 (CS = −1.24, adj. p-value=5.46e-04). Hence, this group contains patients where T2D and other diabetes as well as diseases of the pancreas and disorders of lipoprotein metabolism are not given together. In contrast, I34: nonrheumatic mitral valve disorder is observed more often than expected together with heart failure (I50, CS = 1.83, adj. p-value=0.009) and atrioventricular and left bundle-branch block (I44, CS = 1.53, adj. p-value=0.0018). Interestingly, one of the top genetic signals for symptom group seven maps to MIR8052 (rs6590490, p-value=3.14e-07) that has been associated with pulse pressure (Warren et al., 2017). Comparably, among the top genetic signals for symptom group 4, a group where are large proportion of the patients are characterized by hypotension (I95) and vertigo (R42), are ANLN that has been associated with systolic blood pressure (Parmar et al., 2016).

Glycemic dysregulation

We evaluated five different parameters associated with glycemic dysregulation (glycemic dysregulation, hyperglycemia, check-point detection of fluctuating HbA1c levels, HbA1c level at diabetes onset and amount of HbA1c observations above diagnosis threshold for T1D and T2D [53 mmol/mol]) and found that 2942 patients did not meet any threshold criterion, 2484 met one, 4647 two, 4057 three, 531 four, and 22 met all five criteria. The distribution of HbA1c measurements for T1D and T2D is shown in Figure 1—figure supplement 5. First, we investigated whether there was any difference in mean values of the 20 different biochemical tests (see Material and methods) and subsequently we applied a Kolmogorov-Smirnov test to assess how these distributions differed. We found that the means of 14 of the different biochemical tests were differently distributed between the six groups (adj. p-value<0.01) of patients with different number of dysregulation parameters, and furthermore observed a distinct difference between the not-or-slightly dysregulated patients (groups 0 to 2) and the middle-or-highly dysregulated patients (groups 3 to 5) (Figure 1—figure supplement 6). The group with five parameters showed no significant difference, due to the low number of patients (N = 22). The group with 3–5 parameters showed higher levels of triglyceride and HbA1c, and lower levels of sodium, urine creatinine, C-peptide, hemoglobin, diastolic blood pressure and height. Elevated levels of HbA1c, triglyceride, LDL-cholesterol and cholesterol and lower levels of HDL-cholesterol are known biochemical values associated with glycemic dysregulation and thus verified our findings.

The detection of higher levels of potassium and plasma creatinine as well as the lowered sodium, urine creatinine, hemoglobin levels indicates that these biochemical tests might be used in future prediction of glycemic dysregulation. Glycemic dysregulation is expected to cause renal problems, (identified by elevated plasma creatinine and elevated urine albumin) and hypertension, which is treated with RAS blocking agents (ACE inhibitors and angiotensin two receptor blockers) and diuretic agents, which elevate potassium and lower sodium. The treatment profile of this patient group revealed an enrichment of patients treated by RAS blocking agents in most of the clusters. Based on these observations we considered having at least three of the parameters as the best approximation for a definition of potential glycemic dysregulation.

Using the results from the biochemical analysis, we divided the cohort in two: those with at least three parameters associated with glycemic dysregulation, and those with two or less. In the 71 clusters defined above, 21 had more than 50% patients with at least three parameters (Figure 3, red circles). We found 10 of the 21 clusters in symptom group 3 of which, cluster 5, 24, and 47, were enriched for poor compliance when using the SDC-custom dictionary (adj. p-value=5.9e-03, 1.9e-03 and 2.6e-02, respectively). By further investigating the enrichment of SDC-custom terms (adj. p-values≤0.05) we found that the majority of the 21 clusters had terms related to cardiovascular complications (e.g. beta blocks, ischemia, diuretics and bypass), kidney complications (e.g. nephropathy, edema and albuminuria), metabolic complications (hypoglycemia and insulin chock) and neurologic related disorders (e.g. neuropathy and loss of memory). Furthermore, all the patients in cluster 47 have schizophrenia (N = 76, adj. p-value=8e-141), and behavioral features might therefore account for the glycemic dysregulation. The same could be the case for cluster 24, in which all have epilepsy (N = 108, adj. p-value=7.6e-186).

Genetic characterization of dysregulated patients

To assess if glycemic dysregulation is a diabetic complication or evidence of disease etiology, we further tested whether any SNPs were associated with glycemic dysregulation (n = 2,120). The five top associating signals map to NCKAP5, CLNK, PSD3, KPNA5, and LINC00333 (Supplementary file 5), although not reaching genome-wide significance. Interestingly, two of the genes associated with schizophrenia (LINC00333 [Goes et al., 2015] and NCKAP5 [Draaken et al., 2015]) and PSD3 have also been associated to traits related to urinary and blood metabolite levels, metabolic traits, and triglyceride levels (Raffler et al., 2015; Teslovich et al., 2010; Shin et al., 2014; Rueedi et al., 2014). However, none of the five top ranked genes have been previously linked to glycemic levels or diabetic dysregulation.

Discussion

Previous studies using EHRs in diabetes research have focused on improving clinical decision making (O'Connor et al., 2011), clinical prediction (Miotto et al., 2016), patient management (Cebul et al., 2011), mortality risk (Pantalone et al., 2009; Pantalone et al., 2010), genetic risk factors (Kho et al., 2012), and subgroup identification (Li et al., 2015). Only the study by Miotto et al. (2016) used the different layers of the EHRs, aimed at predictive measures of clinical outcome. A study from the eMERGE consortium extracted phenotypes from EHR narratives by using NLP-based methods (Kho et al., 2011). They used EHR for phenotypic characterization of five main diseases, but a fine-grained analysis of phenotypic characterization within the diseases was not performed. Further, NLP was included only in the phenotypic determination of three of the diseases, not for diabetes determination.

Stratification and subdivision of diabetic cohorts have typically been performed on homogeneous data sets within specific diabetes types such as T1D, T2D, or gestational diabetes (Perry et al., 2012; Ren et al., 2016; Lin et al., 2012; Achenbach et al., 2004). One of the more recent stratification studies of diabetes patients is Li et al. (2015) that identified subtypes of T2Ds of mixed ethnicity using the structured part of EHRs. They detected three distinct subgroups that could be linked to significant SNPs through gene-disease associations in a patient-unspecific manner. Further, elevated HbA1c levels were used to explain one subgroup with microvascular diabetic complications. In contrast to the study by Li et al., we have taken the stratification and characterization several steps further both by investigating a heterogeneous diabetic cohort almost five times as large and obtaining the full comorbidity pattern and symptoms relatedness through mining of the text-narratives using both an ‘exposure-oriented’ and a diagnosis-based dictionary. In addition, we used the biochemical data to produce a severity classification (the five parameters of glycemic dysregulation) and integrating this with both the text-mined and assigned diagnoses, we were able to determine many different, more homogeneous groups of patients with shared symptoms and comorbidities, as well as different levels for glycemic dysregulation.

Another recent diabetes stratification study by Ahlqvist et al. (2018) used a data-driven approach and k-means clustering to subgroup adult-onset diabetes and characterize five subgroups showing differing disease progression and risk of diabetes complications. However, this approach concerned only individuals with type 2 diabetes and a characterization based on six parameters (glutamate decarboxylase antibodies, age at diagnosis, BMI, HbA1c, and homoeostatic model assessment 2 estimates of β-cell function and insulin resistance), and thus clinical narratives, medication, and genetics were not used as we have done in this study.

The text mining approach used in relation to ICD-10 codes was based on level three rather than the more detailed level four since it would increase tremendously the dimensionality of the feature space. While this obviously reflects a less deep phenotyping, for a data set of this size many level four codes would be unique, likely leading to a less stable subsequent clustering and analysis. In fact, our attempt to use the much more fine-grained SNOMED-CT terminology confirmed that a data set needs to be very large for such a fine-grained vocabulary to be useful.

In this work, we deliberately excluded the primary diabetes types without complications, T1D and T2D, and thereby constructed a stratification of the cohort driven solely by comorbidities, complications, other diseases, and symptoms. However, combining different diabetic subtypes can be problematic, since their etiologies differ and disease progression is different across diabetes types, treatment, compliance and lifestyle (Adeghate et al., 2006). Our focus was not to characterize specific comorbidity-related groups within a certain diabetes type, since extensive epidemiological studies of this kind have been done previously. Instead, we focused on the diabetes continuum with the aim of investigating whether it was possible in an unsupervised manner to detect relevant and meaningful diabetic subgroups by comorbidities, symptoms, or level of glycemic dysregulation. Further, we detected novel biochemical and genetic candidates that might relate these to the different cohort subdivisions, such as shared symptom patterns for phenotypically similar patients and the level of glycemic dysregulation. These biochemical and genetic candidates could be potential risk factors for additional complications, especially concerning glycemic dysregulation, that could be verified by further experimental studies. As the cohort is enriched for sicker patients with diabetes melitus complications the features and the overall grouping described would not necessarily be the same if another cohort dominated by prediabetes individuals would have been analyzed.

Despite our focus on the phenotypic variation among diabetes patients, the stratification is restricted by the limited coverage of the genetic data, which lowers the power considerably. We were able to obtain genetic data for 2337 patients, of whom 2125 remained after quality control and stratification. Hence, only 14% of the patients in our final cohort had descriptive genetic information.

By adding biochemical, prescription, and genetic data we observed that the clusters were significantly different from each other on parameters other than comorbidities. By including the text narratives of the EHRs we were able to capture diagnoses that in another context would be considered as a primary diagnosis, for example epilepsy, schizophrenia and cerebral palsy. These diagnoses are not known comorbidities of diabetes but can influence the treatment and management of the diabetes patient. For instance, we observed that all patients in cluster 47 had schizophrenia, which could influence their compliance since the cluster was associated with glycemic dysregulation. We determined this when assessing the level of glycemic dysregulation and found that this cluster indeed showed a high number of patients with at least three parameters for glycemic dysregulation. However, a more in-depth analysis is required to clarify whether the glycemic dysregulation is due to the behavioral effects of schizophrenia, underlying genetic variants, adverse drug reactions due to polypharmacy, or other variables.

Despite our data from both assigned and text-mined diagnoses, misdiagnoses can occur, and we performed a manual inspection of randomly selected EHRs to establish the validity of the data. Furthermore, we observed some patients assigned with different diabetes types, for example first assigned with T1D and later with T2D, and vice versa. Inspecting the biochemical values of GAD65 autoantibodies and comparing them to the primary diagnosis type we found 182 T2D assigned individuals to have GAD65 levels above 10 IU/ml, possibly indicative of T1D or LADA; however, these individuals were not significantly enriched in any cluster. We also observed 621 individuals with GAD65 levels below 10 IU/ml, which is consistent with known late-term effects of T1D (results not shown). An in-depth temporal analysis of these patients with mixed diabetes types could be interesting and integrating biochemical as well as genetic variation data could elucidate which, if any, phenotype might be the most accurate.

In this study, we have used data from a unique cohort of 14,017 patients with diabetes, of which 12,866 had been diagnosed with either T1D or T2D. Integrating the assigned and text-mined ICD-10 and SDC-custom diagnoses, an MCL clustering was carried out which resulted in 172 unique clusters. Of these, 71 had at least 50 patients, which were subsequently divided into groups with shared symptoms. Investigating the complication enrichment and comorbidity patterns in the clusters and symptom groups we detected clusters described by specific disorders such as hypothyroidism, schizophrenia, and functional intestine disorder as well as unique comorbidity interaction patterns both with and without temporal significance. An interesting approach could be to extend the temporal analysis to investigate how disease progression within and between clusters and symptoms groups develops for multiple diagnoses. This could be done with a trajectory-based approach as done recently by Jensen et al. (2014).

Materials and methods

EHR data

Request a detailed protocol

All data originate from the Steno Diabetes Center Copenhagen (SDCC), a specialized diabetes hospital in the Capital Region of Denmark. In Denmark patients with type 1 diabetes (T1D) are followed in hospital outpatient clinics such as SDCC, and the T1D patients studied comprise 35% of all adult patients with T1D in the Capital Region of Denmark. Patients with type 2 diabetes (T2D) are referred from primary care for treatment optimization, typically for a period of six to twelve months. When treatment goals are reached, and they have no diabetic complications, they are referred back to general practice. Patients needing intensive control and treatment, because of micro- or macrovascular complications, are offered life-long follow-up at the SDCC. At any time, approximately 2000 patients with complicated T2D are followed at the SDCC. Generally, the patients registered in the SDCC electronic patient records are representative of Danish patients with T1D and the 10% most complicated patients with T2D (Jørgensen et al., 2016). Moreover, the patient followed at SDCC are comparable to patients followed in all Danish hospital diabetes outpatient clinics in terms of distribution of age and duration of diabetes. The data comprise all communications and contacts recorded at the hospital over a period of 19 years (1993–2012) for 14,017 patients. This includes, primary diagnoses, prescriptions and laboratory tests, 1.2M clinical narrative entries, 420 different types of laboratory tests with 4.15M laboratory measurements and a total number of 440,555 drug prescriptions. On average, each patient had 85 clinical narratives with an average length of 34 words (212 characters). In addition, genetic data from several research projects have been linked to the patients and added to the EHRs.

Text-mining dictionaries, tagging and corpus matches

Request a detailed protocol

An in-house developed framework for mining Danish text was used for the analysis (Roque et al., 2011; Eriksson et al., 2013). The algorithm tags words in the text narratives in a named entity recognition (NER) fashion based on supplied dictionaries. In this study, we used two main dictionaries: The International Classification of Disease version 10 (ICD-10) truncated to level 3 (e.g. E10: Type 1 Diabetes), and a complementary ‘exposome-oriented’ dictionary (SDC-custom). The latter holds terms related to diabetes specific subtypes (e.g. MODY and LADA), complications (e.g. the different severities of neuropathy, retinopathy and nephropathy), treatments and examinations (e.g. gastric bypass, renography, and beta blockers), lifestyle and lifestyle related disorders (e.g. obesity, exercise level, smoking), and compliance. The SDC-custom dictionary was developed in collaboration with physicians at the SDCC (see Supplementary file 6 for a translated and condensed version). The Danish ICD-10 version currently contains roughly 20,000 unique descriptions of clinical concepts, each with a unique ICD-10 code.

The NER used for dictionary matching, in addition, performs lemmatization and de-latinization of tagged words, accounts for language negations or subject negations (e.g. ‘the patient’s mother had retinopathy’), and performs fuzzy matching with a Hamming distance of 1 (e.g. ‘diabtes’ is transformed to its correct spelling ‘diabetes’). A thorough explanation of the algorithm is provided (Simon et al., 2019, manuscript in preparation). Other details, for example on ‘negation scope’, that is the position of negations relative to the negated term in Danish, have been published previously (Thomas et al., 2014).

Running the text-mining algorithm (Simon et al., 2019, manuscript in preparation) on the SDCC corpus with the two dictionaries (ICD-10 and SDC-custom) recognized 1,028,593 entities from the dictionaries in 12,504 patients (80.5% of the entire corpus). None of the remaining patients had any non-trivial match between the dictionaries and EHR narratives. The two dictionaries shared some general terms, for example T1D and T2D; these duplicate matches were removed and 941,087 unique code matches remained. Of these, 267,404 were fuzzy matches representing 4181 unique variants. The variants were manually validated, resulting in removal of 10,952 (4.1%) matches. After removal of negated sentences (n = 255,302) 594,600 code-to-text matches in 12,467 patients were left.

Patient phenotype vectors from assigned and text mined codes

Request a detailed protocol

The structured ICD-10 codes assigned to patients during their contact with SDCC were extracted from the EHRs, along with all ICD-10 codes captured by mining the text parts of the EHRs. The two ICD-10 lists were combined, but to prevent the primary, assigned diabetes types from dominating the patient stratification, diagnosis codes for diabetes without complications (E10 and E109, in total 3740 codes, and E11 and E119, in total 3624 codes) were removed. Approximately 8% of the assigned codes were removed in this way. The list of codes and their frequencies for each patient were transformed using the BM25 weighting scheme (Robertson and Walker, 1994), which scores a code c in patient P, accounting for the code frequency in all patients, frequency of the codes in the patient (document frequency), and number codes in the patient record (document length), see Equation 1.

(1) Score(p,c)=i1nIDF(ci)f(ci,p)(k1+1)f(ci,p)+k1(ab+b|p||pave|)

Here, IDF(c) is the inverse document frequency for the code c computed as

IDF(ci)=logNn(ci)+0.5n(ci)+0.5

With N being the total number of patients and n(c) the number of patients with a given code ci , and the term f(ci,p) is the frequency of code ci in patient p. The number of codes associated with each patient vector, P, is given by the length of the vector, |p|, and the average number of codes in the entire corpus is |pave|. Finally, b and k1 are free parameters that determine to what extent document length is considered (b) and how much the scoring equation resembles a normal TF-IDF (k1), respectively. The value of b was set to 0.75 and does not fully account for the document length (b=1) and k1 was set to 1.2 giving a low resemblance of TF-IDF (k1).

Clustering patients from Cosine similarities

Request a detailed protocol

All patients were clustered using their pairwise cosine similarities calculated from the BM25 transformed code vectors. A cosine similarity ≥ 0.5 was set as a cut-off prior to clustering, to minimize the number of edges in the subsequent patient network. To increase the variance of the cosine similarities, these were scaled from the interval 0.5–1 to 10–100. We wanted to do a network based clustering and therefore used Markov Clustering (MCL) with the inflation parameter set to 1.2 and rest left as default (Van Dongen, 2000). Different inflation parameters were tested and evaluated based on the efficiency, mass fraction, and area fraction parameters.

Grouping clusters in symptom related groups

Request a detailed protocol

We organized the clusters into symptom groups based on the frequency of their symptom codes using ICD-10 chapter XVIII level 1, for example R50-69: General symptoms and signs. We used a Euclidean distance and applied a hierarchical clustering using Ward.D as the agglomeration method since we wanted to expose the hierarchical relationship amongst the clusters. The entire analysis was performed using R (version 3.2.1).

Enrichment analysis of diagnosis codes

Request a detailed protocol

The MCL clusters were tested for ICD-10 and SDC-custom codes found more often than expected, using a binominal test while correcting for sex and birth decade. The metadata such as average age, days at SDCC, and diabetes duration (from the date of diabetes diagnosis until the end of the study) were calculated, and further p-values for each cluster were obtained using a Wilcoxon test against the remaining clusters. In both analyses, p-values were adjusted using Benjamini-Hochberg correction for multiple testing, where a p-value≤0.05 was considered significant.

Comorbidity patterns for diagnosis pairs

Request a detailed protocol

We performed three independent analyzes without considering the clusters by applying Fischer’s exact tests to obtain p-values for all diagnosis pairs within the SDCC corpus: 1) p-values for observing the codes together, 2) p-values for observing diagnosis A prior to diagnosis B, and 3) p-values for observing diagnosis B prior to diagnosis A. P-values from the three different sets were adjusted using Bonferroni correction for multiple testing, and the pairs were subsequently ranked based on these values. To detect whether the pairs were observed more together than expected we applied a comorbidity score as described in Roque et al. (2011). For the temporal pairs, we also applied an ANOVA test to investigate whether any of these pairs were unique for a symptom group. All p-values were corrected for multiple testing, and an adjusted p-value≤0.05 was considered significant.

Robustness of the MCL generated clusters

Request a detailed protocol

To assess quantitatively the stability of the clusters generated, we constructed various diluted and shuffled realizations of the similarity network used as input to the MCL algorithm. We used a reference clustering similar to the clustering presented in Figure 2B (either by including the patients in the 71 clusters or all patients). The diluted versions were generated by randomly deleting edges with a probability of α, whereas the shuffled realizations were created by shuffling edges between nodes (patients) as described earlier (Karrer et al., 2008). The latter produces a network where the number of edges and vertices are unchanged. An α of zero leaves the reference network unchanged, while a value of 1 leads to a complete randomization of the similarity network. Each of these randomizations of the input were repeated five times for various values of α in the range 0–50% and used as input for the MCL algorithm. The resulting clustering’s were then compared to the reference clustering by means of the Variation of Information measure (VI) (Meilă, 2007) and plotted as function of increasing values of α (see Figure 2—figure supplement 5). The figure includes two horizontal lines corresponding to the value that the VI would take if we were to randomly assign 10% and 20% of the vertices to different random clusters, respectively. This analysis showed that the clustering is stable in relation to removing edges, which is evidence that the cosine metric-based cutoff used does not change the overall structure of the clustering. The shuffling is a more impactful randomization, however despite this, we can still shuffle around 10% of the edges and still retrieve 90% of the patients in the groups of the 71 reference clusters.

Quantitative assessment of glycemic dysregulation

Request a detailed protocol

Glycemic dysregulation was assessed for each patient by evaluating five different parameters. The first two parameters were obtained using the SDC-custom code for dysregulation (sdcL03) and the ICD-10 codes for hyperglycemia (R73 and E89). The remaining three were found by analyzing longitudinal measurements for glycated hemoglobin (HbA1c). Due to a large variation in both the number of measurements and their frequency, HbA1c values were pre-processed. We divided the HbA1c measurements for each patient into segments containing a minimum of five values, spanning a time interval of at least three months (equivalent to the functional lifetime of red blood cells). In total 10,112 patients had HbA1c measurements that fulfilled the criteria, and the subsequent analyses were performed on this sub-population.

We performed three analyses on the longitudinal pre-processed HbA1c data for each patient: 1) a Bayesian analysis of change point detection to find potential peaks of HbA1c values in a patient, 2) analysis of mixed effects models to estimate the HbA1c value at diabetes onset, and, 3) analysis of the frequency of values in different HbA1c bins (e.g. general level for diagnosing T1D or T2D, the critical interval for hyperglycemia etc.) to appoint an HbA1c severity score.

Laboratory test data

Request a detailed protocol

The laboratory tests were longitudinal data such as blood pressure measurements and biochemical analyses of blood and urine samples, and each test was assigned a unique identifier using the NPU-terminology, which is the recommended administration and communication measure of laboratory tests in Denmark (Petersen et al., 2012). In our data, several laboratory tests had an SDC identifier, being from local laboratory facilities at SDCC. Both test IDs, NPU and SDC, were analyzed separately, despite sometimes measuring the same biochemical variables.

In total, 420 different physiological tests were performed across 14,847 patients from the entire corpus. Measurements within and between tests were unbalanced with no general system in measurement interval, frequency, or number of patients who had a test taken. Due to this lack of systematic coverage, only tests that were performed on at least 75% of the entire corpus (10,788 patients) were analyzed (26 tests). However, the test for C-peptide (NPU18004) was also included as it was available for 74.9% of the cohort and is widely used to distinguish T1D and T2D. Measurements outside the biological reference interval for a given test, that is HbA1c measurements below 15 mmol/mol and above 184 mmol/mol, were removed, and for each patient the mean, median and standard deviation for each test with continuous values (20 of the 26 tests) were calculated. If the data was not normally distributed for a test we log-transformed it and normalized all values to mean = 0 and SD = 1. All analyses after assigning patients to clusters were performed on the 10,788 patients.

We applied a MANOVA to test if means among the three different patient groups (clusters, symptom groups or patients being dysregulated) were significantly different, and a Kolmogorov–Smirnov test was applied to investigate whether the distribution of the sample means in the patient groups were significantly higher or lower than means in the remaining groups. All p-values were adjusted using Bonferroni correction for multiple testing, and an adjusted p-value≤0.05 was considered significant.

Drug prescription data

Request a detailed protocol

Prescription data was available for 12,147 patients with a total number of 440,555 drug prescriptions. Drug compounds were identified by the ATC classification system, which is divided into groups at five different levels. In this study, we summarized the data using ATC-codes at level three and four: chemical and pharmacological and therapeutical, respectively.

From the initial set of prescriptions, we manually reviewed 104 drugs which did not have an ATC code in the EHR or were mapped to more than one ATC code. In addition to the manual review, pro.medicin (www.pro.medicin.dk, accessed October 2018) was used to map drug names to their corresponding ATC code. The SDCC prescription data and the WHO Collaborating Centre for Drug Statistics Methodology (www.whocc.no, accessed October 2018) were used for crosschecking. We performed Fisher’s exact test to investigate prescription enrichment (3rd level of the ATC classification) in clusters with at least 50 patients. The p-values were adjusted using Benjamini–Hochberg correction for multiple testing, and an adjusted p-value≤0.05 was considered significant.

Genomic characterization

Request a detailed protocol

A total of 2290 patients with T2D and 1028 patients with T1D from SDCC were genotyped separately using the HumanOmniExpress (24v1) array from Illumina as previously described (Charmet et al., 2018; Steinthorsdottir et al., 2014). Genotypes were called using GenomeStudio, and imputed separately using the Haplotype Reference Consortium (HRC) imputation panel (McCarthy et al., 2016). Prior to imputation, the two datasets were filtered to retain only high-quality samples/SNPs (sample call rate ≥98%, no mislabeled sex, no ethnic outliers, heterozygosity within 2 SD from the mean, SNP call rate ≥98%, no monomorphic SNPs, no Hardy–Weinberg disequilibrium outliers). After imputation, SNPs with minor allele frequency (MAF) <0.01, more than 20% missingness, R square less than 0.30, and duplicate SNPs were removed, and the two datasets were merged retaining only variants common to the two sets. After merging, relatedness between all individuals were calculated and close relatives were excluded. Of the 3318 patients, 2337 had EHR information and could be mapped to clusters. In total 2125 patients passed quality control and were taken forward for genomic characterization. Logistic regression was used to test for genetic differences (PLINK 1.90 beta, https://www.cog-genomics.org/1.9) between the different groups of interest (clusters and symptom groups) and linear regression was used to evaluate the SNPs impact on dysregulation. Cases were defined as all individuals in a given cluster/symptom group, and controls as all individuals not belonging to the respective cluster/symptom group. Glycemic dysregulation was defined as a score ranking from 0 (low) to 5 (high) based on five dysregulation parameters (see section on glycemic dysregulation). All analyses were adjusted for age and sex. The test statistics were adjusted for inflation (population stratification) using the three first principal components estimated using the --pca function in PLINK. Genetic associations were defined based on data derived from the EBI GWAS catalog version 1.0.1 (http://www.ebi.ac.uk/gwas/) unless otherwise stated. A p-value less than 5*10–8 was considered genome-wide significant.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
    Disease comorbidity network guides the detection of molecular evidence for the link between colorectal Cancer and obesity
    1. Y Chen
    2. L Li
    3. R Xu
    (2015)
    AMIA Joint Summits on Translational Science Proceedings. AMIA Joint Summits on Translational Science 2015:201–206.
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
    Diabetes: foot ulcers and amputations
    1. DL Hunt
    (2011)
    BMJ Clinical Evidence 2011:0602.
  23. 23
  24. 24
  25. 25
    Meta-Analysis of Genome-Wide association studies for abdominal aortic aneurysm identifies four new Disease-Specific risk loci
    1. GT Jones
    2. G Tromp
    3. H Kuivaniemi
    4. S Gretarsdottir
    5. AF Baas
    6. B Giusti
    7. E Strauss
    8. FN Van't Hof
    9. TR Webb
    10. R Erdman
    11. MD Ritchie
    12. JR Elmore
    13. A Verma
    14. S Pendergrass
    15. IJ Kullo
    16. Z Ye
    17. PL Peissig
    18. O Gottesman
    19. SS Verma
    20. J Malinowski
    21. LJ Rasmussen-Torvik
    22. KM Borthwick
    23. DT Smelser
    24. DR Crosslin
    25. M de Andrade
    26. EJ Ryer
    27. CA McCarty
    28. EP Böttinger
    29. JA Pacheco
    30. DC Crawford
    31. DS Carrell
    32. GS Gerhard
    33. DP Franklin
    34. DJ Carey
    35. VL Phillips
    36. MJ Williams
    37. W Wei
    38. R Blair
    39. AA Hill
    40. TM Vasudevan
    41. DR Lewis
    42. IA Thomson
    43. J Krysa
    44. GB Hill
    45. J Roake
    46. TR Merriman
    47. G Oszkinis
    48. S Galora
    49. C Saracini
    50. R Abbate
    51. R Pulli
    52. C Pratesi
    53. A Saratzis
    54. AR Verissimo
    55. S Bumpstead
    56. SA Badger
    57. RE Clough
    58. G Cockerill
    59. H Hafez
    60. DJ Scott
    61. TS Futers
    62. SP Romaine
    63. K Bridge
    64. KJ Griffin
    65. MA Bailey
    66. A Smith
    67. MM Thompson
    68. FM van Bockxmeer
    69. SE Matthiasson
    70. G Thorleifsson
    71. U Thorsteinsdottir
    72. JD Blankensteijn
    73. JA Teijink
    74. C Wijmenga
    75. J de Graaf
    76. LA Kiemeney
    77. JS Lindholt
    78. A Hughes
    79. DT Bradley
    80. K Stirrups
    81. J Golledge
    82. PE Norman
    83. JT Powell
    84. SE Humphries
    85. SE Hamby
    86. AH Goodall
    87. CP Nelson
    88. N Sakalihasan
    89. A Courtois
    90. RE Ferrell
    91. P Eriksson
    92. L Folkersen
    93. A Franco-Cereceda
    94. JD Eicher
    95. AD Johnson
    96. C Betsholtz
    97. A Ruusalepp
    98. O Franzén
    99. EE Schadt
    100. JL Björkegren
    101. L Lipovich
    102. AM Drolet
    103. EL Verhoeven
    104. CJ Zeebregts
    105. RH Geelkerken
    106. MR van Sambeek
    107. SM van Sterkenburg
    108. JP de Vries
    109. K Stefansson
    110. JR Thompson
    111. PI de Bakker
    112. P Deloukas
    113. RD Sayers
    114. SC Harrison
    115. AM van Rij
    116. NJ Samani
    117. MJ Bown
    (2017)
    Circulation Research 120:341–353.
    https://doi.org/10.1161/CIRCRESAHA.116.308765
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
    Some simple effective approximations to the 2–Poisson Model for Probabilistic Weighted Retrieval
    1. SE Robertson
    2. S Walker
    (1994)
    SIGIR ’94 Proceedings of the 17th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval. pp. 232–241.
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
    Biological, clinical and population relevance of 95 loci for blood lipids
    1. TM Teslovich
    2. K Musunuru
    3. AV Smith
    4. AC Edmondson
    5. IM Stylianou
    6. M Koseki
    7. JP Pirruccello
    8. S Ripatti
    9. DI Chasman
    10. CJ Willer
    11. CT Johansen
    12. SW Fouchier
    13. A Isaacs
    14. GM Peloso
    15. M Barbalic
    16. SL Ricketts
    17. JC Bis
    18. YS Aulchenko
    19. G Thorleifsson
    20. MF Feitosa
    21. J Chambers
    22. M Orho-Melander
    23. O Melander
    24. T Johnson
    25. X Li
    26. X Guo
    27. M Li
    28. Y Shin Cho
    29. M Jin Go
    30. Y Jin Kim
    31. JY Lee
    32. T Park
    33. K Kim
    34. X Sim
    35. R Twee-Hee Ong
    36. DC Croteau-Chonka
    37. LA Lange
    38. JD Smith
    39. K Song
    40. J Hua Zhao
    41. X Yuan
    42. J Luan
    43. C Lamina
    44. A Ziegler
    45. W Zhang
    46. RY Zee
    47. AF Wright
    48. JC Witteman
    49. JF Wilson
    50. G Willemsen
    51. HE Wichmann
    52. JB Whitfield
    53. DM Waterworth
    54. NJ Wareham
    55. G Waeber
    56. P Vollenweider
    57. BF Voight
    58. V Vitart
    59. AG Uitterlinden
    60. M Uda
    61. J Tuomilehto
    62. JR Thompson
    63. T Tanaka
    64. I Surakka
    65. HM Stringham
    66. TD Spector
    67. N Soranzo
    68. JH Smit
    69. J Sinisalo
    70. K Silander
    71. EJ Sijbrands
    72. A Scuteri
    73. J Scott
    74. D Schlessinger
    75. S Sanna
    76. V Salomaa
    77. J Saharinen
    78. C Sabatti
    79. A Ruokonen
    80. I Rudan
    81. LM Rose
    82. R Roberts
    83. M Rieder
    84. BM Psaty
    85. PP Pramstaller
    86. I Pichler
    87. M Perola
    88. BW Penninx
    89. NL Pedersen
    90. C Pattaro
    91. AN Parker
    92. G Pare
    93. BA Oostra
    94. CJ O'Donnell
    95. MS Nieminen
    96. DA Nickerson
    97. GW Montgomery
    98. T Meitinger
    99. R McPherson
    100. MI McCarthy
    101. W McArdle
    102. D Masson
    103. NG Martin
    104. F Marroni
    105. M Mangino
    106. PK Magnusson
    107. G Lucas
    108. R Luben
    109. RJ Loos
    110. ML Lokki
    111. G Lettre
    112. C Langenberg
    113. LJ Launer
    114. EG Lakatta
    115. R Laaksonen
    116. KO Kyvik
    117. F Kronenberg
    118. IR König
    119. KT Khaw
    120. J Kaprio
    121. LM Kaplan
    122. A Johansson
    123. MR Jarvelin
    124. AC Janssens
    125. E Ingelsson
    126. W Igl
    127. G Kees Hovingh
    128. JJ Hottenga
    129. A Hofman
    130. AA Hicks
    131. C Hengstenberg
    132. IM Heid
    133. C Hayward
    134. AS Havulinna
    135. ND Hastie
    136. TB Harris
    137. T Haritunians
    138. AS Hall
    139. U Gyllensten
    140. C Guiducci
    141. LC Groop
    142. E Gonzalez
    143. C Gieger
    144. NB Freimer
    145. L Ferrucci
    146. J Erdmann
    147. P Elliott
    148. KG Ejebe
    149. A Döring
    150. AF Dominiczak
    151. S Demissie
    152. P Deloukas
    153. EJ de Geus
    154. U de Faire
    155. G Crawford
    156. FS Collins
    157. YD Chen
    158. MJ Caulfield
    159. H Campbell
    160. NP Burtt
    161. LL Bonnycastle
    162. DI Boomsma
    163. SM Boekholdt
    164. RN Bergman
    165. I Barroso
    166. S Bandinelli
    167. CM Ballantyne
    168. TL Assimes
    169. T Quertermous
    170. D Altshuler
    171. M Seielstad
    172. TY Wong
    173. ES Tai
    174. AB Feranil
    175. CW Kuzawa
    176. LS Adair
    177. HA Taylor
    178. IB Borecki
    179. SB Gabriel
    180. JG Wilson
    181. H Holm
    182. U Thorsteinsdottir
    183. V Gudnason
    184. RM Krauss
    185. KL Mohlke
    186. JM Ordovas
    187. PB Munroe
    188. JS Kooner
    189. AR Tall
    190. RA Hegele
    191. JJ Kastelein
    192. EE Schadt
    193. JI Rotter
    194. E Boerwinkle
    195. DP Strachan
    196. V Mooser
    197. K Stefansson
    198. MP Reilly
    199. NJ Samani
    200. H Schunkert
    201. LA Cupples
    202. MS Sandhu
    203. PM Ridker
    204. DJ Rader
    205. CM van Duijn
    206. L Peltonen
    207. GR Abecasis
    208. M Boehnke
    209. S Kathiresan
    (2010)
    Nature 466:707–713.
    https://doi.org/10.1038/nature09270
  58. 58
    Negation scope and spelling variation for text-mining of Danish electronic patient records
    1. CE Thomas
    2. P Bjødstrup Jensen
    3. T Werge
    4. S Brunak
    (2014)
    Proceedings of the 5th International Workshop on Health Text Mining and Information Analysis. pp. 64–88.
  59. 59
  60. 60
  61. 61
  62. 62
    Graph Clustering by Flow Simulation
    1. S Van Dongen
    (2000)
    University of Utrecht.
  63. 63
  64. 64
  65. 65
  66. 66

Decision letter

  1. Alfonso Valencia
    Reviewing Editor; Barcelona Supercomputing Center - BSC, Spain
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The paper describes the stratification of a diabetes cohort based on characteristics extracted from the medical records of a rather homogenous population. These characteristics include other diagnoses and lifestyle factors. The clusters of patients, obtained with an unsupervised Markov clustering method, represent characteristic longitudinal glycemic dysregulation patterns that include the temporal order of comorbidities, as well as genetic features, i.e. SNPs close to some of the known diabetes related genes.

The work opens the doors to the study of the molecular basis of the classical etiological subtypes of diabetes at the light of their relations with the clusters based on comorbidity relationships and symptoms.

Decision letter after peer review:

Thank you for submitting your article "Linking glycemic dysregulation in diabetes to symptoms, comorbidities and genetics through EHR data mining" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

The biomedical goal of the paper is to identify subgroups of diabetes patients. The study is based on a large set of medical records from patients with T1D or T2D. The methodology produces 6 large clusters (MCL) combining structure information (ICD-10 level 3 codes) and symptoms extracted from the free text. The analysis shows enrichment in clinically relevant groups and relations with other diseases. The paper includes the specific analysis of glycemic deregulation, its relation with comorbidities and with previously identified genes/SNPs.

Essential revisions:

The topic is relevant and timely. The technical contribution is at the leading front of the analysis of medical records analysing structured and free text information.

On the critical side there are limitations in clustering analysis, as well as in the use of limited ICD third level data. These serious criticisms have to be addressed completely in a revised version. There are also a number of other issues related with the origin of the data (temporal coverage of the cases) that need to be addressed in detail in the text.

Finally, it is not completely clear what are new results and what is validated (and how). Even, if these are common problems in many large scale studies, it is important to address them clearly in the paper. This point is related with the criticisms on the unclear goal of the paper: technical or medical? The comments of the reviewers reflect that given the situation, you will have to think if you want to prioritise the novelty of the approach instead of making medical claims that might be difficult to sustain.

We ask you to address the following;

More details on the key points above:

1) Circularity

The analysis strategy and reasoning appears to be circular and there are not provided measures of independence between codes, factors, elements, etc. that are introduced or measured at each step. In particular, there are used two vocabularies (which obviously expand the number of "codes" assigned to patients, as it should be; first section Results is not novel) that gave approximately 940,000 codes as explained in Materials and methods. However, it is expected that many of these codes are in fact correlated, non-independent, but instead they are all used to generate patient clusters; for example, in several instances it is emphasized that codes from primary diabetes diagnoses are not used in clustering/analyses, but it is well known a priori that many other codes are significantly correlated with the diagnoses (it is therefore not unexpected that patient clusters differentiate diabetes subtypes). Therefore, it would also not be unexpected that clusters differentiate glycemic levels. In fact, the 71 clusters are defined using distances from ICD-10 chapter XVIII symptoms (subsection “Enriched comorbidity and symptom patterns in diabetes patient clusters”, first paragraph), which includes examination of blood parameters. Again, latter the clusters are defined based on comorbidities, which presumably many are linked to the original codes. Another example of this point is the identification of a "schizophrenia" cluster (subsection “Glycemic dysregulation”, last paragraph), which is expected if the corresponding codes are included. I agree that an observation of differential glycemic levels in this clusters may be interesting, but then first it should be evaluated the overall correlation among the thousands of codes used in the study, and also then the study is limited to assess co-occurrences of terms/codes.

One possible way out could be the comparison results of structured vs. structured+unstructured features.

2) Clustering

As overall conclusion, Discussion section includes the statement that the results show "many different and homogenous groups of patients", but it has to be formally demonstrated that the groups are independent, that they are many more than randomly expected by permuting the same cohort, and that they are more homogenous that some value.

Additionally, questioning the clusters stability, would be necessary to provide a sensitivity analysis on the clusters to illustrate the stability.

3) ICD codes

Truncating the codes to three digits is akin to adding noise by losing information. Please justify your rationale behind truncating the codes. It'll be good if this issue could be addressed in the Discussion or one of the limitations of the work.

4) Data limitations and temporal series

First, there might be problems with the time-dependent comorbidities in a given cluster. There are multiple reasons: (1) the length of the observational window before and after considered in a given cluster is not talked about/defined. This can impact the enrichment of symptoms and hence all the downstream analysis; (2) length observational windows is bound to be different in different clusters, which may introduce the bias in computing comorbidities; (3) it is not clear from the analysis if there exist any relation between the observational window and cluster size.

Indeed, it is not clear how temporality is defined. The heatmap presented in Figure 2B does not help either. Specifically, what is the index date, relative to which pre and postcodes were identified? It is quite likely that the index date would vary for patients, therefore, how does that variability was accounted? Please explain this paragraph in more detail and if possible, then provide a schematic diagram representing temporality and how it was considered in this analysis. In addition, in the fourth paragraph of the subsection “Comorbidity pairs and patterns within symptom related clusters”, it appears the time window before and after is likely to bias the estimation of T2D and elevated blood glucose level. More clarity on this will be helpful.

In any case, a better explanation is needed about the length of the observational window before and after used to compute comorbidities. The length of the observational window (before and after in a given cluster) is bound to confound the disease association. Was this length kept similar for all the clusters? If not, then it is very much possible that certain cultures are (just based on cluster size) are likely to be enriched in certain symptoms. Please provide an in-depth analysis on this.

5) Additional data

Please briefly explain the type of genetic data is included and at what level of omics.

The candidate genetic and biochemical markers of glucose regulation could be of interest if replicated. However, could these have been discovered without the clustering e.g. using simpler regression methods?

It is unclear how the clusters segregate in their genetic basis, which could then support the identification of disease subtypes. The question is not as much which genetic variants are associated with a given cluster, but if the effect estimations are statistically different between clusters, which would then favor the existence of subtypes. It is also unclear if the association signals are genome-wide significant in any case.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Linking glycemic dysregulation in diabetes to symptoms, comorbidities and genetics through EHR data mining" for further consideration at eLife. Your revised article has been favorably evaluated by Naama Barkai as the Senior Editor, a Reviewing Editor, and two reviewers.

The manuscript has been improved but there is still some confusion about the significance of the resulting clusters, significance of the selected thresholds and the subsequent biological interpretation, as outlined below by the first reviewer. The final version should clarify the exploratory nature of the work and the limits of the statistical reliability, as pointed out by the second reviewer.

Reviewer #1:

The authors have done a perceptible effort trying to answer the major concerns initially raised, and I personally appreciate the clarifications and additional data provided.

However, I am still not convinced that the main results are robust enough to support the overall message, in particular regarding the title: "Linking glycemic dysregulation.…to symptoms, comorbidities and genetics…".

First, the issue of "circularity". I might have over-looked the original Discussion that mentioned that primary and secondary diabetes diagnoses were excluded, or my initial comment might have been unclear, but I was referring to the expectation that clusters should obviously identify major known diabetes types; in particular if dozens of diagnosis codes (some of them correlated) and extensive text vocabulary is analyzed. The point is not if the clusters identify what it is already known, but if they identify new etiology/subtypes. Of note, the following example is somewhat inaccurate: "… triple negative breast cancer patients share many features with ER-positive patients for example." Indeed TNBC and ER+ cases do not share cell-of-origin, prognosis, tissue-metastasis preference, neither therapeutic approaches.

My reasoning above is then linked to the title issue: the study of glycemic perturbation is presented in one section (“Glycemic dysregulation”) and remains unclear which clusters, symptoms, are robust regarding statistical differences for "glycemic perturbation" from the others. The thresholds of 0-2, 3-5 parameters, and then the clusters (21 with 50% patients 3 parameters…) appear to be arbitrary. The use of sentences with the following grammatical constructions may not help to understand the analyses: "distinct differentiation…dysregulation parameters". Next, glycemic perturbation is linked to genetic risk in the last Results section, but as initially raised, it is unclear what means "top association signals"; are these significant genome-wide, are they statistically different between the clusters/symptoms? Are cluster-risk interactions significant? Honestly, this remains unclear to me and I hope I am not overlooking to any data. I agree with the authors on their final response to the question of "circularity": 1) To which extent does it reproduce the classification made by doctors? 2) Does it allow us to make a more specific subgrouping of patients?

I believe that the lack of specific subgroups with statistically defined differences across the three features included in the title is evident from the Abstract, which does not detail any result: it just includes a general message of what the study has found as a last sentence.

Reviewer #2:

I am satisfied with the revision of the article. I feel that the questions I raised have been addressed. Limitation in the use of 3-digit ICD9 codes remains, but the response from authors and subsequent changes made in the manuscript are appropriate.

This study may serve as a starting point for more in-depth analysis in subsequent analysis.

https://doi.org/10.7554/eLife.44941.sa1

Author response

Essential revisions:

[…]

1) Circularity

The analysis strategy and reasoning appears to be circular and there are not provided measures of independence between codes, factors, elements, etc. that are introduced or measured at each step. In particular, there are used two vocabularies (which obviously expand the number of "codes" assigned to patients, as it should be; first section Results is not novel) that gave approximately 940,000 codes as explained in Materials and methods. However, it is expected that many of these codes are in fact correlated, non-independent, but instead they are all used to generate patient clusters; for example, in several instances it is emphasized that codes from primary diabetes diagnoses are not used in clustering/analyses, but it is well known a priori that many other codes are significantly correlated with the diagnoses (it is therefore not unexpected that patient clusters differentiate diabetes subtypes). Therefore, it would also not be unexpected that clusters differentiate glycemic levels. In fact, the 71 clusters are defined using distances from ICD-10 chapter XVIII symptoms (subsection “Enriched comorbidity and symptom patterns in diabetes patient clusters”, first paragraph), which includes examination of blood parameters. Again, latter the clusters are defined based on comorbidities, which presumably many are linked to the original codes. Another example of this point is the identification of a "schizophrenia" cluster (subsection “Glycemic dysregulation”, last paragraph), which is expected if the corresponding codes are included. I agree that an observation of differential glycemic levels in this clusters may be interesting, but then first it should be evaluated the overall correlation among the thousands of codes used in the study, and also then the study is limited to assess co-occurrences of terms/codes.

There are several separate issues here, some relate to the interdependencies of terms in medical ontologies and some to how such terms can be used to describe subgroups of patients that share features (stratified medicine).

First, the goal in patient stratification is not to find features which are mutually independent, but to identify subgroups in larger patient populations using observable features (e.g. from physical examinations, laboratory tests, genetics or images). However, across subgroups these will typically not be unique or independent feature-wise. The entire medical profession is about combining absent or present features (also those which are mutually exclusive) in ways such that diagnoses, and diseases can be disambiguated. The medical profession is based on the use and recording of these features, and the secondary use of EHR data as in our work is also based on that premise. We do not aim to construct a new ontology with independent features; moreover, such a goal is presumably entirely unrealistic.

Second, patient subgroups are almost always interrelated, triple negative breast cancer patients share many features with ER-positive patients for example. The nested hierarchy we present with its six overall groups and additional subgroups will be interrelated in some way, especially as we go for characterizing all diabetes patients (in a continuum).

The main comment here relates to circularity, and it seems that we have done a bad job in explaining in the Results section why we have left out the primary and secondary diabetes diagnoses A/B. This was however, already included in the Discussion and perhaps the reviewer overlooked that. Contrary to cancer for example, where the ICD-10 diagnoses are quite reliable and highly detailed (lung cancer is not typically confused with liver cancer even if a patient might have both), the primary codes in a systemic, multi-organ disease like diabetes are used extremely broadly and in a fuzzy way, as the knowledge on robust diabetes subtypes and their characteristics in the context of comorbidities is quite limited. This is well-known in the diabetes EHR data mining literature (for example in the work from Vanderbilt). We have now clarified further (Results section “Comorbidity clustering based on text-mined and assigned diagnosis codes”) that we do not want the clustering to be driven by the unreliable, broad, less etiology-relevant codes (E10 and E11) from the ICD endocrinology chapter, but rather by more objectively observed symptoms, other diseases and lifestyle features. It should not come as a surprise, and it is not a sign of circularity, that there is a tendency to segregate into what is today known as type 1 and type 2 diabetes, as these overall subgroups are well established and likely real. It is the further subgrouping and the subgrouping across all diabetes patients that is the result here, not that we rediscover subtypes that have been established for a long time (since the 1920s). We have now made this much clearer in the text. The diabetes diagnoses, e.g. E10 and E11, are obviously assigned by doctors based on the same measurements and observations that we (via the text mining) use as input features. If one were to train a supervised machine learning algorithm to assign diagnoses based on these features, this would obviously be circular. However, this is not what we do. Instead, what we do is to use a non-supervised clustering approach that stratifies the patients in a purely data-driven manner. This allows us to answer two interesting questions: 1) To which extent does it reproduce the classification made by doctors? 2) Does it allow us to make a more specific subgrouping of patients?

One possible way out could be the comparison results of structured vs. structured+unstructured features.

We agree that this is highly relevant and have already in Figure 1 made such a comparison at the ICD-10 level 3. Perhaps the reviewer overlooked that.

2) Clustering

As overall conclusion, Discussion section includes the statement that the results show "many different and homogenous groups of patients", but it has to be formally demonstrated that the groups are independent, that they are many more than randomly expected by permuting the same cohort, and that they are more homogenous that some value.

Additionally, questioning the clusters stability, would be necessary to provide a sensitivity analysis on the clusters to illustrate the stability.

We fully agree that it is highly relevant to include a quantitative, robustness analysis of the clustering. To assess quantitatively the stability of the clusters generated, we constructed various diluted and shuffled realizations of the similarity network used as input to the MCL algorithm. We used a reference clustering similar to the clustering presented in Figure 2B (either by including the patients in the 71 clusters or all patients). The resulting clustering’s were then compared to the reference clustering by means of the Variation of Information measure (VI) (1, 2) and plotted as function of increasing values of α (see Figure 2—figure supplement 5). The figure includes two horizontal lines corresponding to the value that the VI would take if we were to randomly assign 10% and 20% of the vertices to different random clusters, respectively. This analysis showed that the clustering is stable in relation to removing edges, which is evidence that the cosine metric-based cutoff used does not change the overall structure of the clustering. The shuffling is a more impactful randomization, however despite this, we can still shuffle around 10% of the edges and still retrieve 90% of the patients in the groups of the 71 reference clusters. We have added these references:

1) Meilă, 2007. Comparing clusterings—an information based distance. Journal of multivariate analysis, 98(5), pp.873-895.

2) Van Dongen, S., 2000. Performance criteria for graph clustering and Markov cluster experiments. Technical Report, Centre for Mathematics and Computer Science, Amsterdam, The Netherlands.

3) Karrer, Levina, and Newman, 2008. Robustness of community structure in networks. Physical review E, 77(4), p.046119.

3) ICD codes

Truncating the codes to three digits is akin to adding noise by losing information. Please justify your rationale behind truncating the codes. It'll be good if this issue could be addressed in the Discussion or one of the limitations of the work.

It is correct that we in principle and actually lose information when using level 3, but at the same time we would increase tremendously the dimensionality of the feature space if we went to level 4. In that sense, given the size of the data set where many level 4 codes will be unique, we likely would not be able to make a stable subsequent clustering and analysis. Using level 4 would mean that closely related level 4 codes (which would be same level 3 code) would be counted as unrelated. In fact, our attempt to use the much more fine-grained SNOMED-CT terminology confirmed that a data set needs to be very large for such a fine-grained vocabulary to be useful. The conclusion is likely very similar for ICD-10 level 4, in particular when taking the above remarks on the fuzziness of primary diabetes diagnoses into account. We have now elaborated further on the balance between information loss and coverage across patient features in the Discussion (new fifth paragraph).

4) Data limitations and temporal series

First, there might be problems with the time-dependent comorbidities in a given cluster. There are multiple reasons: (1) the length of the observational window before and after considered in a given cluster is not talked about/defined. This can impact the enrichment of symptoms and hence all the downstream analysis; (2) length observational windows is bound to be different in different clusters, which may introduce the bias in computing comorbidities; (3) it is not clear from the analysis if there exist any relation between the observational window and cluster size.

These remarks are all relevant, but we think the reviewer has overlooked Figure 2—figure supplement 1 that describes the “Density of days in contact with SDCC for each cluster”. This additional characterization shows that we have long observational windows for a significant part of the cohort and that no cluster is considerably biased, or mixed time-wise. Given that the Steno Diabetes Center is a specialty clinic for severe diabetes cases where patients are admitted over long periods of time this is not surprising. This was part of the rationale for choosing this cohort.

Indeed, it is not clear how temporality is defined. The heatmap presented in Figure 2B does not help either. Specifically, what is the index date, relative to which pre and postcodes were identified? It is quite likely that the index date would vary for patients, therefore, how does that variability was accounted? Please explain this paragraph in more detail and if possible, then provide a schematic diagram representing temporality and how it was considered in this analysis. In addition, in the fourth paragraph of the subsection “Comorbidity pairs and patterns within symptom related clusters”, it appears the time window before and after is likely to bias the estimation of T2D and elevated blood glucose level. More clarity on this will be helpful.

We assume that the reviewer refers to Figure 4B and Figure 4 in general (and not Figure 2B). The figure is actually only used to describe the order of diagnosis in observed pairs. The actual time-difference analysis is not used elsewhere in the manuscript and we have therefore decided to remove the one sentence that mentions it. Similarly, the few stars that were inserted in the heatmap have been removed.

In any case, a better explanation is needed about the length of the observational window before and after used to compute comorbidities. The length of the observational window (before and after in a given cluster) is bound to confound the disease association. Was this length kept similar for all the clusters? If not, then it is very much possible that certain cultures are (just based on cluster size) are likely to be enriched in certain symptoms. Please provide an in-depth analysis on this.

This was per the comments above included already, or removed.

5) Additional data

Please briefly explain the type of genetic data is included and at what level of omics.

The candidate genetic and biochemical markers of glucose regulation could be of interest if replicated. However, could these have been discovered without the clustering e.g. using simpler regression methods?

Detailed information about the genetic data can be found in the “Genomic characterization” section under Materials and methods. We have inserted a reference to this section in the “Genomic characterization by SNP association of phenotypically determined clusters” section, where we first present results from this analysis to resolve any inaccuracies.

Further, we think the reviewer means glycemic dysregulation and not glucose regulation here. Per the remark below the genetic markers are analyzed without the clustering. We prefer not to go into the question of whether the biochemical markers could be discovered using a simple regression method. When one knows a pattern, or a signal, simpler “rules” or formulas can often be designed subsequently, but we do not think this is the case here.

Finally, we had by mistake left out references to the studies that generated the genetic data. These are now included (Charmet et al., 2018; Steinthorsdottir et al., 2014).

It is unclear how the clusters segregate in their genetic basis, which could then support the identification of disease subtypes. The question is not as much which genetic variants are associated with a given cluster, but if the effect estimations are statistically different between clusters, which would then favor the existence of subtypes. It is also unclear if the association signals are genome-wide significant in any case.

The genetic characterization of dysregulation analysis has been done independent of the clustering. We have realized that the section title “Genetic characterization of dysregulated patients and clusters” could indicate that the analysis is done cluster-wise, which is not the case due to low coverage of genetic data. We have changed the title accordingly to: “Genetic characterization of dysregulated patients”.

We agree that the segregation of genetic variants across clusters could provide additional support for disease subtypes. However, we think that such analysis should be performed in a cohort where the coverage of genetic data in the clusters are near complete. In our current setup, we believe that basing comparisons of effect estimates on varying sample sizes (that is number of patients with genetic data available in each cluster) is strongly under-powered and could lead to spurious results.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The manuscript has been improved but there is still some confusion about the significance of the resulting clusters, significance of the selected thresholds and the subsequent biological interpretation, as outlined below by the first reviewer. The final version should clarify the exploratory nature of the work and the limits of the statistical reliability, as pointed out by the second reviewer.

Reviewer #1:

The authors have done a perceptible effort trying to answer the major concerns initially raised, and I personally appreciate the clarifications and additional data provided.

However, I am still not convinced that the main results are robust enough to support the overall message, in particular regarding the title: "Linking glycemic dysregulation.…to symptoms, comorbidities and genetics…".

In the paper we are exploring how symptoms, comorbidities and genetics are related to glycemic dysregulation, so we would argue that the title is appropriate. In our opinion, the word “linking” is not particularly loaded, or guaranteeing causal relationships and does not exaggerate the merits of the paper. Data linkage is not a term that in the literature is used to indicate causality. We have already in the text emphasized and described its explorative nature. Other titles could obviously be formulated, but we feel it is a bit hard to follow the reviewer’s concern here.

First, the issue of "circularity". I might have over-looked the original Discussion that mentioned that primary and secondary Diabetes diagnoses were excluded, or my initial comment might have been unclear, but I was referring to the expectation that clusters should obviously identify major known Diabetes types; in particular if dozens of diagnosis codes (some of them correlated) and extensive text vocabulary is analyzed. The point is not if the clusters identify what it is already known, but if they identify new etiology/subtypes.

We have not in the paper anywhere stated that the clusters represent diabetes subtypes in terms of etiology. They are descriptive in relation a full continuum diabetes population, and they are presented as a starting point for assessing to what extent the comorbidity/symptom-based subgrouping represents dissimilar etiologies, different or overlapping pathways at the molecular/trajectory levels, including genetic differences. Several recently published subgroupings have been presented as different etiological subtypes even if they presumably to a larger extent represent ethnic differences (when genetics is included in the feature space), or just physiologically similar subgroups that may not share genetics or other mechanistic features. We feel that we have been careful not to follow the same strategy. We present a comorbidity-based subgrouping that is linked to glycemic dysregulation. This subgrouping prompts precision medicine considerations as changes in disease co-occurrences link to differences in glycemic dysregulation and treatment planning. Precision medicine efforts should zoom in on these differences in outcomes, but as the paper is not a multi-omics effort across subgroups it would be unserious to present speculative, mechanistic models and claim them as different diabetes subtypes.

Of note, the following example is somewhat inaccurate: "… triple negative breast cancer patients share many features with ER-positive patients for example." Indeed TNBC and ER+ cases do not share cell-of-origin, prognosis, tissue-metastasis preference, neither therapeutic approaches.

This particular text in the response is referring to phenotypic similarity and, in that regards, TNBC and ER+ indeed share many similarities, both are breast cancers and therefore associated with many of the same types of symptoms to varying degrees. We do not at all disagree with the reviewer, we just highlight that subgroups of patients with diseases that are different also have many similarities, that is the whole point in subgrouping. In our case, we identify patient subgroups of e.g. cancer/diabetes co-occurrences that differ in terms of glycemic regulation. Please note that this text in not in the manuscript file, but only in the response document, and therefore we do not think we need to change our wording. This also links to the remarks about on similarities and dissimilarities in the context of precision medicine above.

My reasoning above is then linked to the title issue: the study of glycemic perturbation is presented in one section (“Glycemic dysregulation”) and remains unclear which clusters, symptoms, are robust regarding statistical differences for "glycemic perturbation" from the others. The thresholds of 0-2, 3-5 parameters, and then the clusters (21 with 50% patients 3 parameters…) appear to be arbitrary. The use of sentences with the following grammatical constructions may not help to understand the analyses: "distinct differentiation…dysregulation parameters".

The comparison between patients fulfilling 0-2 dysregulation parameters with patients fulfilling 3-5 is motivated by the hierarchical clustering of these groups presented in Figure 1—figure supplement 6. Patients in the 0-2 group are more similar to one another than those in the 3-5 group, and vice versa. Next, we chose to examine clusters where the majority of patients are in the “high” dysregulated group and therefore a cutoff of 50% is used. We have now rewritten the text to make this clearer.

Next, glycemic perturbation is linked to genetic risk in the last Results section, but as initially raised, it is unclear what means "top association signals"; are these significant genome-wide, are they statistically different between the clusters/symptoms? Are cluster-risk interactions significant? Honestly, this remains unclear to me and I hope I am not overlooking to any data.

We have rewritten the last part of the Results section to emphasize that the top associating signals do not reach genome-wide significance. We have not performed any comparisons of genetic risk or cluster-risk interactions between the clusters. While the suggestions from the reviewer are highly interesting, unfortunately, the low coverage of genetics for the clusters does not allow us to answer such questions. We have now also added a sentence to the genetics Materials and methods section with the definition of genome-wide significance in terms of p-value.

https://doi.org/10.7554/eLife.44941.sa2

Article and author information

Author details

  1. Isa Kristina Kirk

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    Contributed equally with
    Christian Simon
    Competing interests
    No competing interests declared
  2. Christian Simon

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    Contributed equally with
    Isa Kristina Kirk
    Competing interests
    No competing interests declared
  3. Karina Banasik

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Formal analysis, Methodology, Writing—original draft, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Peter Christoffer Holm

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Formal analysis, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  5. Amalie Dahl Haue

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Formal analysis, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  6. Peter Bjødstrup Jensen

    1. Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    2. Odense Patient Data Explorative Network (OPEN), Odense University Hospital, Odense, Denmark
    Contribution
    Formal analysis, Supervision, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  7. Lars Juhl Jensen

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Data curation, Formal analysis, Supervision, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
  8. Cristina Leal Rodríguez

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  9. Mette Krogh Pedersen

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Methodology
    Competing interests
    No competing interests declared
  10. Robert Eriksson

    Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Methodology
    Competing interests
    No competing interests declared
  11. Henrik Ullits Andersen

    Steno Diabetes Center Copenhagen, Gentofte, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Validation, Writing—review and editing
    Competing interests
    No competing interests declared
  12. Thomas Almdal

    1. Steno Diabetes Center Copenhagen, Gentofte, Denmark
    2. Department of Endocrinology, Rigshospitalet, Copenhagen, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  13. Jette Bork-Jensen

    Novo Nordisk Foundation Center for Basic Metabolic Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  14. Niels Grarup

    Novo Nordisk Foundation Center for Basic Metabolic Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  15. Knut Borch-Johnsen

    Holbæk Hospital, Holbæk, Denmark
    Contribution
    Conceptualization, Resources, Supervision, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
  16. Oluf Pedersen

    1. Steno Diabetes Center Copenhagen, Gentofte, Denmark
    2. Novo Nordisk Foundation Center for Basic Metabolic Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  17. Flemming Pociot

    1. Steno Diabetes Center Copenhagen, Gentofte, Denmark
    2. Department of Clinical Medicine, Herlev-Gentofte Hospital, Herlev, Denmark
    Contribution
    Resources, Data curation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  18. Torben Hansen

    1. Steno Diabetes Center Copenhagen, Gentofte, Denmark
    2. Novo Nordisk Foundation Center for Basic Metabolic Research, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Supervision, Writing—review and editing
    Competing interests
    No competing interests declared
  19. Regine Bergholdt

    Steno Diabetes Center Copenhagen, Gentofte, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Methodology, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
  20. Peter Rossing

    1. Steno Diabetes Center Copenhagen, Gentofte, Denmark
    2. Department of Clinical Medicine, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Supervision, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    peter.rossing@regionh.dk
    Competing interests
    No competing interests declared
  21. Søren Brunak

    1. Novo Nordisk Foundation Center for Protein Research, University of Copenhagen, Copenhagen, Denmark
    2. Center for Biological Sequence Analysis, Department of Bio and Health Informatics, Technical University of Denmark, Lyngby, Denmark
    Contribution
    Conceptualization, Resources, Data curation, Funding acquisition, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    soren.brunak@cpr.ku.dk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0316-5866

Funding

Danish Council for Strategic Research (0603-00321B)

  • Søren Brunak

Innovation Fund Denmark (5153-00002B)

  • Søren Brunak

Novo Nordisk Foundation (NNF14CC0001)

  • Søren Brunak

Novo Nordisk Foundation (NNF17OC0027594)

  • Søren Brunak

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Alfonso Valencia, Barcelona Supercomputing Center - BSC, Spain

Publication history

  1. Received: January 17, 2019
  2. Accepted: November 16, 2019
  3. Version of Record published: December 10, 2019 (version 1)

Copyright

© 2019, Kirk et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 373
    Page views
  • 48
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    Erika E Kuchen et al.
    Research Article
    1. Computational and Systems Biology
    2. Physics of Living Systems
    Eugenia Lyashenko et al.
    Research Article