Diagnosis value of aberrantly expressed microRNA profiles in lung squamous cell carcinoma: a study based on the Cancer Genome Atlas

TCGA dataset and sample information

Up to February 12, 2017, a total of 504 LUSC samples with related information were obtained from the TCGA database. According to the inclusion criteria: (1) histologic diagnosis was LUSC; (2) data of samples were available on expression of genes and characteristics of patients. For the exclusion criteria: (1) first histologic diagnosis was not LUAD; (2) patients suffered from other malignant neoplasms except LUSC; and (3) overall survival more than five years. At last, this study embraced 332 eligible patients. Among these 332 LUSC, related information of LUSC tissue samples was from 293 patients and adjacent non-tumorous information were from 39 subjects. In addition, personal comprehensive sources including related clinical information and RNA expression data on LUSC patients were also obtained from the Data Coordinating Center. Among these 293 cases, 105 LUSC subjects had lymphatic metastases and 188 LUSC subjects had no lymphatic metastases. Moreover, according to the staging system of the Union for International Cancer Control (UICC), well or moderately differentiated LUSC (stage I–II) were 243 cases, and poorly differentiated LUSC (stage III–IV) were the remaining 50 cases. Because the data were provided by the TCGA, there was no need for approval of the Ethics Committee. This study fully met the guideline of the NIH TCGA human subject protection and data access policies.

In addition, 47 LUSC tissues samples were obtained from the Nanjing Chest Hospital of Southeast University. These samples including tumor tissues and adjacent non-tumor tissues were obtained immediately after surgical resection from patients who do not undergo preoperative radiotherapy or chemotherapy. These patients consented in advance and signed informed consent forms. Tissues were stored in RNA later (Ambion, Austin, TX, USA) at −80 °C until further use. This study was approved by the ethics committee of Zhongda Hospital Southeast University.

Differential analysis of expressed miRNAs in LUSC samples

The LUSC RNA expression data (level 3) of patients were collected from the TCGA Data Portal. The TCGA database provided the normalized expression profile data of RNA sequencing including mRNAs and miRNAs by RNASeqV2 and Illumina HiSeq 2000 miRNA sequencing platforms (Illumina Inc., Hayward, CA, USA), respectively.

Then, abnormally expressed miRNAs were compared in level 3 (fold changes >2, p < 0.05), including LUSC tumor samples vs. adjacent non-tumorous lung samples, lymphatic metastases of LUSC samples vs. non-lymphatic metastases of LUSC samples, and I–II stage vs. III–IV stage, respectively. Figure 1 shows the flow chart for bioinformatics analysis.

Correlation analysis between LUSC-specific intersection miRNAs, characteristics of patients and overall survival

To explore the connection between LUSC-specific intersection miRNAs and clinical information, related miRNAs were found, using the comparative analysis of LUSC miRNA sequencing data in TCGA. The characteristics of patients was divided into five parts including gender, tumor grade, TNM stage, lymph node metastasis, and patient outcome assessment. Subsequently, Student t-test was used to evaluate the association between the LUSC-specific miRNAs and personal information using IBM SPSS Statistics Version 21.

Furthermore, to correlate specific intersection miRNAs with patients’ prognosis characteristics, the univariate Cox proportional hazards regression model and Mantel–Haenszel test were used to explore the association between specific intersection miRNAs and LUSC patient survival. At the same time, the overall survival (OS) curves were evaluated. LUSC-specific intersection miRNAs associated with OS were defined two independent subclasses by using hazard ratios (HRs), whose cutoff was that a positive signature showed HR < 1 and hazardous signature had HR > 1.

Total RNA extraction and qRT-PCR verification

Total RNA was extracted tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the instructions. Concentration and integrity of all extracted RNA were evaluated to be satisfactory.

To confirm the reliability and validity of TCGA data, we randomly selected five specific key miRNAs (miR-205-5p, miR-30a-3p, miR-30a-5p, miR-30c-2-3p and miR-30d-5p) of above intersection miRNAs. Then their actual expression levels in 47 pairs of LUSC samples were determined by quantitative real-time polymerase chain reaction (qRT-PCR).

A two-step protocol of reverse transcription reactions: First, 1 μg of RNA samples were pre-denatured (5 min at 65 °C and held at 4 °C). Then the 9 μl mixture including 2 μl 5 × RT buffer, 0.5 μl RT Enzyme Mix, 0.5 μl miRNA-specific stem-loop RT primers and 6 μl RNase-free water were added in the 1 μg pre-denatured RNA (37 °C × 15 min, 98 °C × 5 min and subsequently held at 4 °C).

Real-time PCR was then performed by Thunderbird SYBR qPCR Mix (QPS-201; TOYOBO, Osaka, Japan) according to the manufacturer’s protocol. The PCR reaction components comprised 1 μl cDNA, 5 μl Thunderbird SYBR qPCR Mix, 0.3 μl PCR primers (RiboBio, Guangzhou, China) and 3.4 μl RNase-free water. Then a two-step protocol [95 °C × 1 min; 40 cycles of (95 °C × 15 s, 60 °C × 30 s, 72 °C × 30 s)] was undertaken in on StepOnePlus PCR System (Applied Biosystems, Waltham, MA, USA). Results were normalized to the expression of U6. The results were expressed as mean ± SD. Each sample expression was calculated by the 2^−ΔΔCt method (ΔCt = Ct_miRNAs − Ct_U6 and ΔΔCt = ΔCt_{tumor tissues} − ΔCt_{adjacent non-tumor tissues}). Paired t-test was applied to for comparison between tumor tissues and the adjacent non-tumor lung tissues. In all cases, differences with p < 0.05 were considered to be statistically significant.

ROC curve analysis of specific key miRNAs

According to the outcome of 332 patients, receiver operating characteristic (ROC) curve was also applied to evaluate the diagnosis value of the specific key miRNAs of LUSC.

The prediction of microRNA target genes

It is known that mRNAs are regulated by corresponding miRNAs. So, to find the differential display of targeted mRNAs, samples were quartered similarly as abnormally expressed miRNAs, whose flow chart of bioinformatics analysis is presented in Fig. 1. In addition, interaction mRNA of LUSC-specific intersection miRNAs was predicted by using the Targetscan and the miRbase (Hsu et al., 2014). Last, according to the results of previous two steps, the mRNAs of the microRNA target genes were obtained.

Differential expression of key mRNA function

To understand the potential biological processes and pathways of aberrant expression of intersection mRNAs, bioinformatics resources from database for annotation, visualization, and integrated discovery (DAVID) were used. We were only interested in the significant level (the enrichment score >2 and p < 0.05) of gene ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways to analyze the potential role of these key mRNAs. Moreover, to further explore the relationship and function of intersection mRNAs, the PPI network was analyzed using the protein–protein interaction (PPI) network via STRING (http://string-db.org; Version 10.5).

Significantly specific miRNAs in LUSC

In this study, 332 LUSC patients were selected from TCGA database. Then the expression of 332 LUSC and adjacent normal tissue miRNAs were compared to identify significantly differential miRNAs through a certain criterion (fold change ≥2 or ≤0.5, and p < 0.05) (Romero-Cordoba et al., 2012; Xu et al., 2014; Zhu et al., 2015).

Then 97 LUSC-associated abnormally expressed miRNAs were identified between 293 LUSC tumor tissue samples and 39 adjacent non-tumor tissues. Further, these 97 miRNAs between tumor stage and lymphatic metastasis were analyzed. And aberrantly expressed miRNAs were selected from four groups: (1) I–II stage (non-lymphatic metastases) LUSC tissues and adjacent non-tumorous lung tissue; (2) III–IV stage (non-lymphatic metastases) LUSC tissues and adjacent non-tumorous lung tissue; (3) I–II stage (lymph node metastasis) LUSC tissues and adjacent non-tumorous lung tissues; and (4) III–IV stage (lymph node metastasis) LUSC tissues and adjacent non-tumorous lung tissues. At last, we found 70, 52, 76 and 84 aberrantly expressed miRNAs in these four groups respectively (fold change >2, p < 0.05). To further enhance data reliability for bioinformatics analysis, 42 aberrantly expressed miRNAs in the intersection of these four groups were selected. Among these 42 intersection miRNAs, 17 miRNAs were upregulated and 25 miRNAs were downregulated (Figs. 2 and 3; Table 1).

Association between LUSC-specific miRNAs and characteristics of patients

After analyzing the correlation between expression of the 42 LUSC-specific intersection miRNAs and their related characteristics of patients (race, gender, age, TNM stage, lymphatic metastases, and patient outcome status in the TCGA database), six specific miRNAs were significantly abnormal in clinical information (p < 0.05, Table 2).

Two miRNAs (hsa-miR-629-3p and hsa-miR-511-5p) were aberrantly expressed in gender, four miRNAs (hsa-miR-30c-2-3p, hsa-miR-130b-5p, hsa-miR-629-3p and hsa-miR-130b-3p) were aberrantly expressed in lymphatic metastasis and one miRNA (hsa-miR-30a-3p) were aberrantly expressed in patient outcome assessment, however, none of aberrantly expressed miRNA was discovered in race, age, TNM stage and tumor stage.

Survival analysis

A univariate Cox model was used to identify the 42 LUSC-specific intersection miRNAs correlated with OS. Five of these miRNAs were significantly associated with the overall survival status of 332 LUSC patients (log-rank p < 0.05): three (miR-139-5p, miR-139-3p and miR-326) negatively (p < 0.05) and two (miR-30d-3p and miR-486-5p) positively (p < 0.05) (Fig. 4).

qRT-PCR verification

The actual expression levels of five specific key miRNAs (miR-205-5p, miR-30a-3p, miR-30a-5p, miR-30c-2-3p and miR-30d-5p) were measured by qRT-PCR. The results demonstrated that miR-205-5p was upregulated in LUSC tumor tissues, while miR-30a-3p, miR-30a-5p, miR-30c-2-3p and miR-30d-5p were significantly downregulated in LUSC tumor tissues. Expression levels of these five miRNAs were consistent with the TCGA data (Fig. 5; Table 1).

ROC curve analysis of specific key miRNAs

Receiver operating characteristic curves were also used to assess diagnosis value of the LUSC-specific intersection miRNAs which related with overall survival. The area under ROC curve (AUC) = 0.930, 0.950, 0.961, 0.919 and 0.952 for miR-30d-3p, miR-139-5p, miR-139-3p, miR-326 and miR-486-5p (p < 0.01, Fig. 6), which are the score all higher the cutoff (0.7) and could be considerable biomarkers for early diagnosis of LUSC.

Prediction of target genes

In 474 LUSC samples from TCGA database, 18,633 mRNAs were obtained from the TCGA database. According to the cutoffs that was fold change ≥2 or ≤0.5, and p < 0.05, significantly aberrantly expressed mRNAs were identified between 293 LUSC tumor tissue samples and 39 adjacent non-tumor tissues. In total, 14,689 LUSC-associated abnormally expressed mRNAs were selected to be further compared between tumor stage and lymphatic metastasis. 3,681 aberrantly expressed mRNAs were selected from comparisons of I–II stage (non-lymph node metastasis) LUSC specimens and adjacent non-cancerous lung tissues, 3,848 from comparisons of I–II stage (lymphatic metastasis) with adjacent non-cancerous lung tissues, 3,275 from comparisons of III–IV stage (non-lymph node metastasis) and adjacent non-cancerous tissues, and 3,885 from comparisons of III–IV stage (lymphatic metastasis) and adjacent non-cancerous tissues. As a result, we selected 2,950 aberrantly expressed mRNAs (1,770 downregulated and 1,180 upregulated) from the intersection of the above four groups (Fig. 7).

The next step was to predict the mRNA of LUSC-specific intersection miRNAs. 577 interaction mRNAs were obtained by using the Targetscan and the miRbase.

Based on the above steps, the interaction mRNAs that appeared simultaneously in previous results were selected. It is shown that 2,950 aberrantly expressed mRNAs in the first step contained all 577 interaction mRNAs in the second step, which meant that the 577 mRNAs were the target genes.

Function analysis

Furthermore, functions of LUSC-specific intersection miRNAs were predicted by interaction mRNAs with DAVID bioinformation.

Eighty-one pathways were indicated by KEGG pathway analysis, and 433 GO terms (p < 0.05 and enrichment >2) were identified by analyzing the enrichment of these genes. Moreover, KEGG pathway analysis showed the most significant networks were Staphylococcus aureus infection (path ID: 05150) and cell cycle (path ID: 04110) (Fig. 8). The result suggested that five of the top 20 pathways were tumorous pathways including cell cycle, cell adhesion molecules (CAMs), cytokine–cytokine receptor interaction, Rap1 signaling pathway, and pathways in cancer. There were other tumor-related pathways such as cAMP signaling pathway, metabolic pathways, phagosome, PI3K-Akt signaling pathway, miRNAs in cancer, small cell lung cancer, and hippo signaling pathway in the other significantly differentially altered pathways. The highest enriched GO terms were mitotic cell cycle (GO: 0000278) and GO 0007165 (GO: 0007268) (Fig. 8). The relationship and function of 577 mRNAs were revealed in the PPI network (Fig. 9).