Exploratory electronic health record analysis with ehrapy

Ehrapy: a framework for exploratory and targeted EHR data analysis

The foundation of ehrapy is a robust and scalable implementation-independent data storage backend which is combined with a series of preprocessing and analysis modules. In ehrapy, EHR data is organized as a data matrix where observations are individual patient visits (or patients, in the absence of follow-up visits), and variables represent all measured quantities (Methods). These data matrices are stored together with metadata of observations, variables, and unstructured annotations. By leveraging the open AnnData (annotated data) data structure that implements this design, ehrapy builds upon established, scalable standards and is compatible with analysis and visualization functions provided by the single-cell omics scverse³⁰ ecosystem. Ehrapy thus allows for larger-than-memory analyses by leveraging HDF5-³⁶ or Zarr-based³⁷ AnnData objects. Both formats are platform, framework and programming language independent formats and readers are also available in R, Julia and Javascript³⁸. We additionally provide a dataset module with more than 20 public loadable EHR datasets in AnnData format to kickstart analysis and development with ehrapy.

For standardized analysis of EHR data, it is crucial that these data are encoded and stored in consistent, reusable formats. Thus, ehrapy requires that input data is organized in structured vectors. Readers for common formats such as CSV, OMOP³⁹, or SQL databases are available in ehrapy. Data loaded into AnnData objects can be mapped against hierarchical ontologies such as MONDO⁴⁰, International Statistical Classification of Diseases (ICD)⁴¹, Human Phenotype Ontology⁴², or Phecodes⁴³ (Methods). Clinical keywords of free text notes can be automatically extracted (Methods).

Powered by scanpy, which scales to hundreds of millions of observations⁴⁴ and the machine learning library scikit-learn⁴⁵, ehrapy provides more than 100 composable analysis functions organized in modules from which custom analysis pipelines can be built. In such pipelines, each function directly interacts with the AnnData object and adds all intermediate results to the object for simple access and reuse of information. To facilitate setting up these pipelines, ehrapy guides analysts through a general analysis pipeline (Figure 1) consisting of data preparation, data preprocessing and knowledge inference. At any step of an analysis pipeline, community software packages can be injected without any vendor lock-in. Since ehrapy is built on open standards, it can be purposefully extended to solve new challenges. Ehrapy natively interfaces with the scientific Python ecosystem via Pandas⁴⁶ and Numpy⁴⁷. The development of deep-learning models for EHR data⁴⁸ is further accelerated through compatibility with pathml⁴⁹, a unified framework for whole-slide image analysis in pathology, and scvi-tools⁵⁰ which provides data loaders for loading tensors from AnnData objects into PyTorch⁵¹ or Jax arrays⁵² to facilitate the development of generalizing foundational models for medical artificial intelligence⁵³.

• Download figure
• Open in new tab

Figure 1.

Schematic overview of ehrapy analysis. (A) Heterogeneous health data are first stored in memory as an Annotated Data (AnnData) object with patients or patient visits as observations and variables as columns. Next, the data are first mapped against ontologies and key terms are extracted from free text notes. (B) The EHR data are subject to quality control where impossible or spurious measurements are removed. Subsequently, categorical data are encoded and missing data are imputed. Data from different sources with data distribution shifts are integrated, embedded, clustered and annotated in a patient landscape. (C) Further downstream analyses depend on the question of interest and can include the inference of causal effects and trajectories, survival analysis, or patient stratification.

In the general ehrapy analysis pipeline, EHR data is first examined to detect quality issues by inspecting feature distributions (Methods). Strong outliers determined this way may represent erroneous measurements that could skew any calculations and can therefore be winsorized, clipped, or removed. Automatic quality control further detects visits and features with high missing rates that ehrapy can then impute (Methods). Subsequently, ehrapy’s normalization and encoding functions (Methods) are used to ensure a consistent numerical data representation that enables data integration while regressing out dataset shift effects (Methods). The size and complexity of EHR data requires lower dimensional representations that can subsequently be visualized, clustered and annotated to obtain a patient landscape (Methods). Such annotated groups of patients can be used for statistical comparisons to find differences in features among them to ultimately learn markers of patient states.

As analysis goals can differ between users and datasets, the general ehrapy analysis pipeline is opened up in the final knowledge inference step such that the right modules can be chosen to best answer the question at hand. Ehrapy provides advanced statistical methods for group comparison and linear models together with extensive support for survival analysis (Methods) enabling the discovery of biomarkers. Furthermore, ehrapy provides functions for causal inference to go from statistically determined associations to causal relations (Methods). Disease progression analysis represents a further analysis goal for EHR data that ehrapy also provides support for in the form of trajectory inference methods (Methods). Patient visits in aggregated EHR data can be regarded as snapshots, offering a limited perspective that may not encompass the entirety of a patient’s clinical picture or the progression of their disease. Individual measurements taken at specific time points might not adequately reflect the underlying progression of disease but result from unrelated variation due to, for example, day-to-day differences^54–56. Therefore, disease progression models should rely on analysis of the underlying clinical data rather than of a notion of time as disease progression in an individual patient may not be monotonous in time. Ehrapy’s integration into the larger scverse ecosystem allows for the usage of advanced trajectory inference methods to facilitate novel views on EHR data through compatible tools like Diffusion Pseudotime⁵⁷ or CellRank⁵⁸. We show that these tools can order snapshots to calculate a pseudotime that can adequately reflect the progression of the underlying clinical process without requiring completeness of real-time measurements in all patients. Given a sufficient number of snapshots, ehrapy increases the potential to capture disease progression that is likely not robustly captured within a single EHR record, but rather across several.

Ehrapy enables patient stratification in unspecified pneumonia cases

To demonstrate ehrapy’s capability to analyze heterogeneous datasets from a broad patient set across multiple care units, we applied our novel strategy to the Paediatric Intensive Care (PIC)³³ database. The PIC database is a single-center, bilingual (English and Chinese) database hosting information on children admitted to critical care units at the Children’s Hospital of Zhejiang University School of Medicine in China. It contains 13,499 distinct hospital admissions of 12,881 individual pediatric patients admitted between 2010 and 2018 for which demographics, diagnoses, doctors’ notes, vital signs, laboratory and microbiology tests, medications, fluid balances, and more were collected (Figure 2A, Supplementary Figure 1, Methods). Following missing data imputation and subsequent preprocessing (Figure 2B-C, Supplementary Figure 2, Methods), we generated a Uniform Manifold Approximation and Projection (UMAP) embedding to visualize variation across all patients using ehrapy (Figure 3A). This visualization of a low dimensional patient manifold shows the heterogeneity of the collected data in the PIC database with malformations, perinatal, and respiratory being the most abundant ICD chapters (Figure 3B).

• Download figure
• Open in new tab

Figure 2.

Preprocessing of the Paediatric Intensive Care (PIC) dataset with ehrapy. (A) Heterogeneous data of the PIC database was stored in “data” (matrix that is used for computations) and “observations” (metadata per patient). During quality control further annotations are added to the “variables” (metadata per feature) slot. (B) Preprocessing steps of the PIC dataset. (C) Minimal example of the function calls in the data analysis pipeline that resembles the preprocessing steps in (B) using ehrapy.

• Download figure
• Open in new tab

Figure 3.

PIC dataset overview and annotation of patients diagnosed with unspecified pneumonia. (A) Uniform manifold approximation and projection (UMAP) of all patient visits in the intensive care unit (ICU) with primary discharge diagnosis grouped by ICD chapter. The organ heterogeneity in the malformation chapter prevents very distinct clusters per chapter. (B) The abundance and complexity of respiratory diseases encouraged us to investigate them further. (C) Respiratory categories show the abundance of influenza and pneumonia diagnoses that we investigated further. (D) We observed the “unspecified pneumonia” sub-group which led us to investigate and annotate it further in more detail. (D) The previously “unspecified pneumonia” labeled patients were annotated using several clinical features (Supplementary Figure 3) of which the most important ones are shown in the heatmap (E). (G) Example disease progression of an individual child with pneumonia illustrating pharmacotherapy over time until positive acinetobacter baumannii swab.

The most common respiratory disease categories (Figure 3C) were pneumonia and influenza (n=1229). We focused on pneumonia to apply ehrapy to a challenging, broad-spectrum disease across age. Pneumonia is a prevalent respiratory infection that poses a substantial burden on public health⁵⁹ and is characterized by inflammation of the alveoli and distal airways⁵⁹. Individuals with pre-existing chronic conditions are particularly vulnerable, as are children under the age of five⁶⁰. Pneumonia can be caused by a range of microorganisms, encompassing bacteria, respiratory viruses, and fungi.

From these cases, we selected the age group “youths” (13 months to 18 years of age) for further analysis, thereby addressing a total of 265 patients. Neonates (0 to 28 days) and infants (28 days to 13 months old) were excluded from the analysis as the disease context is significantly different in these age groups due to distinct anatomical and physical conditions. These 265 patients dominated the group of pneumonia cases and were diagnosed with “unspecified pneumonia” by the collectors of the original database (Figure 3D).

Patients were 61% male, had a total of 277 admissions, a mean age at admission of 54 months (median of 38 months) and an average length of stay (LOS) of 15 days, with a median of 7 days. Of these, 152 patients were admitted to the pediatric ICU (PICU), 118 to the general ICU (GICU), four to the surgical ICU (SICU), and three to the cardiac ICU (CICU). Laboratory measurements generally contained between 12% and 14% missing data with serum procalcitonin (PCT) (24.5%) and C-reactive protein (CRP) (16.8%) being notable exceptions. Parameters assigned as “vital signs” contained between 44% and 54% missing values. Stratifying patients with unspecified pneumonia further enables a more nuanced understanding of the disease, potentially facilitating tailored approaches to treatment.

To deepen clinical phenotyping for the disease group “unspecified pneumonia”, we calculated a k-nearest-neighbor graph to cluster patients into subgroups and visualize these in UMAP space, thereby progressing to the next step in ehrapy’s workflow (Figure 1, Methods). Leiden clustering⁶¹ identified four patient groupings with distinct clinical features that we annotated (Figure 3E).

To identify the laboratory values, medications, and pathogens that were most characteristic for these four groups (Figure 3F), we applied t-tests for numerical and g-tests for categorical data using ehrapy (Supplementary Figure 3, Methods). Based on this analysis, we identified patient subgroups with “sepsis-like, “severe pneumonia with co-infection”, “viral pneumonia” and “mild pneumonia” phenotypes.

The “sepsis-like” group of patients (n=28) were characterized by rapid disease progression as exemplified by an increased number of deaths (adj. P≤ 5.04e-03), indication of multiple organ failure such as elevated creatinine (adj. P ≤ 0.01) or reduced albumin levels (adj. P ≤ 2.89e-04) together with increased expression levels and peaks of inflammation markers including PCT (adj. P ≤ 3.01e-02), whole blood cell count, neutrophils, lymphocytes, monocytes, lower platelet counts (adj. P ≤ 6.3e-02) and changes in electrolyte levels, i.e. lower potassium levels (adj. P ≤ 0.09e-02). Patients we associated with the term “severe pneumonia with co-infection” (n=74) were characterized by prolonged ICU stays (adj. P≤ 3.59e-04), organ affection such as higher levels of creatinine (adj. P ≤ 1.10e-04), increased inflammation markers such as peaks of PCT (adj. P≤ 5.06e-05), CRP (adj. P≤ 1.40e-06), and neutrophils (adj. P≤ 8.51e-06), together with valleys of lower platelet count (adj. P ≤ 5.40e-23), detection of bacteria in combination with additional pathogens such as fungals in sputum samples (adj. P≤ 1.67e-02), and increased application of medication including antifungals (adj. P≤ 1.30e-04) and catecholamines (adj. P≤ 2.0e-02)). Patients in the “mild pneumonia” cluster were characterized by positive sputum cultures in the presence of relatively lower inflammation markers such as PCT (adj. P≤ 1.63e-03) and CRP (adj. P≤ 0.03e-01), while receiving antibiotics more frequently (adj. P ≤ 1.00e-05), and additional medications (electrolytes, blood thinners and circulation supporting medications) (adj P≤ 1.00e-05). Finally, patients in the“viral infection” cluster were characterized by shorter lengths of stays (adj. P ≤ 8.00e-06), a lack of non-viral pathogen detection in combination with higher lymphocyte counts (adj. P≤ 0.01 avg) and peaks of PCT (adj. P ≤ 0.03e-02), reduced application of catecholamines (adj. P≤ 5.96e-07), antibiotics (adj. P≤ 8.53e-06), antifungals (adj. P≤ 5.96e-07).

To demonstrate the ability of ehrapy to examine EHR data from different levels of resolution, we additionally reconstructed a case from the “severe pneumonia with co-infection” group (Figure 3G). In this case, the analysis revealed that CRP levels remained elevated despite broad spectrum antibiotic treatment until a positive “Acinetobacter baumanii” result led to a change in medication and a subsequent decrease in CRP and monocyte levels.

Ehrapy empowers extraction of clinical indicators for unspecified pneumonia

Ehrapy’s survival analysis module allowed us to identify clinical indicators of disease or disease stages that could be used as biomarkers through Kaplan-Meier analysis. We found strong variance in overall aspartate aminotransferase (AST), alanine aminotransferase (ALT), and gamma-glutamyl transferase (GGT) levels (Figure 4A) including changes over time (Supplementary Figure 4A-B) in all four “unspecified pneumonia” groups. Routinely used to assess liver function, studies provide evidence that AST, ALT, GGT levels are elevated during respiratory infections^62,63 including severe pneumonia ⁶⁴, and can guide diagnosis and management of pneumonia in children⁶².

• Download figure
• Open in new tab

Figure 4.

Survival analysis of patients diagnosed with unspecified pneumonia. (A) Line plots of major hepatic system lab measurements per group show variance in the measurements per pneumonia group. (B) Kaplan-Meier survival curves demonstrate lower survival for “sepsis-like” and “severe pneumonia with co-infection” groups. (C) Kaplan-Meier survival curves for children with Gamma Glutamyltransferase measurements outside the norm range display lower survival.

We confirmed reduced survival in more severely affected children (“sepsis-like pneumonia” and “severe pneumonia with co-infection”) using Kaplan-Meier curves and a multivariate log-rank test (Figure 4B; P ≤1.09e-18) through ehrapy. To verify the association of this trajectory with altered AST, ALT, GGT expression levels, we further grouped all patients based on liver enzyme reference ranges (Methods, Table 1). By Kaplan-Meier survival analysis, cases with peaks of GGT (P ≤ 1.4e-02), ALT (P ≤ 2.9e-02), and AST (P ≤ 4.8e-04) in “outside the norm” were found to correlate with lower survival in all groups (Figure 4C, Supplementary Figure 4), in line with previous studies^62,63.

Ehrapy allows for medication-class effect quantification of length of stay by causal inference

Pneumonia requires case-specific medications due to its diverse causes. To demonstrate the potential of ehrapy’s causal inference module, we quantified the effect of medication on ICU length of stay to evaluate case-specific administration of medication. In contrast to causal discovery that attempts to find a causal graph reflecting the causal relationships, causal inference is a statistical process used to investigate possible effects when altering a provided system, as represented by a causal graph and observational data (Figure 5A)⁶⁵. This approach allows researchers to identify and quantify the impact of specific interventions or treatments on outcome measures, thereby providing insight for evidence-based decision-making in healthcare. Causal inference relies on datasets incorporating multiple perturbation conditions, commonly termed “interventions”, to accurately quantify effects.

• Download figure
• Open in new tab

Figure 5.

Causal inference of length of stay affected by different medication types. (A) ehrapy’s causal module is based on the strategy of the tool dowhy. Here, EHR data containing treatment, outcome, and measurements, and a causal graph serve as input for causal effect quantification. The process includes the identification of the target estimand based on the causal graph, the estimation of causal effects using various models, and finally refutation where sensitivity analyses and refutation tests are performed to assess the robustness of the results and assumptions. (B) Curated causal graph using age, liver damage, and inflammation markers as disease progression proxies together with medications as interventions to assess the causal effect on length of ICU stay. (C) Determined causal effect strength on length of stay in days of administered medication categories showcases lower expected length of stays in cases where corticosteroids, carbapenems, and antivirals are administered.

We manually constructed a minimal causal graph with ehrapy (Figure 5B) on records of treatment with corticosteroids, carbapenems, penicillins, cephalosporins, antifungal, and antiviral medication (Supplementary figure 5). We used these as “interventions’’, and assumed the medications affect disease progression proxies such as inflammation markers and markers of organ function. The selection of “interventions” is consistent with current treatment standards for bacterial pneumonia and respiratory distress^66,67.

Based on the approach of the tool dowhy⁶⁸ (Figure 5A), ehrapy’s causal module identified the application of corticosteroids, antivirals, and carbapenems to be associated with shorter lengths of stay, in line with current evidence^60,69,70,71. In contrast, penicillins and cephalosporins were associated with longer lengths of stay (Figure 5C), whereas antifungal medication did not strongly influence length of stay (Figure 5c).

Ehrapy facilitates automated extraction of population-scale risk factors in the UK-Biobank

To demonstrate ehrapy’s scalability and to illustrate the advantages of using a unified data management and quality control framework like ehrapy when performing multimodal survival modeling, we modeled myocardial Infarction risk using Cox Proportional-Hazards Models on UK-Biobank (UKB)³⁴ data. Large population cohort studies like the UK-Biobank enable the investigation of common diseases across a wide range of modalities including genomics, metabolomics, proteomics, imaging data, and common clinical variables (Figure 6A-B). From these, we used a publicly available polygenic risk score for coronary heart disease⁷² comprising 6.6 million variants, 80 nuclear magnetic resonance (NMR) spectroscopy-based metabolomics⁷³ features, 81 features derived from retinal optical coherence tomography^74,75, and the Framingham Risk Score⁷⁶ feature set, which includes known clinical predictors such as age, sex, BMI, blood pressure, smoking behavior, and cholesterol levels. We filtered the feature space by excluding features with more than 10% missingness in the respective population and imputed the remaining missing values (Methods). Further, individuals with comorbidities were excluded from the analyses. We selected 29,216 individuals for whom all of these data types were available (Supplementary Figure 6-7, Methods). Myocardial Infarction, as defined by our mapping to the Phecode nomenclature⁴³, was defined as the endpoint (Figure 6C). We modeled the risk for myocardial Infarction one year after either the metabolomic sample was obtained or imaging was performed.

• Download figure
• Open in new tab

Figure 6.

Analysis of myocardial Infarction risk in the UK-Biobank. (A) The UK-Biobank includes 502,359 participants from 22 assessment centers distributed over the United Kingdom. Most participants have genetic (97 %) and general physical measurement data (93 %) like the Framingham Score features available, while the more complex measures like metabolomics, retinal imaging or proteomics data are only available for a much smaller subset. (B) Ehrapy’s quality control workflow shows a distinct cluster of individuals (bottom right) from the Birmingham assessment center in a UMAP embedding of the retinal imaging data, which we further investigated and found to be artifacts of the image acquisition process, and thus ultimately excluded from subsequent analyses. (C) Myocardial Infarctions are recorded for 15 % of the male, and 7 % of the female study population. (D) For every modality and their respective combinations a linear Cox-Proportional hazards model was fitted to determine the prognostic potential of these for the myocardial Infarction phenotype. Commonly used cardiovascular risk factors show expected positive log-hazard ratios for increased blood pressure or total cholesterol, and negative ones for sampling age and female sex. (E) The predictive performance of individual modalities (dark blue) varies from a C-index of 0.63 to 0.77 for genetic and Framingham Risk respectively. Combining all features yields a C-index of 0.81.

Predictive performance for each modality was assessed by fitting Cox proportional hazards (Figure 6C) models on each of the feature sets using ehrapy (Figure 6D). The age of the first occurrence served as the time-to-event; alternatively, date of death or date of the last record in EHR served as censoring times. Models were evaluated using the Concordance Index (C-Index). Here, the combination of multiple modalities successfully improved the predictive performance for coronary heart disease by increasing the C-index from 0.63 (genetic), to 0.76 (genetics, age and sex) and to 0.77 (clinical predictors) with 0.81 (imaging and clinical predictors) for combinations of feature sets (Figure 6E). Our finding is in line with previous observations of complementary effects between different modalities, where a broader ‘major adverse cardiac event’ phenotype was modeled in the UKB with a neural network achieving a C-index of 0.72 using metabolomics, age and sex⁷⁷. By analyzing the predictive potential of diverse data modalities alone and in combination through the same ehrapy analysis pipeline, we directly assessed the overlap of information contained in either one of them. We observe that adding genetic data improves the predictive potential in all cases because it is the only modality independent of the sampling age and only predicts the other modalities to a limited extent⁷⁸. We also demonstrate that the addition of metabolomic data to the combination of clinical predictors, genetic information and imaging-derived measures did not improve predictive power (Figure 6E).

Imaging-based disease severity projection through fate mapping

To demonstrate ehrapy’s ability to handle diverse image data and recover disease stages, we embedded pulmonary imaging data obtained from COVID-19 patients into a lower dimensional space and computationally inferred disease progression trajectories using pseudotemporal ordering. This describes a continuous trajectory or ordering of individual points based on feature similarity⁷⁹, and has been extensively used e.g. in cellular analyses to recover opaque trajectories or lineages⁵⁷. Continuous trajectories enable the map the fate of new patients onto precise states to potentially predict their future condition.

In COVID-19, a highly contagious respiratory illness caused by the coronavirus SARS-CoV-2 with symptoms ranging from mild flu-like symptoms to severe respiratory distress, chest X-ray images were described to exhibit opacities (bilateral patchy, ground glass) and consolidation in association with disease severity⁸⁰.

We used COVID-19 chest X-ray images from the BrixIA⁸¹ dataset consisting of 192 images (Figure 7A) with expert annotations of disease severity. We made use of the scoring procedure of the BrixIA database which is based on six regions that were annotated by radiologists for disease severity ultimately resulting in a global score for disease severity classification (Methods). We embedded raw image features using a pre-trained DenseNet model (Methods) and further processed this embedding into a nearest-neighbors based UMAP space by the use of ehrapy (Figure 7A, Methods). Fate mapping based on imaging information (Methods) determined a severity ordering from mild to critical cases (Figure 7B-D). Images labeled as “normal” are projected to stay within the healthy group illustrating the robustness of our approach. Images of diseased patients were ordered by disease severity highlighting clear trajectories from “normal” to “critical” states despite the heterogeneity of the X-ray images stemming from e.g., different zoom (Figure 7A). To validate this raw image feature-based approach, we applied a similar workflow (Methods) to 790 eye retina images of patients affected by diabetic retinopathy and confirmed ehrapy’s potential to recover disease severity trajectories while differentiating severe and proliferative cases (Supplementary Figure 8). Here, we are able to additionally separate “Severe” from “Proliferate” cases into distinct clusters with unique trajectories.

• Download figure
• Open in new tab

Figure 7.

Recovery of disease severity trajectory in COVID 19 chest x-ray images. (A) Randomly selected example chest X-ray images of the BrixIA dataset highlight variance of BrixIA images. To determine fates, raw images are used as input for machine learning models that project them into a meaningful lower dimensional space. Using this representation, ehrapy determines pseudotime. Both, a nearest-neighbor graph and the determined pseudotime serve as input for the calculation of a transition matrix to enforce directionality. (B) UMAP visualization of the BrixIA dataset embedding shows a clear separation of disease classes. Disease severity was annotated by a senior radiologist of the BrixIA COVID-19 Project. (C) Calculated pseudotime for all images increases with distance to the “Normal” images. (D) Stream projection of fate mapping in UMAP space showcases disease severity trajectory of the COVID19 chest x-ray images.

Implementation of ehrapy

ehrapy has been implemented in the Python programming language and builds upon a number of existing numerical and scientific open-source libraries, specifically, matplotlib⁸⁹, seaborn⁹⁰, NumPy⁴⁷, numba⁹¹, Scipy⁹², scikit-learn⁴⁵, and Pandas⁹³. All functions are grouped into task specific modules whose implementation is complemented with additional dependencies.

Paediatric Intensive Care database analysis

UK Biobank analysis

COVID-19 Chest-x-ray fate determination

Diabetic retinopathy fate determination