MIMIC-II IAC Introduction#
This tutorial explores the Multiparameter Intelligent Monitoring in Intensive Care II (MIMIC-II) Indwelling Arterial Catheters (IAC) dataset, as subset derived from MIMIC-II, the publicly-accessible critical care database. The database was created for the purpose of a case study in the book “Secondary Analysis of Electronic Health Records”, published by Springer in 2016. In particular, the MIMIC-II IAC dataset was used throughout Chapter 16 (Data Analysis) by Raffa J. et al. to investigate the effectiveness of indwelling arterial catheters in hemodynamically stable patients with respiratory failure for mortality outcomes.
More details on the dataset such as all included features and their description can be found here.
import warnings
warnings.filterwarnings("ignore")
from IPython.display import Image
Image(filename="images/MIMIC-II-database-structure.png", width=400)
In this tutorial we want to explore the MIMIC-II IAC dataset using ehrapy to identify patient groups and their associated features.
The major steps of an analysis with ehrapy include:
Preprocessing and quality control (QC)
Dimensionality reduction
Batch effect identification
Clustering
Additional downstream analysis
Before we start with the analysis of the MIMIC-II IAC dataset, we set up our environment including the import of packages and preparation of the dataset.
Environment setup#
Ensure that the latest version of ehrapy is installed. A list of all dependency versions can be found at the end of this tutorial.
import warnings
warnings.filterwarnings("ignore")
import ehrapy as ep
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
MIMIC-II IAC dataset loading#
ehrapy offers several datasets in AnnData format that can be used out of the box. In this tutorial we will use the MIMIC-II IAC dataset with unencoded features. ehrapy’s default encoding is a simple one-hot encoding in this case. More details on encoding can be seen in the next step.
adata = ep.dt.mimic_2(encoded=False)
adata
AnnData object with n_obs × n_vars = 1776 × 46
layers: 'original'
The MIMIC-II dataset has 1776 patients with 46 features.
Now that we have our AnnData file ready, we can start the analysis using ehrapy and the first step will be to preprocess the dataset.
Analysis using ehrapy#
Preprocessing#
ep.ad.infer_feature_types(adata)
! Features 'aline_flg', 'gender_num', 'service_num', 'day_icu_intime_num', 'hour_icu_intime', 'hosp_exp_flg', 'icu_exp_flg', 'day_28_flg', 'censor_flg', 'sepsis_flg', 'chf_flg', 'afib_flg', 'renal_flg', 'liver_flg', 'copd_flg', 'cad_flg', 'stroke_flg', 'mal_flg', 'resp_flg' were detected as categorical features stored numerically.Please verify and correct using `ep.ad.replace_feature_types` if necessary.
Detected feature types for AnnData object with 1776 obs and 46 vars ╠══ 📅 Date features ╠══ 📐 Numerical features ║ ╠══ abg_count ║ ╠══ age ║ ╠══ bmi ║ ╠══ bun_first ║ ╠══ chloride_first ║ ╠══ creatinine_first ║ ╠══ hgb_first ║ ╠══ hospital_los_day ║ ╠══ hr_1st ║ ╠══ icu_los_day ║ ╠══ iv_day_1 ║ ╠══ map_1st ║ ╠══ mort_day_censored ║ ╠══ pco2_first ║ ╠══ platelet_first ║ ╠══ po2_first ║ ╠══ potassium_first ║ ╠══ sapsi_first ║ ╠══ sodium_first ║ ╠══ sofa_first ║ ╠══ spo2_1st ║ ╠══ tco2_first ║ ╠══ temp_1st ║ ╠══ wbc_first ║ ╚══ weight_first ╚══ 🗂️ Categorical features ╠══ afib_flg (2 categories) ╠══ aline_flg (2 categories) ╠══ cad_flg (2 categories) ╠══ censor_flg (2 categories) ╠══ chf_flg (2 categories) ╠══ copd_flg (2 categories) ╠══ day_28_flg (2 categories) ╠══ day_icu_intime (7 categories) ╠══ day_icu_intime_num (7 categories) ╠══ gender_num (2 categories) ╠══ hosp_exp_flg (2 categories) ╠══ hour_icu_intime (24 categories) ╠══ icu_exp_flg (2 categories) ╠══ liver_flg (2 categories) ╠══ mal_flg (2 categories) ╠══ renal_flg (2 categories) ╠══ resp_flg (2 categories) ╠══ sepsis_flg (1 categories) ╠══ service_num (2 categories) ╠══ service_unit (3 categories) ╚══ stroke_flg (2 categories)
Categorical features could either already be stored numerically (e.g., as 0/1 for flags) or as another type such as strings. Such categorical features need an encoding. Here, we identify service_unit and day_icu_intime as categorical features stored non-numerically. We will therefore encode them first with one-hot encoding. This ensures that no ordering is preserved for the respective features. ehrapy also offers other encoding functions.
adata = ep.pp.encode(adata, encodings={"one-hot": ["service_unit", "day_icu_intime"]})
adata
AnnData object with n_obs × n_vars = 1776 × 54
obs: 'service_unit', 'day_icu_intime'
var: 'feature_type', 'unencoded_var_names', 'encoding_mode'
layers: 'original'
ep.ad.feature_type_overview(adata)
Detected feature types for AnnData object with 1776 obs and 54 vars ╠══ 📅 Date features ╠══ 📐 Numerical features ║ ╠══ abg_count ║ ╠══ age ║ ╠══ bmi ║ ╠══ bun_first ║ ╠══ chloride_first ║ ╠══ creatinine_first ║ ╠══ hgb_first ║ ╠══ hospital_los_day ║ ╠══ hr_1st ║ ╠══ icu_los_day ║ ╠══ iv_day_1 ║ ╠══ map_1st ║ ╠══ mort_day_censored ║ ╠══ pco2_first ║ ╠══ platelet_first ║ ╠══ po2_first ║ ╠══ potassium_first ║ ╠══ sapsi_first ║ ╠══ sodium_first ║ ╠══ sofa_first ║ ╠══ spo2_1st ║ ╠══ tco2_first ║ ╠══ temp_1st ║ ╠══ wbc_first ║ ╚══ weight_first ╚══ 🗂️ Categorical features ╠══ afib_flg (2 categories) ╠══ aline_flg (2 categories) ╠══ cad_flg (2 categories) ╠══ censor_flg (2 categories) ╠══ chf_flg (2 categories) ╠══ copd_flg (2 categories) ╠══ day_28_flg (2 categories) ╠══ day_icu_intime (7 categories); one-hot encoded ╠══ day_icu_intime_num (7 categories) ╠══ gender_num (2 categories) ╠══ hosp_exp_flg (2 categories) ╠══ hour_icu_intime (24 categories) ╠══ icu_exp_flg (2 categories) ╠══ liver_flg (2 categories) ╠══ mal_flg (2 categories) ╠══ renal_flg (2 categories) ╠══ resp_flg (2 categories) ╠══ sepsis_flg (1 categories) ╠══ service_num (2 categories) ╠══ service_unit (3 categories); one-hot encoded ╚══ stroke_flg (2 categories)
Quality Control (QC)#
Demographics distribution#
To see if we have strong differences by demographics, we can check these features in a violin plot.
ep.settings.set_figure_params(figsize=(4, 3), dpi=100)
ep.pl.violin(adata, keys=["age"], groupby="service_unit")
Missing values#
ehrapy’s pp.qc_metrics() function will calculate several useful metrics such as the absolute number and percentages of missing values and properties like the mean/median/min/max of all features. The percentage of missing values is important as features with too many missing values should not be included.
obs_metric, var_metrics = ep.pp.qc_metrics(adata)
obs_metric
| missing_values_abs | missing_values_pct | |
|---|---|---|
| 0 | 0 | 0.000000 |
| 1 | 12 | 22.222222 |
| 2 | 0 | 0.000000 |
| 3 | 3 | 5.555556 |
| 4 | 0 | 0.000000 |
| ... | ... | ... |
| 1771 | 1 | 1.851852 |
| 1772 | 1 | 1.851852 |
| 1773 | 3 | 5.555556 |
| 1774 | 1 | 1.851852 |
| 1775 | 1 | 1.851852 |
1776 rows × 2 columns
var_metrics
| missing_values_abs | missing_values_pct | mean | median | standard_deviation | min | max | iqr_outliers | |
|---|---|---|---|---|---|---|---|---|
| ehrapycat_service_unit_FICU | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_service_unit_MICU | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_service_unit_SICU | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Friday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Monday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Saturday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Sunday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Thursday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Tuesday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| ehrapycat_day_icu_intime_Wednesday | 0 | 0.000000 | NaN | NaN | NaN | NaN | NaN | True |
| aline_flg | 0 | 0.000000 | 0.554054 | 1.000000 | 0.497070 | 0.000000 | 1.000000 | False |
| icu_los_day | 0 | 0.000000 | 3.346498 | 2.185000 | 3.355316 | 0.500000 | 28.240000 | True |
| hospital_los_day | 0 | 0.000000 | 8.110923 | 6.000000 | 8.154862 | 1.000000 | 112.000000 | True |
| age | 0 | 0.000000 | 54.379660 | 53.678585 | 21.056923 | 15.180230 | 99.110947 | False |
| gender_num | 1 | 0.056306 | 0.577465 | 1.000000 | 0.493963 | 0.000000 | 1.000000 | False |
| weight_first | 110 | 6.193694 | 80.075948 | 77.000000 | 22.483765 | 30.000000 | 257.600006 | True |
| bmi | 466 | 26.238739 | 27.827316 | 26.324846 | 8.206940 | 12.784877 | 98.797134 | True |
| sapsi_first | 85 | 4.786036 | 14.136606 | 14.000000 | 4.113085 | 3.000000 | 32.000000 | True |
| sofa_first | 6 | 0.337838 | 5.820904 | 6.000000 | 2.334006 | 0.000000 | 17.000000 | True |
| service_num | 0 | 0.000000 | 0.552928 | 1.000000 | 0.497191 | 0.000000 | 1.000000 | False |
| day_icu_intime_num | 0 | 0.000000 | 4.054054 | 4.000000 | 1.993911 | 1.000000 | 7.000000 | False |
| hour_icu_intime | 0 | 0.000000 | 10.585586 | 9.000000 | 7.922733 | 0.000000 | 23.000000 | False |
| hosp_exp_flg | 0 | 0.000000 | 0.137387 | 0.000000 | 0.344256 | 0.000000 | 1.000000 | True |
| icu_exp_flg | 0 | 0.000000 | 0.095721 | 0.000000 | 0.294208 | 0.000000 | 1.000000 | True |
| day_28_flg | 0 | 0.000000 | 0.159347 | 0.000000 | 0.365999 | 0.000000 | 1.000000 | True |
| mort_day_censored | 0 | 0.000000 | 614.329825 | 731.000000 | 402.996046 | 0.000000 | 3094.080078 | True |
| censor_flg | 0 | 0.000000 | 0.720158 | 1.000000 | 0.448922 | 0.000000 | 1.000000 | False |
| sepsis_flg | 0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | False |
| chf_flg | 0 | 0.000000 | 0.119932 | 0.000000 | 0.324883 | 0.000000 | 1.000000 | True |
| afib_flg | 0 | 0.000000 | 0.116554 | 0.000000 | 0.320888 | 0.000000 | 1.000000 | True |
| renal_flg | 0 | 0.000000 | 0.033784 | 0.000000 | 0.180672 | 0.000000 | 1.000000 | True |
| liver_flg | 0 | 0.000000 | 0.055743 | 0.000000 | 0.229425 | 0.000000 | 1.000000 | True |
| copd_flg | 0 | 0.000000 | 0.088401 | 0.000000 | 0.283877 | 0.000000 | 1.000000 | True |
| cad_flg | 0 | 0.000000 | 0.069257 | 0.000000 | 0.253890 | 0.000000 | 1.000000 | True |
| stroke_flg | 0 | 0.000000 | 0.125000 | 0.000000 | 0.330719 | 0.000000 | 1.000000 | True |
| mal_flg | 0 | 0.000000 | 0.144144 | 0.000000 | 0.351236 | 0.000000 | 1.000000 | True |
| resp_flg | 0 | 0.000000 | 0.318131 | 0.000000 | 0.465751 | 0.000000 | 1.000000 | False |
| map_1st | 0 | 0.000000 | 88.246998 | 87.000000 | 17.590711 | 5.000000 | 195.000000 | True |
| hr_1st | 0 | 0.000000 | 87.914977 | 87.000000 | 18.753561 | 30.000000 | 158.000000 | True |
| temp_1st | 3 | 0.168919 | 97.792194 | 98.099998 | 4.539520 | 32.000000 | 104.800003 | True |
| spo2_1st | 0 | 0.000000 | 98.432995 | 100.000000 | 5.510842 | 4.000000 | 100.000000 | True |
| abg_count | 0 | 0.000000 | 5.984797 | 3.000000 | 8.681962 | 0.000000 | 115.000000 | True |
| wbc_first | 8 | 0.450450 | 12.320396 | 11.300000 | 6.597979 | 0.170000 | 109.800003 | True |
| hgb_first | 8 | 0.450450 | 12.551584 | 12.700000 | 2.200953 | 2.000000 | 19.000000 | True |
| platelet_first | 8 | 0.450450 | 246.083145 | 239.000000 | 99.837223 | 7.000000 | 988.000000 | True |
| sodium_first | 5 | 0.281532 | 139.559006 | 140.000000 | 4.724875 | 105.000000 | 165.000000 | True |
| potassium_first | 5 | 0.281532 | 4.107623 | 4.000000 | 0.794499 | 1.900000 | 9.800000 | True |
| tco2_first | 5 | 0.281532 | 24.416657 | 24.000000 | 4.990763 | 2.000000 | 62.000000 | True |
| chloride_first | 5 | 0.281532 | 103.839074 | 104.000000 | 5.732664 | 78.000000 | 133.000000 | True |
| bun_first | 5 | 0.281532 | 19.277809 | 15.000000 | 14.362833 | 2.000000 | 139.000000 | True |
| creatinine_first | 6 | 0.337838 | 1.095706 | 0.900000 | 1.083171 | 0.000000 | 18.299999 | True |
| po2_first | 186 | 10.472973 | 227.623270 | 195.000000 | 144.817841 | 22.000000 | 634.000000 | False |
| pco2_first | 186 | 10.472973 | 43.413836 | 41.000000 | 13.976388 | 8.000000 | 158.000000 | True |
| iv_day_1 | 143 | 8.051802 | 1622.907946 | 1081.529175 | 1676.615567 | 0.000000 | 13910.000000 | True |
All properties will be added to the respective layers. Categorical features can be found in the obs layer, while numerical features are in the var layer of the AnnData object. When inspecting both layers, we see that our QC properties were added for each feature if possible.
adata.obs.head(4)
| service_unit | day_icu_intime | missing_values_abs | missing_values_pct | |
|---|---|---|---|---|
| 0 | SICU | Friday | 0 | 0.000000 |
| 1 | MICU | Saturday | 12 | 22.222222 |
| 2 | MICU | Friday | 0 | 0.000000 |
| 3 | SICU | Saturday | 3 | 5.555556 |
adata.var.tail(4)
| feature_type | unencoded_var_names | encoding_mode | missing_values_abs | missing_values_pct | mean | median | standard_deviation | min | max | iqr_outliers | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| creatinine_first | numeric | creatinine_first | NaN | 6 | 0.337838 | 1.095706 | 0.900000 | 1.083171 | 0.0 | 18.299999 | True |
| po2_first | numeric | po2_first | NaN | 186 | 10.472973 | 227.623270 | 195.000000 | 144.817841 | 22.0 | 634.000000 | False |
| pco2_first | numeric | pco2_first | NaN | 186 | 10.472973 | 43.413836 | 41.000000 | 13.976388 | 8.0 | 158.000000 | True |
| iv_day_1 | numeric | iv_day_1 | NaN | 143 | 8.051802 | 1622.907946 | 1081.529175 | 1676.615567 | 0.0 | 13910.000000 | True |
We can also represent the missing values in a histogram for both obs and var features.
axd = plt.figure(constrained_layout=True, figsize=(8, 4), dpi=100).subplot_mosaic(
"""
AB
"""
)
sns.histplot(
adata.obs["missing_values_pct"], ax=axd["A"], bins=30, color="#54C285"
).set(title="pct of missing values: obs")
sns.histplot(
adata.var["missing_values_pct"], ax=axd["B"], bins=30, color="#1FA6C9"
).set(title="pct of missing values: var")
[Text(0.5, 1.0, 'pct of missing values: var')]
We can also check which features have the highest percentage of missing values in both obs and vars.
adata.obs.loc[
adata.obs["missing_values_pct"] == adata.obs["missing_values_pct"].max(), :
]
| service_unit | day_icu_intime | missing_values_abs | missing_values_pct | |
|---|---|---|---|---|
| 1732 | SICU | Thursday | 14 | 25.925926 |
| 1751 | MICU | Tuesday | 14 | 25.925926 |
adata.var.loc[
adata.var["missing_values_pct"] == adata.var["missing_values_pct"].max(), :
]
| feature_type | unencoded_var_names | encoding_mode | missing_values_abs | missing_values_pct | mean | median | standard_deviation | min | max | iqr_outliers | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| bmi | numeric | bmi | NaN | 466 | 26.238739 | 27.827316 | 26.324846 | 8.20694 | 12.784877 | 98.797134 | True |
Overall, the percentage of missing values in all features is rather low, however, still some features are not complete.
Features with missing values can introduce a bias in the data, making the processing and analysis challenging. To prevent loss of information due to dropping of multiple features, we can fill up the missing values by performing an imputation. Here, we infer the missing values based on the exisitng part of the data.
To perform this efficiently, we suggest to drop features if the percentage of missing values is very high (>60%). In our data, there is no need to drop any feature, since none exceeds more than 27% missing values (BMI, vars).
Missing data imputation#
ehrapy offers many options to impute missing values in an AnnData object.
Here, we use KNN imputation with 5 neighbors (n_neighbors=5, the default value). The KNN algorithm uses proximity to predict the missing values of a feature by finding the k closest neighbors to the missing value and then imputing the missing value based on the non-missing values in the neighborhood.
ehrapy offers two backends for the nearest neighbors search; scikit-learn and faiss. While faiss is faster for large datasets, scikit-learn is robustly reproducible across different machines.
ep.pp.knn_impute(adata, backend="scikit-learn", n_neighbours=5)
After recalcuating the QC metrices, we can check again the percentage of missing values.
ep.pp.qc_metrics(adata)
Show code cell output
( missing_values_abs missing_values_pct
0 0 0.0
1 0 0.0
2 0 0.0
3 0 0.0
4 0 0.0
... ... ...
1771 0 0.0
1772 0 0.0
1773 0 0.0
1774 0 0.0
1775 0 0.0
[1776 rows x 2 columns],
missing_values_abs missing_values_pct \
ehrapycat_service_unit_FICU 0 0.0
ehrapycat_service_unit_MICU 0 0.0
ehrapycat_service_unit_SICU 0 0.0
ehrapycat_day_icu_intime_Friday 0 0.0
ehrapycat_day_icu_intime_Monday 0 0.0
ehrapycat_day_icu_intime_Saturday 0 0.0
ehrapycat_day_icu_intime_Sunday 0 0.0
ehrapycat_day_icu_intime_Thursday 0 0.0
ehrapycat_day_icu_intime_Tuesday 0 0.0
ehrapycat_day_icu_intime_Wednesday 0 0.0
aline_flg 0 0.0
icu_los_day 0 0.0
hospital_los_day 0 0.0
age 0 0.0
gender_num 0 0.0
weight_first 0 0.0
bmi 0 0.0
sapsi_first 0 0.0
sofa_first 0 0.0
service_num 0 0.0
day_icu_intime_num 0 0.0
hour_icu_intime 0 0.0
hosp_exp_flg 0 0.0
icu_exp_flg 0 0.0
day_28_flg 0 0.0
mort_day_censored 0 0.0
censor_flg 0 0.0
sepsis_flg 0 0.0
chf_flg 0 0.0
afib_flg 0 0.0
renal_flg 0 0.0
liver_flg 0 0.0
copd_flg 0 0.0
cad_flg 0 0.0
stroke_flg 0 0.0
mal_flg 0 0.0
resp_flg 0 0.0
map_1st 0 0.0
hr_1st 0 0.0
temp_1st 0 0.0
spo2_1st 0 0.0
abg_count 0 0.0
wbc_first 0 0.0
hgb_first 0 0.0
platelet_first 0 0.0
sodium_first 0 0.0
potassium_first 0 0.0
tco2_first 0 0.0
chloride_first 0 0.0
bun_first 0 0.0
creatinine_first 0 0.0
po2_first 0 0.0
pco2_first 0 0.0
iv_day_1 0 0.0
mean median \
ehrapycat_service_unit_FICU NaN NaN
ehrapycat_service_unit_MICU NaN NaN
ehrapycat_service_unit_SICU NaN NaN
ehrapycat_day_icu_intime_Friday NaN NaN
ehrapycat_day_icu_intime_Monday NaN NaN
ehrapycat_day_icu_intime_Saturday NaN NaN
ehrapycat_day_icu_intime_Sunday NaN NaN
ehrapycat_day_icu_intime_Thursday NaN NaN
ehrapycat_day_icu_intime_Tuesday NaN NaN
ehrapycat_day_icu_intime_Wednesday NaN NaN
aline_flg 0.554054 1.000000
icu_los_day 3.346498 2.185000
hospital_los_day 8.110923 6.000000
age 54.379660 53.678585
gender_num 0.577590 1.000000
weight_first 79.983435 77.070000
bmi 27.690542 26.635570
sapsi_first 14.100113 14.000000
sofa_first 5.819595 6.000000
service_num 0.552928 1.000000
day_icu_intime_num 4.054054 4.000000
hour_icu_intime 10.585586 9.000000
hosp_exp_flg 0.137387 0.000000
icu_exp_flg 0.095721 0.000000
day_28_flg 0.159347 0.000000
mort_day_censored 614.329825 731.000000
censor_flg 0.720158 1.000000
sepsis_flg 0.000000 0.000000
chf_flg 0.119932 0.000000
afib_flg 0.116554 0.000000
renal_flg 0.033784 0.000000
liver_flg 0.055743 0.000000
copd_flg 0.088401 0.000000
cad_flg 0.069257 0.000000
stroke_flg 0.125000 0.000000
mal_flg 0.144144 0.000000
resp_flg 0.318131 0.000000
map_1st 88.246998 87.000000
hr_1st 87.914977 87.000000
temp_1st 97.793153 98.099998
spo2_1st 98.432995 100.000000
abg_count 5.984797 3.000000
wbc_first 12.318232 11.300000
hgb_first 12.554178 12.700000
platelet_first 246.088964 239.000000
sodium_first 139.561374 140.000000
potassium_first 4.107511 4.000000
tco2_first 24.413345 24.000000
chloride_first 103.842905 104.000000
bun_first 19.272072 15.000000
creatinine_first 1.095608 0.900000
po2_first 230.324887 204.000000
pco2_first 43.393581 41.000000
iv_day_1 1613.196694 1150.000000
standard_deviation min \
ehrapycat_service_unit_FICU NaN NaN
ehrapycat_service_unit_MICU NaN NaN
ehrapycat_service_unit_SICU NaN NaN
ehrapycat_day_icu_intime_Friday NaN NaN
ehrapycat_day_icu_intime_Monday NaN NaN
ehrapycat_day_icu_intime_Saturday NaN NaN
ehrapycat_day_icu_intime_Sunday NaN NaN
ehrapycat_day_icu_intime_Thursday NaN NaN
ehrapycat_day_icu_intime_Tuesday NaN NaN
ehrapycat_day_icu_intime_Wednesday NaN NaN
aline_flg 0.497070 0.000000
icu_los_day 3.355316 0.500000
hospital_los_day 8.154862 1.000000
age 21.056923 15.180230
gender_num 0.493852 0.000000
weight_first 21.885150 30.000000
bmi 7.254844 12.784877
sapsi_first 4.140549 3.000000
sofa_first 2.332028 0.000000
service_num 0.497191 0.000000
day_icu_intime_num 1.993911 1.000000
hour_icu_intime 7.922733 0.000000
hosp_exp_flg 0.344256 0.000000
icu_exp_flg 0.294208 0.000000
day_28_flg 0.365999 0.000000
mort_day_censored 402.996046 0.000000
censor_flg 0.448922 0.000000
sepsis_flg 0.000000 0.000000
chf_flg 0.324883 0.000000
afib_flg 0.320888 0.000000
renal_flg 0.180672 0.000000
liver_flg 0.229425 0.000000
copd_flg 0.283877 0.000000
cad_flg 0.253890 0.000000
stroke_flg 0.330719 0.000000
mal_flg 0.351236 0.000000
resp_flg 0.465751 0.000000
map_1st 17.590711 5.000000
hr_1st 18.753561 30.000000
temp_1st 4.535759 32.000000
spo2_1st 5.510842 4.000000
abg_count 8.681962 0.000000
wbc_first 6.585644 0.170000
hgb_first 2.197587 2.000000
platelet_first 99.646368 7.000000
sodium_first 4.719064 105.000000
potassium_first 0.793555 1.900000
tco2_first 4.984581 2.000000
chloride_first 5.725436 78.000000
bun_first 14.344715 2.000000
creatinine_first 1.081579 0.000000
po2_first 139.134923 22.000000
pco2_first 13.331953 8.000000
iv_day_1 1619.062172 0.000000
max iqr_outliers
ehrapycat_service_unit_FICU NaN True
ehrapycat_service_unit_MICU NaN True
ehrapycat_service_unit_SICU NaN True
ehrapycat_day_icu_intime_Friday NaN True
ehrapycat_day_icu_intime_Monday NaN True
ehrapycat_day_icu_intime_Saturday NaN True
ehrapycat_day_icu_intime_Sunday NaN True
ehrapycat_day_icu_intime_Thursday NaN True
ehrapycat_day_icu_intime_Tuesday NaN True
ehrapycat_day_icu_intime_Wednesday NaN True
aline_flg 1.000000 False
icu_los_day 28.240000 True
hospital_los_day 112.000000 True
age 99.110947 False
gender_num 1.000000 False
weight_first 257.600006 True
bmi 98.797134 True
sapsi_first 32.000000 True
sofa_first 17.000000 True
service_num 1.000000 False
day_icu_intime_num 7.000000 False
hour_icu_intime 23.000000 False
hosp_exp_flg 1.000000 True
icu_exp_flg 1.000000 True
day_28_flg 1.000000 True
mort_day_censored 3094.080078 True
censor_flg 1.000000 False
sepsis_flg 0.000000 False
chf_flg 1.000000 True
afib_flg 1.000000 True
renal_flg 1.000000 True
liver_flg 1.000000 True
copd_flg 1.000000 True
cad_flg 1.000000 True
stroke_flg 1.000000 True
mal_flg 1.000000 True
resp_flg 1.000000 False
map_1st 195.000000 True
hr_1st 158.000000 True
temp_1st 104.800003 True
spo2_1st 100.000000 True
abg_count 115.000000 True
wbc_first 109.800003 True
hgb_first 19.000000 True
platelet_first 988.000000 True
sodium_first 165.000000 True
potassium_first 9.800000 True
tco2_first 62.000000 True
chloride_first 133.000000 True
bun_first 139.000000 True
creatinine_first 18.299999 True
po2_first 634.000000 True
pco2_first 158.000000 True
iv_day_1 13910.000000 True )
axd = plt.figure(constrained_layout=True, figsize=(8, 3), dpi=100).subplot_mosaic(
"""
AB
"""
)
sns.histplot(adata.obs["missing_values_pct"], ax=axd["A"], bins=5, color="#54C285").set(
title="pct of missing values: obs", xlim=(0, 30)
)
sns.histplot(adata.var["missing_values_pct"], ax=axd["B"], bins=5, color="#1FA6C9").set(
title="pct of missing values: var", xlim=(0, 30)
)
[Text(0.5, 1.0, 'pct of missing values: var'), (0.0, 30.0)]
Data distribution#
Depending on the measurement and the unit of a measurement the value ranges of features may be huge. Clusterings and differential comparisons especially may be greatly influenced by exceptionally big values.
axd = plt.figure(constrained_layout=True, figsize=(8, 3), dpi=100).subplot_mosaic(
"""
AB
"""
)
sns.histplot(adata.var["min"], ax=axd["A"], bins=30, color="#54C285").set(
title="minimum values"
)
sns.histplot(adata.var["max"], ax=axd["B"], bins=30, color="#1FA6C9").set(
title="maximum values"
)
[Text(0.5, 1.0, 'maximum values')]
Moreover, features which have a very high coefficient of variation can strongly influence dimensionality reduction. However, since the coefficient of variation performs weak with features that have small means, we only select those which have no small mean.
adata.var["coefficient.variation"] = (
adata.var["standard_deviation"] / adata.var["mean"]
) * 100
adata.var.loc[(adata.var["coefficient.variation"] > 50) & (adata.var["mean"] > 50),]
| feature_type | unencoded_var_names | encoding_mode | missing_values_abs | missing_values_pct | mean | median | standard_deviation | min | max | iqr_outliers | coefficient.variation | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mort_day_censored | numeric | mort_day_censored | NaN | 0 | 0.0 | 614.329825 | 731.0 | 402.996046 | 0.0 | 3094.080078 | True | 65.599297 |
| po2_first | numeric | po2_first | NaN | 0 | 0.0 | 230.324887 | 204.0 | 139.134923 | 22.0 | 634.000000 | True | 60.408115 |
| iv_day_1 | numeric | iv_day_1 | NaN | 0 | 0.0 | 1613.196694 | 1150.0 | 1619.062172 | 0.0 | 13910.000000 | True | 100.363593 |
The standard deviations and coefficients of variation of the features iv_day_1 (input fluids by IV on day 1 in mL) and po2_first (first PaO_2 in mmHg) are very high with strong spread between minimum and maximum values. These features require normalization.
Normalization#
ehrapy offers several options to normalize data. While it is possible to normalize all numerical values at once with the same normalization function, normalizing only the features with high spread, here iv_day_1 and po2_first, can be sufficient. Log normalization with an offset of 1 to add pseudocounts seems appropriate.
Note: When features with negative values should be normalized you have to use the pp.offset_negative_values() function prior normalization.
ep.pp.log_norm(adata, vars=["iv_day_1", "po2_first"], offset=1)
AnnData object with n_obs × n_vars = 1776 × 54
obs: 'service_unit', 'day_icu_intime', 'missing_values_abs', 'missing_values_pct'
var: 'feature_type', 'unencoded_var_names', 'encoding_mode', 'missing_values_abs', 'missing_values_pct', 'mean', 'median', 'standard_deviation', 'min', 'max', 'iqr_outliers', 'coefficient.variation'
uns: 'service_unit_colors', 'normalization'
layers: 'original', 'raw_norm'
after normalization we can calculate the QC metrices again and check the distribution.
ep.pp.qc_metrics(adata)
Show code cell output
( missing_values_abs missing_values_pct
0 0 0.0
1 0 0.0
2 0 0.0
3 0 0.0
4 0 0.0
... ... ...
1771 0 0.0
1772 0 0.0
1773 0 0.0
1774 0 0.0
1775 0 0.0
[1776 rows x 2 columns],
missing_values_abs missing_values_pct \
ehrapycat_service_unit_FICU 0 0.0
ehrapycat_service_unit_MICU 0 0.0
ehrapycat_service_unit_SICU 0 0.0
ehrapycat_day_icu_intime_Friday 0 0.0
ehrapycat_day_icu_intime_Monday 0 0.0
ehrapycat_day_icu_intime_Saturday 0 0.0
ehrapycat_day_icu_intime_Sunday 0 0.0
ehrapycat_day_icu_intime_Thursday 0 0.0
ehrapycat_day_icu_intime_Tuesday 0 0.0
ehrapycat_day_icu_intime_Wednesday 0 0.0
aline_flg 0 0.0
icu_los_day 0 0.0
hospital_los_day 0 0.0
age 0 0.0
gender_num 0 0.0
weight_first 0 0.0
bmi 0 0.0
sapsi_first 0 0.0
sofa_first 0 0.0
service_num 0 0.0
day_icu_intime_num 0 0.0
hour_icu_intime 0 0.0
hosp_exp_flg 0 0.0
icu_exp_flg 0 0.0
day_28_flg 0 0.0
mort_day_censored 0 0.0
censor_flg 0 0.0
sepsis_flg 0 0.0
chf_flg 0 0.0
afib_flg 0 0.0
renal_flg 0 0.0
liver_flg 0 0.0
copd_flg 0 0.0
cad_flg 0 0.0
stroke_flg 0 0.0
mal_flg 0 0.0
resp_flg 0 0.0
map_1st 0 0.0
hr_1st 0 0.0
temp_1st 0 0.0
spo2_1st 0 0.0
abg_count 0 0.0
wbc_first 0 0.0
hgb_first 0 0.0
platelet_first 0 0.0
sodium_first 0 0.0
potassium_first 0 0.0
tco2_first 0 0.0
chloride_first 0 0.0
bun_first 0 0.0
creatinine_first 0 0.0
po2_first 0 0.0
pco2_first 0 0.0
iv_day_1 0 0.0
mean median \
ehrapycat_service_unit_FICU NaN NaN
ehrapycat_service_unit_MICU NaN NaN
ehrapycat_service_unit_SICU NaN NaN
ehrapycat_day_icu_intime_Friday NaN NaN
ehrapycat_day_icu_intime_Monday NaN NaN
ehrapycat_day_icu_intime_Saturday NaN NaN
ehrapycat_day_icu_intime_Sunday NaN NaN
ehrapycat_day_icu_intime_Thursday NaN NaN
ehrapycat_day_icu_intime_Tuesday NaN NaN
ehrapycat_day_icu_intime_Wednesday NaN NaN
aline_flg 0.554054 1.000000
icu_los_day 3.346498 2.185000
hospital_los_day 8.110923 6.000000
age 54.379660 53.678585
gender_num 0.577590 1.000000
weight_first 79.983435 77.070000
bmi 27.690542 26.635570
sapsi_first 14.100113 14.000000
sofa_first 5.819595 6.000000
service_num 0.552928 1.000000
day_icu_intime_num 4.054054 4.000000
hour_icu_intime 10.585586 9.000000
hosp_exp_flg 0.137387 0.000000
icu_exp_flg 0.095721 0.000000
day_28_flg 0.159347 0.000000
mort_day_censored 614.329825 731.000000
censor_flg 0.720158 1.000000
sepsis_flg 0.000000 0.000000
chf_flg 0.119932 0.000000
afib_flg 0.116554 0.000000
renal_flg 0.033784 0.000000
liver_flg 0.055743 0.000000
copd_flg 0.088401 0.000000
cad_flg 0.069257 0.000000
stroke_flg 0.125000 0.000000
mal_flg 0.144144 0.000000
resp_flg 0.318131 0.000000
map_1st 88.246998 87.000000
hr_1st 87.914977 87.000000
temp_1st 97.793153 98.099998
spo2_1st 98.432995 100.000000
abg_count 5.984797 3.000000
wbc_first 12.318232 11.300000
hgb_first 12.554178 12.700000
platelet_first 246.088964 239.000000
sodium_first 139.561374 140.000000
potassium_first 4.107511 4.000000
tco2_first 24.413345 24.000000
chloride_first 103.842905 104.000000
bun_first 19.272072 15.000000
creatinine_first 1.095608 0.900000
po2_first 5.236330 5.323010
pco2_first 43.393581 41.000000
iv_day_1 6.566182 7.048387
standard_deviation min \
ehrapycat_service_unit_FICU NaN NaN
ehrapycat_service_unit_MICU NaN NaN
ehrapycat_service_unit_SICU NaN NaN
ehrapycat_day_icu_intime_Friday NaN NaN
ehrapycat_day_icu_intime_Monday NaN NaN
ehrapycat_day_icu_intime_Saturday NaN NaN
ehrapycat_day_icu_intime_Sunday NaN NaN
ehrapycat_day_icu_intime_Thursday NaN NaN
ehrapycat_day_icu_intime_Tuesday NaN NaN
ehrapycat_day_icu_intime_Wednesday NaN NaN
aline_flg 0.497070 0.000000
icu_los_day 3.355316 0.500000
hospital_los_day 8.154862 1.000000
age 21.056923 15.180230
gender_num 0.493852 0.000000
weight_first 21.885150 30.000000
bmi 7.254844 12.784877
sapsi_first 4.140549 3.000000
sofa_first 2.332028 0.000000
service_num 0.497191 0.000000
day_icu_intime_num 1.993911 1.000000
hour_icu_intime 7.922733 0.000000
hosp_exp_flg 0.344256 0.000000
icu_exp_flg 0.294208 0.000000
day_28_flg 0.365999 0.000000
mort_day_censored 402.996046 0.000000
censor_flg 0.448922 0.000000
sepsis_flg 0.000000 0.000000
chf_flg 0.324883 0.000000
afib_flg 0.320888 0.000000
renal_flg 0.180672 0.000000
liver_flg 0.229425 0.000000
copd_flg 0.283877 0.000000
cad_flg 0.253890 0.000000
stroke_flg 0.330719 0.000000
mal_flg 0.351236 0.000000
resp_flg 0.465751 0.000000
map_1st 17.590711 5.000000
hr_1st 18.753561 30.000000
temp_1st 4.535759 32.000000
spo2_1st 5.510842 4.000000
abg_count 8.681962 0.000000
wbc_first 6.585644 0.170000
hgb_first 2.197587 2.000000
platelet_first 99.646368 7.000000
sodium_first 4.719064 105.000000
potassium_first 0.793555 1.900000
tco2_first 4.984581 2.000000
chloride_first 5.725436 78.000000
bun_first 14.344715 2.000000
creatinine_first 1.081579 0.000000
po2_first 0.684555 3.135494
pco2_first 13.331953 8.000000
iv_day_1 1.833505 0.000000
max iqr_outliers
ehrapycat_service_unit_FICU NaN True
ehrapycat_service_unit_MICU NaN True
ehrapycat_service_unit_SICU NaN True
ehrapycat_day_icu_intime_Friday NaN True
ehrapycat_day_icu_intime_Monday NaN True
ehrapycat_day_icu_intime_Saturday NaN True
ehrapycat_day_icu_intime_Sunday NaN True
ehrapycat_day_icu_intime_Thursday NaN True
ehrapycat_day_icu_intime_Tuesday NaN True
ehrapycat_day_icu_intime_Wednesday NaN True
aline_flg 1.000000 False
icu_los_day 28.240000 True
hospital_los_day 112.000000 True
age 99.110947 False
gender_num 1.000000 False
weight_first 257.600006 True
bmi 98.797134 True
sapsi_first 32.000000 True
sofa_first 17.000000 True
service_num 1.000000 False
day_icu_intime_num 7.000000 False
hour_icu_intime 23.000000 False
hosp_exp_flg 1.000000 True
icu_exp_flg 1.000000 True
day_28_flg 1.000000 True
mort_day_censored 3094.080078 True
censor_flg 1.000000 False
sepsis_flg 0.000000 False
chf_flg 1.000000 True
afib_flg 1.000000 True
renal_flg 1.000000 True
liver_flg 1.000000 True
copd_flg 1.000000 True
cad_flg 1.000000 True
stroke_flg 1.000000 True
mal_flg 1.000000 True
resp_flg 1.000000 False
map_1st 195.000000 True
hr_1st 158.000000 True
temp_1st 104.800003 True
spo2_1st 100.000000 True
abg_count 115.000000 True
wbc_first 109.800003 True
hgb_first 19.000000 True
platelet_first 988.000000 True
sodium_first 165.000000 True
potassium_first 9.800000 True
tco2_first 62.000000 True
chloride_first 133.000000 True
bun_first 139.000000 True
creatinine_first 18.299999 True
po2_first 6.453625 True
pco2_first 158.000000 True
iv_day_1 9.540435 True )
adata.var["coefficient.variation"] = (
adata.var["standard_deviation"] / adata.var["mean"]
) * 100
adata.var.loc[(adata.var["coefficient.variation"] > 50) & (adata.var["mean"] > 50),]
| feature_type | unencoded_var_names | encoding_mode | missing_values_abs | missing_values_pct | mean | median | standard_deviation | min | max | iqr_outliers | coefficient.variation | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mort_day_censored | numeric | mort_day_censored | NaN | 0 | 0.0 | 614.329825 | 731.0 | 402.996046 | 0.0 | 3094.080078 | True | 65.599297 |
The strong spread of iv_day_1 and po2_first was succesfully removed. Now that we normalized the influence of these features, we can continue with dimensionality reduction.
Dimensionality reduction reduces the number of features (dimensions) by projecting the data to a lower dimensional latent space retaining as much information as possible. This is very useful for high dimensional data, since it reduces complexity and facilitates visualization.
Dimensionality reduction#
Principle Component Analysis (PCA)#
As a next step, we reduce the dimensionality of the dataset with principal component analysis (PCA).
We can also visualize the principal components with ehrapy using the components argument.
scanpy, which ehrapy uses under the hood, provides many options for computing a PCA. The option randomized with a random state is particularly reproducible across different machines. However, different BLAS/LAPACK backends used on different machines have slight differences in SVD. For exact reproducibility, working on in a containerized environment is essential.
ep.pp.pca(adata, svd_solver="randomized", random_state=42)
ep.pl.pca(adata, color="service_unit", components=["1,2", "3,4"])
ep.pl.pca(adata, color="service_unit", components=["5,6", "7,8"])
To inspect certain PCs further, we can inspect the PC dimensionality loadings which highlight the features that contribute strongest to the selected PC.
ep.pl.pca_loadings(adata, components="1, 2")
ep.pl.pca_loadings(adata, components="3, 4")
Uniform Manifold Approximation and Projection (UMAP)#
The reduced representation can then be used as input for the neighbors graph calculation which serves as the input for advanced embeddings and visualizations like Uniform Manifold Approximation and Projection (UMAP)
ehrapy provides multiple implementations for neighborhood search. By setting transformer="sklearn", brute-force, but robustly reproducible implementation across machines is available. Faster options for the transformers argument are available, too.
ep.pp.neighbors(adata, transformer="sklearn", n_pcs=10)
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
ep.tl.umap(adata)
Checking for Batch effects#
Before exploring the data further, we need to see if we have a batch effect. A batch effect can e.g. arise from different collection units or collection days. To check if our data contains a batch for those feautures, we visualize the service_unit and the day_icu_intime.
ep.settings.set_figure_params(figsize=(6, 5))
ep.pl.umap(
adata,
color=[
"service_unit",
"day_icu_intime",
],
wspace=0.5,
size=20,
title=["Service unit", "Day of ICU admission"],
)
The embeddings suggest that there’s no strong effect by the aforementioned potential confounders.
Selected features on UMAP#
Now we can also highlight other relevant features on the UMAP. Interesting features could be demographics, hospital statistics and lab parameters.
Demographics#
ep.pl.umap(
adata,
color=["gender_num", "age"],
wspace=0.5,
size=20,
title=["sex (1 = male; 0=female)", "age"],
)
Hospital statistics#
ep.pl.umap(
adata,
color=["icu_los_day", "hosp_exp_flg"],
wspace=0.5,
size=20,
cmap=ep.pl.Colormaps.grey_red.value,
title=["length of stay in ICU (days)", "death in hospital (1 = yes, 0 = no)"],
)
Comorbidities#
ep.pl.umap(
adata,
color=["liver_flg", "stroke_flg"],
cmap=ep.pl.Colormaps.grey_red.value,
title=["Liver disease", "Stroke"],
ncols=2,
size=20,
)
Lab parameters#
ep.pl.umap(
adata,
color=["hr_1st", "platelet_first", "po2_first", "pco2_first"],
wspace=0.5,
ncols=2,
size=20,
title=["Heart Rate", "Platelets (K/u)", "PaO2 (mmHg)", "PaCO2 (mmHg)"],
)
Cluster analysis#
To make more sense of the embedding it is often times useful to determine clusters through e.g. community detection as implemented in the Leiden algorithm. Moreover, clustering allows for unbiased detection of features that are changed between clusters and therefore intersting for us.
Cluster identification#
The implementation in ehrapy allows for the setting of a resolution which determines the number of found clusters. It is often times useful to play around with the parameter.
ep.tl.leiden(adata, resolution=0.3, key_added="leiden_0_3")
The leiden algorithm added a key to obs (leiden_0_3) that stores the clusters. These can subsequently be visualized in the UMAP embedding.
adata.obs.head(4)
| service_unit | day_icu_intime | missing_values_abs | missing_values_pct | leiden_0_3 | |
|---|---|---|---|---|---|
| 0 | SICU | Friday | 0 | 0.0 | 0 |
| 1 | MICU | Saturday | 0 | 0.0 | 1 |
| 2 | MICU | Friday | 0 | 0.0 | 2 |
| 3 | SICU | Saturday | 0 | 0.0 | 0 |
ep.pl.umap(adata, color=["leiden_0_3"], title="Leiden 0.3", size=20)
Next, we can explore certain features which are special for certrain clusters and could therefore be used for annotation.
Cluster features#
ep.tl.rank_features_groups(adata, groupby="leiden_0_3")
ep.settings.set_figure_params(figsize=(4, 4), dpi=100)
ep.pl.rank_features_groups(adata, key="rank_features_groups", ncols=2)
We can also get the top features per cluster as a DataFrame.
df = ep.get.rank_features_groups_df(adata, group=["0", "1", "2", "3", "4", "5"])
df = df.loc[(df["logfoldchanges"] > 0) & (df["pvals_adj"] < 0.05),]
E.g. we can check the top marker of cluster 2.
df.loc[df["group"] == "2",]
| group | names | scores | logfoldchanges | pvals | pvals_adj | |
|---|---|---|---|---|---|---|
| 108 | 2 | platelet_first | 29.306358 | 169.371445 | 3.440125e-142 | 1.857667e-140 |
| 109 | 2 | censor_flg | 313.916103 | 1.000000 | 3.062986e-70 | 8.270063e-69 |
| 111 | 2 | day_28_flg | 204.287557 | 1.000000 | 2.422363e-46 | 3.270190e-45 |
| 112 | 2 | hosp_exp_flg | 172.626270 | 1.000000 | 1.975095e-39 | 2.133103e-38 |
| 113 | 2 | mort_day_censored | 12.352182 | inf | 3.438069e-33 | 3.094262e-32 |
| 115 | 2 | icu_exp_flg | 115.279432 | 1.000000 | 6.835411e-27 | 4.613902e-26 |
| 118 | 2 | hgb_first | 8.391182 | 1.317119 | 1.587689e-16 | 7.794108e-16 |
| 119 | 2 | wbc_first | 8.116502 | 3.879998 | 1.420268e-15 | 6.391207e-15 |
| 120 | 2 | afib_flg | 56.636729 | 1.000000 | 5.242345e-14 | 2.177590e-13 |
| 121 | 2 | liver_flg | 44.480150 | 1.000000 | 2.569492e-11 | 9.910898e-11 |
| 122 | 2 | chf_flg | 42.702421 | 1.000000 | 6.373338e-11 | 2.294402e-10 |
| 123 | 2 | copd_flg | 32.675664 | 1.000000 | 1.088921e-08 | 3.675108e-08 |
| 124 | 2 | mal_flg | 31.651167 | 1.000000 | 1.845037e-08 | 5.860706e-08 |
| 125 | 2 | cad_flg | 23.716682 | 1.000000 | 1.116109e-06 | 3.348327e-06 |
| 126 | 2 | resp_flg | 19.761231 | 1.000000 | 8.774456e-06 | 2.493793e-05 |
| 130 | 2 | renal_flg | 17.080558 | 1.000000 | 3.582720e-05 | 8.411603e-05 |
| 132 | 2 | stroke_flg | 15.606439 | 1.000000 | 7.798856e-05 | 1.684553e-04 |
| 133 | 2 | aline_flg | 6.102413 | 1.000000 | 1.349974e-02 | 2.803793e-02 |
| 134 | 2 | hour_icu_intime | 39.942517 | 1.000000 | 1.559863e-02 | 3.119726e-02 |
| 135 | 2 | sodium_first | 2.332679 | 0.826686 | 1.987340e-02 | 3.832726e-02 |
From this table we can also extract the top features in every cluster and highlight those either on the UMAP or as violins plots by cluster.
top_features = df.groupby("group").head(5)
top_features = pd.Series(top_features["names"].unique())
top_features
0 day_28_flg
1 censor_flg
2 hosp_exp_flg
3 icu_exp_flg
4 age
5 mort_day_censored
6 platelet_first
7 sofa_first
8 mal_flg
dtype: object
ep.settings.set_figure_params(figsize=(3.8, 2), dpi=100)
ep.pl.violin(adata, keys=["censor_flg", "mort_day_censored"], groupby="leiden_0_3")
ep.pl.violin(adata, keys=["platelet_first", "age"], groupby="leiden_0_3")
ep.pl.violin(adata, keys=["sapsi_first", "copd_flg"], groupby="leiden_0_3")
ep.pl.violin(adata, keys=["sofa_first", "liver_flg"], groupby="leiden_0_3")
Cluster annotation#
With the knowledge of the cluster features, together with the UMAP plots from above we can annotate the clusters.
adata.obs["annotation"] = "NA"
annotation = {
"0": "liver+/sofa+",
"1": "weight+",
"2": "age+/stroke+/deceased+",
"3": "platelet+",
"4": "age+/malignancy+/copd+/deceased+",
"5": "age+",
}
adata.obs["annotation"] = [
annotation[l] if l in annotation.keys() else l for l in adata.obs["leiden_0_3"]
]
ep.settings.set_figure_params(figsize=(4, 3), dpi=100)
ep.pl.umap(
adata,
color="annotation",
size=20,
palette={
"weight+": "#007742",
"platelet+": "#54C285",
"age+/stroke+/deceased+": "#087A96",
"age+/malignancy+/copd+/deceased+": "#1FA6C9",
"age+": "#F4CC47",
"liver+/sofa+": "#57C8B9",
"platelet+/heart_rate+": "#ABEC7D",
},
)
Additional downstream analysis#
After these basic ehrapy analysis steps, additional downstream analysis can be performed (see also other tutorials).
PAGA#
It might also be of interest to infer trajectories to learn about dynamic processes and stage transitions. ehrapy offers several trajectory inference algorithms for this purpose. One of those is partition-based graph abstraction (PAGA).
ep.tl.paga(adata, groups="leiden_0_3")
ep.pl.paga(
adata,
color=["leiden_0_3", "day_28_flg"],
cmap=ep.pl.Colormaps.grey_red.value,
title=["Leiden 0.3", "Died in less than 28 days"],
)
ep.tl.umap(adata, init_pos="paga")
ep.pl.umap(adata, color=["annotation"])
ep.tl.draw_graph(adata, init_pos="paga")
WARNING: Package 'fa2-modified' is not installed, falling back to layout 'fr'.To use the faster and better ForceAtlas2 layout, install package 'fa2-modified' (`pip install fa2-modified`).
ep.tl.draw_graph(adata, init_pos="paga")
ep.pl.draw_graph(adata, color=["leiden_0_3", "day_28_flg"], legend_loc="on data")
WARNING: Package 'fa2-modified' is not installed, falling back to layout 'fr'.To use the faster and better ForceAtlas2 layout, install package 'fa2-modified' (`pip install fa2-modified`).
Exporting results#
We save all of our computations and our final state into an .h5ad file.
It can then be read again using the :func:ehrapy.io.read function.
ep.io.write("mimic_2.h5ad", adata)
Conclusion#
The MIMIC-II IAC dataset comprises electronic health records (EHR) summarized in 46 features from 1776 individuals. This high dimensional data is not easy to interpret and many interesting and previously unknown features can be overseen when just focusing on selected well-defined features. To overcome this hurdle, we applied ehrapy on the MIMIC-II IAC dataset.
ehrapy is based on the AnnData data structure and scanpy pipeline to allow for efficient analysis. We used the build-in functions to preprocess the data, perform QC with imputation of missing data and reduce the dimensionality, resulting in PCA and UMAP embeddings. After performing all these steps, we explored the data by visualizing multiple features on the UMAP embedding, giving a first glance at the patient structure. To identify patient groups in and unbiased fashion, we clustered our data using the Leiden algorithm resulting in 7 different patient clusters. Calculation of cluster-specific features allowed us to annotate the clusters according to the most prominent markers. We saw a strong difference between patients that deceased, had higher age and severe comorbidities such as a stroke and COPD (clusters 2+3) and those that had milder features such as increased platelets and weight (clusters 0+1). Close to these two clusters were two additional clusters that harbored more severe features such as increased heart rate (cluster 5) and high SOFA score with liver disease (cluster 6), indicating potential patient trajectories. Cluster 4 clustered apart from all the others and consists of patients that deceased several months/years after leaving the ICU.
To explore the patient fate, survival and a case study in more detail, continue with our other tutorials or go back to our tutorial overview page.
References#
Raffa, J. (2016). Clinical data from the MIMIC-II database for a case study on indwelling arterial catheters (version 1.0). PhysioNet. https://doi.org/10.13026/C2NC7F.
Raffa J.D., Ghassemi M., Naumann T., Feng M., Hsu D. (2016) Data Analysis. In: Secondary Analysis of Electronic Health Records. Springer, Cham
Goldberger, A., Amaral, L., Glass, L., Hausdorff, J., Ivanov, P. C., Mark, R., … & Stanley, H. E. (2000). PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals. Circulation [Online]. 101 (23), pp. e215–e220.
McInnes et al., (2018). UMAP: Uniform Manifold Approximation and Projection. Journal of Open Source Software, 3(29), 861, https://doi.org/10.21105/joss.00861
Traag, V.A., Waltman, L. & van Eck, N.J. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep 9, 5233 (2019). https://doi.org/10.1038/s41598-019-41695-z
Wolf, F.A., Hamey, F.K., Plass, M. et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol 20, 59 (2019). https://doi.org/10.1186/s13059-019-1663-x
Package versions#
ep.print_versions()