Genetic alterations associated with 18F-fluorodeoxyglucose positron emission tomography/computed tomography in head and neck squamous cell carcinoma.

BACKGROUND: We investigated the relationship between genetic alterations and 18F-FDG PET/CT findings in head and neck squamous cell carcinoma (HNSC).
METHODS: Using mRNA-sequences of HNSC samples (480 patients) from the Cancer Genome Atlas (TCGA) portal, gene coexpression networks were constructed via a weighted correlation network analysis (WGCNA) algorithm, and their association with the tumor-to-blood signal ratio on 18F-FDG PET/CT data (21 patients) was explored. An elastic-net regression model was developed to estimate the PET tumor-to-blood ratio from the gene networks and to derive an FDG signature score (FDGSS). The FDGSS was evaluated with regard to clinical variables and general mutational profiles, as well as alterations to oncogenic signaling pathways.
FINDINGS: The FDGSS values differed across clinical stages (p = 0.027), HPV-status (p< 0.001), and molecular subtypes of HNSC (p< 0.001). Multivariate Cox regression demonstrated that FDGSS was an independent predictor for overall (p = 0.019) and progression-free survival (p = 0.024). FDGSS positively correlated with total mutation rate (p = 0.016), aneuploidy (p < 0.001), and somatic copy number alteration scores (p < 0.001). CDKN2A in the cell cycle pathway (q = 0.014) and the TP53 gene in the TP53 pathway (q = 0.005) showed significant differences between high and low FDGSS patients.
CONCLUSION: FDGSS based on the gene coexpression network was associated with the mutational landscape of HNSC. 18F-FDG PET/CT is therefore a valuable tool for the in vivo imaging of these cancers, being able to visualize the glucose metabolism of the tumor and allow inferences to be made on the underlying genetic alterations in the tumor.

Introduction

Genetic alterations are the hallmark of cancer and enable malignant proliferation, angiogenesis, invasion, and resistance to cell death [1]. The differential expression of genes in tumor cells affects their phenotype and ultimately their diverse characteristics such as cancer aggressiveness, treatment response, and patient prognosis [2]. As targeted drugs for specific cancer driver genes or molecular pathways have been investigated and gradually introduced into the clinic, assessments of the genetic alteration status of tumors have become more important for therapeutic management [3]. For example, the addition of epidermal growth factor receptor inhibitors to chemotherapy or radiation therapy regimens has yielded survival benefits in randomized phase III trials [4,5]. Moreover, novel approaches targeting other oncogenic signaling pathways for head and neck squamous cell carcinoma (HNSC) have also been widely investigated [6]

Glycolysis, by which cancer cells convert glucose to lactate in an oxygen-independent manner to generate ATP, instead of performing oxidative phosphorylation, is a hallmark of the metabolic phenotype of cancer. It involves intermediates required for biosynthesizing nucleotides and fatty acids, and reducing agents [7]. A body of evidence indicates that the alteration of multiple signaling pathways resulting from oncogene and tumor suppressor gene mutations can be attributed to this metabolic reprogramming [8,9]. Recent studies reported that this metabolic switch reciprocally facilitates the epigenetic regulation of gene expression which in turn contributes to tumorigenesis [10]. 18F-fluorodeoxyglucose positron emission tomography/computed tomography (18F-FDG PET/CT) enables the non-invasive assessment of increased glucose metabolism in a tumor, and is one of the standard imaging modalities for evaluating HNSC [11], with the semi-quantitative parameters of FDG uptake such as maximum standardized uptake value (SUVmax) known to have significant prognostic value [12]. Accordingly, 18F-FDG PET/CT parameters are highly likely to be associated with alterations to genes and oncogenic signaling pathways in HNSC lesions. A comprehensive understanding of the link between 18F-FDG PET/CT and genetic modifications in cancer cells will therefore enhance our understanding of the biology of 18F-FDG PET/CT and cancer metabolism, knowledge of this relationship is currently lacking. Therefore, in this study, we evaluated such associations using HNSC samples from the Cancer Genome Atlas (TCGA) projects.

Materials and methods

mRNA data from the TCGA-HNSC dataset

Level three RNA-sequencing data for HNSC patient samples were downloaded from the publicly-available TCGA data portal (https://portal.gdc.cancer.gov/). These data contained 20 531 genes from 520 primary tumors, obtained using an Illumina HiSeq RNASeqV2 (Illumina, San Diego, CA, USA). The raw read counts were normalized and genes with low signal intensity were filtered out with a quantile cut-off of 0.50 using the ‘TCGAbiolinks’ R package [13]. Clinical data for the patients were obtained via the ‘TCGAbiolinks’ R package and the cBioPortal website (https://www.cbioportal.org/). Two patients were excluded as they did not have clinical follow-up data available for the determination of overall survival (OS) or progression-free survival (PFS). A further 38 patients with a history of prior treatment or with a previous synchronous malignancy were also excluded. A final cohort of 480 patients with 10 165 mRNA transcripts and available clinical data was analyzed.

PET/CT data and analysis

From the Cancer Imaging Archive (TCIA, https://www.cancerimagingarchive.net/) website, we identified 29 patients in our TCGA-HNSC dataset for whom 18F-FDG PET/CT data were also available. Of these 29 patients, three with post-operative PET/CT and one for whom the Digital Imaging and Communications in Medicine (DICOM) files were not feasible for semi-quantitative analysis were excluded. An experienced nuclear medicine physician working in a blind manner extracted the maximum standardized uptake value (SUVmax) of the primary tumor using Metavol software (https://www.metavol.org/). The DICOM files on the TCIA website were acquired using a variety of PET/CT scanners and reconstruction algorithms; therefore, to compensate for multi-scanner issues, cases with a metabolic tumor volume equal to or higher than 4.2 ml according to a threshold of 40% of the SUVmax were retrieved, and the maximum tumor-to-blood ratio (TBRmax) was calculated as the tumor SUVmax divided by the SUVmean of the blood pool [14,15]. Finally, the TBRmax values on 18F-FDG PET/CT were obtained for 21 patients.

In these 21 included patients, FDG was injected at a median dose of 514 MBq (range, 307–784 MBq). PET/CT images were acquired 1–2 h after administration of tracer (median, 77 min). The PET images were heterogeneous in terms of scanner modality, acquisition, and reconstruction protocol, as they were mostly acquired as part of routine care (not as part of a controlled clinical trial). Nevertheless, all the PET/CT cameras (Discovery RX [n = 2], LS [n = 4], ST [n = 6], and STE [n = 7]; GE Healthcare, Biograph Truepoint [n = 2]; Siemens Healthineers) were equipped with ordered subset maximum likelihood (OSEM) iterative reconstruction; although the reconstruction parameters differed across the institutions that supplied the data, most of the iteration, subset, and post-smoothing parameters used were the standard settings for each camera generation (e.g., two iterations, 15 subsets, 6 mm full-width-at-half-maximum [FWHM] Gaussian smoothing). Moreover, most of the scanners were without time-of-flight and point-spread-function modeling, resulting in a voxel size of approximately 4 mm and 7–8 mm FWHM spatial resolution. Therefore, we believe that the spatial resolution and signal-to-noise ratio did not differ substantially according to the different scanners on which the PET/CT data were acquired.

Exploration of TBRmax-associated gene networks by weighted correlation network analysis

To explore the gene networks that correlated with TBRmax, weighted gene correlation networks were generated using our gene expression data after log2 transformation with the ‘WGCNA’ R package [16]. In bioinformatics, weighted correlation network analysis (WGCNA) has been widely applied to find gene modules (clusters of highly correlated genes) and summarize clusters by means of the module eigengene or an intramodular hub gene for relating modules to external traits (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/). A ‘signed network adjacency’ matrix was first constructed from the gene expression data from the full study cohort of 480 patients using a soft thresholding power of ‘6′, which was selected to ensure scale-free topology and provide sufficient node connectivity. The adjacency matrix was then transformed into a ‘topological overlap matrix (TOM)’ to minimize the effects of noise and spurious associations, and the corresponding dissimilarity was calculated as ‘1 – TOM’. Hierarchical clustering of the dissimilarity (1 – TOM) was performed using an ‘average’-linkage method to produce a clustering tree of genes; the branches of the clustered groups of the genes were highly interconnected. The gene network modules were identified by cutting the branches off the clustering tree using the ‘dynamicTreeCut’ R package, setting a cut-height value of ‘0.99′, a deep split of ‘2′, and a minimum module size of ‘50′ [15,16].

For the 21 patients with PET/CT data, module eigengenes, defined as the first principal component scores of the gene expression data in a given gene network module, were calculated. A Spearman correlation test was then performed between the TBRmax and module eigengenes from each gene network module. Modules with p-values < 0.05 were defined as TBRmax-associated gene networks. Hub genes of TBRmax-associated gene networks were identified according to the Spearman correlation coefficient between the module eigengenes and the expression profiles of individual genes in the respective modules [15]. The functional profiles of TBRmax-associated gene networks were retrieved by gene ontology (GO) analyses on biological process, pathway, and molecular function using the ‘TCGAanalyze_EAcomplete’ function in the TCGAbiolinks R package.

Elastic-net model for tumor metabolism prediction

A statistical model using elastic-net regression was built to predict the TBRmax values from the transcriptome data for the 21 study patients with available PET/CT data. The training data were log2 transformed mRNA transcripts for the top 10 hub genes from significantly correlated modules, which were preprocessed by subtracting the means and dividing by the standard deviation (centered and scaled). As the number of predictor variables (top 10 hub genes from the TBRmax-associated gene networks) was relatively large compared with the number of outcome variables (TBRmax values for 21 patients), elastic-net regression, a convex combination of ridge and lasso regularized regression, was chosen to deal with overfitting problems. The elastic-net model was fitted using the ‘glmnet’ and ‘caret’ R packages to simultaneously tune alpha (mixing percentage between ridge and lasso) and lambda (regularization parameter) values via leave-one-out cross-validation (LOOCV). Root mean square error (RMSE) was used as the loss function for the elastic-net model. After training, the model was applied to the full 480-patient TCGA-HNSC dataset. The output of the model was defined as the FDG signature score (FDGSS).

The patients’ FDGSS values were compared in terms of human papillomavirus (HPV)-positivity, molecular subtype, clinical stage, and pathologic grade using either two-sample t-tests or an analysis of variance (ANOVA) test. Univariate and multivariate Cox regression analyses of OS and PFS with the baseline characteristics (age at diagnosis, clinical stage, and histologic grade) were performed. The median FDGSS value was used to stratify the patients into high and low FDGSS groups. Kaplan–Meier curves were generated to compare the OS and PFS between these two groups, and log-rank tests were used to evaluate the statistical significance of any differences.

Landscape of the genetic alterations with regard to the FDGSS

The precompiled Mutation Annotation Format (MAF) files containing somatic variant data for the TCGA-HNSC dataset obtained from the National Cancer Institute were downloaded and analyzed using the ‘TCGAmutations’ and ‘maftools’ R packages, respectively [17]. Among the 507 patients initially retrieved, we identified 462 patients that matched with our mRNA data. The total number of somatic mutations and non-silent mutations per megabase (Mb) were obtained for each sample. Mutant-allele tumor heterogeneity (MATH) scores, defined as the median absolute deviation of the mutant-allele fraction divided by the median mutant-allele fraction, were derived to infer intra-tumor genetic heterogeneity via the methodology proposed by Mroz et al. [18]. Scores for the somatic copy number alteration (SCNA) burden, including aneuploidy scores (the total sum of chromosome arm-level alteration), the fraction altered (the fraction of bases pairs present in the copy number profiles deviating from the baseline ploidy), and the number of segments (the total number of segments present in the copy number profile of each patient), were also derived using previous reports on TCGA pan-cancer analysis [19], [20], [21]. The Pearson correlation test was used to evaluate associations between the FDGSS and mutational scores. The total and non-silent mutation rates and the number of segments were log-transformed because of their right-skewed distributions.

Genetic alterations in the oncogenic signaling pathway were evaluated with regard to the FDGSS. Specifically, 10 canonical signaling pathways and their member genes, which are possible key pathways related to cancer drivers or potential therapeutic targets, were assessed via comprehensive analyses of somatic mutations, copy number alterations, epigenetic silencing, and fusions and structural rearrangements of member genes showing frequent genetic alterations [22]. The alteration status was obtained at both the pathway and gene levels: FDGSS values were compared between samples with altered or intact oncogenic signaling pathways using a two-sample t-test. The frequency of alterations in the member genes of significant pathways was compared between high and low FDGSS groups using Fisher's exact test. All p-values from multiple comparisons were adjusted using false discovery rate (FDR) to compensate for type I errors in the null hypothesis. P-values or FDR-adjusted p-values (q-values) of < 0.05 were defined as statistically significant. All statistical analyses were performed using R software (version 3.6.0; R Foundation for Statistical Computing, Vienna, Austria).

Results

FDGSS estimation based on TBRmax-associated gene networks

The characteristics of the included patients are presented in Table 1. WGCNA analysis identified 14 coexpression gene network modules (Fig. 1a). Spearman correlation tests between module eigengenes and the TBRmax values demonstrated that blue, magenta, and turquoise modules were TBRmax-associated gene networks (Fig. 1b–e; p = 0.013, rho = –0.54 for blue; p = 0.034, rho = –0.47 for magenta; and p = 0.031, rho = 0.47 for turquoise). The functional profiles of the blue, magenta, and turquoise modules were annotated using GO and included “G-protein coupled receptor pathway”, “immune reaction”, and “translation” (Supplementary file: Table S1). We selected the top 10 hub genes from each module to build a prediction model for TBRmax. The correlations between the log2-transformed expression profiles of the hub genes and the FDGSS are shown in Fig. 1f–h.

Table 1

Clinical characteristics of all patients (n = 480) and patients with PET/CT available (n = 21).

Characteristics	TCGA-HNSC (n = 480)	PET/CT (n = 21)
Gender (F:M)	123:357 (25.6%:74.4%)	4/17 (19.0%:81.0%)
Age at diagnosis (y)	60.9 ± 11.7 (19–90)	57.0 ± 14.0 (26–85)
Clinical stage
1	17 (3.6%)	2 (9.5%)
2	88 (18.8%)	4 (19.0%)
3	101 (21.6%)	3 (14.3%)
4 (a/b/c)	261 (246/9/6, 55.9%)	12 (12/0/0, 57.1%)
Tumor site
Hypopharynx	8 (1.7%)	0 (0.0%)
Larynx	106 (22.1%)	5 (23.8%)
Oral cavity	106 (22.1%)	0 (0.0%)
Oropharynx	50 (10.4%)	6 (28.6%)
Overlapping lesion	62 (12.9%)	3 (14.3%)
Tongue	142 (29.7%)	5 (23.8%)
Others*	6 (1/1/3/1, 1.3%)	2 (1/1/0/0, 9.5%)
Pathologic grade†
1	60 (12.9%)	0 (0.0%)
2	285 (61.3%)	13 (61.9%)
3	115 (24.7%)	8 (38.1%)
4	5 (1.1%)	0 (0.0%)
Race‡
American Indian or Alaska native	2 (0.4%)	0 (0.0%)
Asian	11 (2.4%)	0 (0.0%)
Black or African American	46 (9.8%)	3 (14.3%)
White	409 (87.4%)	18 (85.7%)

Data are expressed as number (proportion) or mean ± standard deviation (range).

*Others include palate (not specified), pharynx (not specified), lips, and mandible.

†,‡Data were available for 465 and 468 patients, respectively.

Fig. 1

Gene clustering dendrogram generated using WGCNA. The bars represent corresponding coexpression network module colors consisting of highly interconnected genes in the 480 HNSC study samples (a). Bar plot representing the p-values from correlation tests between the module eigengenes and TBRmax. The dotted line denotes a p-value of 0.05 (b). Scatter plots showing associations between eigengenes of significantly correlated modules and TBRmax (c–e). Bar plots illustrating the correlation coefficients between the log2 transformed expression values of the top 10 hubgenes and TBRmax (f–h).

Elastic-net implementation with LOOCV for TBRmax estimation revealed that the lowest RMSE of the predicted TBRmax was 3.30 at alpha and lambda values of 0.325 and 1.61, respectively (Supplementary file: Fig. S1). The FDGSS, which is predicted via the elastic-net model, showed moderate correlation with TBRmax in the 21 patients (Fig. 2a; p < 0.001, rho = 0.69). The FDGSS was then calculated for all 480 patients from the TCGA-HNSC datasets included in this study. The FDGSS values were significantly higher in the HPV-negative tumors than in the HPV-positive tumors (Fig. 2b; p < 0.001). FDGSS also varied significantly across different clinical stages (Fig. 2c; p = 0.027) and among the four molecular subtypes of HNSC (Fig. 2d; p < 0.001). However, FDGSS did not differ across the pathologic grades of tumor (p = 0.831). Cox regression analyses further demonstrated that FDGSS is an independent predictor of OS (Table 2; adjusted HR, 1.12; p = 0.019) and PFS (adjusted HR, 1.12; p = 0.024). Kaplan–Meier curves indicated that the high-FDGSS group of HNSC cases showed significantly shorter OS and PFS (Fig. 3).

Fig. 2

Scatter plots showing correlations between FDGSS and TBRmax in the 21 study patients with available 18F-FDG PET/CT data (a). Boxplots comparing FDGSS with HPV-positivity (b), clinical stage (c), and molecular subtype (d) of the HNSC tumors.

Table 2

Cox proportional hazard model for overall and progression-free survival outcomes in the study patients.

		Univariate analysis		Multivariate analysis
Outcomes	Variables	HR (95% CI)	p	HR (95% CI)	p
Overall survival	Age	1.02 (1.01–1.04)	<0.001	1.02 (1.01–1.04)	<0.001
	Clinical stage	1.11 (0.94–1.30)	0.209
	Histologic grade	1.01 (0.82–1.24)	0.920
	FDGSS	1.14 (1.03–1.25)	0.009	1.12 (1.02–1.23)	0.019
Progression-free survival	Age	1.01 (0.99–1.02)	0.356
	Clinical stage	1.18 (0.99–1.41	0.061	1.16 (0.97–1.38)	0.100
	Histologic grade	1.08 (0.87–1.35)	0.488
	FDGSS	1.14 (1.03–1.26)	0.009	1.12 (1.02–1.25)	0.024

All variables were regarded as continuous variables.

CI = confidence interval; HR = hazard ratio; NA = not applicable.

Fig. 3

Kaplan–Meier curves comparing overall survival (a) and progression-free survival (b) between patients stratified into high and low groups according to their median FDGSS value.

Genetic and pathway alteration signature of the FDGSS

A general overview of the somatic mutational profiles of our current study samples with high and low FDGSS is illustrated in Supplementary file: Figure S2. The log-transformed total and non-silent mutation rates showed a positive correlation with FDGSS (Fig. 4a–b; r = 0.11, p = 0.016 and r = 0.11, p = 0.023, respectively). The MATH scores were also significantly correlated with FDGSS (Fig. 4c; r = 0.10, p = 0.027). Regarding the scores for SCNA burden, the FDGSS values were significantly correlated with the aneuploidy scores (Fig. 4d; r = 0.19, p < 0.001), log-transformed number of segments (Fig. 4e; r = 0.38, p < 0.001), and the fraction altered (Fig. 4f; r = 0.29, p < 0.001).

Fig. 4

Scatter plots of correlations between FDGSS and the log-transformed total (a) and non-silent mutation rates (b), MATH score (c), aneuploidy score (d), log-transformed number of segments (e), and fraction altered (f).

The FDGSS was next evaluated in relation to alterations in the oncogenic signaling pathways (Fig. 5a). Among the 10 pathways assessed, a significantly higher FDGSS was indicated in samples with altered cell cycle (q = 0.014), TP53 (q = 0.006), and TGF-beta (q = 0.024) pathways than in samples with these pathways intact. The average numbers of alterations in the member genes of each pathway are shown according to high and low FDGSS in Fig. 5b. Specifically, the CDKN2A (cell cycle pathway) and TP53 (TP53 pathway) tumor suppressor genes were differentially altered in patients with a high FDGSS in comparison with cases showing a low FDGSS (Fig. 5c; 73% vs 59%, q = 0.014, and 78% vs 64%, q = 0.005, respectively). Regarding the TGF-beta pathway, the TGFBR2 and SMAD4 genes showed marginally significant alterations according to the unadjusted p-values, but these were not statistically significant following FDR adjustment (q = 0.296 for both). In addition, the MGA gene in the MYC pathway showed a significantly higher frequency of alterations in the high-FDGSS group (12% vs 4%, q = 0.028).

Fig. 5

Clustered boxplot showing a comparison of FDGSS with 10 oncogenic signaling pathway alterations; * and ** represent q-values of < 0.05 and < 0.01, respectively (a). Clustered scatter plot representing the average number of altered genes of each pathway in high and low FDGSS samples (b). Comprehensive heatmap illustrating member gene alterations in cell cycle, TP53, and TGF-beta pathways, in accordance with the FDGSS group. Each gene is described by its pathway and role in tumor development (right of the heatmap). The frequency of alterations in each gene was compared between the high and low FDGSS groups (c).

Discussion

The major findings of our study are as follows: 1) TBRmax is associated with gene networks that function in G-protein coupled receptor signaling, immune responses, and protein translation; 2) the FDGSS, which is calculated on the basis of these gene networks, is significantly associated with the molecular subtype and clinical stage of HNSC and patient survival outcomes; 3) the FDGSS is positively correlated with the mutational burden and genetic heterogeneity of HNSC; and 4) specific alterations in the cell cycle, TP53, and TGF-beta pathways are significantly associated with a high FDGSS.

Glycolysis is promoted by cell membrane receptor tyrosine kinase and G-protein coupled receptors via several intrinsic pathways, and is upregulated by the increased gene transcription and protein translation of glucose transporters and glycolytic enzymes [23,24]. The immune reactions primarily mediated by T-lymphocytes are closely related to the glucose levels in the tumor microenvironment, as glucose is required for the activation, proliferation, and differentiation of T cells [25]. HPV-negative HNSCs exhibit a higher level of glycolysis and a higher SUVmax than the HPV-positive forms [26,27]. A previous study by Ottensmeier et al. indicated that the gene network that negatively correlates with glucose metabolism is highly active in atypical and mesenchymal subtypes of HNSC compared with the basal and classical subtypes [28]. It is also well known that a high SUVmax is indicative of a poorer prognosis in HNSC [12]. These earlier findings support our current results that show an association between FDGSS and clinical variables in our HNSC study population.

Our analyses of HNSC cases also found that FDGSS was positively correlated with the genetic instability of tumor cells, as represented by the total and non-silent mutation rates, SCNA scores, and genetic heterogeneity. The total mutational load is known to be closely related to the cancer incidence rate, its response to treatment, and the patient's prognosis [29,30], and also to the metabolic reprograming characterized by an increase in carbohydrate metabolism, which involves comprehensive energy metabolism pathways including glycolysis, pentose phosphate pathways, gluconeogenesis, and glycogen metabolism [31]. The MATH score is a quantitative measure of intra-tumoral genetic heterogeneity reflecting the fraction of the whole-exome sequence showing mutant alleles, and is associated with the prognosis of HNSC patients [18] and the risk of metastasis in colon cancer [32]. A previous study of TCGA-HNSC samples reported that MATH scores were strongly associated with metabolic parameters on 18F-FDG PET/CT, including the metabolic tumor volume and total lesion glycolysis [33]. Aneuploidy and focal copy number alterations are considered to be different classes of SCNA and arise via different underlying mechanisms [19,34]. An abnormal number of chromosomes is associated with a higher pathologic grade and poor clinical outcomes [35]. Chromosomal aberrations cause a wide range of gene alterations associated with metabolic disruption, a decrease in mitochondrial activity, and an elevated reactive oxygen species level, which in turn contributes to the Warburg effect. Changes in metabolism also increase the error rate in mitosis, causing genetic instability and vice versa [7,36].

The genetic instability that we found to be related to increased FDG uptake may cause dysregulation of canonical oncogenes and signaling pathways, leading to tumorigenesis. We found that a high FDG uptake is related to alterations in the TP53 gene in the TP53 pathway, and to CDKN2A in the cell cycle pathway. TP53 is a key tumor suppressor gene that downregulates glycolytic activity via 1) repression of GLUT1 and GLUT4 expression, 2) upregulation of the TIGAR gene and subsequent lowering of the fructose-2,6-bisphosphate levels in cells, and 3) activation of SCO2 gene expression to promote oxidative phosphorylation [7,8]. The CDKN2A gene encodes the p14ARF protein that is translated via alternative splicing and that stabilizes p53 and inhibits the cell cycle in a TP53-dependent manner [37]. Alterations in TP53 and CDKN2A thus co-occur, as shown in Figure S2. The PI3K/AKT/mTOR pathway is known to play a pivotal role in the Warburg effect via the activation of glucose transporter and hexokinase. However, we found no significant association between FDGSS and PI3K pathway alteration in our current experiments; we speculate that the much lower frequency of alterations in genes of this pathway (including AKT1–3 [0–2%], MTOR [0%], and PTEN [5%]) in our study population may have limited the statistical power of our observations.

This study had several limitations of note. First, our data includes only a small number of patients with available 18F-FDG PET/CT, and these had heterogeneous tumor origins and HPV subtypes. Our FDGSS calculations may have been subject to confounders or overfitted to the small sample because of the lack of external validation, even though we demonstrated that the clinical and prognostic significance of FDGSS in the whole population was comparable to that of the TBRmax reported in previous studies. The TCGA-HNSC dataset is the only publicly-available data source that allows for comprehensive analyses of 18F-FDG PET/CT, gene expression, and specific gene and pathway alterations in HNSC. We hope our preliminary findings in humans will be further examined in large clinical studies in the near future. Second, the scale of the FDGSSwas different from, and smaller than, that of TBRmax. Hence, our FDGSS data should not be directly interpreted or applied in the same manner as TBRmax. This limitation may have also decreased the difference between tumors with high and low FDG uptake, and may have therefore caused an underestimation of the statistical significance of our comparison between FDGSS and genetic alterations determined by Pearson's correlations and two-sample t-tests. Third, the PET/CT images from the TCIA data portal were obtained using various acquisition protocols, devices, and reconstruction algorithms, all of which can affect SUVmax. We selected patients with a measurable tumor size to minimize the impact that the resolution differences between the PET scanners would have on SUVmax, and divided SUVmax by the mean blood pool activity to compensate for differences in imaging protocols. However, this will not have completely removed the variability. We speculate that this is one possible explanation for the moderate degree of correlation between FDGSS and TBRmax.

Conclusion

Our predicted TBRmax values (FDGSS), based on the gene coexpression network and elastic-net model, were significantly associated with the clinical characteristics and survival of TCGA-HNSC patients. FDGSS showed positive correlations with total and non-silent mutation rates, genetic heterogeneity, and copy number alterations. A high FDGSS was related to a higher rate of alteration in the cell cycle pathway via the CDKN2A suppressor gene, and in the TP53 oncogenic signaling pathway via the TP53 tumor suppressor gene. The findings from our present comprehensive analysis may support the biologic significance of FDG uptake in HNSC from the perspective of gene and pathway alterations.

Declaration of Competing Interest

The authors declare no conflicts of interest in relation to this study.