Skip to content
changlabtwPublic

About

cross-ethnic insulin resistance prediction

Topics

Resources

Stars

0 stars

Watchers

0 watching

Forks

Repository files navigation

Cross-ethnic machine learning identifies insulin resistance in non-diabetic adults and its DNA methylation correlates

DOI License: MIT

Analysis code for a gradient boosting model that predicts insulin resistance in non-diabetic adults from measurements a routine blood panel already contains — no fasting insulin assay required. The model is trained on pooled US (NHANES) and Korean (KNHANES) cohorts, interpreted with SHAP, and applied to the Taiwan Biobank cohort, whose predicted groups are then compared by DNA methylation.

Insulin resistance is defined by HOMA-IR > 2.5, where HOMA-IR = fasting insulin (µU/mL) × fasting glucose (mg/dL) / 405. Participants with diabetes are excluded, so the task is detecting the pre-diabetic state in people not yet diagnosed.

Study design

Non-diabetic adults from NHANES and KNHANES, both with measured HOMA-IR, train and test the classifiers. Transfer between the two cohorts is assessed in both directions before they are pooled. The model trained on the pooled set is then applied to a Taiwan Biobank cohort that has no fasting-insulin measurement and therefore no observed label, and the subset of those participants with methylation array data is compared by predicted IR status.

%%{init: {'flowchart': {'nodeSpacing': 90, 'rankSpacing': 55, 'wrappingWidth': 340}}}%%
flowchart TB
    subgraph COHORTS[" "]
        direction LR
        NHANES["NHANES<br/>Non-diabetic adults<br/>n = 11,660"]
        KNHANES["KNHANES<br/>Non-diabetic adults<br/>n = 15,138"]
        NHANES <-.->|bidirectional transfer| KNHANES
    end
    POOLED["Pooled training set<br/>n = 26,798<br/>Gradient-boosting classifiers<br/>IR defined as HOMA-IR &gt; 2.5"]
    TWB["Taiwan Biobank<br/>n = 92,734, no fasting-insulin measurement<br/>IR status predicted by the 20-feature CatBoost model"]
    MET["Methylation subset<br/>n = 1,199 with EPIC array data<br/>Compared by predicted IR status"]

    COHORTS --> POOLED
    POOLED --> TWB
    TWB --> MET

    classDef cohort fill:#DEEBF7,stroke:#0072B2,stroke-width:1.5px,color:#14213D
    classDef pooled fill:#FDF1DC,stroke:#E69F00,stroke-width:1.5px,color:#14213D
    classDef external fill:#DCF1EA,stroke:#009E73,stroke-width:1.5px,color:#14213D
    classDef omics fill:#F7E6F0,stroke:#CC79A7,stroke-width:1.5px,color:#14213D

    class NHANES,KNHANES cohort
    class POOLED pooled
    class TWB external
    class MET omics
    style COHORTS fill:none,stroke:#999999,stroke-width:1px,stroke-dasharray: 5 4
Loading

The same figure is available as TikZ in output/study_design.tex for the manuscript.

Results

Cohort Raw records Analysed Insulin resistant
NHANES 1999–2012 (US) 71,916 11,660 5,183 (44.5%)
KNHANES 2019–2021 (Korea) 22,559 15,138 4,226 (27.9%)
Pooled training set — 26,798 9,409 (35.1%)
Taiwan Biobank (external) 196,963 92,734 19,596 (21.1%), predicted

Taiwan Biobank measures no fasting insulin, so it carries no observed label and no accuracy can be computed on it. It is used to ask whether the model's predictions behave like real labels.

Pooled model, 30% hold-out of 26,798 participants:

Model Features AUC Accuracy Sensitivity Specificity NPV PR-AUC
CatBoost 241 0.880 0.799 0.778 0.810 0.871 0.816
Soft voting (3 models) 241 0.880 0.799 0.782 0.809 0.873 0.817
CatBoost, top 20 by SHAP 20 0.880 0.799 0.785 0.807 0.874 0.815

Reducing 241 engineered features to the twenty most influential costs nothing measurable, which is the practical result: the deployable model needs twenty inputs.

SHAP summary plot of the twenty highest-ranked features

What the model relies on. Each point is one of the 8,040 test-set participants, positioned by how much that feature pushed their prediction and coloured by the feature's value. Liver enzymes rank among the strongest predictors, which is the finding that ties insulin resistance to fatty liver disease.

Uncertainty. Resampling the 8,040-participant test set 1,000 times with replacement puts a 95% interval around each figure above. For the deployed twenty-feature model: AUC 0.880 (0.872–0.887), sensitivity 0.785 (0.769–0.799), NPV 0.874 (0.864–0.882). Every model's AUC interval overlaps every other's, which is the quantitative form of the sentence above. The interval covers sampling variation in the test set alone — not a different split, a different hyperparameter search or a different cohort. Full table in output/final_model_ci.xlsx.

Calibration. The model ranks participants well and scores them too high. The calibration slope is 1.02, so the shape of the relationship is right, but the intercept is −0.574: mean predicted probability is 0.43 against an observed prevalence of 0.351, and all ten deciles fall below the diagonal (Brier score 0.139; see Calibration_CatBoost.png). This is the intended cost of the class weighting — positives are up-weighted by the negative-to-positive ratio precisely to buy sensitivity, and that inflates the probabilities. Two consequences are worth stating plainly. The sensitivity, specificity and NPV above are empirical counts at a 0.5 threshold and are unaffected by it. But the raw predicted probability should not be read as an absolute risk without recalibration. Taiwan Biobank carries no observed label, so no calibration can be computed there.

Calibration curves for the 241-feature and top-20 CatBoost models

Ranks well, scores too high. Each point is a decile of predicted probability plotted against how often those participants were actually insulin resistant; a perfectly calibrated model sits on the diagonal. Both models track it in shape but sit consistently below it.

Cross-cohort transfer. Each of the four classifiers is scored on the whole cohort it never saw, in both directions, so that the direction of transfer is not confounded with the choice of algorithm. Discrimination barely moves with either: AUC is 0.861–0.862 sending NHANES to KNHANES and 0.860–0.862 sending KNHANES to NHANES, against within-cohort ceilings of 0.868 and 0.878.

The operating point is what changes, and it inverts with direction. A NHANES-trained model reaching KNHANES is conservative — sensitivity around 0.62 against specificity around 0.88 — while a KNHANES-trained model reaching NHANES is the mirror image, sensitivity around 0.89 against specificity around 0.61. All four algorithms agree in both directions, so this is a property of the shift between the cohorts rather than of any one model, and it is the argument for pooling: neither single-population model carries its threshold across.

Differential methylation. The 1,199 Taiwan Biobank participants with array data were profiled on the Illumina Infinium MethylationEPIC BeadChip, not the earlier 450K array: its manifest carries 866,895 probes, of which the 863,904 cg sites are annotated on GRCh37/hg19 — the assembly the downstream gene-structure lookups use. Probes are restricted to those cg sites — dropping the 2,932 ch non-CpG and 59 rs control probes — then to the autosomes, so that a group difference cannot simply reflect the sex imbalance between the groups. A reading is kept only when its detection p-value is below 0.001, and a probe only when every participant has a usable reading for it, so that no test runs on a varying subset. That cascade is 866,895 → 863,904 → 844,316 → 332,284 probes tested.

Of those, 22 differ between the predicted groups at q < 0.05 and |log2FC| > 0.2, spanning 12 named genes including COL25A1 (cg22266749), PSMA6 (cg02987832) and ERV3-1 (cg06513015). The manifest annotates those probes with 13 symbols, but UBE2QP1 is an alias of UBE2Q2P1, the official symbol. Two caveats bound how far they should be read: the grouping variable is a model prediction rather than a measurement, and no cross-reactive or SNP-overlap probe filter is applied — the manifest's SNP_ID, SNP_DISTANCE and SNP_MINORALLELEFREQUENCY columns are carried but unused, and published cross-reactive probe lists are not consulted. These 22 probes are therefore hypothesis-generating candidates for follow-up, not a filtered final set.

Volcano plot of 332,284 CpG sites tested between the predicted groups

All 332,284 tested probes at once. Effect size runs left to right, confidence bottom to top, and the dashed lines mark the two thresholds a probe must pass. The 22 that pass both are the highlighted points.

Repository layout

configs/          Tuned hyperparameters, cohort file inventory, probe annotation, paths
data/             Only the directories are tracked; their contents never are
  raw/            Raw cohort files, placed here by you (see Data availability)
  processed/      Intermediate tables written by notebook 01
output/           Result tables and figures (committed)
models/           Final CatBoost models and manifest.csv
notebooks/        01–08, one per step of the analysis
src/
  data/           Cohort readers: nhanes.py, knhanes.py, twb.py, io.py
  features/       Interaction generation and feature selection
  models/         Training, evaluation, SHAP, regression
  methylation/    Extraction, differential testing, annotation, enrichment
  viz/            Figures and tables
tests/            Unit tests, and a check that this README matches the outputs
run_all.sh        Runs the notebooks in dependency order
requirements.txt  Pinned dependencies
environment.yml   The conda environment the results were produced in
CITATION.cff      Machine-readable citation metadata, read by Zenodo
README.md         This file
LICENSE           MIT, for the code only

Setup

Requires Python 3.12.8. The pinned versions in requirements.txt are part of the reproduction contract rather than a lower bound: the boosting libraries are seeded only by their own defaults, so the reported metrics are guaranteed reachable under this exact set and not necessarily under another.

python -m venv .venv && source .venv/bin/activate
pip install -r requirements.txt

The reported results were produced under conda, so environment.yml recreates that environment exactly if you prefer it:

conda env create -f environment.yml && conda activate crossir

Place the raw cohort files under data/raw/ following the layout in configs/cohorts.yaml. To read them from somewhere else, create configs/paths.local.yaml overriding only the keys you need:

raw_data_root: /path/to/your/data/raw

Reproducing the analysis

Run ./run_all.sh, which executes the notebooks in the order they depend on each other. Each reads from data/processed/ and writes to output/, so a run leaves every table and figure in the repository regenerated.

# Notebook Produces
01 01_data_preparation.ipynb Cleaned cohort tables, HOMA-IR label, 241 engineered features
02 02_descriptive_statistics.ipynb stats.xlsx and stats_original.xlsx, boxplots, correlation matrix, HOMA-IR distributions
03 03_homair_regression.ipynb Continuous HOMA-IR regression across six models
04 04_feature_ablation.ipynb 9 / 17 / 57 / 241 feature sets compared
05 05_cross_ethnic.ipynb Within-cohort hold-out, and every model scored on the other cohort
06 06_final_model_and_shap.ipynb Pooled model, SHAP (default and colour-blind palettes), top-20 retrain, bootstrap intervals, calibration, panel-labelled figures, models/
07 07_twb_validation.ipynb Taiwan Biobank scoring, stats_twb.xlsx and stats_twb_original.xlsx, TG/HDL-C comparison, Q-Q plot
08 08_differential_methylation.ipynb Volcano plot, met-point.xlsx, candidate.csv, enrichment

The orderings that matter are the data dependencies: 01 first, since everything reads the tables it writes, and then 06 → 07 → 08, because 07 scores Taiwan Biobank with the model 06 saves and 08 splits participants by the label 07 predicts. Notebooks 02–05 need only 01 and can run in any order.

Notebook 08 takes about twelve minutes, most of it 1.3 million statistical tests; everything else runs in seconds to a few minutes.

pytest tests/ runs without any cohort data. Alongside unit tests of the label definition, the feature generator and the multiple-testing correction, it re-derives every number quoted in this README from the committed files in output/ and models/, so the two cannot drift apart.

Figure colours. Five figures use colour to carry meaning in a way that fails for the common forms of colour vision deficiency, and each has a _color_blind.png sibling drawn from the Okabe–Ito palette: correlation_matrix (a red-to-blue diverging scale), volcano_plot_of_DNA_methylation_analysis (red points against a green threshold line), ORA (a categorical palette pairing red with green), and HOMA-IR and qqplot_NHANES+KNHANES_vs_TWB (red reference lines). SHAP_summary_plot_new has one for the same reason. The remaining figures already encode safely and have no sibling: ConfusionMatrix_CatBoost is a single-hue sequential map, ROC_PR_CatBoost is orange on blue, and the three boxplots — boxplot_of_variable_by_IR_and_Race, boxplot_of_variable_by_IR_and_Race(Add TWB) and boxplot_of_TG_HDL_C_by_IR_and_Datasets — separate IR− from IR+ with blue and orange. In every pair the data, axes and ordering are identical — only the palette differs.

Calibration_CatBoost.png is drawn from the same palette from the start — blue and orange, with the diagonal in neutral grey and the two models additionally separated by marker shape — so it needs no sibling either. The three figures shown above are the colour-blind-safe versions, for the same reason the siblings exist at all.

Panel labels. The two halves of the discrimination figure are written as separate files, so ROC_PR_CatBoost_a.png and ConfusionMatrix_CatBoost_b.png are copies carrying the panel letters (a) and (b) for use in a manuscript. The letter sits in a strip added above the plot rather than drawn over it, so the unlabelled originals and the labelled copies are pixel-identical everywhere the figure itself is drawn.

Data availability

Cohort Access
NHANES (1999–2012) Public. US CDC / National Center for Health Statistics: https://wwwn.cdc.gov/nchs/nhanes/
KNHANES (2019–2021) Public, subject to the agency's terms of use. Korea Disease Control and Prevention Agency: https://knhanes.kdca.go.kr/
Taiwan Biobank Restricted. Available to qualified researchers by application to Academia Sinica: https://www.twbiobank.org.tw/

This study's use of Taiwan Biobank data was approved under TWBR11308-07, and the research protocol under IRB approval NCCU-REC-201904-I017. No individual-level data from any cohort is distributed with this repository; data/raw/ and data/processed/ are excluded from version control.

The methylation analysis cannot be re-run from its raw input here. Notebook 08 begins from per-participant methylation tables rather than from the vendor array archive, which is restricted-access Taiwan Biobank data and is not distributed with this repository. The extraction step is kept in src/methylation/extract.py as the specification of how those tables were produced, for researchers with their own approved Taiwan Biobank methylation access.

Those tables hold already-normalised beta values: the archive supplies one *_nor.txt export per participant with AVG_Beta and Detection Pval columns, and this pipeline applies no normalisation, background correction or batch adjustment of its own. Which normalisation the data provider ran is therefore outside this repository — anyone reproducing or extending the analysis should take that detail from the Taiwan Biobank methylation documentation accompanying their own data release, since it is not recoverable from the exports themselves.

Citation

If you use this software or its results, please cite the accompanying paper and this repository. CITATION.cff carries the metadata in machine-readable form, which GitHub's Cite this repository button will render as BibTeX or APA.

Every version is archived on Zenodo under 10.5281/zenodo.23018637, which always resolves to the most recent one. To cite the exact code behind the paper, use the v1.0.0 DOI, 10.5281/zenodo.23018638.

License

MIT — see LICENSE. The cohort data is not covered by this licence and remains subject to each provider's own terms.

About

cross-ethnic insulin resistance prediction

Topics

Resources

Stars

0 stars

Watchers

0 watching

Forks

Releases

Packages

Contributors

Languages