Immuno-Genomic Landscape of Cutaneous Melanoma Identifies a Prognostic and Immunotherapeutically Relevant Gene Signature

doi:10.21203/rs.3.rs-617208/v1

Download PDF

Primary research

Immuno-Genomic Landscape of Cutaneous Melanoma Identifies a Prognostic and Immunotherapeutically Relevant Gene Signature

https://doi.org/10.21203/rs.3.rs-617208/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

The tumorigenesis of Skin cutaneous melanoma (SKCM) is still a mystery. Our study conducted a comprehensive analysis of the immune cell infiltration in the TME of SKCM. Based on the differential expression genes in the cluster grouped by the immune infiltration status, a set of hub genes related to the clinical prognosis of SKCM and tumor immune infiltration were explored.

Methods

We analyzed the immune cell infiltration in two independent cohorts, and then assessed the relationship between the internal pattern of immune cell infiltration and SKCM characteristics, including clinicopathological features, potential biological pathways and gene mutations. We further divided the three clusters of differential genes into two groups with different unique biological processes. The Signature gene-A gene set was mainly manifested as exon skipping (ES) in SKCM patients, while the Signature gene-B gene set has no obvious alternative splicing form. Subsequently, we not only analyzed the genetic variation of the two signatures, but also constructed a ceRNA regulatory network..LASSO Cox regression was utilized to find the immune infiltration signature and the risk score of SKCM.

Result

We finally obtained 13 Hub genes, and calculated the risk score based on the coefficient of each gene to further explore the impact of the high and low-risk score on the biologically related functions and prognosis of SKCM patients.The correlation between the risk score and the clinicopathological characteristics of SKCM patients indicated that the low risk score was associated with TMECluster-A classification (P <0.001) and metastatic SKCM (P <0.001). We finally obtained 13 Hub genes which showed different prognostic effects in pan-cancers. The IHC staining results showed that Ube2L6, SRPX2, IFIT2 were higher expression while CLEC4E, END3, KIR2DL4 were lower expression in 25 melanoma specimens.

Conclusion

We performed a comprehensive assessment of SKCM's immune environment and constructed a set of unprecedented immune signatures related to the immune landscape (EDN3、CLEC4E、SRPX2、KIR2DL4、UBE2L6、IFIT2), which are correlated with the different prognosis and drug response of SKCM. The immune gene signature we constructed can be used as a robust prognostic biomarker of SKCM and a predictor of immunotherapy effect.

Cancer Biology

Skin cutaneous melanoma

Prognostic biomarker

Tumor microenvironment

Gene signature

Immuno-genomic landscape

Clinicopathological characteristics

Skin cutaneous melanoma (SKCM) occurs most commonly in skin cancer, and its morbidity, mortality and disease burden are constantly rising(1). Tumor staging is significantly related to the prognosis of melanoma, so early diagnosis and treatment can bring a good overall survival (OS) rate. In the past few years, more and more people have realized that the balance condition of the immune system are determinants of the incidence and recurrence of cancer(2, 3). Preventing immune destruction has become the current mainstream view in cancer research(4). SKCM has become one of the most lethal human malignancies. As the deadliest form of skin cancer, it makes up almost 75% of skin cancer lethality and is common in the world. To make matters worse, the death rate is still increasing despite of public health campaigns in many countries. It is difficult and costly to take measures once transferred because of its perplexed etiology, such as multiple gene mutations. Despite the five-year survival of early-stage SKCM patients exceeds 95%(5, 6), the reported survival time for advanced stage melanoma barely exceeded one year(7). Today, patients with primary melanoma really require surgical resection as the first choice. However, advanced melanoma is highly aggressive, making it insensitive to radiotherapy and chemotherapy (8). There is still no effective systemic therapy to treat advanced SKCM(7). The extensive crosstalk between tumor-infiltrating immune system and tumor matrix in TME is crucial to tumor progression. Tumor-infiltrating immune cells (TIICs) exhibit tumor-promoting effects according to tumor types. In the vast majority of cancers, CD8 + T cells play crucial roles in TME, inhibiting the proliferation and invasion of malignant cells. T cell mediated immune responses to melanoma antigens have been extensively documented(9, 10). More and more evidences indicate that the intratumoral immune landscape has become an important factor affecting patient survival and immunotherapy response. And immunotherapy and targeted therapy can benefit patients with melanoma(11). However, it has been reported that ICIs block different immunomodulatory pathways, leading to an increase in CD4 + and CD8 + T cells and a decrease in regulatory T cells(12). A large number of recent studies have shown that immune system dysfunction plays a key role in the development of SKCM(13, 14), leaving many with an unmet need. However, there are still approximately 50% of patients do not derive benefit even with immune checkpoint inhibitors (ICI) (15, 16).

Considering the good effect of immunotherapy in SKCM patients, there is still an urgent need to establish the characteristics of SKCM patients based on a set of immune-specific marker genes. So far, the overall picture of immune cell infiltration in the SKCM tumor microenvironment has not been yet clarified. Therefore, it is urgent to find new prognostic biomarkers and treatment targets of SKCM. Therefore, it is possible to effectively distinguish benefiting groups by predicting the efficacy of immunotherapy, and predicting the prognosis of SKCM patients.

In this study, we integrated the TCGA-SKCM independent cohort and validated the predictive model in four other independent cohorts (GSE8401, GSE35640, GSE15605 and GSE46517) from GEO datasets, to develop and verify a new set of personalized immune signature model. In addition, immunohistochemistry (IHC) staining and qRT-PCR were used to verify the difference in the expression level of 6-gene signature in 25 frozen SKCM tissue samples and normal tissues in the The First People’s Hospital of Foshan. We also analyzed the clinical and pathological characteristics of all existing SKCM patients, somatic cell copy number variation (CNV), tumor mutation burden (tumor mutation burden, TMB), and gene variable splicing. We used a variety of bioinformatics methods to estimate the abundance of immune cells in the TME of SKCM patients, and analyzed the correlation between the genomic characteristics of the immune landscape and the pathological characteristics and prognosis of SKCM. In the end, we discovered a unique set of immune gene signatures that have not been reported. These gene signatures that affect the immune landscape may be closely related to the different prognosis and treatment response of SKCM. On this basis, we have established an unprecedented set of immune signatures, which can be used as robust prognostic biomarkers and predictors of immunotherapy effects for SKCM patients.

2.1 Acquisition and preprocessing of SKCM expression datasets

Data from two publicly available datasets were incorporated into our study. RNA-seq data were extracted from 471 patients along with clinical features from the TCGA-SKCM cohort,as well as 1 healthy sample. Clinical data were downloaded by the University of California Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/). Somatic copy number variation (CNV) data and tumor mutation burden (TMB) were also downloaded through the Genomic Data Commons (GDC) Data Portal(17). Data on alternative splicing (AS) events were downloaded from the TCGASpliceSeq database. “RCircos”(18) package was used to draw a map of genome-wide CNV analysis of SKCM patients was performed in 23 pairs of chromosomes. Somatic mutation data were downloaded in mutation annotation format (MAF) form, and maftools package(19) was used for visualization. We obtained SKCM microarray data from two GeneChips (GPL570[HG-U133_Plus_2]Affymetrix Human Genome U133 Plus 2.0 Array and GPL96-57554 [HG-U133A] Affymetrix Human Genome U133A Array) in GSE84014(20)༌GSE356405(21)༌GSE156056(22) and GSE46517⁽²³⁾ of Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo) as the validation sets. The limma(24) package in R was used for gene expression normalization. The sva(25) package (3.20) was used to correct for plate batch effects.

2.2 tumor-infiltrating immune cells (TIICs) analysis and clustering of samples

CIBERORT(26), an algorithm that quantifies the proportion of TIICs with 547 signature genes for inferring representation in bulk tumor transcriptomes, was used to evaluate the cell components of TIICs. Proportions of 22 types of TIICs for the cohort from TCGA were estimated. The cluster analysis was performed using the unsupervised hierarchical clustering method, an algorithmic approach, that groups individuals with similar observations based on Euclidean distance method. Three immune microenvironment subtypes were defined via ConsensuClusterPlus package(27). The procedure was repeated 1000 times to make the stratification more stability.

2.3 Differential gene expression pattern clustering analysis

In order to determine the genes associated with TME cell infiltration patterns, we divided patients into three different types of patterns according to TME, which were defined as TMECluster-A, B, and C. Using limma package was used to analyze the DEGs between these three groups of SKCM patients, and the significantly different DEGs were defined as genes with an absolute log value of fold-change > 1, and FDR < 0.05. The overlapping DEGs among three group sets were further analyzed as specific genes by the VennDiagram(28) package. The k-means clustering algorithm, which is an unsupervised clustering method was used to cluster these specific genes into meaningful groups in GEO datasets according to the expression of specific genes in TMEcluster-ABC. Meanwhile, based on changes in gene expression, specific genes were divided into two groups: Signature genes-A and B.

2.4 RNA sequence expression analysis in GEPIA

GEPIA (http://gepia.cancer-pku.cn/index.html) is an online data processing webpage with RNA sequence expression information of 9736 tumors and 8587 normal tissue samples(29). In our study, differential expression of hub genes between tumor and normal tissues, correlation analysis, and the survival predication were explored using the GEPIA database. We use student’s t test for analyzing the correlation between the expression and clinicopathologic feature. P value < 0.05 indicates statistically significance.

2.5 Gene set enrichment analysis (GSEA)

Gene Ontology (GO) analysis was performed to illustrate the unique biological signifcance of gene expression signature genes. Gene function were categorized in three series: cellular components (CC), molecular functions (MF) and biological processes (BP). Crucial pathways were identifed by KEGG pathway. GO annotation were visualized with the R package clusterProfiler(30). GSEA(31) is a calculation method to identify the potential biological mechanisms between two biological states. GSEA was conducted in the molecular signatures data-base (MSigDB) (http://software.broad institute.org/gsea/index.jsp), which provided hallmark gene sets to predict biological processes between normal and SKCM samples.

2.6 Evaluation of patient biological characteristics

We further analyzed the correlation between different groups and some biologically related processes. A set of gene sets to store genes related to certain biological processes were collected, including immune checkpoints, antigen processing, CD8 + T cell, epithelial-mesenchymal transition (EMT) markers, such as EMT1, EMT2 and EMT3, angiogenesis characteristics,pan-fibroblast TGF-β response characteristics, WNT characteristics, DNA damage repair, mismatch repair, nucleotide excision repair, DNA replication and antigen processing and presentation.

2.7 Establishment of Immune characteristic model and clinical prediction model

We adopted a two-step method to establish a signature-based risk score. We constructed a prognostic model. Univariate Cox regression was used to analyze the impact of signature genes on prognosis of SKCM patients according to the cutoff with P < 0.05. We sought to validate the potential prognostic hub genes to enhance risk profiling after surgery and to define new targets in prediction models for SKCM patients. In order to prove the significance of risk score combined with clinicopathological characteristics for personalized evaluation of patient prognosis, we first tested the expression correlation of hub gene in different tumors in the TCGA database, and the predictive ability of risk score for the prognosis of patients with different tumors. A least absolute shrinkage and selection operator (LASSO) regression is a regularized approach commonly used for high-dimensional predictors selection. A system of risk score was established by LASSO Cox proportional hazards model to identify gene signature to predict the SKCM overall survival. A predictive score was developed by using the weighted sum of genes, with the coefficients of LASSO regularization.

riskS𝑐𝑜𝑟𝑒 = ∑𝐶𝑜𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑡 (ℎ𝑢𝑏 𝑔𝑒𝑛𝑒𝑖) ∗ 𝑚𝑅𝑁𝐴 𝐸𝑥𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛 (ℎ𝑢𝑏 𝑔𝑒𝑛𝑒𝑖)

Factors related with overall survival (OS) and clinical pathological characteristics were evaluated respectively by univariate and multivariate analysis with Cox and logistic regression. Variables with a value of p < 0.05 in the multivariate analysis were included in the prognostic model. The performance and discriminative ability was assessed using Harrell’s concordance index. Nomograms were constructed to predict the 3-year, 5-year, and 10-year survival rates of SKCM patients based on predictive models with identified prognostic factors. Calibration is the predictions from the nomogram compare to the observed outcomes.

2.7 Statistical analysis

All statistical analysis was applied by R version 3.6.2. Differences in continuous variables between two groups were estimated by independent Student's t test, and the differences between non-normally distributed variables were performed using the Mann–Whitney U rank sum test. (Wilcoxon rank sum test). Pearson’s χ2 test or Fisher’s exact test were used for categorical variables. used to calculate the correlation coefficient between different genes. P-values and Hazard ratios were obtained from univariate Cox proportional-hazards regression models using the R package survival(32). and the log-rank test is used to evaluate the significance of the difference in survival time between the two groups. Receiver operating characteristic(ROC) curve analysis was conducted using the “pROC” package(33) to evaluate prognostic capabilities of different risk models and the time-dependent AUC value. Univariate and multivariate Cox analyses were used to determine independent prognostic factors. All statistical tests were two- sided and statistical significance was set at p < 0.05.

2.8 Immunohistochemistry (IHC)

2.8.1 Patients tissue specimens

Tumor tissues were collected from 25 melanoma patients of the Han Chinese group from 2015 to 2020. Patient-informed consent was obtained and approved by the The First People's Hospital of Foshan Subject Review Board.

2.8.2 Immunohistochemistry (IHC) staining

Paraffin embedded tissues were sectioned at 4µm for IHC. Antigen retrieval was performed by incubating the samples in a citrate buffer (pH 6.0) for 15 min at 100℃ in a microwave oven and naturally cooled to room temperature. After blocked with a mixture of methanol and 0.75% hydrogen peroxide, sections were incubated overnight with appropriate dilatation of primary antibodies (CLEC4E, 1:500, Sigma; END3, 1:500, Sangon Biotech; IFIT2, 1:500, proteintech; SRPX2, 1:600, proteintech; Ube2L6, 1:500, Abcam; KIR2DL4, 1:500, Abcam) followed by incubation with a secondary antibody either conjugated with HRP (goat anti-rabbit, 1:500, Cell Signaling Technology). Sections were washed three times with PBS for incubation with AEC (ZSGB-BIO).

2.8.3 Histologic scoring and analysis

All specimens were examined by the cross-product (H score) of the percentage of tumor cell staining at each of the 3 staining intensities. Intensity of immunopositivity was scored as follows: None, 0; weak, 1; moderate, 2; and strong, 3. For example, a particular tumor may have 50% cell staining at intensity = 1 and 50% of cell staining at intensity = 3, for a combined H score of 200 [(50×1) + (50×3) = 200], which yielded a range from 0 to 300. The final score was graded by H score as follows: Low, H score 0-100; Moderate, H score 101–200; and High, H score 201–300.

3.1 Immune infiltration analysis related to SKCM patients

We obtained the gene expression data of SKCM patients from TCGA database, and 22 different immune cell infiltrations in each sample were analyzed through the CIBERSORT algorithm (Fig. 1A). The TME cell network depicts the interactions between tumor immune cells (Supplementary Fig. 1). Comprehensive status of cell lineage and its impact on overall survival of patients with SKCM were also analysed (Fig. 1B). The TME cell network demonstrated that Cell cluster − B, C and D were positively correlated with each other, including T cells CD8, B cells naive, Plasma cells, T cells CD4 naive, NK cells activated, Monocytes, Dendritic cells resting, T cells regulatory (Tregs), Monocytes and T cells CD4 memory activated. Moreover, there was also a significant positive correlation among immune cells in Cell cluster-A, such as gamma delta T cells, NK cells resting, Neutrophils, Mast cells activated, Macrophages M0, Eosinophils and Dendritic cells activated. Meanwhile, Macrophages M0, T cells CD4 memory activated, Monocytes were performed with a significant reverse correlation. NK cells resting and NK cells activated showed a same significant negative correlation. We also performed an exploratory analysis to measure the survival benefit and potential risks. The results showed that T cells CD4, T cells CD8, NK cells activated, Tregs, Dendritic cells and Macrophages M1 associated with shortened OS while NK cells resting, Neutrophils, Macrophages M0 and Dendritic cells activated associated with prolonged OS. In order to construct the best number of clusters and classifierion, we used the ConsensusClusterPlus package to assess the stability of the clustering structure and divided SKCM patients into TMECluster-A, B and C. Unsupervised hierarchical clustering was used to analyze the normalized immune cell fractions. Heatmap showed the correlation between the infiltration abundance of 22 type of immune cells and the immune scores in the three groups (Fig. 1C). The results showed that patients with TMECluster-A had higher immune scores, while patients with TMECluster-B and C were mainly related to tumor purity and stromal scores. In addition, PCA analysis results show that based on the expression data of SKCM patients, the three TMECluster groups can be significantly distinguished (Fig. 1D). The survival analysis related to the TME phenotype showed that TMEcluster-C (N = 92) was associated with better prognosis (log-rank test, P < 0.001) (Fig. 1E).

Subsequently, we further selected and analyzed the gene expression profile data and immune cell infiltration abundance of SKCM samples in the GEO database. Also use the ConsensusClusterPlus package to evaluate the stability of clustering (Fig. 2A). Subsequently, the heatmap showed that patients with TMECluster-A showed higher immune scores (Fig. 2B), which is consistent with the results in the TCGA database. This was in line with previous studies reported that melanoma is a highly immune dependent malignant tumor.

3.2 Construction of TME-related signature gene in SKCM patients

TME has a crucial impact on tumor epigenetics, metastasis, and immune escape. We used limma package to analyze the 486 differentially expressed genes among the three groups of TMECluster in the TCGA database to determine the potential biological characteristics of different TME phenotypes. In the same way, 330 DEGs were analyzed in the GEO datasets. Unsupervised clustering using these DEGs provided three distinct clusters: GeneCluster-A,B and C (Fig. 3B). Simultaneously, according to the differential expression of DEGs in different clusters, it was further divided into Signaure gene-A and B. GO enrichment analysis showed that Signature gene-A and B exhibit different unique biological processes. Signature gene-A was associated with the overexpression of immune-activated genes (Fig. 3C), while Signature gene-B showed up-regulation of genes related to stromal and transmembrane receptors (Fig. 3D). There were significant differences in the infiltration of TME cells (Fig. 3E) and the enrichment of related biological pathways (Fig. 3F) between the three GeneClusters, which was consistent with the previous functional enrichment results.

3.3 CeRNA regulatory network construction and signature gene expression in SKAM patients

We further analyzed the genetic variants of Signature gene-A and B, including single nucleotide polymorphisms (SNPs), copy number variations (CNVs), and alternative splicing patterns. The differences in the frequency of mutations in SKCM driver genes were explored between the two clusters. SNP analysis shows that there were obvious mutations in both gene sets. Somatic SPTA1 mutations were more frequent in TMEcluster-A (Fig. 4A). And TNN was the most frequently observed somatic mutations in TMEcluster-B (Fig. 4B). Moreover, the study of CNV change frequency showed that CNV changes are common in two groups, and most of them were concentrated on the amplification of copy numbers. Then, we showed the positions of these CNV changes on the chromosome in Fig. 4C and D. In the alternative splicing analysis, the Signature gene-A was mainly manifested as exon skipping (ES) in SKCM patients (Fig. 4E), while the Signature gene-B had no obvious alternative splicing forms (Fig. 4F). In addition, based on the gene expression of Signature gene-A and B, we screened out LncRNAs that may be related to immunity with correlation coefficients > 0.4 and P < 0.05. Then, the starBase database (http://starbase.sysu.edu.cn/) was used to obtain targeted differentially expressed miRNAs to construct the regulatory ceRNA network of mRNA-miRNA-LncRNA interaction (Fig. 4G).

3.4 Construction of immune-related prognostic gene signature

To better predict the impact of immune characteristics on the prognosis of patients, we have constructed a new prognostic-related risk scoring system. The Signature gene-A and B genes were incorporated into univariate Cox analyses, and 233 genes related to prognosis were obtained (P < 0.05). The LASSO Cox analyses was further used for dimensionality reduction and model construction, and finally a total of 16 hub genes were included in the risk scoring model (Fig. 5A). A risk score (RS) formula was established by including individual normalized gene expression values weighted by their LASSO Cox coefficients. We calculated the risk score for each SKCM patient according to locking of coefficients in each gene signature. The high/low-risk groups were divided according to the median risk score. Kaplan-Meier analysis showed that patients with high-risk scores had a relatively poor prognosis (Fig. 5B). The time-dependent ROC curve analysis also showed that the risk score has a good predictive ability on OS (Fig. 5C). The area under the curve (AUC) of 1-year, 3-year, and 5-year OS wre 0.713, 0.694, and 0.734, respectively. In addition, the distribution of each patient's risk score, survival status, and gene expression map are shown in Fig. 5D.

3.5 Gene Set Enrichment Analysis (GSEA)

Subsequently, we analyzed the impact of high and low risk groups on the biologically relevant functions of SKCM patients. GSEA analysis showed that pathways related to metabolism and oxidative phosphorylation were mainly enriched in the high-risk group (Fig. 6A and C), while pathways related to immune response, including cytokine signaling pathway, JAK-STAT signaling pathway, and natural killer cell-mediated cytotoxicity were significantly enriched in low-risk patients (Fig. 6B and D). At the same time, the expression of TME cells (Fig. 6E), as well as some other pathways, such as angiogenesis, mismatch-related features, and stromal-related features (Fig. 6F), there were significant differences between the high and low-risk groups of SKCM patients (P < 0.05).

3.6 Correlation analysis of the risk score and clinicopathological characteristics

We assessed the correlation between the risk score and the clinicopathological characteristics of SKCM. The analysis results showed that the low-risk score was correlated with TMECluster-A classification (P < 0.001; Fig. 7A) and metastatic SKCM (P < 0.001; Fig. 7B); at the same time, it was correlated with gender (P = 0.031; Fig. 7E), Pathological stage (P = 0.0027; Fig. 7F) and T stage (P < 0.001; Fig. 7G) were significantly correlated. There was no significant correlation between the risk score and TMB, age, M stage and N stage (P > 0.05; Fig. 7C, D, H and I).

3.7 Application value of the risk score in pan-cancer

Correlation analysis showed that the expression of 13 hub genes in different tumors has a significant correlation (Fig. 8A). The heterogeneity between tumors caused the risk score have different prognostic effects on different cancers (Fig. 8B). Univariate and multivariate Cox analysis showed that the risk score is an independent risk factor for predicting the prognosis of SKCM patients (Table 3; Fig. 8C). We included the index of P < 0.05 in the multivariate Cox model to construct a nomogram to predict the OS of 1, 2, 3-year survival probability in SKCM patients.The C-indexes (0.732 (95% CI: 0.697–0.767)) of the nomogram was used to calculate the discriminative ability of the nomogram, showing a high degree of discrimination. Moreover, the calibration showed a great agreement between the 1-year, 2-year, and 3-year OS estimates with the actual observation values of SKCM patients by comparing the nomogram (Fig. 8E).

3.8 The prognostic value of hub genes in SKCM patients

We focused on the potential prognostic value of EDN3, CLEC4E, SRPX2, KIR2DL4, UBE2L6, IFIT2 as an immunoscore for melanoma patients.In the group of novel hub gene, low levels of UBE2L6、KIR2DL4、IFIT2 and CLEC4E expression were observed in tumor cells, while SRPX2 and EDN3 were highly expressed in SKCM. However, only UBE2L6 and IFIT2 have significant expression differences in the TCGA-SKCM cohort. This may require a larger sample size to verify the differential expression of these molecules in SKCM. The open online tool GEPIA was used to analyze the prognostic values of these novel hub gene. We found low expression of UBE2L6, KIR2DL4, IFIT2 and CLEC4E were significantly associated with poor prognosis in SKCM.

2.9 Ube2L6, SRPX2, IFIT2 were higher expression in IHC melanoma specimens

We evaluated these gene expressions in melanoma tissues. The IHC staining results showed that Ube2L6, SRPX2, IFIT2 were higher expression while CLEC4E, END3, KIR2DL4 were lower expression in 25 melanoma specimens.(Fig. 9)

The incidence of SKCM continues to increase every year. Tumor metastasis contributes greatly to extremely poor prognosis of SKCM patients. Although a series of treatment methods such as radiotherapy, chemotherapy, targeted therapy and immunotherapy have improved the survival rate of advanced patients(34). Especially the application of immunotherapy in SKCM has made great breakthroughs(35). However, there are still many patients exhibit resistance to novel frontier of cancer immunotherapy, and failing to obtain the expected curative effect. Therefore, there is still an urgent need to determine more therapeutic targets and prognostic biomarkers. This study conducted a comprehensive study on the immune cell infiltration, tumor immune subtypes and the prognostic value of immune cells in SKCM patients. At present, some studies have tried to combine multiple molecular events to improve the sensitivity of diagnosis and prognosis. Previous studies have also constructed prognostic markers based on immune-related genes to estimate gene expression heterogeneity(36). Therefore, the prognostic characteristics we propose here are reliable in predicting the survival of SKCM patients. More importantly, we also performed independent microarray data on the four SKCM data sets in the GEO database and 25 cases from First People's Hospital of Foshan to verify the prognostic efficacy. Cross-platform research increases the reliability of our novel prognosis model. We used the CIBERSORT algorithm to comprehensively evaluate the immune cell infiltration status of the SKCM cohort transcriptome data in the TGCA database, obtained 22 different immune cell infiltration abundances in each patient. We depicted the tumor immune cell interaction network, including the overall situation of the cell lineage and its impact on the overall survival of SKCM patients is reviewed and verified in the GEO database SKCM patient cohort. After evaluating the stability of clustering, we hierarchically clustered tumor samples based on the essential differences in immune cell infiltration patterns. The PCA analysis based on the expression profile data of SKCM patients can distinguish the three TMECluster groups well. We next analyzed the immune-related scores of the three groups, tumor purity and matrix scores. The results showed that patients with TMECluster-A had higher immune scores, while patients with TMECluster-B and C were mainly related to tumor purity and matrix scores. And in the GEO database GSE8401, GSE35640, GSE15605 and GSE46517 data sets verified the gene expression profile data and immune cell infiltration of SKCM patients, and obtained consistent clustering results. In order to have a new understanding of the relationship between the above grouping and TME phenotype, we conducted a series of survival analyses and found that TMEcluster-C (n = 92) is associated with a better prognosis (log-rank test, P < 0.001). In order to explore the potential biological characteristics of different TME phenotypes, we performed a difference analysis on the three TMEClusters in the TCGA database, and displayed the intersection of 486 differentially expressed genes (DEGs) with a Venn diagram. Based on the above differential genes, we finally obtained 330 repeated DEGs from the four data sets in the GEO database. We use unsupervised clustering to divide SKCM patients in the GEO database into three different subtypes based on the expression of DEGs, namely GeneCluster-A, B and C. At the same time, according to the expression of DEGs in different groups, they are divided into Signature gene-A and B. Subsequent GO enrichment analysis showed that Signature gene-A and B showed different unique biological processes. Signature gene-A involves the overexpression of immune-activated genes, and the gene set overexpressing Signature gene-B is mainly manifested in the upregulation of genes related to matrix and transmembrane receptors. This suggests that the Signature gene-A gene set may have an impact on the immunophenotypic landscape of TME. In addition, the analysis results show that there are significant differences in the expression of TME cells and the enrichment of some related biological pathways among the three GeneCluster groups. We downloaded the copy number variation (CNV) of the patient's somatic cells, analyzed the somatic mutation data and tumor mutation burden (TMB) of each patient, and used it to explore the signature gene-A and B. Single nucleotide polymorphisms (SNPs) and copy number variations (CNVs) in SKAM patients. And obtain the gene variable splicing information of SKCM patients in TCGA from the TCGASpliceseq website. The results of SNP analysis showed that there were obvious mutations in both gene sets, and most of the changes in CNV were concentrated on the amplification of copy number. The variable splicing analysis showed that the Signature gene-A gene set was mainly manifested as exon skipping (ES) in SKCM patients, while the Signature gene-B gene set had no obvious form of variable splicing. This result highlights the influence of the Signature gene-A gene set in the immune landscape. In addition, according to the expression of Signature gene-A and B gene sets, we screened out the more relevant LncRNAs as immune-related LncRNAs, and constructed a ceRNA regulatory network of mRNA-miRNA-LncRNA interaction. In order to better predict the impact of immune characteristics on the prognosis of patients, we have constructed a new prognostic-related risk scoring system. We incorporated Signature gene-A and B gene sets into the single-factor Cox model, and obtained 233 prognostic-related genes. The subsequent LASSO Cox analysis was further used for dimensionality reduction and model construction, and finally 13 hub genes were obtained. Among them, the risk score corresponding to each sample is calculated based on the correlation coefficient of hub genes. Survival analysis shows that patients with high-risk scores have a relatively poor prognosis. ROC curve analysis. The score has a good predictive ability for OS. The area under the curve (AUC) of 1-year, 3-year, and 5-year OS are 0.713, 0.694, and 0.734, respectively. We further used GSEA analysis to evaluate the correlation between high and low risk groups and biological characteristics. The results showed that pathways such as metabolism-related pathways and oxidative phosphorylation were mainly enriched in the high-risk group; while immune response-related pathways, including cytokine signaling pathway, JAK-STAT signaling pathway, and natural killer cell-mediated cytotoxicity, etc. were significantly enriched in low-risk patients. This was consistent with the reports of SKCM(37). At the same time, the expression of TME cells, as well as some other pathways, such as angiogenesis, mismatch-related characteristics, and matrix-related characteristics, are significantly different between the high and low-risk groups. Finally, we assessed the correlation between the risk score and the clinical classification and pathological characteristics of SKCM patients, and found that low risk scores are often associated with patients' TMECluster-A classification (P < 0.001) and metastatic SKCM (P < 0.001). In different tumors, 13 hub genes showed different prognostic effects. Univariate and multivariate Cox analysis results show that risk score can be used as an independent risk factor for predicting the prognosis of SKCM patients. We selected statistically significant indicators in the multivariate Cox model to construct a nomogram to predict the OS of patients with SKCM. C-Index indicates that the model has a high degree of discrimination (0.732 (0.697–0.767)). Validation of the calibration curve exhibited good concordance between the estimated values and actual probability.

Ubiquitin/ISG15-conjugating enzyme E2 L6 (Ube2L6) has been reported to promotes insulin resistance, and hepatic steatosis and is related to the cisplatin resistance.(38, 39) In lung cancer, SRPX2 is invasive by upregulating FAK/SRC/ERK pathway and also could lead to pancreatic cancer drug -resistance in PI3K/AKT. Also, NFATc3/SRPX2 axis participant in human embryonic stem cell differentiation, indicating SRPX2 play a critical role in pan-cancers.(40–42) Interferon-induced protein with tetratricopeptide repeats (IFIT) protein may be a new therapeutic target for cancer therapy and can be used as a prognostic marker for cancers like glioblastoma and pancreatic cancer.(43) In previous studies, CLEC4E is generally related to the tuberculosis and could be used to constrain tuberculosis through autophagy against drug-resistant.(44, 45) Killer cell immunoglobulin-like receptor 2DL4 (KIR2DL4) is expressed by NK cells and some reports suggested it to be an intervene for Cancer Immunotherapy.(46) In our IHC, some of those relationships are verified in SKCM patients. There were not many studies about the relationship between SKCM and those hub genes, so our study innovatively came up with their underlying relationships and find the novel important signature.

In short, we have developed and verified an unprecedented set of effective prognostic markers based on immune infiltrating cells, which have certain potential application value in predicting the clinical prognosis of patients with SKCM and the benefit of immunotherapy. This study provides a systematic view of the immune-related characteristics in SKCM and suggests their good prognostic performance.

Ethics approval and consent to participate

Patient-informed consent was obtained and approved by the The First People's Hospital of Foshan Subject Review Board.

Consent for publication

All authors are aware of and agree to the content of the paper and their being listed as a co-author of the paper.

Acknowledgments

Not applicable.

Authors' contributions

Ronghua Yang, Yidan Sun and Tianqi Chen participated in original draft preparation; Xiaobing Pi, Shijian Xiang collected and analysed the data; and Runxing Gao, Sitong Zhou and Xiaodong Chen reviewed and edited the manuscript. Data curation, Ronghua Yang and Yidan Sun; Visualization, Tianqi Chen, All authors have read and agreed to the published version of the manuscript. All authors have read and approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (Nos. 82002913 and 81772136) , Foundation of Foshan City (Nos. FS0AA-KJ218-1301-0034, 2018AB003411) and Special fund of Foshan Summit plan(Nos.2019C002, 2019D008, 2019A006 and 2020A015).

Data availability

Supplementary Figure 1 could be obtained in the Supplementary Material for correlation analysis of 22 tumor infiltrating immune cells in TCGA database.

Competing interests statement

The authors declare no potential conflflict of interest such as employment, consultancies, stock ownership, honoraria, paid expert testimony, patent applications/registrations, and grants or other funding. No author is a member of the editorial board of this journal.

RL S, KD M, Id AJ, - O, Cancer statistics. 2018. - CA Cancer J Clin 2018 Jan;68(1):7–30 doi: 103322/caac21442 Epub 2018 Jan 4. (– 1542–4863 (Electronic)):- 7–30. - eng.
Ott PA, Hu Z, Keskin DB, Shukla SA, Sun J, Bozym DJ, et al. An immunogenic personal neoantigen vaccine for patients with melanoma. Nature. 2017 Jul 13;547(7662):217 – 21. PubMed PMID: 28678778. Pubmed Central PMCID: PMC5577644. Epub 2017/07/06. eng.
Keskin DB, Anandappa AJ, Sun J, Tirosh I, Mathewson ND, Li S, et al. Neoantigen vaccine generates intratumoral T cell responses in phase Ib glioblastoma trial. Nature. 2019 Jan;565(7738):234–9. PubMed PMID: 30568305. Pubmed Central PMCID: PMC6546179. Epub 2018/12/21. eng.
AJ G, AM N, SV CLL, W B. F, D K, et al. - The prognostic landscape of genes and infiltrating immune cells across human. - Nat Med 2015 Aug;21(8):938–945 doi: 101038/nm3909 Epub 2015 Jul 20. (– 1546-170X (Electronic)):- 938 – 45. - eng..
- RSEWOBAJ. Cancer statistics, 2011: the impact of eliminating socioeconomic and racial. - CA Cancer J Clin 2011 Jul-Aug;61(4):212 – 36 doi: 103322/caac20121 Epub 2011 Jun. (– 1542–4863 (Electronic)):- 212 – 36. - eng.
JF T, RA S, RF K. - Cutaneous melanoma in the era of molecular profiling. - Lancet 2009 Aug 1;374(9687):362-5 doi: 101016/S0140-6736(09)61397-0. (– 1474-547X (Electronic)):- 362-5. - eng.
LA F, SD C, MJ K, RM A. - Toward a molecular classification of melanoma. - J Clin Oncol 2007 Apr 20;25(12):1606-20 doi: 101200/JCO2006060442. (– 1527–7755 (Electronic)):- 1606-20. - eng.
AG G. D G. - Strategies for early melanoma detection: Approaches to the patient with nevi. - J Am Acad Dermatol. 2009 May;60(5):719–35. quiz 736-8 doi:. (– 1097–6787 (Electronic)):- 719 – 35; quiz 36 – 8. - eng.
NG BMAG, CMS, AL N, JR. S, et al. - Clonal analysis of cytotoxic and regulatory T cell responses against human melanoma. - J Exp Med 1989 Jun 1;169(6):1961-76 doi: 101084/jem16961961. (– 0022-1007 (Print)):- 1961-76. - eng.
Tian QE, Li HD, Yan M, Cai HL, Tan QY, Zhang WY. Astragalus polysaccharides can regulate cytokine and P-glycoprotein expression in H22 tumor-bearing mice. World journal of gastroenterology. 2012 Dec 21;18(47):7079-86. PubMed PMID: 23323011. Pubmed Central PMCID: PMC3531697. Epub 2013/01/17. eng.
DL B, MJ RM, AC R, LP HBS. R, et al. - Long-term outcomes of a phase I study of agonist CD40 antibody and CTLA-4 blockade. - Oncoimmunology 2018 Aug 20;7(10):e1468956 doi: 101080/2162402X20181468956. (– 2162–4011 (Print)):- e1468956. - eng.
MA C, - WMHYJPA PD-1 and CTLA-4 combination blockade expands infiltrating T cells and reduces. - Proc Natl Acad Sci U S A 2010 Mar 2;107(9):4275-80 doi: 101073/pnas0915174107. (– 1091–6490 (Electronic)):- 4275-80. - eng.
RA FA, P S, M R, SW MRM. M, et al. - Tumor-infiltrating lymphocyte grade is an independent predictor of sentinel lymph. - J Clin Oncol 2012 Jul 20;30(21):2678-83 doi: 101200/JCO2011378539 Epub 2012. (– 1527–7755 (Electronic)):- 2678-83. - eng.
SR S, LE M, CC S, KA SC, DP H. D, et al. - Prognostic value of B cells in cutaneous melanoma. - Genome Med 2019 May 28;11(1):36 doi: 101186/s13073-019-0647-5. (– 1756-994X (Electronic)):- 36. - eng.
FS H, SJ OD, DF M, RW W, JB JAS. H, et al., inventors- Improved survival with ipilimumab in patients with metastatic melanoma.
SL T, FS H, SN JRB, DC G. S, DF M, et al., inventors- Safety, activity, and immune correlates of anti-PD-1 antibody in cancer.
RL G, AP H, HE VF, DR V, WA L. K, et al. - Toward a Shared Vision for Cancer Genomic Data. - N Engl J Med 2016 Sep 22;375(12):1109-12 doi: 101056/NEJMp1607591. (– 1533–4406 (Electronic)):- 1109-12. - eng.
Zhang X. The clinical observation of radix astragali polysaccharide inoculation fluid union DC-CIK cure Lacking in vital energy syndrome breast cancer patient. FUJIAN UNIVERSITY OF TRADITIONAL CHINESE MEDICINE. 2011:1–34. Epub 2011-06-1.
Darvin P, Toor SM, Sasidharan Nair V, Elkord E. Immune checkpoint inhibitors: recent progress and potential biomarkers. Experimental & molecular medicine. 2018 Dec 13;50(12):1–11. PubMed PMID: 30546008. Pubmed Central PMCID: PMC6292890. Epub 2018/12/14. eng.
Sanmamed MF, Chen L. A Paradigm Shift in Cancer Immunotherapy: From Enhancement to Normalization. Cell. 2018 Oct 4;175(2):313 – 26. PubMed PMID: 30290139. Pubmed Central PMCID: PMC6538253. Epub 2018/10/06. eng.
Sanmamed MF, Chen L. A Paradigm Shift in Cancer Immunotherapy: From Enhancement to Normalization. Cell. 2019 Jan 24;176(3):677. PubMed PMID: 30682374. Epub 2019/01/27. eng.
Kantoff PW, Higano CS, Shore ND, Berger ER, Small EJ, Penson DF, et al. Sipuleucel-T immunotherapy for castration-resistant prostate cancer. The New England journal of medicine. 2010 Jul 29;363(5):411–22. PubMed PMID: 20818862. Epub 2010/09/08. eng.
Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017 Jan 18;541(7637):321 – 30. PubMed PMID: 28102259. Epub 2017/01/20. eng.
ME R, CW BPDWYH. L, W S, et al. - limma powers differential expression analyses for RNA-sequencing and microarray. - Nucleic Acids Res 2015 Apr 20;43(7):e47 doi: 101093/nar/gkv007 Epub 2015 Jan 20. (– 1362–4962 (Electronic)):- e47. - eng.
C SCDY, Id LYL. O, X S, et al. - Identification of crucial genes in abdominal aortic aneurysm by WGCNA. - PeerJ 2019 Oct 8;7:e7873 doi: 107717/peerj7873 eCollection 2019. (– 2167–8359 (Print)):- e7873. - eng.
AM N, CL L, MR G, Id O, AJ G, Id O, et al. - Robust enumeration of cell subsets from tissue expression profiles. - Nat Methods 2015 May;12(5):453-7 doi: 101038/nmeth3337 Epub 2015 Mar 30. (– 1548–7105 (Electronic)):- 453-7. - eng.
Monti S, Tamayo P, Mesirov JP, Golub TR. Consensus Clustering: A Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data. Mach Learn. 2003;52(1):91–118.
Azarov I, Peskov K, Helmlinger G, Kosinsky Y. Role of T Cell-To-Dendritic Cell Chemoattraction in T Cell Priming Initiation in the Lymph Node: An Agent-Based Modeling Study. Frontiers in immunology. 2019;10:1289. PubMed PMID: 31244840. Pubmed Central PMCID: PMC6579912. Epub 2019/06/28. eng.
Shao P, Zhao LH, Zhi C, Pan JP. Regulation on maturation and function of dendritic cells by Astragalus mongholicus polysaccharides. Int Immunopharmacol. 2006 Jul;6(7):1161–6. PubMed PMID: 16714220. Epub 2006/05/23. eng.
Yu G, Wang LG, Han Y, He Q. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. Omics A Journal of Integrative Biology. 2012;16(5):284–7.
Ding Z, Li Q, Zhang R, Xie L, Shu Y, Gao S, et al. Personalized neoantigen pulsed dendritic cell vaccine for advanced lung cancer. Signal transduction and targeted therapy. 2021 Jan 20;6(1):26. PubMed PMID: 33473101. Pubmed Central PMCID: PMC7817684. Epub 2021/01/22. eng.
Wculek SK, Amores-Iniesta J, Conde-Garrosa R, Khouili SC, Melero I, Sancho D. Effective cancer immunotherapy by natural mouse conventional type-1 dendritic cells bearing dead tumor antigen. Journal for immunotherapy of cancer. 2019 Apr 8;7(1):100. PubMed PMID: 30961656. Pubmed Central PMCID: PMC6454603. Epub 2019/04/10. eng.
X R, N T, N AH. T, F L, JC S, et al., inventors- pROC: an open-source package for R and S + to analyze and compare ROC curves.
NB W, Id O. B W, M G, M R, K T, et al. - Tumor microenvironment-derived S100A8/A9 is a novel prognostic biomarker for. - J Immunother Cancer 2019 Dec 5;7(1):343 doi: 101186/s40425-019-0828-1. (– 2051 – 1426 (Electronic)):- 343. - eng.
TC M, O H, DC S, TM B, JS W, AJ O, et al. - Epacadostat Plus Pembrolizumab in Patients With Advanced Solid Tumors: Phase I. - J Clin Oncol 2018 Nov 10;36(32):3223–3230 doi: 101200/JCO2018789602 Epub 2018. (– 1527–7755 (Electronic)):- 3223-30. - eng.
Chai C, Ma Q, Zhou L, Niu J, Wang T, Shi Y, et al. Astragalus polysaccharide can enhance the anti-tumor effect of mouse macrophages and natural killer cells. Chinese Journal of Microbiology Immunology. 2019;39(4):292–7.
Chai W, He X, Zhu J, Lv C, Zhao H, Lv A, et al. Immunoregulatory effects of astragalus polysaccharides on myeloid derived suppressor cells in B16-F10 tumor bearing mice. Chinese Journal of Basic Medicine in Traditional Chinese Medicine. 2012;18(1):63–5. Epub 2011-06-11.
Wei W, Li Y, Li Y, Li D. Adipose-specific knockout of ubiquitin-conjugating enzyme E2L6 (Ube2l6) reduces diet-induced obesity, insulin resistance, and hepatic steatosis. Journal of pharmacological sciences. 2021 Apr;145(4):327–34. PubMed PMID: 33712284. Epub 2021/03/14. eng.
Murakami M, Izumi H, Kurita T, Koi C, Morimoto Y, Yoshino K. UBE2L6 is Involved in Cisplatin Resistance by Regulating the Transcription of ABCB6. Anti-cancer agents in medicinal chemistry. 2020;20(12):1487–96. PubMed PMID: 32329696. Epub 2020/04/25. eng.
Chen H, Zeng Y, Shao M, Zhao H, Fang Z, Gu J, et al. Calcineurin A gamma and NFATc3/SRPX2 axis contribute to human embryonic stem cell differentiation. Journal of cellular physiology. 2021 Jan 3. PubMed PMID: 33393109. Epub 2021/01/05. eng.
Gao Z, Wu J, Wu X, Zheng J, Ou Y. SRPX2 boosts pancreatic cancer chemoresistance by activating PI3K/AKT axis. Open medicine (Warsaw. Poland). 2020;15(1):1072–82. PubMed PMID: 33336063. Pubmed Central PMCID: PMC7718643. Epub 2020/12/19. eng.
Li X, Liu J, Sun H, Zou Y, Chen J, Chen Y, et al. SRPX2 promotes cell proliferation and invasion via activating FAK/SRC/ERK pathway in non-small cell lung cancer. Acta biochimica Polonica. 2020 Jun 18;67(2):165 – 72. PubMed PMID: 32550700. Epub 2020/06/19. eng.
Pidugu VK, Pidugu HB, Wu MM, Liu CJ, Lee TC. Emerging Functions of Human IFIT Proteins in Cancer. Frontiers in molecular biosciences. 2019;6:148. PubMed PMID: 31921891. Pubmed Central PMCID: PMC6930875. Epub 2020/01/11. eng..
Kabuye D, Chu Y, Lao W, Jin G, Kang H. Association between CLEC4E gene polymorphism of mincle and pulmonary tuberculosis infection in a northern Chinese population. Gene. 2019 Aug 20;710:24 – 9. PubMed PMID: 31075410. Epub 2019/05/11. eng.
Pahari S, Negi S, Aqdas M, Arnett E, Schlesinger LS, Agrewala JN. Induction of autophagy through CLEC4E in combination with TLR4: an innovative strategy to restrict the survival of Mycobacterium tuberculosis. Autophagy. 2020 Jun;16(6):1021–43. PubMed PMID: 31462144. Pubmed Central PMCID: PMC7469444. Epub 2019/08/30. eng.
Attia JVD, Dessens CE, van de Water R, Houvast RD, Kuppen PJK, Krijgsman D. The Molecular and Functional Characteristics of HLA-G and the Interaction with Its Receptors: Where to Intervene for Cancer Immunotherapy? International journal of molecular sciences. 2020 Nov 17;21(22). PubMed PMID: 33213057. Pubmed Central PMCID: PMC7698525. Epub 2020/11/21. eng.

SupplementFigure101.tif

Download PDF

Version 1

posted

You are reading this latest preprint version

Immuno-Genomic Landscape of Cutaneous Melanoma Identifies a Prognostic and Immunotherapeutically Relevant Gene Signature

Status:

Version 1

Abstract

Figures

1. Introduction

2. Methods And Materials

2.1 Acquisition and preprocessing of SKCM expression datasets

2.2 tumor-infiltrating immune cells (TIICs) analysis and clustering of samples

2.3 Differential gene expression pattern clustering analysis

2.4 RNA sequence expression analysis in GEPIA

2.5 Gene set enrichment analysis (GSEA)

2.6 Evaluation of patient biological characteristics

2.7 Establishment of Immune characteristic model and clinical prediction model

2.7 Statistical analysis

2.8 Immunohistochemistry (IHC)

2.8.1 Patients tissue specimens

2.8.2 Immunohistochemistry (IHC) staining

2.8.3 Histologic scoring and analysis

3. Results

3.1 Immune infiltration analysis related to SKCM patients

3.2 Construction of TME-related signature gene in SKCM patients

3.3 CeRNA regulatory network construction and signature gene expression in SKAM patients

3.4 Construction of immune-related prognostic gene signature

3.5 Gene Set Enrichment Analysis (GSEA)

3.6 Correlation analysis of the risk score and clinicopathological characteristics

3.7 Application value of the risk score in pan-cancer

3.8 The prognostic value of hub genes in SKCM patients

2.9 Ube2L6, SRPX2, IFIT2 were higher expression in IHC melanoma specimens

Discussion

Conclusion

Declarations

References

Supplementary Files

Status:

Version 1