Bioinformatics profiling integrating a four immune-related long non-coding RNAs signature as a prognostic model for papillary renal cell carcinoma

Acquisition of IR-lncRs

Transcriptome data and clinical data of patients with pRCC were fetched from TCGA database. Next, transcriptome data was processed to convert the data ensembl ID into gene names. Following that, transcriptome data were divided into lncRNA and mRNA. From the Immune system process M13664 and Immune response M19817 of Molecular Signatures Database, we identified 331 pRCC IRGs, of which 17 lncRNAs were validated to be the IR-lncRs by correlation analysis.

The relevance of IR-lncRs and prognosis

Based on COX Regression model, we then identified 5 IR-lncRs which were associated with prognosis (sIRlncRs), such as AP001267.3, SNHG16, AC021054.1, AC026471.3 and ADAMTS9-AS1. The relationships between these sIRlncRs and prognosis were clearly illustrated in the forest map (Figure 1).

Clinicopathologic characteristics of the high-risk group and the low-risk group

The top 4 sIRlncRs (AP001267.3, AC026471.3, SNHG16 and ADAMTS9-AS1) among the 5 sIRlncRs were included to establish the risk evaluating model, by which the pRCC samples were divided into the high-risk group and the low-risk group based on the intermediate risk score (Figure 2A). The mortality rate constantly increased with the higher risk score (Figure 2B). And with the increase of risk score, the expression levels of AC026471.3 and SNHG16 were elevated, while AP001267.3 and ADAMTS9-AS1 expressed decreasingly (Figure 2C). The survival of the high-risk group was significantly shorter than that of the low-risk group (Figure 3).

The clinical application of the IRRS and the relationships between the IRRS and clinicopathologic indicators

To investigate the relevance of the sIRlncRs and clinicopathological features of pRCC, we analyzed the correlation between the risk score and the clinical and demographic characteristics, such as age, gender, stage, T-stage, N-stage and M-stage. Under the IRRS, the scores of older patients (Figure 4A), female patients (Figure 4B), patients with advanced stage (Figure 4C), advanced T-stage (Figure 4D), advanced M-stage (Figure 4E) and advanced N-stage (Figure 4F) were significantly increased. The above results elucidate some clinical and demographic characteristics that are sensitive to the IRRS and further corroborate the clinicopathological application value of the model.

We also analyzed the relationships between the compositions of the IRRS and the aforementioned tumor characteristics which all illustrated in Table 1. The expression levels of SNHG16 and AC026471.3 were higher in female patients (Figure 5A–5B). However, the expression levels of AP001267.3 and ADAMTS9-AS1 were lower in female patients (Figure 5C–5D).

We further found the expression levels of AP001267.3 and ADAMTS9-AS1 were gradually decreased in the more advanced stage, T-stage, M-stage and N-stage, while the expression levels of AC026471.3 and SNHG16 were enhanced (Figure 6A–6D). Detection of the roles of different sIRlncRs in indicating tumor characteristics provides insight into the further discovery of biomarkers.

We then conducted the independent risk analysis, the results showed age, stage, N-stage, M-stage and risk score were significantly correlated with OS in univariate analysis (P<0.05). But in the multivariate analysis, M-stage and risk score showed more remarkable correlations with OS (Table 2). The ROC curve represented the accuracy of the model. The AUC of risk score was 0.958 (Figure 7). These results suggested the risk score was an independent prognostic factor.

Analysis of the immune status of the high-risk and low-risk population

We employed the PCA to detect the different distribution patterns between the low-risk group and the high-risk group by the immune gene sets and the genome-wide expression sets. In the IRGs set, the low-risk group and the high-risk group were observably separated with the lower immune scores in the low-risk group (Figure 8A). While we didn’t detect the significant the separation of the immune scores on the basis of the genome-wide expression profiles (Figure 8B). The results of GSEA further proved the functional annotation, with the more active immune-related responses and processes in the high-risk group (Figure 8C and 8D).

In order to verify whether the immune genome accurately reflects the status of the tumor immune microenvironment, we analyzed the relationships between the sIRlncRs and immune cell infiltration (Figure 9A–9F). We found that only B cell showed the most significant relationship with sIRlncRs (Figure 9A).

SNHG16 and ADAMTS9-AS1 were respectively high-expression and low-expression in pRCC patients especially in female with advanced T-stages

We previously explored the roles of SNHG16 and ADAMTS9-AS1 on predicting prognosis by bioinformatics methods (Supplementary Figure 1). To further verify the relationships between sIRlncRs and the clinicopathologic features, as well as discover the roles of sIRlncRs on indicating clinical prognosis, we detected the expression levels of SNHG16 and ADAMTS9-AS1 in carcinoma and adjacent tissues of pRCC patients with different genders and T-stages. Given age was insignificant associated with the five sIRlncRs, it was removed in the verification process. As illustrated in Figure 10, compared with adjacent tissues, the higher expression of SNHG16 (Figure 10A) was detected in carcinoma tissues, but ADAMTS9-AS1 (Figure 10B) showed the reversed expression level. Besides, the expression of SNHG16 (Figure 10C) and ADAMTS9-AS1 (Figure 10D) were significantly more and lower respectively in T3 and 4 stages compared with that in earlier T-stages. Female patients were detected to obtain the higher expression levels of SNHG16 (Figure 10E) and the lower expression levels of ADAMTS9-AS1 (Figure 10F) than male.

Ethics statement

Informed consent forms have been signed by all patients before this study. The research protocol has been approved by the Ethics Committee of The First Affiliated Hospital of Chongqing Medical University and is based on the ethical principles of medical research involving human subjects in the Helsinki Declaration.

Clinical renal samples

A total of 22 clinical samples of pRCC and corresponding adjacent tissues were collected from patients who received treatment in the First Affiliated Hospital of Chongqing Medical University from October 2018 to December 2019. The collected tissue samples were immediately frozen in liquid nitrogen and then stored at –80°C until RNA extraction.

Download and pretreatment of data

Transcriptome RNA-sequencing of pRCC samples were downloaded from the The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/), that contained data from 32 non-tumor tissues and 289 pRCC tissues. The clinical data about these patients was also downloaded and extracted (the OS of patients≤30 days were excluded because these patients probably died for unpredictable factors such as infection and hemorrhage). Raw data was collected to do further analyses. These data were currently updated in November 11, 2019. RNA-seq results were combined into a matrix file by a merge script in the Perl language (http://www.perl.org/). The Ensembl IDs of genes was converted into a matrix of gene symbols by the Ensembl database (http://asia.ensembl.org/index.html).

Immune-related long non-coding RNA acquisition

The Molecular Signatures Database v4.0 (Immune system process M13664, Immune response M19817, http://www.broadinstitute.org/gsea/msigdb/index.jsp) was used to explore the immune-related gene participating in the immune process. Immune related gene was used to establish the immune score of pRCC gene by GSEA. The correlation between the immune score and the expression of lncRNA in pRCC patients was analyzed by Pearson correlation analysis. IR-lncRs was identified by a criterion of |r|>0.8 and P<0.001.

Survival-related IR-lncRs

IR-lncRs associated with clinical outcomes in pRCC patients were regarded as survival-related IR-lncRs (sIRlncRs). Univariate COX analysis was used to screen sIRlncRs (p<0.01). Hazard ratio (HR) was used to specified sIRlncRs into protective and deleterious portion. These sIRlncRs were selected for follow-up study.

Establish the immune-related risk score model (IRRS)

To verify the reliability, sIRlncRs were submitted for the multivariate analysis, while the integrated sIRlncRs were still used as an independent prognostic indicator to develop the IRRS (p<0.05) (Table 3). In order to explore the heterogeneous clinical prognostic outcomes, based on the differential expression of sIRlncRs, we performed a risk score model to divide pRCC patients into the high-risk group and the low-risk group according to the median risk score which was chosen as the cutoff point. IRRS was established based on the expression data multiplied by Cox regression coefficients. The formula was as followed, [Expression level of AP001267.3 * (-1.213353)] + [Expression level of SNHG16 * (0.581303)] + [Expression level of AC026471.3* (0.779272)] + [Expression level of ADAMTS9-AS1* (-0.233562)]. The values of IRRS were employed to evaluate various subtypes of pRCC patients. To further investigate the relevance of the sIRlncRs and clinicopathological features of pRCC, we analyzed the relationship between the IRRS and clinicopathologic characteristics, of which the “TNM staging method” is the most common way to describe the tumor status. The division of “T-stage” pRCC was based on the maximum diameter of tumors and extent of tumor invasion, with the bigger tumors and more extensive invasive statues in the more advanced T stages. “N-stage” reflects the lymph node metastasis conditions with more metastatic lymph nodes in the more advanced N stages. “M-stage” is distinguished according to whether the tumor exhibits distant metastasis, and advanced M stages usually represent poor tumor conditions. In addition, “stage” is a comprehensive method combining T-stage, N-stage and M-stage to divide patients with pRCC into I, II, III and IV stage.

Bioinformatics analysis

The survival ROC curve was employed to verify the prognostic performance through the survival ROC package of the R software. The survival time, survival status and risk scores of patients with pRCC were used to predict the prognosis over a 5-year period, then the ROC curve was drawn and the value of AUC was calculated. Abscissa was the false positive rate and ordinate represented the true positive rate. Principal component analysis (PCA) was displayed to demonstrate the expression levels of pRCC samples and gene set enrichment analysis (GSEA) was used to detect the different functional phenotypes between the low-risk group and high-risk group. The tumor infiltrating immunocytes were evaluated through the TIMMER database. The levels of immune infiltration in pRCC patients were downloaded, and the relationship between IRRS and immune cell infiltration was calculated.

Real-time quantitative PCR

Triazole (Invitrogen) was used to extract total RNA from tissues and cell lines under various experimental conditions according to the manufacturer's instructions. cDNA Synthesis Kit (Osaka, Japan of TaKaRa) combining with RNA (1μg) was utilized to reverse transcribed cDNA. The reaction steps were as follows: 37°C for 15min, 85°C for 5s, and quantitative polymerase chain reaction (qPCR) was performed on an ABI 7500 real-time PCR system (Applied Biosystems) using SYBR-Green method (TaKaRa). The reaction cycle conditions were performed (95°C 30s, followed by 40 cycles of 95°C for 5 s and 60°C for 34 s). Relative expression level of lncRNAs normalized to β-actin was calculated by the 2^−ΔCt method. The primer sequences are shown in Table 4. Three replicate assays were performed for each cDNA sample.

Statistical analysis

Univariate Cox regression analysis and Pearson correlation analysis were used to identify the target IR-lncRs. Kaplan-Meier curve was used to evaluate the OS between low-risk group and high-risk group. Univariate and multivariate Cox regression analysis were used to verify the independent prognostic factors for pRCC patients. All statistical analysis was conducted using SPSS21.0 software (SPSS Inc, Chicago, IL) and GraphPad Prism5 (GraphPad Software Inc, La Jolla, CA). Varieties in clinical parameters were determined using independent t-tests. P<0.05 was considered significantly statistical difference.