Integrated analysis of two-lncRNA signature as a potential prognostic biomarker in cervical cancer: a study based on public database

The Cancer Genome Atlas database and samples

RNA sequencing data from 307 cases of cervical squamous cell carcinoma (CESC) was downloaded from the TCGA database (up to May 26, 2017, https://portal.gdc.cancer.gov/), including individuals with clinical information. The following patient or lncRNA exclusion criteria were applied: 1) no complete tissue samples were present for data analysis; 2) histologic diagnosis was not CESC; 3) other malignancies apart from CESC were present; 4) The differentially expressed lncRNAs which were present in no more than 10% of the samples were eliminated. Overall, 289 CC patients were included in the present study providing 289 CC tissue samples and three adjacent normal tissue samples. The clinical data included staging information and outcomes of CESC patients. Processing this data did not violate the requirements of TCGA’s human protection and data access policies. (http://cancergenome.nih.gov/publications/publicationguidelines).

In addition, 31 CC patient tissue specimens (primary solid tumor and solid normal tissue) were collected from the Zhongda Hospital (Nanjing, China) of Southeast University, between 2016 and 2017, for qRT-PCR analysis. The tissues were snap-frozen in RNAlater (Ambion, Austin, TX, USA) after surgical resection and stored immediately in liquid nitrogen for subsequent total RNA extraction and analysis. These 31 patients (aged 23–64 years) were diagnosed with CC based on the histopathology and clinical history. All patients signed informed consent and this study was approved by the Zhongda Hospital Southeast University ethics committee.

Identification of differentially expressed lncRNAs

The TCGA database provides normalized count data for RNA sequencing through the RNASeqV2 system, including lncRNA and mRNA expression profiles. The CESC level 3 lncRNA sequencing raw data were obtained through Illumina HiSeq 2000 RNA sequencing platforms (Illumina Inc., Hayward, CA, USA). Data were already normalized by the TCGA. The data of tumor CC tissues were already normalized to the adjacent normal tissues. RNA-Seq expression level read counts are normalized using two related methods: The Fragments per Kilobase of transcript per Million mapped reads (FPKM) and The upper quartile FPKM (FPKM-UQ). The differentially expressed RNAs included those upregulated and downregulated with fold changes >2 and <0.5, respectively, and adjusted false discovery rate at P < 0.05. To detect the differential expression of lncRNAs, samples were divided into three groups based on FIGO stage I, stage II and stages III–IV. The intersection of the lncRNAs was selected for further analysis. In addition, sequencing results of differentially expressed lncRNAs that showed no change in more than 10% of all samples were eliminated.

Identification of the specific lncRNA prognostic signature

The expression level of each differentially expressed lncRNA was transformed by log2 to calculate the risk score. Subsequently, we further analyzed clinical features related to lncRNAs in CC. The univariate Cox proportional hazard model was used to analyze the effects of risk score and clinical features on the OS of CC patients (Bair & Tibshirani, 2004; Gao et al., 2016). The multivariate Cox proportional hazard model was used to evaluate the prognostic value of these OS-related lncRNAs. Based on the previously reported risk score model (Zhao & Sun, 2007), we constructed a prognosis-related risk score based on lncRNA expression levels. The formula is: Riskscore = explncRNA1 × β lncRNA1 + explncRNA2 × β lncRNA2+…+explncRNAn × β lncRNAn, where exp represents expression level and β represents the regression coefficient from the multivariate Cox regression model (Zeng et al., 2017). We used the median as the cutoff point in risk score. The CC patients were divided into high- and low-risk score group, respectively (Zhou et al., 2016). The receiver operating characteristic (ROC) curve analysis within 5 years was used to calculate the predictive value of the risk score for time-dependent outcomes (Heagerty, Lumley & Pepe, 2000).

qRT-PCR verification of specific lncRNA expression in CC tissues and ROC curve analysis

We used qRT-PCR to analyze actual expression levels and validate the accuracy and reliability of the two lncRNAs in 31 newly diagnosed CC patients. All results were normalized to the control reference gene GAPDH. Following the standardized manufacturer’s protocol, total RNA was isolated from tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA purity was determined using a NanoDrop 2000 spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). RT reactions and qRT-PCR were both conducted according to the manufacturer’s protocol using the Reverse Transcription System and qPCR Master Mix kit (Promega, Madison, WI, USA), respectively. Step One Plus^TM PCR System (Applied Biosystems, Foster City, CA, USA) was used to detect the expression levels of lncRNAs. All the primers were produced by Generay Biotech Co., Ltd. (Shanghai, China). qRT-PCR results were calculated by the 2^−ΔΔCt method using the formulas ΔCt = Ct_ΔlncRNAs – Ct_ΔGAPDH and ΔΔCt = ΔCt_{tumor tissues} – ΔCt_{adjacent non-tumor tissues} (Wang, Chen & Liu, 2015). The ROC curve was used to evaluate the diagnostic value of the expression levels of the two lncRNAs. A P value of <0.05 was considered statistically significant.

Cell culture and qRT-PCR verification of expression of two lncRNAs in cervical cancer cell lines

Three human CC cell lines (Hela, SiHa and C33-a), and an immortalized human cervical epithelial cell line (H8) were purchased from Shanghai Chuyu Biological Co., Ltd. The cells were cultured in 5% CO₂ humidified atmosphere at 37 °C, using high-glucose Dulbecco’s modified Eagle’s medium (HyClone, South Logan, UT, USA) supplemented with 10% fetal bovine serum, 100 U/mL penicillin and 100 mg/mL streptomycin. qRT-PCR and total RNA extraction conditions, reagents and methods are the same as for tissues.

Statistical analysis

Statistical analysis was performed by IBM SPSS Version 24.0. The final results are shown as means ± SD. Student’s t-tests were used to compare two groups of sequencing data. In all cases, P value of <0.05 was considered statistically significant. Fold change was used to analyze the statistical significance of results. Clinical parameters and risk scores were screened and compared using univariate and multivariate Cox regression models. In addition, based on the expression of lncRNA levels in CC patients, we performed ROC curve and area under curve (AUC) analyses to judge diagnostic value. R language was used as the main tool for generation of ROC curves (R Core Team, 2017).

Identification of significantly differentially expressed lncRNAs and clinical characteristics of patients

There were 289 CC primary solid tumors and three solid normal tissues with clinical patient information obtained from TCGA which were included in the present study. The average age of patients was 47.31 ± 13.62 years. The OS time was 1045.70 ± 67.07 days, 71 patients died. Significant differentially expressed lncRNAs were identified based on the criteria of fold change >2 and <0.5 at P value <0.05. According to inclusion–exclusion criteria, we obtained 49 differentially expressed lncRNAs included in the intersection of the three groups’ analyzed results (Fig. 1A; Table 1). The available clinical features from the TCGA database are shown in Table 2.

Construction of an lncRNA signature significantly associated with prognostic features

The 49 differentially expressed specific lncRNAs were further analyzed. We analyzed race, age, body mass index, human papillomavirus, tobacco use, clinical stage, tumor-node-metastasis staging system, tumor status and menopause in the TCGA database. In the univariate Cox proportional hazard model, four of the 49 differential expressed lncRNAs had important prognostic value (P < 0.05; Table 3). Further analysis using multivariate Cox regression showed that only two lncRNAs were important and independent biomarkers for OS in CC patients: ILF3-AS1 and RASA4CP (P < 0.05; Table 3). These two specific lncRNAs were positively correlated with OS (log-rank P < 0.05; Figs. 1B and 1C). Using ROC curve analysis, the AUC values for ILF3-AS1 and RASA4CP were calculated as 0.963 and 0.988, respectively (P < 0.05; Fig. S1). The AUC of these two lncRNAs together was 0.991, which was higher than that of each lncRNA taken alone (P < 0.05). We next built a risk score for predicting the prognostic value. The formula was Risk score = expILF3-AS1 × (−0.703) + expRASA4CP × (−0.576). The 289 patients were divided into low-risk (n = 144) and high-risk (n = 145) groups (Fig. 2). The survival time of patients in the low-risk score group was 1152.48 ± 104.27 days compared to 939.66 ± 83.98 days in the high-risk score group. The risk score largely predicted 5-year survival of CC patients, as the AUC upon ROC curve analysis was 0.607 (Fig. 3A). Furthermore, Kaplan–Meier curves showed that low-risk group was positively correlated with OS, and the survival time of patients in the low-risk group was longer than that of the high-risk group (P < 0.001; Fig. 3B).

The prognostic value of the two-lncRNA signature and other clinical features

Based on the results of univariate Cox proportional hazard model analysis, some of the clinical features may predict poorer survival outcomes of CC patients (Table 2). In addition, we further analyzed by the multivariate Cox proportional hazard model tumor status (P < 0.001) and the risk score (P = 0.002), which were determined to be two independent prognostic factors of CC. The Kaplan–Meier curves of the aforementioned clinical features are shown in Fig. S2. The results showed that clinical stage (P = 0.001), T stage (P < 0.001), N stage (P < 0.001), M stage (P = 0.011) and tumor status (P < 0.001) were related to OS. We also evaluated the relationship between the risk score based on the differentially expressed two-lncRNA signature and the clinical features, and the risk score showed prognostic value for predicting the status of T stage (AUC = 0.607, P = 0.028; Fig. S3). The expression levels of two differentially expressed lncRNAs in the low and high score groups in TCGA are shown in Fig. 3C. The results revealed that the expression level of lncRNA RASA4CP was significantly different between the low-risk and high-risk groups (P < 0.05).

qRT-PCR verification of the expression level of two lncRNAs in tissues and ROC curve analysis

To confirm the expression levels of ILF3-AS1 and RASA4CP in CC, we analyzed their actual expression levels in 31 pairs of newly diagnosed CC clinical samples by qRT-PCR (Table S1). The qRT-PCR results confirmed that the expression of the ILF3-AS1 and RASA4CP was downregulated, consistent with the TCGA results (Fig. 4A). The actual expression levels of the two lncRNAs in the low-risk score and high-risk score groups in the clinical sample are shown in Fig. 4B. To evaluate the specific diagnostic potential of ILF3-AS1 and RASA4CP, ROC curve analysis was performed. AUC values were 0.704 and 0.685 for ILF3-AS1 and RASA4CP, respectively, (P < 0.05; Figs. 4C and 4D). This indicates that each lncRNA could be an important biomarker for diagnosis of CC. However, the AUC value for both lncRNAs was 0.723, higher than that any single lncRNA alone (P < 0.05; Fig. 4E), suggesting improved diagnostic efficiency using the two-lncRNA signature.

qRT-PCR verification of the two lncRNAs’ expression in cervical cancer cell lines

To detect ILF3-AS1 and RASA4CP expression in CC, we further examined gene expression pattern from three CC cell lines by qRT-PCR. Compared to the immortalized cervical epithelial cell line H8, expression of ILF3-AS1 and RASA4CP was downregulated in the human CC cell lines HeLa and SiHa, respectively. The expression pattern of the two lncRNAs is consistent with the TCGA results (Fig. 4F).