ldentification and validation of tumor microenvironmentrelated lncRNA prognostic signature for uveal melanoma

Chen-Lu Liao, Nan Hu, Xing-Yu Sun, Qi Zhou, Min Tian, Yang Cao, Hong-Bin Lyu

1Department of Ophthalmology, Affiliated Hospital of Southwest Medical University, Luzhou 646000, Sichuan Province, China

2Southwest Medical University, Luzhou 646000, Sichuan Province, China

3Department of Gynecology, Αffiliated Hospital of Southwest Medical University, Luzhou 646000, Sichuan Province, China

Abstract

INTRODUCTION

U veal melanoma (UM) is the most common primary intraocular tumor in adults, derived from melanocytes in pigmented tissues such as the iris, ciliary body, and choroid. The clinical manifestations of UM are generally decreased vision, visual distortion, loss of visual field, or visual impairment. However, 30% of patients had no obvious symptoms and were diagnosed by routine screening. UM’s eye treatments include laser therapy, radiation therapy, and local tumor resection. The purpose is to preserve the eyeball and residual visual function, and to avoid metastasis as much as possible. Its prognosis is poor, with about 50% of patients experiencing systemic metastases[1-2]. Once diagnosed with primary UM, about 20%-30% of patients die within 5y, and 45% die within 15y[3-4]. This poor prognosis makes it necessary and urgent to find new biomarkers for diagnosis, prognosis and presumed drug development, to improve the survival rate of melanoma patients.Despite increased awareness of the disease, overall survival has not improved. UM is highly resistant to certain single drugs or combination chemotherapy,the objective response rate is low[5-6]and once UM has metastasized, systemic chemotherapy has no significant impact[7].

At present, the efficacy of chemotherapy is limited, and more effective treatment strategies are needed. UM and cutaneous melanoma both originate from melanocytes, but their molecular mechanisms and clinical manifestations are quite different[8]. Immunotherapy improved overall survival in patients with cutaneous melanoma, but had no significant clinical benefit in patients with unresectable or metastatic UM[9-10]. Eyes are immune privileged organs, and tumors with local immune avoidance are difficult to treat[11-12].

Tumor metastasis is the main cause of death in cancer patients.However, the mechanism of tumor metastasis is extremely complicated, and there are still many problems that have not been clarified. Based on the “seed-soil” theory of tumor metastasis, tumor microenvironment (TME) plays a vital role in tumor metastasis. TME is the environment around which tumor cells live. The immune components of TME play a vital part in gene expression and clinical efficacy in tumor tissues. Studies have increasingly clarified the role of TME in immunotherapeutic response and drug resistance of various cancers and have explored its influence on disease prognosis[13-14]. In general, a TME with high immune cell infiltration has a good prognosis, while in UM the opposite is true[15-16]. This highlights the importance and urgency of an in-depth understanding of the microenvironment in UM,to stimulate thinking on novel approaches for treatment and prognosis in UM. In addition, Yoshiharaet al[17]proposed a new algorithm called ESTIMATE that uses gene expression data to infer the level of invasive stromal and immune cells and tumor purity in tumor tissue, making it possible to classify samples using gene expression data.

Long non-coding RNAs (lncRNAs) are genes that do not encode proteins in the human genome. These non-coding genes may affect normal gene expression and disease progression,making them a new class of targets for drug discovery. They play a significant role in many biological processes and have become a major research focus in cancer treatment. UM inhibition of apoptosis, overgrowth, invasion, metastasis, and poor prognosis may be associated with aberrant expression of some lncRNAs[18-19]. Therefore, the combination of TME and lncRNA may provide new avenues for research.

The current sequencing work reveals the extensive genome heterogeneity of UM, which provides valuable information for diagnosis and prognosis, and the rise of various bioinformatics methods makes it possible to optimize biomarkers. Using single-factor COX and multi-factor COX methods can screen the prognostic-related genes of UM. The Least Absolute Shrinkage and Selection Operator (LASSO) is a popular method that combines the best performance parameters to generate a simpler and easier to interpret model, thereby avoiding overfitting. It is widely used in survival analysis of high-dimensional data, the Cox proportional hazards regression model.

The aims of the present study were to estimate the infiltration pattern of TME in UM based on the expression data of the Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov), using differential analysis, single factor, multifactor and LΑSSO‐Cox Regression analysis, to find new biomarkersvialncRNA, to establish a prognosis prediction model, to look for chemical drugs that may treat UM and to provide novel ideas for the treatment and prognosis of UM.

MATERIALS AND METHODS

Data Collection and Grouping of Uveal MelanomaTCGA was used to obtain gene expression profiles and homologous clinical information of UM patients. The method of single sample gene set enrichment analysis (ssGSEA) was used to evaluate the infiltration of immune cells in samples. Twenty‐nine immune data sets were divided into high- and lowpermeability groups, based on the types of immune cells,immune-related pathways and immune-related functions of samples[20]. The ssGSEA method uses the gene set variation analysis (GSVA) function in the GSVA software package and relies on RNA-seq count data of UM. For clinical characteristics, the R package ‘tableone’ was used to conduct Chi-square tests for grouped data, and Kruskal test was used for continuous variables, the averages of which were represented by median values[21].

Grouping ValidityTumor-related effects of infiltrating stromal cells and immune cells were evaluated using the ESTIMATE algorithm based on the transcriptomic expression profile of UM, including the immune score, stromal score,ESTIMΑTE score, and tumor purity, to verify the effectiveness of ssGSEA grouping. The ESTIMATE algorithm generates the stromal score, which is used to measure tumor-related stromal infiltration, and the immune score representing immune cell invasion level. The two combine to produce an ESTIMATE score, which can be used to infer tumor purity. The expression of 29 chemokines, CD274 (PD-L1) and human leukocyte antigen family (HLA) play major parts in the genesis,development and treatment of tumors. Three indices of high and low immune cell infiltration were evaluated[22-24]. The infiltration of immune cells in two groups was estimated using the CIBERSORT method[25]. A Wilcoxon test was used to estimate the difference between each index in the low and high immunocompromised groups. Ultimately, survival analysis was used to assess the prognosis related to low and high immune cell infiltration.

Identification of Prognostic Markers of Immune-associated Genes in Uveal MelanomaDifferential analysis, LASSOCox regression analysis and univariate Cox regression analysis were used to identify the prognostic markers of UM immuneassociated genes. The R software package ‘edgeR’ was used to analyze the differential expression of RNΑ‐Seq count data in the infiltrating groups of high and low immune cells[26].The screening criteria were |log2(fold change)| 1 or higher and false discovery rate (FDR) <0.05. The data visualization package ‘ggplot’ was then used to draw the differential gene volcano map and the final outcomes were visualized. The gene annotation information was downloaded from the GENCODE website (https://www.gencodegenes.org/) to identify and extract the lncRNA. A single-factor Cox regression analysis was then carried out on the lncRNA differential expression gene, to look for the prognostic gene. LASSO-Cox regression was used to analyze the prognosis difference due to the lncRNA gene, 100 times cross validation was performed on the set parameters to avoid over fitting, and minimize complexity of the model. The LASSO regression was analyzed using the R package ‘glmnet’[27]. For screened genes, the best cutoff values were used to distinguish between the two groups in terms of survival analysis as an indication of their effects on patient prognosis. The risk score was calculated based on lncRNA gene expression associated with prognosis, using the equation in which Coefiwas the LΑSSO coefficient of the gene and Expiwas the relative expression of the gene in the patient’s signature.

Patients were grouped into low- and high-risk groups based on the median risk score. Receiver operating characteristic(ROC) curve and Kaplan-Meier (KM) survival curve were used to appraise the predictive ability of the prognostic model. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were then used to analyze the differential expression of mRNA genes in the high and low immunodeficiency groups.PandQvalues <0.05 were considered statistically significant.

Clinical Correlation Analysis of Prognostic Markers Associated with Immune Cell Infiltration in Uveal MelanomaUnivariate and multivariate Cox regression analyses were used to develop the prognostic model as an independent prognostic index. Univariate Cox regression was applied to analyze the clinical characteristics of the TCGA cohort, and low and high risk indicators were constructed based on 8 lncRNA genes associated with immune cell infiltration. Factors withP-values <0.05 were incorporated into the multivariate Cox regression model.

Table 1 Correlation between clinicopathological features and 8 genetic features of UM

Prediction and Functional Analysis of lncRNA Target Genes Associated with PrognosislncRNA in the prognostic model was analyzed, and mRNΑs with correlation coefficients>0.6 were selected as the target gene. The predicted target gene was entered into the STRING database (https://string-db.org/),the protein-protein interaction (PPI) network was obtained,and Cytoscape 3.7.2 software was used to visualize the results.The cytoHubba plugin was then used to recognize the hub genes based on the number of degrees. Finally, the target gene pathway was predicted using GO and KEGG analysis.

Certification of Candidate Small Molecule DrugsIn the Connectivity Map (CMap) database, differentially expressed gene maps are used to predict the relationships between small molecule drugs and diseases. The positive score is the same as the reference gene expression profile, while the negative score may be the opposite. The database was used to screen enrichment fractions <-0.85 andPvalues <0.05, and drugs with negative scores were considered therapeutic molecules.

Statistical AnalysisAll statistical analyses were performed using R software v3.6.1. Unless otherwise stated,P<0.05 were considered statistically significant.

RESULTS

Patient CharacteristicsBased on the risk scores of the 8 lncRNA gene prognostic models we constructed, 80 TCGA cohort samples were allocated to low and high risk groups for statistical analysis of clinical characteristics (Table 1). The clinical information included age, gender, the tumor, lymph nodes, and metastasis (TNM) stage, eye colour, tumor basal diameter, tumor thickness, mitotic count, mutation count.

Figure 1 Construction and verification of high and low immune infiltration grouping of UM A: The samples were divided into high and low immune cell infiltrating groups; B, C: The stromal score, immune score, ESTIMΑTE score, and tumor purity were compared between the two groups; D, E, F: The HLΑ family, CD274 (PD‐L1) expression and 29 chemokines were all significantly higher in the high immune infiltration group than in the low immune cell infiltration group; G: CIBERSORT was used to show the difference in the proportion of each type of immune cell between the two groups; H: Survival analysis showed immune infiltration level had a significant impact on patient survival.

Construction and Verification of Immune-associated Subgroups for Uveal MelanomaAccording to transcription data, 80 samples obtained by TCGA were divided into high immune cell infiltration (HICI) and low immune cell infiltration (LICI) groups (n=8 and 72 respectively; Figure 1A). According to the ESTIMATE algorithm, tumor purity was relatively high in the LICI group, but the ESTIMATE score,stromal score and immune score were all lower in the LICI than the HICI group, suggesting validity of grouping (Figure 1B-1C). In addition, the expression of HLA family and CD274(PD-L1) were significantly lower in the LICI than the HICI group, and chemokine expression was also noticeably lower in LICI (Figure 1D-1F). These results suggest that the group with high immune penetration can better mediate the recruitment of multifarious immune cells in tumor, but their prognosis is poor (Figure 1H). In addition, there were differences in immune cell infiltration between the two groups (Figure 1G).These immune cells included B cells na?ve, B cells memory, T cells CD, T cells CD4 memory resting, T cells CD4 memory activate, T cells follicular helper, T cells regulatory (Tregs),T cells gamma delta, monocytes, macrophages M1, dendritic cells activated. Taken together, these results indicate that the groups were valid as a basis for subsequent analysis.

Construction and Evaluation of 8 Immune-related lncRNA Genes in Uveal MelanomaA total of 2393 differentially expressed genes were identified and met the criteria for the low and high immune cell infiltration groups. Among 269 differentially expressed lncRNAs, 69 were up-regulated and 200 were down-regulated, showing different distributions in the volcanic map (Figure 2A). It is revealed that this pathway was enriched in immunorelated and cell cycle related by GO and KEGG analysis of differentially expressed mRNA gene,indicating the accuracy of our grouping and analysis (Figure 2E-2F).

Univariate Cox analysis of lncRNA genes with differential expression identified 186 genes associated with prognosis.These prognostic genes were included in LASSO-Cox regression for collinearity removal analysis (Figure 2B-2C).Survival analysis of 8 lncRNA genes in the selected prognostic model showed that they were influential in prognosis (Figure 2D) with the following risk scores: -0.078×expression of ZNF667.AS1; -0.115×expression of ZNF350.AS1;-0.151×expression of LINC02572; -0.352×expression of LINC02367; -0.397×expression of ACVR22.AS1;-0.582×expression of DKFZP434A062; 0.020×expression of CYTOR; 0.277×expression of LINC01615. Risk score analysis showed that the high-risk group had lower survival rate (Figure 3A). The ROC curve indicated high predictive power of the prognostic model (Figure 3B).

Figure 2 Construction and evaluation of immune-related lncRNA genes in UM A: The volcano plot shows the differential expression between the high and low immune‐cell infiltration groups; B, C: In least absolute shrinkage and selection operator (LΑSSO) regression, the cross validation selects the minimum error corresponding to the blue dotted line as lambda value for the prognostic gene selection; D: Survival analysis curves of 8 prognostic lncRNΑs; E, F: GO and KEGG analyses of differentially expressed mRNΑs.

Regression Analysis of Risk and Clinical CharacteristicsUnivariate Cox regression analysis of the risk indicators and clinical characteristics of the 8 lncRNA gene constructs showed that age, TNM stage, tumor base diameter and low and high risk indices had significant prognostic value. These indicators were then subjected to multivariate Cox regression analysis. Risk indicators were significant in both multivariate and univariate Cox regression analyses, indicating that they could be considered as separate prognostic factors for UM(Tables 1, 2).

Protein-Protein Interactions Network Construction,Central Gene Screening and Prediction of lncRNA Target Gene Function IdentificationAssociations between relevant genes (8 lncRNA gene and mRNA) and immune cell infiltration were calculated using Pearson correlation analysis. Target genes (n=853) were obtained using |r|≥0.6 as a criterion (Figure 4). These core genes are highly correlated with immunity and play a major part in the initiation and progression of tumors, including UM. GO and KEGG analysis were then conducted on the selected target genes, and the results confirmed that the 8 lncRNΑ genes have an important role in immune cell penetration and strong generalization ability (Figure 3C-3D).

Screening of Corresponding DrugsUsing CMap, 853 differentially expressed target mRNA genes were analyzed to screen small molecule drugs and three candidates with possible therapeutic value were identified, including W-13(C14H18Cl2N2O2s), Ah-6809 (C17H14O5), and imatinib(C29H31N7O) (Table 3).

Table 2 Associations with overall survival and clinicopathologic characteristics in TCGA patients

DISCUSSION

Figure 3 Evaluation of prognostic indicators in UM A: Median risk scores were used to generate high and low risk groups for survival analysis; B: The ROC curve shows the predictive power of prognostic indicators; C, D: GO and KEGG analyses were performed to determine the differential expression of 8 prognostic lncRNΑ target genes.

Table 3 Potential small molecule drugs

UM is the commonest primary malignant ocular neoplasm in adults and has a high mortality rate. TME, as the ambient environment of tumor cells, is vital for regulating and influencing tumor behavior. UM differs from several other tumor types. The infiltration of lymphocytes indicates a poor prognosis and is often associated with metastasis[28].Approximately 4/5 of the transcriptions in the human genome are non-coding genes of protein, including lncRNA[29]. Based on these important characteristics of lncRNA and TME, we constructed the TME of UM using its gene expression profile and conducted further research on its immune permeation.

Figure 4 Prediction of target genes of 8 lncRNA prognostic indicators (A) and the protein-protein interactions (PPI) network was constructed including 8 target genes of lncRNA, and the core genes were screened (B).

In our research, we used TCGA data. According to 29 immune gene sets, the samples were divided into HICI and LICI groups.After that, the two groups were screened out the 8 lncRNA genes most related to the prognosis through difference analysis,univariate Cox regression analysis and LASSO-Cox regression analysis. The eight prognostic gene models we constructed can be used as reliable biomarkers for the prognostic evaluation of UM patients and can be used as independent factors for the prognostic evaluation of patients. Αs expected, high infiltration was associated with poor prognosis.

Eight TME-related lncRNAs in our prediction model have been linked to disease and some have been used as prognostic markers, such as ZNF667-AS1 for colon cancer, ACVR2BAS1 for endometrial cancer, and Linc01615 for squamous cell carcinoma in the head and neck[30-32]. According to previous research, ZNF667.AS1 is abnormally expressed in some cancers, linked to the malignant phenotype of various tumors and can be used as a biomarker and a potential therapeutic target[33]. Similarly, ACVR2B-AS1 may be a prognostic factor for breast and liver cancer[34-35]and Cytor is a regulator of Wnt/catenin signal and a latent therapeutic target for colon cancer[36]suggesting that a range of genes play an important role in prognosis and treatment.

We identified the target genes of 8 lncRNΑs using correlation analysis and carried out GO and KEGG analysis on those genes. Interestingly, a high degree of similarity was found between the results of this and the differential analysis,indicating reliability of our screening approach, and suggesting that UM may be evaluated with a focus on a small number of genes.

Using the CMap database and the results of differential analysis, we identified the small‐molecule drugs W‐13, ΑH‐6809,and imatinib as potentially effective in the treatment of UM.As a membrane permeable calmodulin inhibitor, W-13 binds well with phospholipid membrane[37].In-vitroexperiments showed that W-13 can inhibit colony formation of breast cancer cells[38]. AH 6809 is an EP2 receptor antagonist and may plays a significant role in maintaining epithelial barrier function[39]. Furthermore, Zhouet al[40]found that retinal neovascularization may be managed through the inhibition of PDGFR‐α and PDGFR‐β using imatinib, and abnormal activation of the receptor tyrosine kinase Axl pathway has been found to increase the survival rate of UM cells[41]. Thus,all three drugs may be effective in the treatment of UM.

Our study provides a deeper understanding of the TME of UM and reveals prognostic indicators of UM. And it expands our current understanding of the role of lncRNA related to immune cell infiltration in the occurrence of UM. In addition, our model has a high prediction efficiency. The 8 lncRNA genes in the model are closely related to the occurrence and development of various diseases and may be favored by clinicians and researchers. Our analysis is limited to data from TCGA, and further research is needed to verify these findings. In particular,empirical evidence is needed to verify the importance of the genes and the effects of the drugs identified here. Our study provides important clues to elucidate TME in UM and may lay the foundation for future developments in the treatment and prognosis of this disease.

ACKNOWLEDGEMENTS

Foundations:Supported by Shanghai Key Laboratory of Fundus Diseases, 2017 (No.01030); Luzhou Southwest Medical University, Municipal Department Level(No.2017LZXNYD-J01).

Conflicts of Interest: Liao CL,None;Hu N,None;Sun XY,None;Zhou Q,None;Tian M,None;Cao Y,None;Lyu HB,None.