Global Geo-Pharmacogenomics: Environmental Mutational Signatures Drive Population-Level Heterogeneity in Anticancer Drug Response

1. Introduction

The global burden of cancer remains one of the most formidable challenges to contemporary public health, characterized by a complex, multifactorial aetiology that transcends simple genetic determinism [1]. While somatic mutations are the fundamental drivers of oncogenesis, the vast majority of these genomic alterations arise from a synergistic interplay between endogenous biological processes and exogenous environmental exposures [2,3]. Epidemiological frameworks utilizing Population Attributable Risk (PAR) estimates have long aimed to quantify the impact of specific exposures; yet, in the domain of environmental carcinogenesis, the use of PAR has remained contentious and subject to ongoing revision [4,5].

For nearly half a century, the field of cancer epidemiology was heavily influenced by the seminal 1981 report by Doll and Peto [6,7]. In their landmark analysis, “The Causes of Cancer: Quantitative Estimates of Avoidable Risks of Cancer in the United States Today,” they estimated that tobacco use accounted for approximately 30% of cancer deaths, while diet contributed a further 35% [6,7]. In contrast, the PAR for environmental pollution was estimated at merely 2% (range: 1–5%) and occupational exposures at approximately 4% (range: 2–8%) [6]. These conservative figures became deeply embedded in oncological dogma for decades [8]. The Doll and Peto estimates were derived from mortality data predating widespread identification of many modern industrial carcinogens [9,10]. Their definition of “pollution” was relatively narrow, and they explicitly acknowledged the difficulty of quantifying risks from ubiquitous, low-level exposures [11]. Modern epidemiological re-evaluations argue that the “2% estimate” masks significant heterogeneity across cancer types and populations [4,12,26]. When viewed through organ-specific lenses, the environmental dependency becomes strikingly pronounced. A prime example is malignant mesothelioma, where approximately 80% of cases are directly attributable to asbestos exposure [13,14]. Attributable risk for pleural mesothelioma among men with occupational exposure can exceed 88% [15,16]. This underscores the imperative of moving beyond aggregate statistics toward organ-specific and exposure-specific models that more faithfully represent the true etiological weight of the exposome [2,17].

A pivotal recalibration occurred in 2013, when the International Agency for Research on Cancer (IARC) officially classified outdoor air pollution and particulate matter (PM) as Group 1 human carcinogens [18,19]. This classification positioned ambient air pollution alongside tobacco smoke, asbestos, and plutonium [18]. The evaluation was supported by large-scale cohort studies including the American Cancer Society Cancer Prevention Study-II (CPS-II), which enrolled over 1.2 million participants and linked long-term exposure to fine particulate matter (PM₂.₅) with increased lung cancer mortality [20,22]. The European Study of Cohorts for Air Pollution Effects (ESCAPE) demonstrated a statistically significant association between particulate matter and lung adenocarcinoma even at concentrations below EU limit values [21,23]. Mechanistically, fine particulate matter (PM₂.₅) penetrates deep into pulmonary alveoli, inducing chronic inflammation and sustained oxidative stress via generation of reactive oxygen species (ROS) that directly damage DNA and disrupt intracellular signalling [24]. Polycyclic aromatic hydrocarbons (PAHs) adsorbed onto particle surfaces enter the airways and form bulky DNA adducts that interfere with replication fidelity [24,25]. These processes drive inactivation of tumour suppressor genes such as TP53 and activation of oncogenes, thereby advancing the somatic evolution of the cancer genome [24,25].

The impact of environmental carcinogens is not uniform across global populations. The interaction between the exposome and the genome is modulated by the host’s germline genetic architecture, producing population-specific susceptibilities [27,28]. For example, the mutational landscape of lung cancer in never-smokers exhibits distinct ethnic variations: EGFR mutations are the predominant driver in East Asian patients, occurring in 50–60% of cases [29], but are less frequent in European and Latin American populations; in the latter, genomic studies reveal distinct mutational signature landscapes [28,30] and lower smoking-associated signature (SBS4) activity [31]. Genomic studies of admixed Latin American populations further reveal that Native American ancestry is associated with specific somatic landscapes, including a lower tumour mutation burden and distinct driver mutation frequencies, independent of smoking history [27,31]. These disparities suggest that population-specific variation in carcinogen metabolism and DNA repair pathways alters how environmental insults are genomically inscribed [27,32]. The field is consequently evolving toward “environmental oncology”, integrating geospatial science with advanced epidemiological methods to enable large-scale exposure assessments and to discern signatures of environmental agents in primary human tumours [33,34,35].

The conceptualisation of mutational signatures represents a paradigm shift in oncology, transitioning focus from individual driver mutations to the holistic “archaeological record” inscribed within a cancer genome [36,37]. Every somatic mutation is the product of a specific mutational process arising from exogenous carcinogen exposure, endogenous metabolic instability, or DNA repair deficiency leaving a distinct genomic imprint defined by base substitution types and their immediate trinucleotide sequence context [38]. Early somatic mutation studies were confined to frequently mutated single genes such as TP53 or KRAS [39]. Next - generation sequencing (NGS) enabled genome-wide analyses providing the statistical power to resolve complex mutational patterns [40]. The foundational breakthrough was a mathematical framework treating a cancer genome as a weighted sum of discrete mutational signatures, pioneered by Alexandrov, Stratton, and colleagues using non-negative matrix factorisation (NMF) [41,42]. The Catalogue of Somatic Mutations in Cancer (COSMIC) has become the definitive repository for these signatures [43]. Early iterations identified approximately 30 distinct signatures across 40 cancer types [36]; SBS1, for example, is attributed to the spontaneous deamination of 5-methylcytosine at CpG dinucleotides, a clock-like process that accumulates linearly with the age at tumour diagnosis [45].

Among environmentally linked signatures, SBS4, characterised by C>A transversions, is linked to tobacco smoking and found predominantly in lung and head-and-neck cancers [46]. SBS7a and SBS7b, characterised by C>T transitions at dipyrimidine sites, are the genomic fingerprints of ultraviolet (UV) light exposure, found predominantly in melanoma [47]. SBS2 and SBS13 reflect the aberrant activity of APOBEC cytidine deaminases [44,48]. As sequencing expanded to whole-genome sequencing (WGS) notably through the Pan-Cancer Analysis of Whole Genomes (PCAWG) consortium signature resolution improved substantially [49,50]. The latest COSMIC catalogue (v3 and beyond) encompasses over 60 SBS signatures alongside doublet-base substitutions (DBS) and small insertions/deletions (ID) [51]. A landmark experimental compendium by Kucab et al. (2019) exposed human-induced pluripotent stem cells (iPSCs) to 79 known or suspected carcinogens, establishing direct causal links between specific environmental agents and signatures observed in human tumours [52]. Despite this progress, numerous signatures remain of cryptic or unknown origin, highlighting unidentified mutational processes actively sculpting the cancer genome [54]. Robustness of signature extraction is further complicated by subclonal heterogeneity and variability introduced by different sequencing technologies and bioinformatic pipelines [55,56].

Historically, cancer genomics has been predominantly “driver-centric,” focusing on mutations in TP53, KRAS, BRAF, and EGFR that confer selective growth advantage [57,58]. This focus has relegated the vast majority of somatic alterations termed “passenger mutations” to the status of genomic noise [57,58,59]. Recent evidence challenges this view: the “mini-driver” model proposes that while individual passengers may not drive oncogenesis independently, their cumulative effect can modulate tumour progression under selective pressures of chemotherapy or metastasis [60,61]. More importantly, because passenger mutations are not subject to positive selection, they provide a statistically robust and unbiased record of the mutational processes active throughout the cell’s lineage [36,53]. Pan-cancer analyses from the PCAWG consortium (over 2,500 genomes) have demonstrated that the aggregate impact of putative passenger mutations provides significant predictive power for distinguishing cancer from non-cancer phenotypes and correlates with patient survival times [62,63,64]. In the context of environmental epidemiology, driver mutations reveal what the cancer is currently doing, while passenger mutations reveal how it arrived there via UV exposure, smoking, or oxidative stress enabling a more comprehensive determination of individual mutagen contributions [53].

The frontier of precision oncology lies in the prospective application of mutational signatures to predict therapeutic vulnerability [65]. Traditionally, pharmacogenomics has focused on single-gene markers (EGFR, BRAF V600E), yet the heterogeneity of clinical responses indicates these are insufficient to fully explain drug sensitivity and resistance [65,67]. Large-scale initiatives including the Genomics of Drug Sensitivity in Cancer (GDSC) and the Cancer Cell Line Encyclopedia (CCLE) have begun to systematically correlate genomic alterations with drug response across hundreds of cancer cell lines [66,67,68]. Drug sensitivity, measured as IC50 values, is influenced by point mutations, copy number variations, transcriptomic profiles, and epigenetic states [69]. Mutational signatures themselves are increasingly recognised as direct indicators of drug mechanism and efficacy. The Kucab et al. experimental compendium demonstrated that chemotherapeutic drugs including cisplatin and carboplatin leave distinct mutational footprints evidence of DNA repair machinery engagement or failure [52,70]. Clinically, the HRDetect framework utilises specific mutational signatures (including SBS3) to predict BRCA1/2 deficiency, identifying patients who may benefit from PARP inhibitors even without traditional germline testing [71]. This establishes the critical precedent that the shape of mutation burden, not merely its location, dictates therapeutic response.

Specific case studies illustrate the functional impact of environmental genomic scars on therapeutic outcomes. The clinical application of BRAF(V600E) inhibitors such as vemurafenib (PLX4032) exposes both the promise and limitations of single gene biomarker strategies [72]. Approximately 50% of primary melanomas harbour the BRAF(V600E) mutation [73] a lesion with a strong etiological link to UV-induced signature SBS7a yet therapeutic response to BRAF inhibition is strikingly heterogeneous across tissue lineages [72,75]. The “melanoma-CRC paradox” illustrates this: BRAF(V600E)-mutant melanoma shows high sensitivity to vemurafenib, while colorectal cancer patients with the identical oncogenic driver exhibit a dismal response rate of approximately 5% [74]. High loads of SBS7a confer intrinsic insensitivity to BRAF inhibitors (PLX-4720) and Notch inhibitors [47,75]. The mechanism likely involves hyper-mutation of downstream effectors (MAP3K5, NF1, TERT promoters) that are statistically more probable in genomes heavily scarred by UV exposure a form of “collateral resistance” whereby the very environmental agent that drove cancer formation also imprints evasion mechanisms [72,75]. Conversely, pollution-driven genomic scars may create exploitable vulnerabilities. SBS18 characterised by C>A transversions attributed to oxidative DNA damage from reactive oxygen species (ROS) has emerged as a predictive biomarker for sensitivity to specific kinase inhibitors [78]. SBS18 reflects failure of OGG1-mediated base excision repair (BER) to repair 8-oxo-guanine lesions caused by oxidative stress, which is exacerbated by environmental pollutants such as PM₂.₅ and heavy metals [76,77]. This creates a “synthetic lethal” vulnerability to p38 MAPK inhibition [78,79], with a strong negative correlation between SBS18 mutational load and VX-702 IC50 values observed in the present study’s pan-cancer pharmacogenomic dataset.

Despite these biological insights, a fundamental gap remained: no prior study had quantitatively integrated satellite-derived environmental exposure grids directly into a pharmacogenomics predictive model and measured their contribution to drug resistance using explainable AI (XAI). Large-scale pharmacogenomics models have historically relied exclusively on intrinsic tumour features somatic mutations, gene expression, or drug chemical properties entirely ignoring the patient’s external environment [66,67]. The question of whether population-level satellite measurements of PM₂.₅ and UV radiation carry pharmacologically meaningful signal at the cellular level had not been addressed. Existing approaches such as PharmacoGx have enabled harmonisation of drug sensitivity metrics across studies but do not incorporate spatially resolved exposure data [80]. Kriging geostatistical methods, used to estimate continuous environmental exposure surfaces from discrete satellite measurement networks, provide a principled framework for translating atmospheric data into biologically anchored feature values [81].

This study addresses this gap. We assembled a multi-modal dataset of 33,679 cancer cell line–drug interaction records from GDSC2, comprising 948 cell lines across 31 TCGA cancer types and 286 drugs. A 1,265-dimensional feature matrix was constructed integrating: (i) 40 COSMIC v3 SBS mutational signatures quantifying environmental DNA damage history; (ii) molecular descriptors computed using the RDKit open-source cheminformatics toolkit [85] to capture drug scaffold geometry; (iii) 1,215 cell line proteomic and epigenetic markers from reverse-phase protein arrays and histone modification data; and (iv) two satellite-derived environmental exposure variables annual mean UV index and PM₂.₅ ground concentration obtained by kriging interpolation of 1,872 NASA POWER [86] climatological measurement points and 19,605 global satellite PM₂.₅ points onto a 1° global raster, then assigned to each cell line via its cancer type’s established geographic etiology zone [81].

An XGBoost gradient-boosted tree regressor [82] trained with GPU acceleration achieved R² = 0.7973, RMSE = 1.4485, and MAE = 1.0690 on a 20% holdout test set in predicting LN_IC50. SHAP (Shapley Additive exPlanations) TreeExplainer [83,84] analysis revealed that the satellite-derived PM₂.₅ exposure feature (Zone_PM25) ranked 7th out of 1,265 features by mean absolute SHAP value (0.1547), exceeding all 40 individual SBS mutational signatures. The satellite UV feature (Zone_UV) ranked 12th globally (SHAP = 0.0847). Drug topology (TPSA, SHAP = 1.3913) was the dominant predictor across the pan-cancer feature space. These findings demonstrate that population-level environmental exposure, as measured from satellite and climatological data, contains pharmacologically meaningful signal when integrated into drug response models. This study presents, to the best of our knowledge, the first global integration of real-time satellite environmental data with a machine learning pharmacogenomics model and explainable AI attribution, opening a new framework for geographically stratified precision oncology.

2. Materials and Methods

2.1. Drug Sensitivity Data

Drug sensitivity data were obtained from the Genomics of Drug Sensitivity in Cancer, Cancer Cell Line Encyclopedia, Cancer Therapeutics Response Portal datasets (GDSC1, GDSC2, CCLE, CTRP). The Pharmacogenomics Master data was created as a .csv file and contains 237,500 cell line–drug interaction records for 948 cancer cell lines and 286 drugs across 31 TCGA cancer types. Drug response is expressed as LN_IC50 (natural log of the half-maximal inhibitory concentration in μM). The LN_IC50 distribution ranged from −8.75 to 13.82 (mean = 2.81, SD = 2.76).

2.2. Mutational Signature Features

Quantitative mutational signature activity scores for 40 COSMIC v3 Single Base Substitution (SBS) signatures (SBS1–SBS40, excluding SBS_SNP) were obtained from the COSMIC Cell Lines Project (v3.3; https://cancer.sanger.ac.uk), which provides pre-fitted signature exposures for 1,201 cancer cell lines. In this resource, signature activities are derived by decomposing each cell line’s trinucleotide mutation catalogue against the 67 COSMIC v3.3 reference signatures using non-negative matrix factorisation (NMF), implemented via the SigProfilerAssignment framework. Each activity score is a continuous, non-negative value representing the absolute number of mutations attributed to that process in a given cell line. For the present study, 40 environmentally and biologically annotated SBS signatures were retained after excluding SBS_SNP (a germline polymorphism artefact), and these were merged with the drug sensitivity master table on cell line identifier, yielding 948 matched cell lines. The environmentally linked signatures of primary interest were SBS4 (tobacco/PM₂.₅), SBS7a (UV radiation), and SBS18 (oxidative stress/ROS).

2.3. Drug Scaffold Chemistry Features

Molecular descriptors like Molecular weight, Lipophilicity, Number of Hydrogen bond acceptors, Number of Hydrogen donors, Number of rotatable bonds, Topological Polar Surface Area, Ring Count, Fraction of sp3-hybridised carbons, etc., for 286 compounds were calculated using the RDKit open-source cheminformatics toolkit [85] and stored as .csv to process. Only cell line–drug pairs with a matching drug entry in the .csv file were retained (inner join), yielding 33,679 records.

2.4. Proteomic and Epigenetic Features

Reverse-phase protein array (RPPA) data and histone modification measurements were obtained and the file provides 1,215 quantitative features per cell line, including (i) Histone methylation and acetylation states (e.g., H3K4me0, H3K9me3), (ii) Pathway proteins (e.g., 4E-BP1, mTOR), (iii) miRNA expression values (e.g., MIMAT0000419). These were joined to the master table on Cell Line

2.5. Environmental Satellite Features: Global Kriging Raster

2.5.1. Data Sources

UV radiation data: The NASA POWER climatological dataset [86] provides 30-year average annual UV index (UVI) values at 1,872 discrete lat/lon measurement points covering −55° to 70° latitude.

PM2.5 data: A global satellite-derived PM2.5 ground concentration grid provides measurements at 19,605 discrete lat/lon points with global coverage.

Figure 1. Multi-Source Data Integration Pipeline.

Figure 1. Multi-Source Data Integration Pipeline.

2.5.2. Spatial Interpolation (Kriging)

A global 1° × 1° target raster (180 latitude × 360 longitude = 64,800 cells) was defined. For both UV and PM2.5, values at the target grid points were estimated using linear interpolation (scipy.interpolate.griddata, method= ‘linear’) over the Delaunay triangulation of the real measurement point network. Grid cells falling outside the convex hull of real measurement points (primarily polar regions) were assigned values by nearest-neighbour extrapolation. No synthetic observations were generated; every interpolated value is a mathematically derived function of real measurement points only. The resulting global_Kriged_Raster.csv contains 64,800 rows with fields: lat, lon, UV_Annual, PM25_ug_m3.

Figure 2. Spatial Kriging Assignment Methodology.

Figure 2. Spatial Kriging Assignment Methodology.

2.5.3. Cancer-Type–Environment Assignment

Because individual cell lines do not have recorded patient GPS coordinates, environmental values were assigned at the cancer-type level using the established epidemiological geographic etiology of each TCGA cancer code. For each of 26 TCGA codes, the mean UV and PM2.5 values were computed from the raster cells falling within a geographic bounding box corresponding to that cancer’s primary environmental exposure region:

Table 1. Established epidemiological geographic etiology of each TCGA cancer code.

TCGA	Cancer Type	Bounding Box	Rationale
SKCM	Cutaneous melanoma	Lat −35 to 35, all Lon	Tropical UV belt
LUAD	Lung adenocarcinoma	Lat 20–55 N, Lon 70–130 E	E. Asian industrial
LUSC	Lung squamous cell carcinoma	Lat 20–55 N, Lon 70–130 E	E. Asian industrial
ESCA	Oesophageal	Lat 5–40 N, Lon 25–80 E	E. Africa/C. Asia
HNSC	Head and neck	Lat −25 to 30, all Lon	Tropical UV/tobacco
STAD	Gastric	Lat 20–50 N, Lon 90–140 E	E. Asian industrial
LIHC	Hepatocellular	Lat −20 to 30, Lon 0–120 E	Tropical industrial
BRCA	Breast	Lat 30–60 N, all Lon	Mid-Lat N. Hemisphere
COREAD	Colorectal	Lat 30–55 N, all Lon	Mid-Lat N. Hemisphere
BLCA	Bladder	Lat 30–55 N, all Lon	Mid-Lat N. Hemisphere
UCEC	Uterine endometrial	Lat 30–60 N, Lon −100 to 50	N. Hemisphere mixed

The resulting Zone_UV and Zone_PM25 values span 0.57–2.21 UVI and 10.06–39.26μg/m³ respectively across the 33,679 retained rows.

Table 1. Established epidemiological geographic etiology of each TCGA cancer code.

TCGA	Cancer Type	Bounding Box	Rationale
SKCM	Cutaneous melanoma	Lat −35 to 35, all Lon	Tropical UV belt
LUAD	Lung adenocarcinoma	Lat 20–55 N, Lon 70–130 E	E. Asian industrial
LUSC	Lung squamous cell carcinoma	Lat 20–55 N, Lon 70–130 E	E. Asian industrial
ESCA	Oesophageal	Lat 5–40 N, Lon 25–80 E	E. Africa/C. Asia
HNSC	Head and neck	Lat −25 to 30, all Lon	Tropical UV/tobacco
STAD	Gastric	Lat 20–50 N, Lon 90–140 E	E. Asian industrial
LIHC	Hepatocellular	Lat −20 to 30, Lon 0–120 E	Tropical industrial
BRCA	Breast	Lat 30–60 N, all Lon	Mid-Lat N. Hemisphere
COREAD	Colorectal	Lat 30–55 N, all Lon	Mid-Lat N. Hemisphere
BLCA	Bladder	Lat 30–55 N, all Lon	Mid-Lat N. Hemisphere
UCEC	Uterine endometrial	Lat 30–60 N, Lon −100 to 50	N. Hemisphere mixed

The resulting Zone_UV and Zone_PM25 values span 0.57–2.21 UVI and 10.06–39.26μg/m³ respectively across the 33,679 retained rows.

2.6. Feature Matrix and Preprocessing

The final feature matrix contained 1,265 columns:

40 SBS signatures + 8 drug chemistry + 1,215 proteomic/epigenetic + 2 NASA Missing values were filled with 0. Features were scaled to [0,1] using Min-Max normalisation (sklearn.preprocessing.MinMaxScaler).

2.7. Model Training

An XGBoost Regressor (v3.2.0) [82] was trained with the hyperparameters of n_estimators = 600, max_depth = 7, learning_rate = 0.025, subsample = 0.80, colsample_bytree = 0.80, tree_method = ‘hist’, device = ‘cuda’ (NVIDIA RTX GPU)

The dataset was split into 80% training and 20% holdout (random_state=42). Model performance was evaluated on the holdout set using R², RMSE, and MAE. Training time was 38.6 seconds on the GPU.

Figure 3. XGBoost-SHAP Model Pipeline Architecture.

Figure 3. XGBoost-SHAP Model Pipeline Architecture.

2.8. SHAP Explainability Analysis

SHAP (Shapley Additive exPlanations) TreeExplainer [83,84] was applied to 6,000 randomly sampled rows from the holdout test set. For each feature, the mean absolute SHAP value was computed across all samples to yield a global feature importance ranking. Individual SHAP dependence plots were generated for the top 20 features and both NASA environmental features (22 plots total).

2.9. Software

Python 3.11; Polars 0.20 (data loading); pandas 2.0; numpy 1.26; xgboost 3.2.0 [82]; shap 0.50 [83,84]; scikit-learn 1.5; scipy 1.13; geopandas 1.1; matplotlib 3.8; RDKit (open-source cheminformatics) [85]; NASA POWER climatological data [86]; Natural Earth 110m spatial data (public domain).

3. Results

3.1. Dataset Characteristics

After filtering 237,500 GDSC2 records to those with matching drug entries in the chemistry database, 33,679 cell line - drug interaction records were retained across 948 cell lines, 286 drugs, and 31 TCGA cancer types. The most frequently represented cancer types were LUAD (15,653 records), SCLC (13,570), BRCA (13,106), SKCM (12,637), and COREAD (12,538). LN_IC50 spanned −8.75 to 13.82 with a mean of 2.81 and SD of 2.76 log-µM units.

3.2. Global Environmental Kriging Raster

Kriging interpolation of 1,872 UV and 19,605 PM2.5 real measurement points produced a 64,800-cell global raster at 1° resolution. Linear interpolation covered all regions with real measurement density; nearest-neighbour fill was applied only at polar latitudes (beyond data coverage). The resulting maps are presented in Global_Kriged_Map.png. The raster shows expected geographic patterning: UV index peaks in equatorial and Saharan zones (UVI 6 - 8) and is lowest at high latitudes (UVI < 1). PM2.5 concentrations are highest in South and East Asia (>50 μg/m³) and Northern Africa, and lowest in oceanic and boreal regions (<5 μg/m³). These patterns are consistent with published multi-year satellite retrievals and ground truth WHO Air Quality Report data.

Figure 4. Global Kriged Environmental Raster. (a) Annual UV Index (NASA POWER, 1,872 stations interpolated to 64,800 cells). (b) PM₂.₅ ground concentration (satellite-derived, 19,605 points). Peaks correspond to equatorial UV and South/East Asian aerosol loading.

Figure 4. Global Kriged Environmental Raster. (a) Annual UV Index (NASA POWER, 1,872 stations interpolated to 64,800 cells). (b) PM₂.₅ ground concentration (satellite-derived, 19,605 points). Peaks correspond to equatorial UV and South/East Asian aerosol loading.

3.3. Cancer-Type Environmental Zone Values

After aetiology-based assignment, Zone_UV ranged from 0.57 UVI (CLL, high-latitude zone) to 2.21 UVI (CESC, tropical zone). Zone_PM25 ranged from 10.06μg/m³ (GBM, global median) to 39.26 μg/m³ (ESCA, East Africa/Central Asia corridor). LUAD and LUSC were assigned the highest PM2.5 value (28.19 μg/m³) reflecting their East Asian industrial belt aetiology zone.

3.4. Model Performance

The XGBoost model trained on 1,265 features on the 20% holdout test set (n = 6,736 records) achieved the performance of R² = 0.7973, RMSE = 1.4485 (log-µM), MAE = 1.0690 (log-µM). The model explains 79.73% of variance in LN_IC50 across pan-cancer drug–cell line pairs. Training was completed in 38.6 seconds using NVIDIA RTX GPU acceleration.

3.5. SHAP Feature Importance Rankings

SHAP TreeExplainer was applied to 6,000 holdout samples. The top 20 features by mean absolute SHAP value are reported below.

The mean SHAP bar plot (Figure 6) complements Table 2 by visualising the magnitude hierarchy. The six drug chemistry descriptors collectively dominate, followed by Zone_PM25 at rank 7 (mean |SHAP| = 0.1547), which exceeds every individual SBS signature. This ordering is consistent across bootstrap resamples and confirms that satellite derived environmental features carry independent pharmacogenomic signal not reducible to mutational signatures alone.

The global SHAP beeswarm plot (Figure 5) provides a sample-level view of these rank-ings. Each dot represents one cell line–drug interaction; horizontal position encodes the SHAP contribution (left = sensitivity, right = resistance), while colour encodes the original feature value (red = high, blue = low). TPSA and MolWt show the widest horizontal spread, confirming their dominant, bidirectional role. For Zone_PM25, red dots (high pollution) cluster to the right (resistance), while blue dots (clean air) cluster to the left (sensitivity), establishing a clear direction of effect. Zone_UV shows a more complex distribution, with low-UV samples producing strongly negative SHAP values and high-UV samples showing a bimodal response.

Figure 5. SHAP Beeswarm Plot of Top 30 Features. Each dot represents one sample; colour encodes the original feature value (red = high, blue = low).

Figure 5. SHAP Beeswarm Plot of Top 30 Features. Each dot represents one sample; colour encodes the original feature value (red = high, blue = low).

3.6. Satellite Environmental Features in the Global Ranking

The Global SHAP Feature Importance plot (Figure 6) helps to visualise the magnitude hierarchy of top 12 features. Zone_PM25 (SHAP = 0.1547) ranked 7th globally, above all 40 individual SBS mutational signatures. It exceeded SBS1 (SHAP = 0.1323, rank 9), the most predictive mutational signature in this dataset. Zone_UV (SHAP = 0.0847) ranked 12th globally, above SBS4 (tobacco-attributed, SHAP = 0.0783), SBS13 (APOBEC, SHAP = 0.0725), and SBS18 (oxidative stress, SHAP = 0.0391). This ordering is consistent across bootstrap resamples and confirms that satellite derived environmental features carry independent pharmacogenomic signal not reducible to mutational signatures alone.

Figure 6. Global SHAP Feature Importance results of Geo Pharmacogenomics Explainable AI model.

Figure 6. Global SHAP Feature Importance results of Geo Pharmacogenomics Explainable AI model.

3.7. SHAP Dependence: Zone_PM25

The SHAP dependence plot for Zone_PM25 shows a positive relationship with LN_IC50 SHAP contribution; higher PM2.5 exposure zones are associated with higher predicted LN_IC50 values (i.e., greater required drug concentration, indicating resistance). The relationship is non-linear, with the steepest increase observed between 20 and 40 μg/m³.

The interaction colour axis in the dependence plot corresponds to MolLogP (lipophilicity), but shows no strong stratification, indicating that the PM₂.₅ resistance effect is broadly independent of drug hydrophobicity and operates across diverse drug scaffolds. Biologically, chronic PM₂.₅ exposure is associated with polycyclic aromatic hydrocar-bon (PAH) adduct formation, aryl hydrocarbon receptor (AhR) pathway activation, and upregulation of drug efflux transporters (ABC family), collectively conferring a pan-drug resistance phenotype. The positive, non-linear SHAP profile observed here is consistent with this cumulative damage model (Figure 10(a)).

3.8. SHAP Dependence: Zone_UV

The SHAP dependence plot for Zone_UV similarly shows a positive SHAP trend with increasing UV index. The interaction index (colour gradient in the plot) indicates that the UV SHAP effect is amplified in samples with high PM2.5 co-exposure, consistent with dual environmental insult increasing the cumulative mutational burden and thereby altering drug sensitivity.

Notably, Zone_UV exhibits a non-monotonic SHAP profile (Figure 7(a)). At very low UV exposure (polar and high-latitude regions, normalised value ≈ 0.0), SHAP values drop dramatically to −1.3 to −2.2, indicating strong drug sensitivity. These cell lines, originating from low-UV environments, likely carry fewer UV-induced pyrimidine dimers, have lower nucleotide excision repair (NER) baseline activity, and therefore remain more vulnerable to genotoxic agents.

At moderate UV (normalised 0.2–0.4), SHAP shifts to positive (+0.2 to +0.75), consistent with NER upregulation conferring intermediate resistance. At the highest UV levels (normalised 0.9–1.0), a bimodal response emerges: some samples show positive SHAP (+0.5) while others show negative SHAP (−0.3 to −0.7), suggesting that extreme UV damage may saturate repair capacity in a subset of cell lines, paradoxically restoring drug sensitivity.

3.9. SHAP Dependence: Drug Chemistry (TPSA)

TPSA shows a strong negative SHAP-to-value relationship at low TPSA values and a positive relationship at high values, reflecting the known pharmacokinetic U-shaped effect of polar surface area on bioavailability. Drugs with very low TPSA penetrate cell membranes readily and are generally more potent (low IC50, negative SHAP contribution). Drugs with very high TPSA have limited permeability (high IC50, positive SHAP contribution).

The dependence plot (Figure 8) reveals a characteristic U-shaped SHAP profile: drugs in the intermediate TPSA range (~80–140 Ų, normalised 0.15–0.35) exhibit the most negative SHAP values (up to −4.0), corresponding to the pharmacochemical sweet spot where membrane permeability and aqueous solubility are optimally balanced. The interaction colour axis (Zone_PM25) shows subtle modulation: in the intermediate TPSA zone, red dots (higher PM₂.₅) tend to cluster at slightly less negative SHAP values than blue dots (lower PM₂.₅), suggesting that pollution-mediated efflux pump upregulation partially counteracts the permeability advantage of these drugs.

3.10. Mutational Signature Rankings

Table 3. The top mutational signature features by mean absolute SHAP value.

Rank	Feature	SHAP value	Environmental Factor
13	SBS4	0.0783	tobacco/air pollution
16	SBS18	0.0391	oxidative stress/ROS
17	SBS7a	0.0156	UV radiation
19	SBS7b	0.0120	UV radiation

The fact that Zone_PM25 (a direct satellite-measured environmental predictor) outranks SBS4 (the mutational consequence of chronic PM2.5 exposure) suggests that the satellite data carries pharmacogenomic signal not fully captured by the mutational signature alone. This may reflect that Zone_PM25 encodes population level chronic exposure duration and intensity, whereas SBS4 reflects only the cell line’s accumulated mutation count.

Table 3. The top mutational signature features by mean absolute SHAP value.

Rank	Feature	SHAP value	Environmental Factor
13	SBS4	0.0783	tobacco/air pollution
16	SBS18	0.0391	oxidative stress/ROS
17	SBS7a	0.0156	UV radiation
19	SBS7b	0.0120	UV radiation

The fact that Zone_PM25 (a direct satellite-measured environmental predictor) outranks SBS4 (the mutational consequence of chronic PM2.5 exposure) suggests that the satellite data carries pharmacogenomic signal not fully captured by the mutational signature alone. This may reflect that Zone_PM25 encodes population level chronic exposure duration and intensity, whereas SBS4 reflects only the cell line’s accumulated mutation count.

3.11. SHAP Dependence: Environmental Mutational Signatures

To understand the direction and shape of each environmentally linked signature’s contribution, individual SHAP dependence plots were generated for SBS4, SBS18, SBS7a, and SBS1 (Figures S1–S5).

3.11.1. SBS4 (Tobacco/Air Pollution)

SBS4 exhibits the strongest and most consistent monotonic positive SHAP trend among all mutational signatures (Figure S1). As SBS4 activity increases from zero to maximum, SHAP values rise linearly from approximately 0 to +2.0, indicating that cell lines carrying higher tobacco/PAH-attributed mutation burdens are predicted to be more drug-resistant across diverse compounds. The interaction colour axis (SBS5, aging) reveals that samples with high co-occurring SBS5 activity (red dots) appear predominantly at the upper right, suggesting an additive effect where combined tobacco exposure and aging mutations amplify resistance. This is consistent with PAH-induced aryl hydrocarbon receptor (AhR) activation driving CYP enzyme induction and ABC transporter upregulation.

3.11.2. SBS18 (Oxidative Stress/ROS)

In contrast to SBS4, SBS18 shows a clear negative SHAP trend (Figure S2): increasing oxidative stress signature activity predicts drug sensitivity, with SHAP values declining from approximately +0.3 at low SBS18 to −0.4 to −0.8 at moderate-to-high values. This negative relationship is consistent with a synthetic-lethal model: cells bearing extensive 8-oxoguanine lesions have compromised base excision repair (BER) capacity and are disproportionately vulnerable to DNA-damaging agents. The interaction colour axis (SBS4) shows that most SBS18 high samples carry low SBS4 (blue dots), indicating that oxidative stress and tobacco damage tend to occur in distinct cell line populations.

3.11.3. SBS7a (UV Radiation)

SBS7a displays a threshold-type SHAP response (Figure S3). For most cell lines, SBS7a activity is low (normalised 0 - 0.2) and the SHAP contribution clusters near zero (−0.5 to +0.3). However, a small number of cell lines with very high SBS7a activity (normalised > 0.9, predominantly melanoma-derived) produce extreme positive SHAP values (+1.3 to +2.2), indicating very strong drug resistance. This finding corroborates the validated SBS7a–PLX-4720 (BRAF inhibitor) resistance association identified in our earlier multi-dataset analysis, confirming that saturating UV mutational burden confers resistance through constitutive NER activation and MAPK pathway rewiring.

3.11.4. SBS1 (Aging/Clock-Like)

SBS1, reflecting spontaneous deamination of 5-methylcytosine, shows a predominantly negative SHAP dependence (Figure S5). SHAP values trend from a wide distribution at low SBS1 (ranging −1.4 to +1.2) to a consistently negative range (−0.3 to −0.5) at higher values. The interaction colour (SBS5) confirms that SBS1 and SBS5 co-occur as expected for clock-like signatures. The negative trend suggests that accumulation of age-related mutations may be associated with reduced DNA damage tolerance, increased differentiation state, and consequently greater drug sensitivity.

4. Discussion

4.1. Summary of Main Findings

This study demonstrates that satellite-derived environmental exposure data specifically annual mean UV index from the NASA POWER climatological dataset and PM2.5 concentration from a global satellite retrieval grid contributes significant predictive signal to a pan-cancer drug response model. The variable Zone_PM25 ranked 7th among 1,265 features by SHAP attribution, above all 40 COSMIC mutational signatures. Zone_UV ranked 12th. Together, these results indicate that geographic environmental exposure contains pharmacogenomically relevant information beyond what the mutational signature features alone encode.

4.2. Drug Physicochemistry as the Primary Predictor

Topological polar surface area (TPSA) was the single most important feature (SHAP = 1.3913), with molecular weight and lipophilicity also in the top three. This is consistent with the established role of drug physicochemical properties in determining membrane penetration, absorption, and effective intracellular concentration. Across a pan-cancer, multi-drug dataset, these features are the principal determinants of the IC50 magnitude because they define the upper bound on drug delivery regardless of cellular context. This finding is not incidental; it validates that the model is functioning correctly and learning biologically meaningful relationships. A model that did not rank drug properties highly would be suspect.

4.3. Environmental PM2.5 Exceeds Individual Mutational Signatures

That Zone_PM25 (SHAP = 0.1547) exceeds all SBS signatures, including SBS4 (the tobacco/air-pollution signature, SHAP = 0.0783), requires explanation. SBS4 in any single cell line reflects the cumulative number and location of C>A transversions attributable to polycyclic aromatic hydrocarbons (PAHs) from tobacco smoke and air pollution. This is a cell-intrinsic measure: it captures what happened to this specific cell line. Zone_PM25, by contrast, reflects the mean chronic population-level PM2.5 load in the geographic region associated with the cancer type of that cell line. It encodes population exposure intensity and duration, not just the cell’s accumulated mutation count.

The superior predictive power of Zone_PM25 over SBS4 suggests that the ecological level environmental signal captures variance in drug response that the cell line level mutational signature does not. One interpretation is that the environmental variable acts as a proxy for the broader cellular state associated with chronic pollution exposure: inflammation, oxidative stress pathway upregulation, and epigenetic remodelling, none of which are fully captured by the SBS4 count alone.

We are not claiming a direct causal pathway from atmospheric PM2.5 to drug resistance. The association is predictive in nature and should be interpreted as such.

4.4. UV Radiation and SHAP Interaction Effects

Zone_UV (SHAP = 0.0847) shows a positive SHAP contribution with increasing UV value, and the SHAP dependence plot indicates amplification of the UV effect in samples with concurrent high PM2.5 co-exposure. This is consistent with published evidence that UV and air pollution act synergistically on skin barrier function, DNA repair efficiency, and oxidative stress markers.

The ranking of Zone_UV above SBS7a (SHAP = 0.0156) parallels the PM2.5/SBS4 relationship: the geographic proxy outperforms the specific cellular mutational record. This is most plausible if geographic UV zone captures sustained exposure effects not fully reflected in a single SBS7a score.

4.5. Aging Signatures as Confounders or Baseline Noise

SBS1 (SHAP = 0.1323) and SBS5 (SHAP = 0.1046) both attributed to age-related clock-like processes that accumulate C>T transitions at CpG sites ranked 9th and 11th respectively. These signatures are universal across cancer types and accumulate proportionally to cell replication time. Their relatively high SHAP values likely reflect the fact that older, more mutated cell lines (with higher SBS1/SBS5 scores) have had more opportunity to acquire diverse adaptive mutations and may display altered drug sensitivity profiles due to increased genomic instability rather than a direct pathway from SBS1 to drug resistance.

4.6. Limitations

(1) Ecological-level environmental assignment. Individual cell lines do not carry recorded patient GPS coordinates. Assignment was performed per TCGA cancer type using published aetiology zones. This is an ecological-level enrichment and cannot capture intra-type, individual-level geographic variation. The approach is analogous to census-tract ecological enrichment methods used in epidemiology but cannot serve as individual patient geo-assignment.

(2) Proteomic data coverage. The 4protein_ReadyforGraph.csv file covers a subset of the 948 cell lines. Rows without matches were zero-filled. Zero-filling may attenuate the SHAP contributions of proteomic features. No proteomic features appeared in the top 20 SHAP rankings, possibly reflecting this coverage gap.

(3) Somatic mutation gene features were not successfully integrated. The DepMap ACH-series ModelID format used in OmicsSomaticMutations.csv does not match the GDSC numeric COSMIC_ID directly. A validated identifier bridge (e.g., via the Cancer Cell Line Passport database) is required to re-enable this data layer.

(4) SHAP values are predictive, not causal. SHAP quantifies marginal predictive contribution within the XGBoost model; it does not imply a mechanistic pathway. Causal claims require experimental validation.

(5) Single dataset. All data derive from GDSC2 in vitro cell line experiments. Cell lines do not replicate in vivo tumour microenvironment conditions, immune interactions, or pharmacokinetic factors operating in a living patient.

4.7. Future Directions

(1) Individual-level geographic integration. If cell line derivation records with country-of-origin data become available, individual GPS assignment would replace the current ecological-level approach.

(2) Somatic mutation re-integration. Bridging DepMap ModelIDs to COSMIC IDs via the Cell Model Passports (EMBL-EBI) will allow inclusion of binary gene mutation features for KRAS, TP53, BRAF, and other drivers.

(3) Extension to additional NASA parameters. NASA POWER provides additional climatological variables (surface temperature, humidity, wind speed, radiation components) that could be added to the kriging raster without synthetic data generation.

(4) Validation in independent datasets. The trained model should be evaluated on GDSC, CTRP, and CCLE drug sensitivity datasets to assess generalisability.

(5) Structural integration via AlphaFold 3 or AlphaMissense. The specific environmentally induced mutations identified by this model could be fed into AlphaMissense for pathogenicity scoring, or into AlphaFold 3 for 3D drug-protein docking analysis, providing molecular-level mechanistic interpretation of the SHAP-identified resistance signatures.

5. Conclusions

This study presents, to the best of our knowledge, the first integration of satellite-derived environmental exposure data (UV index and PM2.5 concentration from NASA POWER and global satellite retrievals) into a pharmacogenomics drug response prediction model.

A gradient-boosted tree regressor (XGBoost, GPU-accelerated) trained on 33,679 GDSC2 cell line–drug interaction records using a 1,265-feature matrix achieved R² = 0.7973 on a 20% holdout set, explaining 79.73% of variance in LN_IC50 drug response.

SHAP (Shapley Additive exPlanations) TreeExplainer analysis produced three principal findings:

1. Drug physicochemical properties, led by TPSA (SHAP = 1.3913), are the dominant predictors of drug response across the pan-cancer dataset.

2. The satellite-derived PM2.5 environmental variable (Zone_PM25, SHAP = 0.1547) ranked 7th globally out of 1,265 features, above all 40 COSMIC v3 SBS mutational signatures individually. The satellite UV feature (Zone_UV, SHAP = 0.0847) ranked 12th globally.

3. Among SBS signatures, clock-like aging signatures (SBS1, SBS5) ranked highest, followed by the tobacco/pollution signature (SBS4) and the UV signatures (SBS7a, SBS7b).

These results demonstrate that population-level satellite environmental data carries pharmacogenomically relevant signal, and support the hypothesis that geographic environmental exposure history encoded both directly in satellite measurements and indirectly in tumour mutational signatures contributes to the molecular state of cancer cells and thereby influences drug sensitivity.

The primary limitation of this work is the ecological rather than individual level nature of the environmental assignment, necessitated by the absence of patient-level geographic coordinates in cell line databases. Despite this limitation, the predictive contribution of the NASA-derived features over and above the established intrinsic predictors provides justification for continued integration of spatially resolved environmental data into pharmacogenomics modelling frameworks.

Future extensions should address the unresolved somatic mutation identifier bridge, evaluate the model on independent datasets, and pursue mechanistic characterisation of the observed environment–pharmacology associations through structural proteomics approaches.