Bioinformatics profiling intergrating a three immune-related long non-coding RNA signature as a prognostic model for clear cell renal cell carcinoma

Background: Immune related long non-coding RNAs (IRlncRs) plentiful in immune cells and immune microenvironment (IME) are potential in evaluating prognosis and assessing the effects of immunotherapy. A complete and meaningful IRlncRs analysis based on abundant clear cell renal cell carcinoma (ccRCC) gene samples from The Cancer Genome Atlas (TCGA) will provide insight in this field. Methods: Based on the TCGA dataset, we integrated the expression profiles of IRlncRs and overall survival (OS) in the 611 ccRCC patients. The immune score of each sample was calculated based on the expression level of immune-related genes and used to identify the most meaningful IRlncRs. Survival-related IRlncRs (sIRlncRs) was estimated by calculating the algorithm of difference and COX regression analysis in ccRCC patients. Based on the median immune-related risk score (IRRS) developed from the screened sIRlncRs, the high-risk and low-risk components were distinguished. Functional annotation was detected by gene set enrichment analysis (GSEA) and principal component analysis (PCA), and the immune composition and purity of the tumor was evaluated by microenvironment cell population records. The expression levels of three sIRlncRs were verified in various tissues and cell lines. Results: A total of 39 IRlncRs were collected by Pearson correlation analyses among immune score and the lncRNA expression. A total of 7 sIRlncRs were significantly associated with the clinical outcomes of ccRCC patients. Three sIRlncRs (ATP1A1-AS1, IL10RB-DT and MELTF-AS1) with the most significant prognostic values were enrolled to build the IRRS model in which the OS of in the high-risk group was shorter than that in the low-risk group. The IRRS was identified as an independent prognosis factor and correlated with the OS. The high-risk group and low-risk group illustrated different distributions in PCA and different immune status in GSEA. Besides, we found the more significant expression in certain ccRCC cell lines and tumor tissues of ccRCC patients compared with the HK-2 and adjacent tissues respectively. Additionally, the expression levels of lncR-MELTF-AS1 and ATP1A1-AS1 in the carcinoma and adjacent tissues of ccRCC patients with various T-stages. We found the more significant expression in certain ccRCC cell lines and tumor tissues of ccRCC patients compared with the HK-2 and adjacent tissues respectively. Additionally, the expression levels of lncR-MELTF-AS1 and IL10RB-DT were remarkably enhanced along the more advanced T-stages, but the lncR-ATP1A1-AS1 showed the inverse gradient. The verification results suggest the specific sIRlncRs illustrating significant differences in variable risk factors probably indicate the progression and prognosis of ccRCC, and they are ergastic to be the personalized

Therefore, an increasing number of immunotherapy drugs including PD-1/PD-L1 blocking agents have been approved in the treatment of ccRCC and inhibiting immune checkpoint shown satisfied results. 4 However, parts of patients remain response poorly, emerged resistance or progression. 5 In addition to immunotherapy, other alternative therapies, such as vascular endothelial growth factor-tyrosine kinase inhibitors, also shown their efficiency and limitations including drug resistance. 6,7 Therefore, series of researchers are studying the underlying mechanisms, of which tumor immune escaping was regarded as one of the most probable reasons.
Consisting of immune cells, stromal cells together with other molecules, tumor microenvironment (TME), as a crucial regulator of gene expression, participated in the oncogenesis, development and prognosis. 8,9 Although immune-related genes (IRGs) including immune-related long non-coding RNAs (IRlncRs) in ccRCC have been explored recently, and some markers plentiful in immune cells and IME are potential in assessing and predicting the sensitivity and efficacy of immunotherapy, their practical effects on prediction of prognosis and therapeutic potential remain problematic. 10,11 Therefore, forecast the progression and prognosis of ccRCC through several novel and sensitive biomarkers might provide the more personalized guideline and more appropriate therapeutic schedule.
As a brand new class of transcripts absent of potential coding proteins, long non-coding RNAs (lncRNAs) have been demonstrated to be critically involved in the tumorigenesis, tumor progression and tumor immune response. 12 Additionally, lncRNAs significantly affect the tumor immune process including antigen exposure and recognition, as well as immune infiltration. 13 Multiple researchers identified certain dysregulated lncRNAs participating in the occurrence and progression of various tumors. LncR-SNHG1 appears to enhance the immune escaping ability in breast cancer. 14 In the context of lung cancer, lncR-NKILA was proved to display the apoptosis-promoting roles on cytotoxic T lymphocytes so that inhibited tumor immunity. 15 Besides, lncR-DILC abundant in ccRCC was identified to suppress tumor progression by stabilizing PTEN. 16 Based on the overwhelming strength and significant roles of lncRNAs on tumor immunity, their potentials of predicting progression and prognosis are being widely investigated. LncR-ZNF180-2 was examined to express significantly in advanced RCC. 17 LncR-MCM3AP-AS1 in hepatocellular carcinoma was closely associated with prognosis. 18 Although a class of studies have identified some differentially expressed lncRNAs in various tumors, the generally predictive roles of IRlncR in ccRCC remains unclear. Discovering some promising prognostic markers and investigating the underlying molecular mechanisms are eagerly anticipated.
Therefore, we designed the study to give an insight into the clinical potency of IME and IRlncRs on prognosis estimation of ccRCC. We extracted series of IRlncRs in IME forecasting poor prognosis in patients with ccRCC, and combining their clinicopathologic characteristics, we further evaluated the connections with overall survival (OS). The results shine light on clarifying the approaches and underlying mechanisms of IRlncRs on ccRCC and establish a more personalized precision predicting model for immunotherapy.

Methods
Human renal cell lines ccRCC tissues and normal adjacent tissues were collected from 50 patients admitted to the First Affiliated Hospital of Chongqing Medical University and diagnosed with ccRCC. Human normal renal cell lines HK-2 and renal cancer cell lines were purchased from the American Type Culture Collection (Manassas, Virginia, USA). Cells were cultured by DMEM and 1640 medium, which were supplemented with 10% fetal bovine serum (FBS), 100 u/ml penicillin and 100 mg/ml streptomycin (Gibco, Gaithersburg, MD, USA). Incubate cells at 37 °C in 5% CO 2 . The medium was changed every 2-3 days.

Data Download And Pretreatment
A series of transcriptome RNA-sequencing data of ccRCC samples were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov/), that contained data from 72 non-tumor tissues and 539 ccRCC. Clinical data about these patients was also downloaded and extracted (the OS of patients ≤ 30 days were excluded because these patients probably died of unpredictable factors such as infection and hemorrhage). Raw data was prepared to do further analyses. These data were current updated in November 11, 2019. RNA-seq results were combined into a matrix file using a merge script in the Perl language(http://www.perl.org/). Next, the Ensembl database (http://asia.ensembl.org/index.html) was used to convert the Ensembl IDs of genes into a matrix of gene symbols.

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 specified immune-related gene participating in the immune process. Immune related gene was used to establish the immune score of ccRCC gene by GSEA. Pearson correlation analysis was used to analyse the correlation between immune score and the expression of lncRNA in ccRCC patients. IRLNRS was identified by a standard of |r|>0.8 and P < 0.001.

Survival-related IRlncRs (sIRlncRs)
IRlncRs associated with clinical outcomes in ccRCC patients were identified as sIRlncRs. sIRlncRs were selected by univariate COX analysis using R software survival packages(p < 0.01). Hazard ratio (HR) was used to clarify sIRlncRs into protective and deleterious portion. These sIRlncRs were specified for subsequent research.

Create The Immune-related Risk Score (IRRS) Model
In order to verify the reliability of the sIRlncRs, sIRlncRs has been submitted for multivariate analysis, while the integrated sIRlncRs were still used as an independent prognostic indicator to develop the IRRS model. To distinguish the heterogeneous clinical prognostic outcomes, based on the differential expression of sIRlncRs, we performed a IRRS model to divide ccRCC patients into high-risk group and low-risk group. IRRS model was based on expression data multiplied by Cox regression coefficients.  Table 1. Three assays per cDNA sample. Table 1 The primer sequences of IL10RB-DT, MELTF-AS1 and ATP1A1-AS1 Note: F primer: forward primer; R primer: reverse primer.

Statistical analysis
Univariate Cox regression analysis and Pearson correlation analysis were used to identify the target IRlncRs. Kaplan-Meier curve was used to evaluate the OS between high-risk group and low-risk group.
Univariate and multivariate Cox regression analysis were displayed to verify the independent prognostic factors for ccRCC 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.

Acquisition of IRlncRs
Transcriptome data and clinical data were fetched from TCGA database, and then 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 ccRCC IRGs, of which 39 lncRNAs were validated to be the IRlncRs by correlation analysis.
3. Clinicopathologic characteristics of the low-risk group and the high-risk group The top three sIRlncRs (ATP1A1-AS1 IL10RB-DT and MELTF-AS1) among the 7 sIRlncRs were included to establish the IRRS model, by which the ccRCC samples were divided into the high-risk group and the low-risk group (Fig. 2A). The intermediate IRRS was regarded as the critical value. The mortality rate constantly increased with the higher IRRS (Fig. 2B). And with the increase of IRRS, the expression levels of IL10RB-DT and MELTF-AS1 were elevated, while the ATP1A1-AS1 expressed decreasingly (Fig. 2C). The survival of the high-risk group was significantly shorter than that of the low-risk group (Fig. 3).

The Relationships Between Irrs And Clinicopathologic Indicators
To further investigate the relevance of the sIRlncRs and clinicopathological features of ccRCC, we

Analysis of the immune status of the low and high-risk population
Based on the immune gene sets and the genome-wide expression sets, we employed the principal component analysis (PCA) to detect the different distribution patterns between the low-risk group and the high-risk group. In the IRGs set, the low-risk group and the high-risk group were separated into two parts of which the low-risk group having the lower immune scores than the high-risk group (Fig. 6A). While we didn't detect the significant separation of the IRRS on the basis of the genomewide expression profiles (Fig. 6B). The results of GSEA further proved the functional annotation, with the more immune-related responses and processes in the high-risk group (Fig. 6C and D). tissues showed the much higher expression levels of lncR-MELTF-AS1 and IL10RB-DT than that of adjacent tissues, besides, the expression levels were more significant in T3 and 4 stages compared with that in earlier T-stages (Fig. 7B). However, the lncR-ATP1A1-AS1 illustrated the inverse gradient ( Fig. 7B).

Discussion
The increasing cognitions about the importance of immunization activities in tumorigenesis, progression and prognosis offer more inspirations in explorations about tumor immunology, such as immunotherapy and immune-related predicted markers. Although innovated and multimodal therapeutic strategies including targeted therapy and immunotherapy have provided more novel options and prolonged survival of plenty of patients with ccRCC, the curative effects of another part of patients remain unsatisfactory. 2,6,19 Multiplied studies revealed the individual variation at the genetic level should be responsible for the phenomenon. 20,21 Therefore, more explorations and understanding of ccRCC characterized by variant molecular heterogeneous, and further identifying more sensitive and effective immune prognostic indicators are above rubies.
The larger numbers of researches highlighting the vital importance of immune cells infiltration motivate us to investigate the potential of IRGs in distinguishing expression in different clinicopathologic status. 9,22,23 Although certain immune-related miRNAs, mRNAs and protein indicators have been identified to forecast treatment outcomes of ccRCC and some risk models based on differentially expressed IRGs have been established to evaluate survival, increasing numbers of researches are constantly reporting that without the mission encoding proteins, lncRNAs obtain more specificity on indicating tumor actual condition than other types of markers. 24 Besides, as the main regulating factors, IRlncRs have been detected to participate in multiplied Immunization responses and activities. 25 Given the inherent strengths of lncRNAs on cancer biological processes and the remarkable immune correlation of IRlncRs, exploring their values on predicting prognosis in ccRCC patients is eagerly anticipated.
In the present study, 611 patients with ccRCC were included in a genome-wide analysis for lncRNAs, combining with 311 IRGs screened in Molecular Signatures Database v4.0 (Immune system process M13664, Immune response M19817), 39 IRlncRs were identified eventually. We further detected the relation between the prognosis of ccRCC patients and the expression levels of the 39 IRlncRs, of which 7 sIRlncRs indicated the significant correlation with OS. Utilizing multivariate Cox and risk scoring model, we further identified 3 IRlncRs to establish a risk evaluating model which was available to distinguish ccRCC patients into the high-risk group and low-risk group with obviously differences of OS. Due to the molecular heterogeneity, the accuracy and sensitivity of the present clinical risk factors in predicting the survival of ccRCC patients remain unsatisfying, we further validated the predicting value of the three sIRlncRs by multivariate analysis. We found the three sIRlncRs were independent of traditional risk factors and molecular characteristics. We also used the principal component analysis (PCA) method to study the different distribution patterns between the low-risk group and the high-risk group based on the genome-wide IRlncRs and IRGs expression set. According to the IRGs set, the low-risk group and the high-risk group tended to be divided into two parts, with the low-risk group having lower immune scores than the high-risk group. When PCA was performed on the basis of genome-wide expression sIRlncRs and IRGs profiles, the immune statues of these groups shown no significant separation. GSEA was employed to further verify the functional annotation, and we found the more abundant immune-related responses and processes in the high-risk group.
Therefore, immune-related scores are bound up with the immune status of ccRCC, with higher scores indicating the poor prognosis. These findings suggest that the risk evaluating scores based on the three sIRlncRs can contribute to identify the high-risk patients from patients with the same clinical or molecular characteristics, thereby realize individualized and appropriate therapeutic strategy. suggest the specific sIRlncRs illustrating significant differences in variable risk factors probably indicate the progression and prognosis of ccRCC, and they are ergastic to be the personalized molecular biomarkers to assess the infiltration of immunocyte and predict treatment outcomes.
Although we validated the effects of some sIRlncRs on predicting prognosis and further verified the expression level of lncR-MELTF-AS1, IL10RB-DT and ATP1A1-AS1 in certain ccRCC cell lines and tumor tissues, there are still some limitations in the study. Firstly, in the evaluating model, the assessments of proteomics, metabolomics and immunogenomics shouldn't be ignored to complete the more comprehensive analysis. Then, the clinical application values of these sIRlncRs remain undefined.
Thirdly, in addition to lncR-MELTF-AS1, IL10RB-DT and ATP1A1-AS1, other sIRlncRs enrolled in the IRRS model also need to be detected.

Conclusion
In the present manuscript, we comprehensively analyzed and verified the roles of the IRlncRs on indicating clinical outcomes of ccRCC. Certain sIRlncRs with remark clinical relevance were identified to be valuable for the latent monitoring and prognosis values for ccRCC patients. Our results develop a reliable and referable risk evaluating model to assess prognosis of ccRCC and provide the novel insight into immunological researches and treatment strategies.