Radiomic biomarkers of tumor immune biology and immunotherapy response.

Immunotherapies are leading to improved outcomes for many cancers, including those with devastating prognoses. As therapies like immune checkpoint inhibitors (ICI) become a mainstay in treatment regimens, many concurrent challenges have arisen - for instance, delineating clinical responders from non-responders. Predicting response has proven to be difficult given a lack of consistent and accurate biomarkers, heterogeneity of the tumor microenvironment (TME), and a poor understanding of resistance mechanisms. For the most part, imaging data have remained an untapped, yet abundant, resource to address these challenges. In recent years, quantitative image analyses have highlighted the utility of medical imaging in predicting tumor phenotypes, prognosis, and therapeutic response. These studies have been fueled by an explosion of resources in high-throughput mining of image features (i.e. radiomics) and artificial intelligence. In this review, we highlight current progress in radiomics to understand tumor immune biology and predict clinical responses to immunotherapies. We also discuss limitations in these studies and future directions for the field, particularly if high-dimensional imaging data are to play a larger role in precision medicine.

Introduction

With improvements in high-throughput computing, tumor imaging data are playing an increasingly important role in aiding the comprehension of tumor biology [1] and guiding decision-support for precision medicine [2], [3], [4]. Radiomics refers to the process of quantifying patient-specific tissue characteristics from images to produce high-dimensional minable features, or image biomarkers. These features can span from simple descriptive quantitative information, such as geometry (morphologic properties like sphericity) and intensity (descriptive statistics like mean and variation), to complex higher order relationships like texture (e.g. correlations between image voxels indicating properties like heterogeneity) [3], [4], [5], [6]. Such images are typically obtained from CT, MRI, and PET, and can be further processed with filters or transforms (e.g. Gabor, wavelet) prior to feature extraction. The underlying hypothesis is that certain aspects of tumor biology, such as vascularization, degree of inflammation, spatial organization and heterogeneity, define architectural features that can be captured by high-resolution imaging and quantified mathematically.

Radiomic data can be powerful resources for clinical oncology for several reasons. First, imaging studies are non-invasive. Second, almost all cancer patients undergo imaging during their diagnosis and treatment course, making these types of data more universally available. Third, imaging studies can capture the full spatial organization of tumors, thereby overcoming the critical limitations of sampling from biopsies. Last but not least, imaging data have the potential to be used longitudinally to track and provide biomarkers of treatment response and resistance phenotypes.

For more than a decade, studies have highlighted the potential for imaging data to proxy molecular tumor phenotypes, a field that has been commonly referred to as imaging genomics or radiogenomics [4], [7], [8]. Some of the earliest work in the field focused on gene expression correlates of semantic (i.e. clinician-defined) tumor features like presence of edema and necrosis [9]. More recently, advances in computing have enabled high throughput as well as automated quantitative feature extraction. Rather than being limited to a handful of features, radiogenomic correlations can now be performed across hundreds of image characteristics [3], [10], [11], [12]. A larger array of features that define more complex tumor architectures may bolster the analysis of complex and heterogenous areas of tumor biology, such as the tumor immune microenvironment, that have eluded our understanding thus far.

Understanding the tumor immune microenvironment is paramount given its critical role in determining tumor response to anti-cancer therapies like immune checkpoint inhibitors (ICIs). Inflammation of the tumor microenvironment (TME) can be pro-tumor or anti-tumor depending on the immune context [13]. Chronic inflammation can promote tumor growth by enabling autocrine signaling loops that drive cell proliferation, increased angiogenesis, immune evasion, and metastasis [14]. Both innate and adaptive immune cells interact with each other and the tumor to regulate inflammation, ultimately dictating the course of tumor progression. The balance of anti-tumor immune cells, such as cytotoxic T cells and NK cells, with cells that suppress the anti-tumor immune response, particularly regulatory T cells (Tregs), myeloid derived suppressor cells (MDSCs), and tumor associated macrophages (TAMs), also governs the extent to which cancers respond to certain therapies. While immune-suppressive cells like TAMs undesirably mitigate the response to chemotherapy [15] and ICIs [16], tumor-infiltrating lymphocytes (TILs) are positive predictors of treatment response following chemoradiation [17], [18], [19] and ICIs [20], [21], [22]. The presence of certain surface markers and gene expression profiles also determines the tumor immune state. For instance, high expression of PD-L1 by cells in the TME associates with better response to PD-1 blockade, whereas tumors without expression of PD-L1 respond poorly [21], [22], [23], [24].

Last but not least, tumor response to immunotherapy is also complicated by incidences of partial increase in volume followed by response, known as pseudoprogression. The prevalence of pseudoprogression may be as high as 20% [25], which can result in inappropriate discontinuation of therapy. Identifying biomarkers to distinguish pseudoprogression from progression or necrosis is critical to providing proper care. As evidenced by these gaps in our knowledge of how to effectively administer immunotherapy, finding robust biomarkers that represent the tumor immune state is a priority. Recent studies suggest that imaging informatics may be poised to address these challenges.

Overview of radiogenomic studies on tumor immune biology and immunotherapy response

Based on a comprehensive review of the literature (search methodology in Fig. 1), we identified 54 studies in which imaging features were associated with tumor immune phenotypes (Table 1 and Supplemental Table 1) or response to immunotherapy (Table 2 and Supplemental Table 2). The following types of cancers were investigated: lung [3], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44] (20), brain [45], [46], [47], [48], [49], [50], [51], [52], [53], [54] (10), breast [55], [56], [57], [58], [59], [60], [61], [62], [63] (9), liver [9], [41], [64], [65], [66], [67] (6), skin [44], [68], [69], [70], [71] (5), head and neck [41], [72], [73], [74] (4), along with single studies focused on several other cancer types [75], [76], [77]. The central questions collectively addressed by these studies were 1) whether imaging features could consistently and robustly represent underlying immune biology and 2) whether radiomic biomarkers and radiogenomic models could predict immunotherapy response.

Fig. 1

Search criteria for studies reporting radiogenomic associations or associations between imaging features and response to immunotherapy.

Table 1

Overview of key studies reporting radiogenomic associations with tumor immune phenotypes.

Reference	Tumor type	Modality	Primary study (ortraining set)	Validation	Notable immune associations	Features and feature classes assessed	Summary of findings	Statistics
3	Lung (NSCLC)	CT	89	0	Leukocyte activation and regulation of immune system process	440 radiomic features, including intensity, shape, texture and wavelet features	Lymphocyte activation, leukocyte activation, and regulation of immune system process positively associated with statistics total energy. Lymphocyte activation positively associated with shape compactness.	FDR < 0.2 and NES > 0 from GSEA of ranked correlation of features with GO gene sets
26	Lung (NSCLC)	PET/CT	25	147 (63external cohort, 84 validation)	Antigen presentation and processing, immune response, NFkB signaling	14 features related to SUV including intensity, distribution, and spatial metrics	Positive associations include pSUV mean with cell cycle and immune response; pSUV PCA2 with antigen presentation and processing; pSUV max with NFkB on network analysis; Multivariate-pSUV with cell/antigen processing, immune response.	FDR < 0.05, enrichment of gene sets from GeneSigDb, DAVID, MSigDb, and Reactome
32	Lung (NSCLC)	CT	89	60	Type 2 helper T-cell (Th2) signature	239 radiomic features used in machine learning models to predict tumor immune microenvironment (TIME) signaturescomputed via GSVA	Type 2 helper T-cell (Th2) expression signature positive correlation with skewness, kurtosis, variance, and informational measure of correlation(IMC).	AUC 0.684 (test), p = 0.027 for linear discriminant model
35	Lung (NSCLC)	PET/CT	57	0	PD-L1 expression by IHC	SUVmax, MTV, TLG; radiomic featuresincluding size, shape, first-order, and second-order features	Low coarseness and higher GLZLM_ZLNU associated with PD-L1 (medium/high vs low).	p = 0.025 (coarseness), p = 0.035 (GLZLM_ZLNU), Kruskal-Wallis test
36	Lung (NSCLC)	PET/CT	374	0	PD-L1	SUVmax, SUVmean, primary (-P) and combined (-C) MTV and TLG	SUVmax positively correlated with PD-L1 on multivariable analysis; TLG-P and TLG-C onunivariate analysis.	p < 0.01 for TLG-P/C, univariate logistic regression; p < 0.001 for SUVmax, multivariate logistic regression
37	Lung (NSCLC)	PET/CT	55	0	CD8, PD-1 TILexpression by IHC	SUVmax, SUVmean	SUVmax/mean positively correlated with CD8, PD- 1 (but not PD-L1, CD68).	p = 0.027 (SUVmax/mean) for CD8; p = 0.017 (SUVmax), p = 0.009 (SUVmean) for PD-1
38	Lung (NSCLC)	CT	262	89	Gene modules enriched for immune pathways, NFkB activation	636 radiomic features, bi-clustering used to establish modules of radiomic- pathway coherency in training set, which identified 13 modules in validation set	Three modules M2, M9, M12 (quantified textural entropy and dispersion image intensity values) associated with overall survival were enriched for immune system. For M10, shape compactness and sphericity predicts NFkB activation.	AUC 0.66 (p = 0.003) for M10 feature prediction; FDR < 0.05 for all reported module associations
39	Lung (NSCLC)	CT	114	176	PD-L1 and CD3expression by IHC	490 features; final model of 4 features (mean, standard deviation, and uniformity as primary features; GLCM_homogeneity as secondary feature) clustered based on PD-L1/CD3 expression	Inferences of associations based on model clustering (but not shown explicitly in univariate analysis): low PD-L1/high CD3 (Cluster D) associated with low mean, uniformity, GLCM_homogeneity and high SD.	Multinomial regression (p = 0.01 training; p < 0.001 validation)
40	Lung	PET/CT	263	0	PD-L1 expression by IHC	SUVmax	SUVmax positively associated with PD-L1.SUVmax predicts PD-L1 positivity.	p < 0.001, Spearman correlation; AUC 0.797, p < 0.0001 for logistic regression prediction model
41	Mixed (HNSCC, NSCLC, HCC, BLCA)	CT	135	219 (119TCGA, 100Gustave Roussy)	TIL density by IHC, CD8 expression signature	84 imaging features for machine-learning trained on CD8 expression; final elastic- net model included 5 features: (-) coefficient for tumor min value, tumor GLRLM_SRHGE; (+) coefficient for ring GLRLM_SRLGE, ring GLRLM_LGRE, ring GLRLM_LRLGE	Radiomics score positively correlated with TIL density and predicts CD8 expression signature.	AUC 0.74 (training), 0.67 (TCGA), and 0.76 (Gustave Roussy) for score prediction of CD8 signaturep = 0.00022, Spearman correlation for score and TIL density
45	Glioma (GBM)	MRI	55	0	Module comprised of genes in IL4, T-cell differentiation and proliferation	79 features per ROI (3 ROIs), including: necrotic edge sharpness, minor axis length, radial distance signal, skewness, median, mean, min	Module 20 (enriched for IL4 and T cell differentiation/proliferation) positively associated with blurry (vs sharp) edge of tumor necrosis.Correlations between several quantitative features and pathways in supplemental heatmaps (see study).	FDR < 0.05, correlations with pathways from KEGG
47	Glioma (GBM)	MRI	91	0	Inflammatory and immune response pathways	Primary features: contrast enhancement (CE), edema (ED), volume (TV), bulk (TB), necrosis (NE)Feature ratios: NE/TV, CE/TV, ED/TV, TB/TV, NE/CE, CE/TB	64 pathways associated with primary features or tumor-volume normalized features. Tumor bulk and necrosis anti-correlated with immune system/response. NE/CE anti-correlated with immune system and NFkB. CE/TB correlated with NFkB and immune response.	FDR < 0.05, GSEA of ranked correlations with GO gene sets
48	Glioma (GBM)	MRI	50	0	Module enriched for dendritic cell biology and adaptive immunity	ADC mean, standard deviation, skewness, kurtosis, and entropy	Negative correlation between mean ADC and module 5 immune gene module (including genes related to dendritic cell biology and adaptiveimmunity).	p = 0.001, Spearman correlation
50	Glioma (GBM)	MRI	35	34 (internal cohort)	CD3 T cell infiltration by IHC	86 radiomic features; 6 used in model (histogram kurtosis, NGTDM contrast, GLSZM small zone size emphasis, GLSZM low gray-level zone emphasis, GLSZM high gray-level zone emphasis, GLSZM small zone high gray emphasis)	Best single predictor for CD3 (T cells) was GLSZM small zone high gray emphasis (AUC 0.79); full 6- feature model performed best.	AUC 0.847 (validation, full model) p = 0.009, Spearman correlation (between predictionand calculated CD3)
51	Glioma (anaplastic)	MRI	91	0	Inflammatory andimmune response pathways	T1-weighted contrast enhancement	Immunity-associated pathways enriched in contrast-enhanced tumors, including immune system, NFkB, T cell activation.	FDR < 0.01 for listed enrichments from GO gene sets
52	Glioma (GBM)	MRI	60	0	Myeloid and lymphoid cell surface markers from RNA	13 features (variations of ADC, nCBV, volume, necrosis)	Tumor-associated macrophages (CSF1R), MDSCs (CD33), and helper T cells (CD4) positively correlated with nCBV. MDSCs (CD49d) and T cells (CD3e) anti-correlated with ADC mean.	p < 0.05, Pearson correlation
53	Glioma (low grade)	MRI	47	84	Inflammatory and immune response pathways	431 radiomic features; 9 used in model: 1 first order, 1 texture, 7 wavelet transform features (see study for model coefficients)	High-risk score enriched for Antigen processing/presentation and NFkB. Specifically, GLRLM run length nonuniformity HHL, GLRLM run percentage HLH, and median HLH positively associated with immune response and NFkB (and innate immune) signaling. GLRLM short run low gray level emphasis LLL positively associated with antigen processing/presentation.	Ontology using DAVID pathways on top 200 genes (p < 0.05, Pearson correlation) for each feature
55	Breast (TNBC)	MRI	112	0	TIL level by H&E	BI-RADS and computed features including shape, margin, internal enhancement characteristics, tumor kinetics (initial and delayed patterns), ADC, and tumor roundness	Tumors in the high-TIL group had a more round shape (vs irregular), circumscribed margin, homogeneous enhancement, and absence of multifocality.	p < 0.0001 (shape), p < 0.0001 (margin), p = 0.0003 (enhancement), p = 0.023 (focality), chi-squared test

Legend: NSCLC: non-small cell lung cancer, HNSCC: head and neck squamous cell carcinoma, HCC: hepatocellular carcinoma, BLCA: bladder urothelial carcinoma, GBM: glioblastoma, TNBC: triple- negative breast cancer, TCGA: the cancer genome atlas, TIL: tumor infiltration lymphocytes, IHC: immunohistochemistry, H&E: hematoxylin and eosin, Treg: regulatory T cell, PD-L1: programmed death-ligand 1, CD: cluster of differentiation, MDSC: myeloid derived suppressor cell, NFkB: nuclear factor kappa B, TNF: tumor necrosis factor, GSVA: gene set variation analysis, SUV: standardized uptake value, MTV: metabolic tumor volume, TLG: total lesion glycolysis, GLCM: gray level co-occurrence matrix, GLRLM: gray level run length matrix, SRHGE: short run high gray-level

emphasis, SRLGE: short run low gray-level emphasis, LGRE: low gray-level run emphasis, LRLGE: long run low gray-level emphasis, ROI: region of interest, ADC: apparent diffusion coefficient, NGTDM: neighborhood gray tone difference matrix, GLSZM: gray level size zone matrix, CBV: cerebral blood volume, BI-RADS: breast imaging-reporting and data system, BPE: background parenchymal enhancement, HGRE: high gray-level run emphasis, LRHGE: long run high gray-level emphasis, SRE: short run emphasis, MCC: maximum correlation coefficient, IMC: information measure of correlation, ZLNU: zone length non-uniformity, GSEA: gene set enrichment analysis, DAVID: database for annotation, visualization, and integrated discovery, KEGG: Kyoto encyclopedia of genes and genomes, MWU: Mann-Whitney U

Table 2

Overview of key studies of imaging and radiomic predictors of response to immunotherapies.

Reference	Tumor type	Modality	Primary study (ortraining set)	Validation	Treatment	Clinical endpoint	Features and feature classes assessed	Individual feature findings	Model findings and statistics
27	Lung (NSCLC)	CT	228	0	Single agent anti- PD-1/anti-PD- L1 with/without anti-CTLA4	Time to progression or hyperprogression at 2 months	600 radiomic features used to develop TTP (time-to-progression) and HPD (hyperprogressive disease) models	Univariate analysis showed Tumor 3D laws E5L5E5 significant (p < 0.05) in TTP model; tumor border NGTDM (Neighboring Gray Tone Difference Matrix) strength significant (p < 0.01) in HPD model.	TTP 4-feature model: AUC 0.717 TTP clinical-radiomic model: AUC 0.804HPD 1-feature model: AUC 0.674HPD clinical-radiomic model: AUC 0.865
30	Lung (NSCLC)	PET/CT	109	0	Anti-PD-1/PD- L1	Response based on RECIST 1.1	SUVmax, TMTV	High TMTV, but not SUVmax, associated with shorter OS and absence of disease clinical benefit on multivariate Cox analysis.	p = 0.004, OS; p = 0.045, DCB
31	Lung (NSCLC)	PET/CT	99	95 (47 test, 48 prospective)	Anti-PD-1/PD- L1	Response based on RECIST 1.1	790 radiomic features (including PET, CT, and KLD) and filtered to obtain set of 8 features in multiparametric radiomic signature (mpRS)	Multiparametric radiomic signature (mpRS) predicts PFS, OS, and DCB. DCB predicted by higher textural heterogeneity and convexity and lower SUVmean and HU.	mpRS for DCB: AUC 0.86 (training),0.83 (test), 0.81 (prospective validation)mpRS for PFS/OS: p < 0.01 for all cohorts, log-rank test
33	Lung (NSCLC)	CT	35	24 (external cohort)	Anti-PD-1	PFS and OS	Radiomic features consisting of morphological, histogram and texture parameters (GLCM) used to develop nomogram-based texture score	High volume, entropy, higher GLCM-entropy, higher GLCM-dissimilarity, and lower GLCM- correlation all associated with worseprognosis.	Texture score for low vs high risk predicts PFS (p = 0.03) and OS (p = 0.04) in validation set, log-rank test
35	Lung (NSCLC)	PET/CT	57	0	Anti-PD-1/PD- L1	Response based on RECIST 1.1	SUVmax, MTV, TLG; radiomic features including size, shape, first- order, and second-order features (see study for specific associations)	Higher MTV and TLG associated with stable/progressive disease. High tumor volume, TLG, heterogeneity (e.g. skewness and kurtosis, and six textural features) associated with progressive disease.	p < 0.05 for all associations, Kruskal- Wallis test (see study for specific p- values)
41	Mixed (HNSCC, NSCLC, HCC, BLCA)	CT	135	137 (external cohort)	Anti-PD-1/PD- L1	Response based on RECIST 1.1	84 imaging features for machine- learning trained on CD8 expression; 8 variables were used for the final elastic- net model including 5 features: (-) coefficient for tumor min value, tumor GLRLM_SRHGE; (+) coefficient for ring GLRLM_SRLGE, ring GLRLM_LGRE, ring GLRLM_LRLGE	High radiomics score associated with DCB and better OS.	p = 0.013 (DCB) and p = 0.0081 (OS) byCox model; p = 0.0022 on multivariate analysis
43	Lung (NSCLC)	CT	50	89 (62 cohort 1,27 cohort 2)	Anti-PD-1/PD- L1	Response based on RECIST 1.1	99 texture features x5 statistics computed for each + 24 shape features; 8 most stable and discriminative DelRADx (delta radiomics) features determined after 6–8 weeks of ICI used to train linear discriminant analysis (LDA) classifier (see study for model coefficients)	All 8 top features in DRS (DelRADx risk score) model associated with response. High DRS associated with worse OS. High intranodular mean Haralick entropy and skewness of Haralick correlation and high perinodular skewness of Laws associated with worse OS on univariate analysis.	DRS for response prediction: AUC0.88 (training), 0.85 (validation 1), 0.81 (validation 2)Response associations for individual features in DRS: p < 0.001, Wilcoxon rank-sum testOS associations: p < 0.05, Cox model
44	NSCLC, Melanoma (metastatic)	CT	133	70 (internal cohort)	Anti-PD-1	Response based on RECIST 1.1	5865 radiomic features based on image and transforms; 10 features selected using unsupervised feature selection and machine learning (see study for model coefficients)	Radiomic biomarker performed well for NSCLC and poorly for melanoma at predicting lesion level response, but better as aggregate score for predicting OS in both cancer types. Responding lesions had more irregular patterns (Wavelet HLH GLSZM Zone Entropy), compactness, sphericity (low SVR). For melanoma, response was positively associated with heterogeneity (GLCM Difference Entropy).	For lesion-level response: p < 0.05, Kenward–Roger testFor lung mets: AUC 0.83, p < 0.001, MWU testFor nodal mets: AUC 0.78, p < 0.001, MWU testAggregate OS: p < 0.01, MWU test
68	Melanoma (brain mets)	MRI	88	17	Anti-CTLA4 with/without anti- PD-1/anti-PD- L1	RANO-BM criteria	21 radiomic features, including first and second order texture features, as well as Gabor, Sobel, and Laplacian of Gaussian (LoG) edge features; OS analysis using univariate Cox regression	Higher LoG mean/SD, GLCM entropy, Gabor mean most associated with better OS on univariate analysis. LoG mean positively associated with OS in validation set.	p = 0.001 (univariate analysis, training); p = 0.003 (validation)
69	Melanoma (metastatic)	PET/CT	112	0	Anti-PD-1 or dual anti-PD- 1/anti-CTLA4	True progression vs pseudoprogression based on RECIST 1.1	172 radiomic features per lesion including shape, intensity, and texture; logistic regression models based on blood LDH/S100, volume, radiomics, and Delta radiomics	Delta radiomics models perform better than single time point models. Blood-radiomics model combining LDH level at TP1 (3 mo) and relative change of CT coarseness between TP1 and TP0 (baseline) performed best; higher LDH and larger decrease in CT coarseness indicated lower chance of pseudoprogression.	AUCs: 0.68 (PET-based), 0.69 (CT-based), 0.79 (Delta fractal dimension model), 0.78 (Delta coarseness model),0.82 (blood-radiomics model)
70	Melanoma (mucosal and cutaneous)	PET/CT	56	0	Anti-CTLA4 with/without anti- PD-1/anti-PD- L1	Response based on RECIST 1.1	SUVmax/mean, TMTV, TLG, BLR (bone marrow-to-liver ratio)	For Muc-M: high SUVmax/mean associated with worse OS.For Cut-M: TMTV, TLG and BLR associated with worse OS; TMTV and BLR with worse PFS and disease control.	Muc-M: p = 0.02 (SUVmax), p = 0.03(SUVmean), univariate Cox Cut-M OS: p = 0.009 (TMTV/TLG),p = 0.04 (BLR), multivariate Cox Cut-M PFS: p = 0.004 (TMTV), p = 0.02(BLR), multivariate Cox
71	Melanoma	PET/CT	90	110	Anti-PD-1, Anti- CTLA4, or dual therapy	Not discussed	SUVmax, MTV, SLR (spleen-liver- ratio)	SLR independently associated with survival (high SLR with short OS) after Ipilimumab (but not anti-PD-1); validated in external cohort; lowest MTV quintile associated with better OS.	SLR: p = 0.008 (PFS), p = 0.0002 (OS),Cox modelMTV: p = 0.03 (OS, as continuous variable), Cox model
76	Urothelial (metastatic)	CT	42	21	Anti-PD-1/PD- L1	Response based on RECIST 1.1	49 radiomic features; 5 features in final logistic regression model: (+) coefficient for SD, Entropy, Inverse difference moment, GLRLM_HGRE; (-) coefficient for Cluster tendency; combined model included visceral organ involvement	Low signature value associated with improved disease control and survival; visceral organ involvement associated with poor response; combined model performed best (AUCs shown to the right). Risk of disease progression based on combined model associated with worse PFS/OS.	AUCs (training): 0.87 (response), 0.77 (disease control)AUCs (validation): 0.87 (response),0.88 (disease control)PFS (validation): log-rank p = 0.044 OS (validation): log-rank p = 0.035

Legend: NSCLC: non-small cell lung cancer, HNSCC: head and neck squamous cell carcinoma, HCC: hepatocellular carcinoma, BLCA: bladder urothelial carcinoma, PD-L1: programmed death- ligand 1, PD-1: programmed cell death protein 1, CTLA4: cytotoxic T-lymphocyte-associated protein 4, RECIST: response evaluation criteria in solid tumors, PFS: progression free survival, OS: overall survival, DCB: durable clinical benefit, RANO-BM: response assessment in neuro-oncology brain metastases, SUV: standardized uptake value, (T)MTV: (total) metabolic tumor volume, KLD: Kullback-Leibler divergence, TLG: total lesion glycolysis, GLCM: gray level co-occurrence matrix, GLRLM: gray level run length matrix, SRHGE: short run high gray-level emphasis, SRLGE: short run low gray-level emphasis, LGRE: low gray-level run emphasis, LRLGE: long run low gray-level emphasis, LDH: lactate dehydrogenase, HGRE: high gray-level run emphasis, GLSZM: gray level size zone matrix, MWU: Mann-Whitney U.

The most frequent radiogenomic associations involved gene or pathway analysis, staining for immune cell surface markers, or identification of tumor infiltrating lymphocytes (TILs). Identifying imaging biomarkers of response to immunotherapy was an objective in 18 studies. All except one study (on dendritic cell therapy) focused on immune checkpoint blockade. Table 1, Table 2 summarize key findings, focusing on studies with cohort size 50 or that included validation studies. Additional studies are detailed in Supplemental Tables 1 and 2. Highlighted further are features associated with immune phenotypes in more than one cancer type or by more than one study within a cancer type based on MRI/CT (Fig. 2, Fig. 3) and PET (Fig. 4).

Fig. 2

Heatmap depicting associations between lower order radiomic features on CT/MRI and either immune phenotypes or response to immunotherapy. Included features were reported to have significant associations in > 1 study. Immune: NK: natural killer, TIL: tumor infiltrating lymphocyte, TLR: toll-like receptor, CTLA4: cytotoxic T-lymphocyte-associated protein 4, IL: interleukin, MDSC: myeloid derived suppressor cell, CD: cluster of differentiation, NFKB: necrosis factor kappa B, PD-L1: programmed death-ligand 1, ICI: immune checkpoint inhibitor; Features: ADC: apparent diffusion coefficient, SD: standard deviation.

Fig. 3

Heatmap depicting associations between higher order radiomic features on CT/MRI and either immune phenotypes or response to immunotherapy. Included features were reported to have significant associations in > 1 study. Immune: TNF: tumor necrosis factor, NK: natural killer, TIL: tumor infiltrating lymphocyte, TLR: toll-like receptor, TGFB: transforming growth factor beta, CTLA4: cytotoxic T-lymphocyte-associated protein 4, CD: cluster of differentiation, NFKB: necrosis factor kappa B, PD-L1: programmed death-ligand 1, PD1: programmed cell death protein 1, ICI: immune checkpoint inhibitor; Features: CoLIAGe: co-occurrence of local anisotropic gradient orientations, SD: standard deviation, GLCM: gray level co-occurrence matrix, IMC: information measure of correlation, GLRLM: gray level run length matrix, HGRE: high gray level run emphasis, SRHGE: short run high gray level emphasis, GLSZM: gray level size zone matrix, NGTDM: neighborhood gray tone difference matrix.

Fig. 4

Heatmap depicting associations between imaging features on PET and either immune phenotypes or response to immunotherapy. Included features were reported to have significant associations in > 1 study. Immune: PD-L1: programmed death-ligand 1, PD1: programmed cell death protein 1, PFS: progression free survival, OS: overall survival, CD: cluster of differentiation, ICI: immune checkpoint inhibitor; Features: MTV: metabolic tumor volume, SUV: standardized uptake value, TLG: total lesion glycolysis.

Lung cancer

Numerous immune phenotypes have been associated with radiological features from PET and CT studies of lung cancer. These phenotypes include pathway expression, histological assessment of T cell infiltration and exhaustion (PD-1/PD-L1), as well as response to immune checkpoint inhibition (ICI). In particular, both pre-treatment and post-treatment predictors of ICI response were identified in tumor metabolic studies.

Among tumor and TME related phenotypes, signaling pathways like leukocyte antigen presentation and NF-kB have been shown to associate with both CT radiomic and PET-FDG uptake features [26], [38]. PET and CT studies have also identified associations with T cell infiltration. Notably, two key studies reported that heterogeneous texture features were positively associated with infiltration: of non-exhausted T cells by Tang et al. [39] and type 2 helper T cells by Yoon et al. [32], both of which were validated in external cohorts. However, another large mixed cancer study (including NSCLC) by Sun et al. found that homogeneous tumors with heterogeneous peripheries, rather, were positively associated with T cell infiltration [41]. Given that lung cancers comprised only a quarter of the tumors in this mixed analysis, the result suggests that texture heterogeneity associations may be dependent on cancer type.

Based on PET studies, high SUV features and total lesion glycolysis (TLG) were consistently associated with higher T cell infiltration and PD1/PD-L1 expression [29], [30], [34], [36], [37], [40]. Even though PD1/PD-L1 are often considered signs of T cell exhaustion, high uptake may not be an indicator of a broader immune-inhibitory TME. For instance, high SUV features did not correlate with CD68 (TAM) expression [37]. Overall, heterogeneity textures on CT along with FDG-avidity may be useful for guiding ICI use in lung cancer by predicting the presence of T cells as well as PD-L1 status.

Pre-treatment image features have been associated with response to ICI. At baseline, lower values of PET whole body SUVmax [29], MTV (metabolic tumor volume) [30], [35], TLG [35], and MMVR (metabolic to morphologic volume ratio) [28] were associated with either better tumor response. Notably, Mu et al. were the first to employ a multiparametric radiomic texture model combining PET and CT features to predict tumor response and survival [31]. Along with others, they found low complexity CT features that associated with better survival outcomes, including greater tumor roundness (low convex hull density) and compactness [31], [33], [44], as well as lower volume and entropy [33], [44]. Texture features that associated with improved outcomes included heterogeneity (short run low gray emphasis, short zone emphasis) and irregularity (higher zone entropy and lower gray-level co-occurrence matrix [GLCM] homogeneity) [31], [44], which is concordant with the association between heterogeneous textures and T cell infiltration.

Restaging imaging features taken after ICI treatment were also predictive of response in smaller prospective studies [34], [42], in which patients were re-imaged after single agent PD-1 inhibition. Lower values for SUVmax, MTV, and TLG were associated with better PFS when measured as early as 1 month after treatment [42]. Notably, TLG was independently associated with PFS and OS on multivariate analysis in both studies, whereas histological biomarkers like PD-L1, CD4, and CD8 did not associate with response.

Highlighting the importance of interval feature analysis, Khorrami et al. recently were the first to show that a model built on “delta” radiomics features determined before and after 2–3 cycles of anti-PD-1/PD-L1 was predictive of OS [43]. In particular, these results suggest that changes in uptake features may be predictive of outcomes after ICI, possibly independent of baseline tumor markers. Therefore, more effort should be made to conduct similar prospective work.

Glioma

To date, immunotherapy clinical trials have indicated that the majority of glioblastomas (GBM) respond poorly [78]. This is unsurprising, as GBMs have less T cell infiltration compared to most solid tumors [79], and their TME is generally immunosuppressive [78]. Therefore, there still exists a clinical demand for better biomarkers of response and, in general, a better understanding of the immune TME.

A variety of MRI modalities have been applied toward radiogenomic immune associations. Associated features not only included semantic, first order, and texture features, but also diffusion, perfusion, and kinetic measurements. While the majority of associations were isolated to single studies, there were some commonalities as well as discrepancies. Semantic high contrast enhancement was noted to be positively associated with inflammatory signals from toll-like receptors in GBM and cytokines in low grade glioma (LGG) [46], [51]. In general, volume was negatively associated with a variety of adaptive immune signals [49], [52], which is consistent with larger, more aggressive tumors harboring immunosuppressive environments. Moreover, coarse and heterogeneous textures associated positively with adaptive immunity, including CD3 expression (by RNA and IHC) in GBM [50] and cytokine and IFNy signaling in LGG [51]. However, certain features like semantic necrosis had ambiguous associations. Necrosis was reported to be both positively and negatively associated with inflammatory immune signals from pathways involving cytokines and NF-kB [38], [46].

Significant radiogenomic immune associations have also been determined via diffusion-weighted imaging (DWI) analysis. Mean apparent diffusion coefficient (ADC) correlated negatively with expression of a gene module in which the most enriched pathway was “dendritic cell maturation” [48], as well as markers of both myeloid and lymphoid immune cells [52]. Given that ADC correlates inversely with cellularity in glioma [80], these results suggest that multiple types of immune cells may contribute to the observed hypercellularity in low-ADC tumors.

Changes in DWI parameters may also provide value in predicting immunotherapy response. Recently, Cuccarini et al. found that up to 4 months after immunotherapy via dendritic cell vaccine, an increase in NK cell response was associated with a decrease in normalized minimum ADC (rADCmin) [54]. Interestingly, higher pre-vaccination rADCmin, which is consistent with higher baseline permissiveness to immune cell infiltration, was associated with better PFS and OS [54]. While increased baseline diffusion may be prognostic, ΔADC appears to be the relevant interval biomarker for predicting immunotherapy response. Overall, studies support the conclusion that lower ADC and higher CBV values are associated with immune cell infiltration but are not definitive of pro- versus anti-tumor composition.

Breast cancer

Breast cancer has been a challenging arena for immunotherapy, which has demonstrated limited efficacy in clinical trials to date [81]. Only between 4 and 16% of breast tumors have been reported as highly enriched for lymphocytes [81], which are associated with improved prognosis [82]. Furthermore, the degree of immune infiltration differs significantly depending on breast cancer subtype [83]. This emphasizes the need for adequate patient selection and improved biomarkers that capture tumor immune heterogeneity.

To date, most of the radiogenomic immune associations in breast cancer have been based on MRI features and histological assessment of tumor infiltrating lymphocytes (TILs). Within triple negative breast cancer (TNBC), Ku et al. found that tumors with > 50% TILs by H&E tended to have rounder shapes, a circumscribed margin, and internal homogeneous contrast enhancement [55]. A later study by Celebi et al. confirmed the same shape and enhancement associations [63]. While many of the first order and morphology feature associations were consistent in both TNBC and broader breast cancer cohorts, certain features may be subtype specific. In contrast to the lack of ADC associations in TNBC [55], for less aggressive breast cancer cohorts, mean ADC was positively associated with TIL levels [61], [63].

Breast cancer peritumoral features were also informative of immune phenotypes. Wu et al. reported an association between textural heterogeneity (higher GLCM information measure of correlation) and tumor necrosis factor (TNF) signaling pathways, a finding which they validated in TCGA cohorts, the sole validation effort in breast cancer to date [59]. Peritumoral heterogeneity may also provide valuable information about immune infiltration [58], with TIL density associating with top peritumoral Gabor-filtered features within 0–3 mm of HER2 + tumors [58].

Generally speaking, radiogenomic associations with TIL levels appeared to be consistent across several studies. More importantly, they may be independently valuable for predicting TIL infiltration. For example, the associated textural heterogeneity features identified by Wu et al. were orthogonal to other verified TIL correlates like cytolytic index (a metric of expression of perforin and granzyme A) [60], [84]. Inclusion of cytolytic index into a composite model improved prediction of TIL density as well as prognosis in TNBC cohorts [60].

Liver cancer

There have been a limited number of studies reporting radiogenomic immune associations in liver cancer, with limited overlap between associated features. In 2007, Kuo and colleagues were first to identify associations between immune system and antigen presentation pathways and semantic image features [9]. Liao et al. in 2019 were next to use CT radiomics to evaluate CD8 + T cell infiltration. Within their radiomic model, three GLCM entropy texture features (indicating randomness) and four gray level run features (with both short and long runs) were positively associated with CD8 + T cell infiltration [66], with a higher model score predicting better survival. They further validated this model in an external cohort. More recently, Chen et al. developed and validated a radiomic model of homogeneous texture features to evaluate tumor T cell infiltration [65].

Given the high perfusion and oxygenation of the liver, studies have also used multiparametric MRI (mpMRI) techniques to capture pharmacokinetic features. A prospective study identified several parameters, including oxygenation, blood flow, and diffusion, associated with immune markers [64]. In particular, the presence of TAMs and T cells correlated with lower cellularity (high ADC) and better tissue perfusion (lower relaxation rates), which is consistent with more vascularized and edematous regions of hepatic lesions allowing for increased immune infiltration. In a more recent study, the same group incorporated texture features into their analysis and found several features associated with gene expression of CTLA4 and PD1 [67], though the biological significance of mpMRI textures is understudied.

So far, no radiogenomic correlates of ICI response have been identified for HCC.

Melanoma

Melanoma has been an exciting area of investigation for radiogenomic predictors of response to ICIs, particularly given the success of such therapies in the metastatic setting. FDG uptake features have been consistently associated with ICI treatment outcomes, as previously reviewed [85]. For instance, lower baseline features, such as SUVmax and MTV, were associated with better OS across different types of metastatic melanoma [70], [71]. However, radiomic studies on this topic are a recent development. In patients with brain metastases, less complex lesion borders on MRI (higher mean tumor edge values after Laplacian of Gaussian transform) were associated with better OS [68]. Heterogeneous textures on CT were also predictive of better response to treatment for metastatic lesions [44]. Furthermore, changes at 3 months after ICI to more homogeneous textures (higher delta CT coarseness and lower delta CT fractal dimension) were associated with pseudoprogression [69]. This again emphasizes the importance of longitudinal assessment of radiomic associations with tumor response and underscores the need for prospective studies.

Head and neck cancer

Associations with immune phenotypes have been observed in limited studies via both MR and CT imaging of head and neck cancers. ADC was negatively correlated with CD3 levels by IHC in a small cohort study of oropharyngeal squamous cell carcinoma [72]. Associations have also been reported for HNSCC based on gene expression data. In the largest radiogenomic study to date, Zhu et al. reported that CT features such as volume and textural homogeneity were negatively associated with multiple innate and adaptive immune pathways, whereas the same pathways were positively associated with convexity, sphericity, and heterogeneous texture features [73]. In another recent PET radiomics study, immune system pathways were associated with both shape and texture features [74].

Developing a radiomic model for TILs has been challenging in HNSCC. Sun et al. built a CT-based radiomic model on expression of a CD8 signature in > 10 cancer types, with a predominance of lung and head and neck cancers [41]. In a validation cohort of TCGA tumors, this radiomic score correlated with TILs for several individual cancer types, but not for HNSCC, possibly due to their heterogeneous immune status [41]. While checkpoint inhibition has shown promising results in clinical trials to date [86], [87], the diverse and complicated immune landscape of HNSCC [88] may require larger cohorts to stratify analyses by tumor subtypes when performing predictive modeling of immunotherapy response.

Current challenges and future directions

In recent years, there has been enormous interest in radiomic modeling of tumor immune biology and immunotherapy response. Imaging features have been associated with similar immune phenotypes across studies and cancer types (Fig. 2, Fig. 3, Fig. 4). These associations not only involved first order, shape, and size features but higher order texture features as well. For example, sphericity (roundness) and sharper borders positively associated with immune activity both at a gene expression level, and on histology [55], [63], and predicted response to ICI. Of note, while first-order features recapitulated well across cancers (Fig. 2), texture feature associations did not demonstrate similar consistency. For example, measures of heterogeneity (GLCM information measure of correlation 1, GLRLM run length nonuniformity, GLSZM gray level nonuniformity) associated negatively with expression of immune gene signatures in head and neck cancer [73], whereas the opposite association appeared to be true of breast cancer and glioma [53], [62] (Fig. 3). While this may be due to tumor type-specific biological differences when measuring more complex, textural qualities of tumor images, it also raises concerns about the reproducibility and biological relevance of the calculated textures themselves. In general, PET features associations were consistent between studies both within and between cancers (Fig. 4). FDG uptake features like SUVmax associated positively with expression of PD-L1 in esophageal [77] and lung cancer [29], [30], [34], [36], [40]; SUVmax, MTV and TLG associated with ICI response and survival outcomes in melanoma [70], [71] and lung cancer [30], [34], [35], [42]. Furthermore, many of the radiomic features and models discussed were independent predictive variables for response to immunotherapy, and outperformed established biomarkers like PD-L1 [42], [43]. Overall, radiogenomics studies have yielded many promising results. However, there are numerous issues that the field faces. Whether radiogenomic models will be effective independently or only as an adjunct to histology and genomic profiling will depend on how these challenges, as discussed below, are met in the future.

Tumor cohort size and validation

Among the major limiting factors for radiogenomic associations and modeling are the size of tumor cohorts and the lack of validation. Only 11 of 42 studies of radiogenomic associations and 8 of 18 studies of imaging feature associations with immunotherapy response had validation cohorts (Fig. 5). Median cohort sizes were 56 (IQR: 42.25–113.5) and 84 (IQR: 58.5–161.5) for primary and validation cohorts respectively for radiogenomic studies. Median cohort sizes were 56.5 (IQR: 36.75–106.5) and 79.5 (IQR: 23.25–98.75) for primary and validation cohorts respectively for immunotherapy response studies. Smaller study sizes and lack of external validation increased the risk that many of the reported associations were either false positives or relevant only to the training cohort (poor generalizability). The lack of overlap between most associations analyzed in these studies may be a consequence not only of the limited number of studies specifically focused on immune phenotypes but also the lack of effort to reproduce findings of existing studies (both internally and by outside groups). Issues of reproducibility may improve over time with better access to data via shared resources like the Imaging Data Commons as well as a community effort to make validation a standard part of radiogenomic studies.

Fig. 5

Size distributions of primary and validation cohorts for studies reporting radiogenomic associations and associations between imaging features and immunotherapy response.

Biomarker selection

The types of immune phenotypes associated with radiomic features should also be scrutinized carefully. While many of these associations will help elucidate the biological basis for otherwise complex imaging phenotypes, it is unclear how to interpret certain associations. For instance, the outputs of several studies were radiogenomic heatmaps, in which modules represented associations between subsets of features and genes. In these cases, it was uncertain which modular features were most robust and driving the associations, hence more likely to be validated. Furthermore, not all immune phenotypes have a pre-defined clinical significance in terms of prognosis and response prediction. Even widely validated biomarkers like PD-L1 do not always predict response to immune checkpoint inhibition [89]. Therefore, while radiomic models can predict TIL status and expression of important immune genes and proteins, they might be better used as direct biomarkers [42]. Rather than being limited to use as proxies of existing molecular markers, the future of radiogenomic modeling may lie in its longitudinal applications.

Currently, few biomarkers exist to measure the temporal response of tumors to immunotherapy. Recent studies have demonstrated the importance of radiomic measurements before and after treatment in predicting response to immunotherapy, specifically for checkpoint inhibition and dendritic cell therapy [34], [42], [43], [54]. Changes to FDG uptake features, rADCmin, and perinodular Gabor were shown to be associated with response and prognosis, and ought to be validated in future studies. Delta radiomic features have also provided new insights into pseudoprogression [69]. Moreover, interval radiomic biomarkers have unexplored utility for understanding therapeutic resistance. In tumors that are recalcitrant to therapy, changes to the tumor and surrounding TME may underlie potential mechanisms of resistance. Few studies have captured the biology of these changes as they are occurring, as this would typically entail repeated invasive procedures like tumor biopsies. Non-invasive imaging taken sequentially during the course of treatment may help identify radiomic biomarkers of resistance early on, thereby allowing for real-time decision support and changes to clinical management.

Feature reproducibility

Feature reproducibility has been an ongoing concern in radiomics. For semantic features like “necrosis” and “edema,” variability of subjective assessment can result in opposing observations. For example in GBM, “necrosis” was reported to be both positively [46] and negatively [38] correlated with expression of immune pathways. The automation of radiomic feature extraction bypasses concerns of clinician-reader bias but is not without issues. Variation in imaging acquisition parameters, and pre-processing techniques have shown a significant effect on imaging feature calculations, not only affecting their reproducibility, but also making the features inconsistent within a single dataset. Where appropriate, methods for intensity standardization, such as through referencing of healthy tissue, or more robust methods based on statistical learning, should be considered [90].

Subsequently, after image acquisition and processing, there are two major sources of bias that can lead to lack of uniform reproducibility in radiomic feature extraction. The first is implementation of the mathematical feature definitions to quantify information within ROIs. Concerns have been raised regarding the actual implementation even when the same feature definitions were used [91], which has been evident amongst open source frameworks [92], [93], [94], [95], [96]. Efforts like MITK Phenotyping [97] and the Image Biomarker Standardization Initiative [91] have been developed to standardize test data and centralize these platforms, showing promise for unified feature extraction. The second, and arguably larger source of radiomic features bias, arises from the manual tumor segmentations. While manual contours are often performed by individuals with expert domain knowledge and are often seen as gold-standards for radiomic analysis, they are nonetheless prone to inter- and intra- observer variability [98]. Moreover, tumors can be notoriously difficult to contour due to their unclear borders. This bias should be kept in consideration during radiogenomic analysis interpretation [99], [100]. Ongoing efforts aim to auto-segment the tumor may alleviated these issues [101], but there are concerns over the accuracy of these techniques.

Recently, there has been growing support for using deep learning approaches that mitigate many of the aforementioned biases and reproducibility issues [102]. Deep learning allows the algorithm to define its own features instead of relying on pre-defined features, which can obviate the need for manual segmentation. However, a major draw-back of these deep learning methods is their need for copious training data. Additionally, the underlying meaning behind these self-defined features derived from deep learning algorithms are unclear and the subject of active investigation [103], [104]. There have been several approaches that attempt to bridge the gap between traditional radiomics and deep learning [105], and these may be important studies for the future of radiogenomic approaches that integrate deep learning algorithms into their workflows.

Data analysis and modeling

Current use and integration of large scale “omics” data is primarily retrospective and aimed at hypothesis generation. This underscores the importance of taking a rigorous approach to exploratory analysis and model generation. Some of the goals are to reduce user bias and also to mitigate issues of multiple hypothesis testing (including p-hacking) and overfitting. A study by Chalkidou et al. demonstrated an alarming number of published radiomic studies at the time had high type 1 error probabilities and did not reach statistical significance [106]. Therefore, appropriate statistical methodology, such as correcting for multiple comparisons, is crucial in radiomic analysis. Feature selection is also a key issue for model generation, as radiomics often leads to the creation of high-dimensional feature spaces that can lead to overfitting. Therefore, proper dimensionality reduction techniques should be employed [107].

An increasing number of radiomic and radiogenomic studies are utilizing machine learning in their workflows. While machine learning avoids certain errors common to conventional statistics, it is prone to its own pitfalls. This is particularly true when considering relatively small data sets, on the order of hundreds of training samples or smaller, which forms the bulk of current imaging studies. Improper use of training datasets and data leakage from training to evaluation datasets are factors that have led to overly optimistic results in many studies. It is crucial to use hold-out data sets for evaluation of models and to prevent any data set leakage. Finally, it is highly encouraged to use validation sets of data that are independent from the training set and ideally from multi-institutional sources so as to maximize model generalizability.

As described previously, there is considerable optimism regarding the advent of deep learning in quantitative imaging. Deep learning has made significant strides in the field of genomics and radiomics independently, so it is natural to believe an artificial intelligence framework rooted in deep learning will play a considerable role in radiogenomics as well. However, the key issues of interpretability and overall model transparency have been at the forefront of whether these techniques are suitable for clinical implementation. Fortunately, the deep learning interpretability issue is becoming an increasingly studied problem domain [108]. Solutions can include the analysis of the algorithm itself [109], as well as utilization of novel algorithmic structures that inherently lend themselves to higher levels of interpretability [110]. Therefore, future radiogenomic studies that implement deep learning in their analysis should keep these interpretability considerations in mind.

Building composite models and data integration

One of the key takeaways from current radiogenomic studies is the value of composite models, which combine radiomic models with other covariates, such as clinical and molecular features. Currently, composite models utilizing imaging data are limited by data completeness. Certain efforts have been directed at standardizing multi-omic and imaging data, including development of public data repositories through the support of the NCI Moonshot program. The Imaging Analysis Working Group has compiled high-resolution hematoxylin and eosin (H&E) imaging to allow for quantization of lymphocytic infiltration patterns across 13 of 33 cancer types in TCGA [111]. Progress in machine learning in digital tumor pathology are being accompanied by advances in multiplexed IHC as well as immunostaining approaches like mass cytometry (CyTOF), which provide in-depth spatial characterization of tumor immune composition.

Building better composite models requires integrating the next generation of multi-omic data. In areas where genomic and transcriptomic data do not sufficiently capture tumor biology by linking genotype to tumor phenotype, these next generation “omic” approaches frequently provide better insight [112], [113], [114], [115]. The integration of these large, disparate data sources will necessitate better software and workflow management. Platforms like MultiAssayExperiment have been developed with the express purpose of providing data objects and structures for this type of integrative analysis [116].

Data transparency, reporting and best practices

To foster better study design and data transparency, there has been a recent trend toward pre-registering studies in databases such as the Open Science Framework (OSF) to promote data use integrity and encourage standard practices [117]. Additionally, the utilization of curated public repositories for open-sourcing of computational analysis tools (e.g. GitHub) and datasets (e.g. NCI Cancer Research Data Commons) [118] will further foster transparency and reproducibility of radiogenomic studies.

Radiogenomics has much to gain in consolidating and standardizing methods of analysis in an effort to accurately compare studies, as has been previously described for quantitative imaging biomarkers [119]. Success in this arena will depend on developing good practice guidelines. Some of the factors that ought to be considered as part of good practices have been discussed previously and agglomerated into a radiomics quality score (RQS) [1], [4]. The reporting of radiogenomic studies aimed at model development and validation can benefit from broader guidelines, such as those recommended by the TRIPOD group [120].

Moving forward, the utility of radiomic features and radiogenomic models with regards to tumor immunology will continue to be twofold: 1) predicting response to immunotherapy and 2) comprehending and modeling immune biology. Based on these considerations and the major challenges discussed above, our recommendations for best practice guidelines for future studies are summarized in (Table 3) [1], [120], [121], [122].

Table 3

Recommendations for conducting and reporting studies that investigate radiogenomic associations with tumor immune phenotypes.

Process	Considerations	Recommendations
Study design	Study registration	Pre-register studies in databases such as the Open Science Framework (OSF)
	Cohort selection	Focus on specific molecular subtypes or subclasses of cancers may enable more accurate radiogenomic modelsMeta-analysis of multiple cohorts can be used to achieve more generalizable models
	Study design	Prospective study design to enable longitudinal feature assessment may be ideal for generating models to predict immunotherapy response and identify biomarkers of resistanceFor retrospective study design, statistical and modeling approaches should be decided a priori
Evaluating molecular data	Tumor and TME gene expression data procurement and processing	RNA-seq for assessing gene expression, refer to Conesa et al. 2016 for a review of good data practices [121]RNA-seq may be eventually supplanted by single-cell RNA-seq, which can improve the ability todistinguish tumor versus immune cell gene expression
	Pathway and immune infiltration analysis	Software like Gene Set Enrichment Analysis (GSEA), Ingenuity Pathway Analysis, DAVID, Metascape are standard for pathway enrichment analysisApproaches including single sample GSEA (ssGSEA), CIBERSORT, and Immunoscore useful for more specific quantification of types of tumor immune cell infiltration
	Cell markers by IHC	Specific staining of cell surface markers remains the gold standard for quantifying immune cell infiltrationTo increase staining throughput, consider using tissue microarrays and multiplexed IHC
	Quantifying TILs by H&E	H&E allows for good quantitation of TILs, but is often subject to clinician-reader biasBest clinical practices are outlined in Salgado et al. 2015 [122]
Image acquisition, processing, and extraction	Image acquisition parameters	Use standardized acquisition parameters
	Image pre-processing	Normalize voxel intensities of images, particularly MRI, to more accurately and reproducibly extractradiomic features
	Feature definition and extraction	Use feature standardization platforms, such as MITK Phenotyping and the Image BiomarkerStandardization Initiative
	Tumor segmentation	Use multiple independent observers if segmenting manually or consider semi-automatic/automaticapproaches to maximize reproducibility
	Deep learning	Utilize algorithm visualization methodology, such as saliency maps, to increaseinterpretability/explainability/transparency
Modeling and data analysis	Feature selection	Reduce feature dimensionality such as through regression modeling (e.g. LASSO Cox, Elastic Net) or using intra-class feature similarity measures (e.g. intra-class correlation coeffcient) to preventoverfitting and improve feature reliability
	Model design	Best performing models for predicting prognosis and immunotherapy response are likely achieved by combining radiogenomics models with other covariates into composite modelsCorrect for multiple hypothesis testing where appropriate
	Machine learning	Use hold-out data sets for evaluation of models and to prevent any data leakage from training to evaluation setsValidate on data that are independent from the training set and ideally from multi-institutional sources
Data transparency and reporting	Public data and code repositories	Share code in open-source repositories like GitHubShare imaging data in public repositories like the Imaging Data Commons (IDC)
	Radiomics quality score (RQS)	Report RQS score (out of 36) developed by Sanduleanu et al. 2018 [1]
	Study reporting checklists	Use of TRIPOD 22-item checklist for model development and validation [120]

Legend: DAVID: database for annotation, visualization, and integrated discovery, IHC: immunohistochemistry, TIL: tumor infiltrating lymphocyte, H&E: hematoxylin and eosin, MITK: medical imaging interaction toolkit, LASSO: least absolute shrinkage and selection operator, TRIPOD: Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis.

Conclusions

The goal of developing more robust prediction models from radiomics and other tumor properties is to ultimately provide better personalized care for patients. As more targeted therapies and immunotherapies emerge from clinical trials, the need for strategies to augment clinical decision-making will increase. Precision medicine stands to benefit tremendously from more accurate characterization of the 3D tumor immune microenvironment, identifying prognostic biomarkers for immunotherapy, and providing clinical guidance in the setting of immunotherapy resistance. Based on studies to date, there is growing evidence that tumor immune states can be characterized by radiomic features alone or in combination with other molecular and histological features. More so, baseline imaging features as well as changes in features appear to provide a new source of predictive power for immunotherapy response. However, current radiomic findings are generally not mature enough to serve as surrogate predictors of immune biology, as few radiogenomic associations have been thoroughly validated. Furthermore, while studies can recapitulate immune associations with low complexity features, associations with texture features are more equivocal. Overall, radiogenomic modeling applied to tumor immune biology is at its nascent stages. Development of clinically relevant radiogenomic models will need to be accompanied by a roadmap for technical and biological validation, as well as thoughtful integration of other biological (“omic”) and clinical data to develop more accurate, composite models.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.