Construction and validation of a novel prognostic signature of microRNAs in lung adenocarcinoma

Data and sources

The preprocessed LUAD mature miRNA expression profiles in the Cancer Genome Atlas (TCGA) database displayed as log2 converted reads per million (log2(RPM + 1)) were downloaded, along with the clinical information of related samples including survival follow-up data, age, gender, tumor node metastasis (TNM) stage, and tumor recurrence. The miRNA expression data from the LUAD samples were based on IlluminaHiSeq_miRNASeq, followed by quantile normalization and log2-scale transformation by R software. All the enrolled samples were 2:1 random sampling without replacement by using sample R package. The randomizing process was blind to the demographic information or clinical information, which could ensure that the training and validation cohorts are independent to clarify the key points we focused on. We use R (http://www.R-project.org) software to make all statistical analyses. A two-sided P < 0.05 was considered statistically significant.

Herein, the predictive model was initially constructed in the training set, and the validation group was employed to provide the unbiased evaluation of the final model in the training set.

In addition, 12 cases of fresh LUAD tissue samples and paired adjacent cancer tissue samples were enrolled from the Department of Thoracic Surgery, the First Affiliated Hospital of Anhui Medical University. All the 24 surgically resected samples were immersed in 10% formalin. Informed consents were signed by all patients and our study was approved by the Ethics Committee of the First Affiliated Hospital of Anhui Medical University (anyiyifuyuanlunshen-Quik-PJ-2019-09-09) (anyiyifuyuanlunshen-Quik-PJ-2020-05-08). For the patients collected from the TCGA project, the ethics and policies were announced on their Website (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga/history/policies).

Establishment and validation of a miRNA signature with the least absolute shrinkage and selection operator (LASSO) regression model

The prognostic significance of every single miRNA was preliminarily assessed by Kaplan–Meier (K-M) survival analysis, followed by a log-rank test to select all the potential miRNAs. The candidate miRNA was considered to be significantly associated with the OS of LUAD patients when the P- value was less than 0.05 (Xu et al., 2017). A group of 1,000 training matrices based on the LUAD patients among the training set was subsequently obtained after the resample model inclusion probability (RMIP). Followed, with the R package “glmnet” (2.0–10), the LASSO bagging cox regression test, the classical and modified method in Cox regression analysis of high dimensional data was performed with 10-fold cross-validation on all the OS related matrices (Tibshirani, 1997). The tuning parameter λ was determined by the standard error (1-SE), which introduced the penalties to the model construction process to avoid over-fitting. The larger λ was, the more the coefficients were shrunk toward zero (Goeman, 2010). In this study, the recommended λ value with minimal mean squared error plus one standard deviation was selected in each Lasso cox regression test. A group of 1,000 matrices was generated after resampling the data points 1,000 times of the training group, and the LASSO Cox regression test was performed. A list of miRNAs, which displayed more than 500 bootstraps (2/3 of all) during the down-sampling test, was eventually picked up as prognostic biomarkers. Finally, an individual’s risk score model for each patient was constructed for predicting prognosis of LUAD patients by including the expression level of each optimal prognostic miRNA, weighted by their estimated regression coefficients of LASSO Bagging Cox regression model as follows:

Risk Score (patient) =∑ i coefficient miRNA Þexpression miRNAiðÞ.

The risk score of each LSCC patient was calculated according to the prognostic predicting formula. Afterward, the LUAD samples with assigned risk scores were classified into a high-risk or low-risk group by using the median as the cutoff point. To confirm the predictive ability of this model, K-M survival curves with the log-rank test in the training and validation sets were also conducted to distinguish the survival difference between the high-risk group and low-risk group. A two-tailed P-value < 0.05 in K-M analyses was considered as statistically significant. In addition, to graphically show the diagnostic value of the classifier system, receiver operating characteristic (ROC) analysis was adopted to compare the sensitivity and specificity of survival prediction. Youden’s index = max [1- (sensitivity + specificity)] (Hilden & Glasziou, 1996). Moreover, the area under the ROC curve (AUC) for evaluating discriminatory ability was calculated as well.

Network and annotation analysis of target genes

We obtained the target genes from TargetScan, miRTarBase, and miRDB (Wong & Wang, 2015; Chou et al., 2018; Agarwal et al., 2015), respectively. To obtain more convincing results, we overlapped the target genes from the three websites and generated a new target gene cluster. Cytoscape software (v3.5.1, San Diego, La Jolla, California 92093, USA) was used for the visualization of the interactive network to assist the understanding miRNA function and regulatory mechanisms. Besides, Hallmark gene set and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to further study the functional roles of the target genes during LUAD progression using the R package “clusterProfiler” (Q-value < 0.05) (Yu et al., 2012).

Immunohistochemistry (IHC)

We randomly chose the SOX9, RASA1, CEP55 and MAP4K4, the target genes of miR-582, for IHC validation, which was performed according to the previous description (Zhang et al., 2019a; Zhang et al., 2019b). The t-test was used to determine the protein expression differences between LUAD tissue specimens and paired adjacent cancer tissue specimens. A two-tailed P-value < 0.05 was considered as statistically significant. The SOX9, RASA1, CEP55 and MAP4K4 antibodies (catalog # GR3253329-3; GR85944-25; GR234851-12; 09000421) were purchased from Abcam (Cambridge, UK; http://www.abcam.com) and Proteintech (Chicago, IL, USA; http://www.ptglab.com).

Construction of a prognostic signature based on miRNAs

The miRNA expression profile and clinical information of 518 LUAD cases were downloaded from the TCGA databases. OS information was available in 499 cases, which were further retained for downstream analysis. Among LUAD patients, the samples were randomly divided into two groups as the training set (329 samples) and validation set (170 samples).

The K-M survival analysis was employed for the preliminary screening of potential OS related miRNAs (Fig. S1). As a result, we selected eight miRNAs based on P < 0.05, which indicated significant correlation with the LUAD patients’ OS. Afterwards, we rebuilt 1,000 resample LUAD training matrices based on the training cohort patients. The RMIP of each OS-related miRNA was listed and ordered from the most frequent one to the least miRNA (Fig. 1A). According to the LASSO cox regression analysis on performance of all the OS related matrices with 10-fold cross validation, six miRNAs were acquired, including miR-1468 (HR = 0.69, 95% CI [0.560–0.870], P = 0.001, co-ef = −0.366), miR-299 (HR = 1.48, 95% CI [1.220–1.800], P <0.001, co-ef = 0.394), miR-4709 (HR = 0.720, 95% CI [0.590–0.890], P = 0.002, co-ef = −0.322), miR-5571 (HR = 0.570, 95% CI [0.370–0.860], P = 0.008, co-ef = −0.570), miR-582 (HR = 1.18, 95% CI [1.050–1.320], P = 0.005, co-ef = 0.166), and miR-3653 (HR = 0.87, 95 % CI [0.730–1.030], P = 0.106, co-ef = −0.143) (Fig. 1B). A risk score model for each patient was constructed to predict prognosis of LUAD patients by including the expression level of each optimal prognostic miRNA, weighted by their estimated regression coefficients of LASSO Bagging Cox regression model as = −0.365 * miR-1468 + 0.394 * miR-299 -0.322 * miR-4709 –0.570 * miR-5571 + 0.166 * miR-582 - 0.143 * miR-3653.

Validation and evaluation of the miRNAs signature

Using the median risk score as the cut-off value, the LUAD patients were assigned into high-risk or low-risk group. The K-M curve with the log-rank test demonstrated significantly different survival between high-risk and low-risk groups in both training (P < 0.0001, Fig. 2B) and validation groups (P = 0.0085, Fig. 2D). To evaluate the classification potential of this model, the ROC analysis was performed and AUC was calculated in both training and validation cohorts. As shown in Figs. 2A and 2C, the AUC of ROC was 0.72 and 0.70 in the training group and validation group, respectively. Besides, considering other partitions to training and validation groups can both increase and decrease the reliability of the classifier, we repeated the step of different training and validation groups multiple times and efficiently calculated the means/variations/confidence intervals for the sensitivity, specificity and AUC to validate the reliability of the classifier as the Table S3. These results proved that the classifier has moderate discriminatory power and potential utility in determining the LUAD patients with low or high risk as presenting high and stable of AUC and sensitivity.

Moreover, nomograms were constructed in the training, validation, and overall groups after incorporating the clinical characteristics (age, sex, TNM stage, and tumor recurrence) (Figs. 3A–3C). The results of ROC analyses identified this six-miRNAs classifier as the most accurate contributor to prognosis compared to other clinical factors, and the nomogram could precisely and steadily judge the OS rate of LUAD.

Subgroup analysis

Age, TNM stage, and tumor recurrence are always important clinical factors influencing the prognosis of LUAD. Therefore, stratification analyses based on TNM stage, tumor recurrence, and other clinical variables were performed to validate whether this classifier had functional discriminatory capacity among different clinical status (Fig. 4). According to the TNM staging system, LUAD patients have categorized two groups (stage I/II, stage III/IV). Within Stage I/II, the OS related to the six miRNAs signature prognostic value remains significantly different (stage I/II: P < 0.001), but not in Stage III/IV subgroup (stage III/IV: P = 0.061). When adjusted by age, sex, and new tumor event after initial treatment, our risk score model still exhibited the reliability and general applicability for distinguishing each group among different clinical factors.

Functional annotation analysis and IHC Validation

The target genes of miR-1468, miR-299, miR-4709, miR-5571, and miR-582 were predicted by the online prediction tool like TargetScan, miRTarBase, and miRDB. The miRNA-miRNA regulatory network was shown in Fig. 5. Additionally, we performed functional annotation analysis for Hallmark gene set and KEGG pathways (Fig. 6), demonstrating that target genes of different miRNAs were significantly enriched in several functional pathways related to cancer development, such as hallmark UV response, Wnt signaling pathway and mTOR signaling pathway.

Moreover, to confirm the protein expression pattern of the target genes in the LUAD tissue, the expression of SOX9, RASA1, CEP55, and MAP4K4 was detected in LUAD clinical specimens by IHC. The results showed that the LUAD patients mostly had high SOX9, RASA1, CEP55 and MAP4K4 protein expression (p < 0.05) (Figs. S2–S5).