Predictive models for stage and risk classification in head and neck squamous cell carcinoma (HNSCC)

Data retrieval and pre-processing

The normalized reads of miRNAs (Illumina HiSeq 2000 platform, miRgene-level, Normalized), mRNAs (llumina HiSeq 2000 platform, Gene-level) expression were retrieved from the TCGA (https://cancergenome.nih.gov). A total of 528 samples (484 primary tumors and 44 Normal Adjacent Tissue (NAT)) with associated clinical information such as TNM stage, status, gender, and age etc. (Table S1) were selected whereas the samples with unknown stage information were removed. Finally, the data contained 27, 74, 78 and 274 samples from stage I–IV respectively (total 453 samples). The stages I and II were treated as early stage (total 101 samples) and stage III and IV as later stage (total 352 samples). Further, miRNAs and mRNAs expression profiles were normalized using the limma-voom (Law et al., 2014) library package in R3.3.5.

Differential miRNA and mRNAs expression analysis

Differential expression (DE) analysis of miRNAs and mRNAs were performed by the limma-voom using student’s t-test. The mRNAs and miRNAs were considered differentially expressed if (|log₂FC ≥ 1| and p-value < 0.05). The volcano plot and heat map were generated using Enhancevolcano (Blighe, 2018) and Complexheatmap (Gu, Eils & Schlesner, 2016) library in R3.3.5.

Data scaling

In the current study the data was scaled using Z-scaling. The following equation was used for the scaling of data set. (1)${\displaystyle x^{\prime }=\frac{x-\overline{x}}{\sigma }.}$ Here, $x-\overline{x}$, represented the expression variance. While σ is the standard deviation of miRNAs expression. The x′ is the normalized miRNAs expression. The ‘Z-score’ method represent that data scaled to a mean of zero and a standard deviation of 1.

ML model building for stage classification

The “caret” library of R package (Kuhn, 2012) that contains several machine learning, complex regression and classification methods was utilized to develop various ML algorithm based classification models. The machine learning techniques were employed for building the models for patient classification into clinical stages using miRNA expression profiles. In this study, four widely used machine learning (ML) methods namely RF, svmR, AdaBoost, and avNNet were used.

The random forest (RF) method commonly uses non-linear regression models and has been applied in a variety of computational studies. It is a simple, interpretive, and flexible method which allows for a large number of predictor variables. Additionally, it also predicts well in the small sample sizes and high genetic heterogeneity (Alexopoulos, 2010). It comprises building trees from bootstrapped training data in such a way that each split is based on a random sample of m variables from a full set of n attributes (Liaw & Wiener, 2002). (2)${\displaystyle IG\left(n\right)=1-\sum _{i=1}^j{\left(pi\right)}^2.}$

Here IG is Gini impurity of a node n is 1- the sum of overall predictor as j of the fraction of each pi square

Support Vector Machine (SVM) is popular ML technique in medical research for molecular function prediction, genomics variation, histopathology image classifications and subtyping of diseases (Bianchi et al., 2011). It is a non-probabilistic algorithm classifier based on the hyperplane line to maximize the maximum margin to separate two classes based on the two vectors points to achieve the best fit classification (Huang et al., 2018). (3)${\displaystyle \int \left(x\right)={\beta }_0+\sum _{i\in S}{\alpha }_{i}K\left(x,x_i\right).}$

Here, function (x) is the kernel function which calculate the similarity between training and predicted x, x_i and α_i represent the parameter corresponding to each training and predictor variable. β₀ isthe constant.

The AdaBoost is an algorithm for producing strong classifiers from the weak classifier. It is based on the penalized weighted matrix of each instance and voting them according to their weight scheme. It is a popular machine learning method for both balance and unbalanced data sets (Zemel & Pitassi, 2001). The AdaBoost equation is following (4)${\displaystyle H\left(x\right)=\mathrm{sign}\left(\sum _{t=1}^T{\alpha }_t{h}_t\left(x\right)\right).}$

Where, h_t(x) is the out of weak classifier of t for input x and α_t is weight given to model. The α_t = 0.5∗ln(1 − E)∕E. Here, E is based on the error rate of model.

The avNNeT is a kind of neural network capable of learning nearly infinite number of mapping functions. The neural network (NN) behaves like a natural human neuron. Every predictor variables (inputs) connects to an output response variable (output), either directly or through backward and forward propagation of single and several of hidden nodes to calculate units (neurons). The depth of an NN corresponds to the number of hidden layers. The calculation performed in NN are performed in a hidden layer called deep network (Huang, Zhu & Siew, 2004). The NN is commonly used for calculation of gene expression, CNV and clinical data for prognosis prediction of complex diseases such as cancer (Mirza et al., 2019). (5)${\displaystyle a_j^l=\sigma \left(\sum _k{w}_{jk}^l{a}_k^{l-1}+b_j^l\right).}$

Here, $a_j^l$ of the jth neuron in the lth layer is related to the activation in the (l − 1)th layer. The σ is vectorising function. The W^l is weight matrix of each layer, l. The jth row and kth column is $w_{jk}^l$ . Also, for each layer l as define j a bias vector $b_j^l$.

The GBM algorithm is used to convert weak learners into strong learners. It is the tree-decision boosting based classification model where each of instances assigned an equal weight at the start. After the examination of first steps, then increase the weight of those observation that are difficult to classify and lower the weight for those that are easy to classify. This is the iterative process to optimized the best fit for classification (Friedman, 2001). The general equation used for the GBM is given below (6)${\displaystyle h_m\left(x\right)=\sum _{j=1}^{j_m}{b}_{jm}1R_{jm}\left(x\right).}$

Here, m-th step fit the decision tree hm (x). jm number leaves in each steps. The b_jm is the value predicted in the region of R_jm.

Data balancing: A challenge in biomedical studies is working with imbalanced data sets i.e., unequal number of normal and disease samples. Unbalanced predictive variable ratio do not meet the assumptions of the machine learning models and its predict biased. Therefore, SMOTE algorithm (Chawla et al., 2002) was used to balance the data using the ‘DMwR’ library in R.3.3.5 package. The SMOTE is a popular algorithm to deal with unbalanced data (Lusa, 2013; Ramezankhani et al., 2016). There are two types of sampling strategies commonly employed in machine learning viz. oversampling and under-sampling for unbiased classification. In under-sampling, the samples are reduced based on k nearest neighbor (kNN) clustering centroid distance while in oversampling minority classes are amplified to balance both predictors. In our data set, the number of late stage patient samples was reduced so as to match the number of early stage patient samples (N = 101) (Fig. S1A) using the under-sampling method and in oversampling (Fig. S1B), early stage data was amplified (N = 352) to balanced early and later-stage samples.

Thereafter, the datasets (under and over sampled) were systematically distributed into training (70%) and test set (30%) separately. In under-sampled data, 141 and 61 samples were used as training and test set respectively whereas in oversampled data 478 and 206 samples were used as training and test set respectively. First, the model was trained using only training set the test set was kept aside for independent external cross-validation. A 10-fold repeated internal cross-validation with 10-time iterations to randomized data set was done to avoid generation of over/under fitting model. The cost functions were optimized (100 to 3,000 with 100 steps per iteration) to achieve accurate classification. The performance of the generated model was examined based on sensitivity, specificity, accuracy, and AUC (ROC), precision and MCC by external independent data. Equations are given below: (7)${\displaystyle \text{Sensitivity}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}$ (8)${\displaystyle \text{Specificity}=\frac{\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}}}$ (9)${\displaystyle \text{Precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}}$ (10)${\displaystyle \mathrm{MCC}=\frac{\mathrm{TP}\times \mathrm{TN}-\mathrm{FP}\times \mathrm{FN}}{\sqrt{\left(\mathrm{TP}+\mathrm{FP}\right)\left(\mathrm{TP}+\mathrm{FN}\right)\left(\mathrm{TN}+\mathrm{FP}\right)\left(\mathrm{TN}+\mathrm{FN}\right)}}.}$

Where TP is true positive; TN is true negative; FP false positive; and FN is false negative, MCC is Matthews correlation coefficient.

Patient’s risk assessment based on signature miRNA expression

Further, we performed L-1 Penalized Estimation in the Cox Proportional Hazards Model using optimal cross-validated likelihood to identify prognostic markers based on the survival status (alive/dead), survival time and normalized miRNAs expression (Acharya et al., 2017). The Kaplan–Meier curve and the Log-rank method were used to estimate the significance of median miRNA expression on the survival of the patients. It evaluates the independent survival of the patient between the two risk groups (low/high). The patient’s clinical parameters such as age, survival status, gender, metastasis, TNM stages, and median expression of significant miRNAs were evaluated for their association with patient’s survival using cox-regression for a univariate and multivariate analysis.

The risk score for individual patients was calculated using L1-penalized (LASSO) likelihood Cox-regression (Goeman, 2010). As earlier, the data was partitioned into training (70%) and test (30%) to assess the performance of the generated models. The prognostic value of individual miRNAs was first calculated by the L1-penalised Cox univariate proportional hazard regression in R tool using a ‘penalized’ package. The significant miRNAs (p-value < 0.05) were used in subsequent Cox multivariate analysis. The prognostic model was built by the linear combination of the expression level of significant miRNAs multiplied with the Cox regression coefficient (β). The standard formula was as follows: (11)${\displaystyle \text{Risk Score}=X_1{\beta }_1+X_2{\beta }_2+\cdots +X_n\beta n+C.}$

Using the median risk score the patients were divided into high- and low-risk groups. The time-dependent receiver-operating characteristic (ROC) curve was created using “survivalROC” package in R to evaluate the performance (specificity and sensitivity) of the generated miRNA prognostic signature (Chawla et al., 2002).

The mRNAs target identification for identified miRNAs

The mRNAs targets of the identified biomarkers (diagnostic and prognostic) were predicted using the IPA (Ingenuity Pathway Analysis). The IPA predict miRNAs target based on three evidence (i) experimental (ii) high confidence (Kozomara & Griffiths-Jones, 2014) and (iii) low confidence. The experimental and high confidence interactions are curated by experts. In the current studies only experimentally validated and high confidence miRNA-mRNA interactions were used. The miRNAs can have a broad range of mRNAs targets and many of them may not be relevant to the disease under consideration. Therefore, the predicted target mRNAs were considered only if they are also significantly differentially expressed in the HNSCC. As a result, 2 sets (set 1 with 230 and set 2 with 262 mRNAs) were finally selected for diagnosis and prognostic miRNA targets. The miRNA-mRNAs interactions were further visualized for better understanding using Cytoscape3.7 (Shannon et al., 2003).

Pathway enrichment analysis

The previously identified genes (230) of set and 262 were further used for the identification of pathways and other molecular functions, network and toxicity enrichment analysis in the IPA (Ingenuity Pathway analysis). The IPA can analyze enriched biological pathways, disease and functions, upstream and downstream interactions, molecular function and toxicity for a given gene set. It also gives the statistical significance for each enriched process.

Gene ontology (GO) enrichment

The Gene Ontology is a vocabulary developed to compare and identify genes based on the processes in which they are involved. The GO is defined into three categories viz. Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). In the current studies, the enrichment analysis was performed using PANTHER (Mi et al., 2019) (http://www.pantherdb.org/geneListAnalysis.do) and enrichR (https://maayanlab.cloud/Enrichr/) with p-value < 0.001 for identification of significant GO terms.

Differentially expressed microRNAs (DE-miRNAs)

There were 171 miRNAs differentially expressed (62 over and 109 under-expressed) in early and 174 miRNAs differentially expressed (73 over and 101 under-expressed) in late stages (Figs. 2A & 2B) (Tables S2A & S2B). The top 5 over and under-expressed miRNAs are given in Figs. 2C and 2D. The heat map of differentially expressed miRNAs in early stage shows a visual distinction between the expression of miRNAs in normal and cancer tissues (Figs. S2A & S3B).

An analysis of the common and stage specific miRNAs revealed that there are 148 miRNAs common to early and late stages, whereas 23 and 26 miRNAs are specific to early and late stages respectively (Fig. 3A) (Table S3). Similar to the DE-miRNAs analysis, the DE-mRNAs of 101 early stage with 44 normal adjacent tissues (NAT) samples were resulted in the 3831 differentially expressed mRNAs (2407 over and 1424 under-expressed) with (|log2FC ≥ 1| and p-value < 0.05) (Fig. S3).

Identification of miRNA signature for classification of early and late stages

An important objective of this study was to identify the miRNAs signature for the early and late (TNM stage) using the miRNAs expression profiles. The signature miRNAs were identified using the ensemble machine learning feature selection method (Random Forest (RF), Support Vector Machine Radial Kernel (svmR), Adaptive Boost (AdaBoost), averaged Neural Network (avNNet), and Gradient Boosting Machine (GBM)). A detailed method for feature selection by minimum Redundancy Maximum Relevance (mRMR) algorithm is given in Fig. 3B. The top 20 predictor miRNAs were identified using each of the machine learning methods (RF, svmR, AdaBoost, avNNet and GBM) and the common among them were chosen for further model building and evaluation (Fig. 3C, Table S4). The ensemble ML based mRMR feature selection method reduce the bias that might be introduced by one or two methods. As a result, a 6 miRNA (let-7c, miR29c, miR-96, miR-133a, miR-196a, miR-519a) with average importance ranging from 90.03 to 78.27 for classification were selected for further studies. The mean importance of these miRNAs is shown in Fig. 3D. Among them, the miR-96 and miR-196a are over-expressed while the rest (miR-133b, miR-29c, miR-519a and let-7c) are under-expressed in normal vs primary tumor samples (Fig. S4). It is interesting to note that miR-519a is significantly differentially expressed in early stage only.

The performance of generated models for stage classification

The 6 miRNAs were selected by all five machine learning methods. The model was built using a training set and an independent test set (70:30) ratio was used for model prediction. The Synthetic Minority Oversampling Technique (SMOTE) algorithm was used for data balancing. In the under sampling total 101 (early and late stage) samples were used for internal model building and performance was evaluated using accuracy, sensitivity specificity, area under the curve (AUC), precision call and MCC of RF, svmR, AdaBoost, avNNet, and GBM. The accuracy of these five model was found to be 0.79 ± 0.03, 0.71 ± 0.04, 0.76 ± 0.03, 0.76 ± 0.03 0.79 ± 0.03, 0.71 ± 0.04, 0.76 ± 0.03, 0.76 ± 0.03, and 0.64 ± 0.042 for under-sampling methods. Similarly, the accuracy for RF, svmR, AdaBoost, avNNet and GBM was 0.87 ± 0.021, 0.86 ± 0.032, 0.87 ± 0.035, 0.85 ± 0.032, and 0.80 ± 0.038 respectively for oversampled data (early and late = 352 samples) The details of the internal top 10 models of ML is given in the Table S5. The detailed performance of each ML algorithm is given in Figs. 4A & 4B.

The final model performance was evaluated by the independent test set using accuracy, sensitivity, specificity, AUC, precision, and MCC. The accuracy for RF, svmR, AdaBoost, avNNet, GBM was found 0.87, 0.87, 0.87, 0.75, and 0.84 for test set generated using oversampled data. Whereas the accuracy was 0.79, 0.76 0.78, 0.68, 0.77 for the above five methods for the test set generated using under-sampled data . The MCC (0.65, 0.62, 0.64, 0.59, and 0.68 was for RF, svmR, AdaBoost, avNNet and GBM) gave better model performance over other confusion matrix categories. The robustness of model was evaluated using the area under curve (AUC) of Receiver Operating Characteristics (ROC) given in Figs. 4C & 4D for over and under sampled data. The RF and AdaBoost found the subsequently higher prediction accuracy in the under and over sampling methods and avNNet was the least prediction accuracy. The overall classification result showed, the generated 6 miRNAs signature can stratify the HNSCC patients on clinical TNM stages with reasonable accuracy based on expression profile.

Risk score calculation

To assess the prognostic performance of the generated model the time-dependent ROC curve was considered for risk score (low and high) calculation. The overall 5-year risk score was calculated of the normalized expression with the coefficient of the LASSO cox-coefficient. As stated in methodology, the risk score can be used to assess the risk for a patient based on the median value expression. The risk score was calculated using the risk score equation as given below and plotted using Kaplan-Meir plot Fig. 5. Kaplan-Meir plot for each of the signature miRNAs is given in Fig. S5. Risk Score = (0.022 × let − 7_expression) + (−1.996 × miR-96_expression) + (6.593 × miR-9_expression) + (1.361 × miR-143_expression) + (−2.938 × miR-379_expression) + (−1.249 × miR-545_expression) + (−1.969 × miR-658_expression) + (−1.446 × miR-3926_expression) − 1.505.

The patient samples were systematically divided into training (368 samples) and test sets (160 samples). The 5 risk score (high and low) prediction was evaluated by AUC. The AUC for the training set and test set was found to be ∼0.86 and ∼0.79, respectively which indicated good performance of the identified signature (Figs. 6A & 6B). This time-dependent risk prediction analysis clearly showed that these miRNAs can be used for risk assessment of HNSCC patients.

Univariate and multivariate analysis for prognostic evaluation

The identified miRNA signature was evaluated against various clinical variables to check their effect on survival using the standard Cox-proportional hazards regression (Table 1). It indicates that advancing pathological stages may have a significant adverse impact on the survival of the patients. The Cox univariate analysis showed that pathological stage, T stage, N stage, metastasis and signature miRNAs were significantly associated with the prognosis. While in multivariate analysis pathological T stage (p-value = 0.017) and signature mRNA (p-value = 0.023) showed significant p-values associated with the poor prognosis (Table S6).

The microRNA-mRNA interactions

It was imperative to understand the role of the identified miRNA signatures in different biological processes. A mapping of the mRNA targets of the identified signatures (6 miRNAs for stage classification and 8 for prognosis) resulted in the identification of 1532 and 1446 mRNAs, respectively. Further, the analysis was done to select those mRNAs that get differentially expressed in various stages of HNSCC. It was found that for the identified diagnostic and prognostic miRNA signatures, there are 2 set of mRNA’ target (set 1 with 230 and set 2 with 262) mRNA that were differentially expressed respectively in the HNSCC patients (Tables S7 & S8).

Biological pathways, diseases functions and network enrichment analysis

As indicated earlier, the pathways enrichment analysis was performed by the Ingenuity Pathway Analysis (IPA) separately for diagnostic targets (set 1) and prognostic targets (set 2).

Analysis of set 1: The top enriched pathways for diagnostic targets were (i) GP6 signalling (ii) Neuroinflammation signalling (iii) Intrinsic Prothrombin Activation Pathway (iv) granulocyte adhesion and diapedisis (Fig. 7). Other enriched pathways are given in Table S9. The KEGG pathways analysis also revealed that majority of the genes were associated with pathways such as cytokine-cytokine receptor interaction, small cell lung cancer, transcriptional misregulation in cancer, ECM-receptor, and glutathione metabolism (Table S10). While WikiPathways enrichment showed, nuclear receptors meta-pathways (ID: WP2882), TGF-beta signalling pathways (WP366), PI3K-Akt-mTOR-signalling pathway (WP4172) (Table S11). The pathway analysis revealed that these miRNAs targets are highly enriched with cancer associated pathways. The top disease and disorders associated with the set 1 are organismal injury and abnormalities, cancer, connective tissue disorder, and skeletal muscle disorder. Top molecular functions and cellular functions were linked with cellular movement, cellular assembly and organization, cellular function and maintenance, cell death & survival and cell development (Fig. 8) (Table S12). The top upstream regulators were miR-29c, TGFB1, estrogen receptor, histone h4, and CREB which are associated with cell growth and proliferation. The miR-29, a potent tumor suppressor, is found under-expressed in this study. Interestingly, the top analysis of the upstream regulatory genes showed that miR-29 directly targets 11 genes (9 over-expressed and 2 under-expressed) associated with head and neck cancer (Fig. S6). A thorough look at these mRNAs revealed that many of them are transcription factors (total 31) and part of cancer gene panels (total 19) (Sondka et al., 2018). These mRNAs were used in the subsequent functional analysis step. It was imperative to note that many of the transcription factors and cancer censes gene are targeted by the let-7c which is well known for its role in a variety of cancers (Fig. S7). The toxicity enrichment analysis revealed increased levels of alkaline phosphatase (ALP). Interestingly, the serum ALP level have been found to be significantly higher in OSCC patients (Acharya et al., 2017).

Analysis of set 2: The pathway enrichment analysis by IPA showed that the top enriched pathways are neuroinflammation signalling, small cell lung cancer, MAPK signalling, and cell cycle: G2/M DNA damage checkpoints (Fig. 9) (Table S13). Top prognostics signature miRNAs are directly targeting several key cancer regulatory pathways such as p53 signalling, ATM signalling, EMT pathways etc., (Fig. 10). The KEGG pathways analysis also revealed that majority of the genes were associated with cancer progression such as protein digestion and absorption, ECM-receptor, transcriptional misregulation in cancer, and focal adhesion. In case of WikiPathways, the ECM and membrane receptors (ID: WP291), miR-509-3p alteration of YAP1/ECM axis (WP3967), and focal adhesion-PI3K-Akt-mTOR-signalling pathway (WP3932) were found enriched (Table S14).

The top regulatory effect networks were MITF, CD3, YAP1, let7, TGF β1, and TNF. MITF is mitochondrial transcription factor genes and have essential roles in cell differentiation, proliferation and survival of cell (Bertolotto et al., 1998). It is also important co-activator of TGF β (Nijman et al., 2006). CD3 is cluster of differentiation 3 activates in early oral premalignant (Öhman et al., 2015). The Hippo-YAP1 pathway is oncogenic in head and neck cancer (Santos-de Frutos, Segrelles & Lorz, 2019). TGF β1 dysregulation is very common in several cancers including HNSCC. It is key regulator of epithelial cell proliferation, growth factor and angiogenesis (Rothenberg & Ellisen, 2012). The dysregulation of TNF genes are involved in almost all types of cancer. The TNF plays an important role in cytokine production and have critical role in immune regulation (Alsahafi et al., 2019).

The top diseases and disorders found to be associated with the genes of this set were cancer, organismal injury & abnormalities, and hematological disorders. Additionally, top associated network functions were associated with cell death and survival, cellular development, cellular growth, and proliferation (Table S15). The miRNA-mRNA interaction analysis showed that there are 16 mRNAs which are part of cosmic consensus cancer gene panel and there are 20 transcription factors among the targets. Six mRNAs were common to transcription factors and cancer gene panels (Fig. S8). Overall enrichment analysis revealed that these miRNAs were directly or indirectly associated with the tumorigenesis and poor prognosis.

Gene ontology enrichment analysis

Set 1: The most enriched biological processes were associated with extracellular matrix proteins organization (GO:0030198), regulation of natural kill cell chemotaxis (GO:2000501) and response to alcohol (GO:0097305). While top molecular functions enriched genes were mRNA 3′-UTR binding (GO:0003730), chemokine activity (GO:0008009) and transforming growth factor beta-activated receptor activity (GO:0005024). The cellular component such as messenger ribonucleoprotein complex (GO:1990124), serine/threonine protein kinase complex (GO:1902554), gamma-tubulin large complex (GO:0000931) were found majority of genes involved in these process (Fig. S9).

Set 2: Many of the enriched GO processes are similar to set 1. However, there are notable differences as well. The most enriched genes were involved in biological process such as collagen fibril organization (GO:0030199), negative regulation of myeloid cell differentiation (GO:0045638), and negative regulation of cell differentiation (GO:0045596). In the molecular functions enriched genes were platelet-derived growth factor binding (GO:0048407), transcriptional activator activity, RNA polymerase II transcription regulatory region sequence-specific binding (GO:0001228), cell proliferation (GO:0008283), metabolic process (GO:0008152), biological adhesion (GO:0022610). The endoplasmic reticulum lumen (GO:0005788), integral component of plasma membrane (GO:0005887), and intermediate filament (GO:0005882) were most enriched cellular component terms (Fig. S10).