Single cell sequencing analysis constructed the N7-methylguanosine (m7G)-related prognostic signature in uveal melanoma

TCGA data downloading and processing

UVM RNA-seq data were downloaded from the TCGA (https://portal.gdc.cancer.gov/) database. Data inclusion criteria were: (1) patients with prior pathological diagnosis of Uveal melanoma and (2) patients with detailed clinical information, including survival time, clinical stages, survival status, and so on. All data were processed and integrated using Perl software.

WCGNA analysis for phenotype-associated modules and genes identification

Weighted gene co-expression network analysis (WCGNA) is an integrative method used for identifying modules and genes that are most significantly related to the sample phenotypes. A WGCNA analysis was performed using the WGCNA R package (version 3.2.2) to assess the relationship between clusters of co-expressed genes and the tumor-related phenotypes.

Establishment of an m⁷G-related gene signature for prognosis

The method for establishing a prognostic signature has been intensively described in our previous articles. The least absolute shrinkage and selection operator (LASSO) regression was performed to construct a prognostic signature. The risk score can be illustrated as follows:

$$\begin{array}{c}\text{Risk}\text{score}={\sum }_{i=1}^n{\beta }_i\ast gene\text{\_}expression\\ ({\beta }_i:coefficientofgenei;\\ gene\text{\_}expressio{n}_i:expressionofgenei)\end{array}$$

The patients were grouped by risk scores, which divided them into high-risk and low-risk groups. Survival analysis was performed accordingly using package Survival (version 3.2).

Single-cell sequencing and data processing

Single-cell sequencing data of Uveal melanoma GSE139829 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo). The data were filtered using R software and the Seurat package (version 4.1.1) to eliminate cells that do not meet the data quality standards. The number of highly variable genes was set at 3000. These 11 samples were integrated using "Harmony" (version 0.1.0). Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) method was used to reduce the dimension of data. Clustering of cells was conducted using the K-Nearest Neighbor (KNN) method with resolution=1.0. Cells were annotated using cell markers obtained from previous literature. Finally, the percentage of m⁷G-related genes in each cell was calculated by importing necroptosis genes using the “PercentageFeatureSet” function.

Tumor immune infiltration analysis

In this study, we employed multiple methods for tumor immune infiltration. Including TIMER, QUANTISEQ, and so on [23, 24]. Immune checkpoint analysis was also performed to examine the immunological differences between the high-risk and the low-risk groups.

Nomogram construction and ROC analysis

The nomogram of the risk score and relevant clinical information was depicted using the package Survival (version 3.2) and package rms (version 6.3).

Cell culture, transfection, and qRT-PCR assay

MUM-2B, C918 cells were cultured in the RPMI1640 medium and OCM-1 cells in DMEM medium (Gibco, Shanghai, China) supplemented with 10% FBS (AuxGene, Gold Coast, Australia), transfected with siRNA (Ribobio, Nanjing, China) according to the manufacturer’s protocol. The sequences of the siRNAs we used are listed in Supplementary Table 1. Total cellular RNAs were extracted using Trizol™ Reagent (Life Technologies, CA, USA), and the reverse transcription and qRT-PCR assay were conducted using a reverse transcription kit and SYBR (Vazyme, Nanjing, China) according to the manufacturer’s protocol.

Cell proliferation assay

5-ethynyl-2’-deoxyuridine (EdU) incorporation assay was performed to examine the cell proliferation capacity of the UVM cell lines. After incubation with EdU for 2 hours, cells were fixed with 4% Paraformaldehyde and stained with BeyoClick™ EdU Cell Proliferation Kit (Beyotime, Nanjing, China) according to the manufacturer’s protocol. A clustering Formation assay was also used to examine the proliferation ability of the cell lines. After transfection, the cells were seeded into a 6-well culture plate (Corning, NY, USA) and placed in a cell culture incubator for 10 days. The cells were then rinsed with PBS, fixed with methanol, stained with crystal violet, and eventually photographed using a digital camera (Canon).

Scratch-wound healing assay

The scratch-wound assay is a simple, reproducible method commonly used to measure basic cell migration parameters. Cells were allowed to grow to confluence before a thin “wound” was introduced by scratching with a 10μL pipette tip. The cells subsequently began to migrate from the wound edge to fill the wound space. The wound size 24 hours after the introduction of the wound was recorded and quantified to assess the migration capacity of the cell lines.

Transwell assay

Successfully transfected UVM cell lines were seeded into the upper well for 24 h and allowed to invade through the transwell plate (Corning, NY, USA). The cells on the inserts were fixed with methanol, stained with crystal violet, and counted under a light microscope to assess the migration capacity of the cells.

Statistical analysis

In vitro, experimental data were processed using the SPSS software (version 26.0). All data were illustrated employing mean ± SD. One-way analysis of variance (ANOVA) and Turkey’s multiple comparisons of the means were employed in the analysis of multiple sets of data unless otherwise specified.

Availability of data and materials

TCGA and GEO databases.

Single-cell analysis of the data

11 samples of UVM were included in our analysis. The batch effects between the 11 samples were negligible, and thus suitable for subsequent analysis (Figure 2A). After quality control and clustering analysis of the data, all cells were divided into 67 subclusters (Figure 2B). Based on the expression of the signature genes, all cells were annotated into five cell types: The B cells, Endothelial cells, Monocytes, macrophage cells, photoreceptor cells, Plasma cells, and T cells (Figure 2C, 2D), as well as the Tumor cells. As shown in Figure 2E, the activation of the m⁷G phenotype varied in different cell types. Using differential analysis, the genes with adjusted p -value<0.05 was set, and, | avg_log₂FC |> 0.3 were selected and 1153 genes closely associated with the m7G phenotype were defined.

Figure 2. Single cell sequencing analysis. (A) No significant batch effects were found in the 11 samples in the single-cell dataset. (B) All cells were divided into 67 subclusters. (C, D) Cell annotation. (E) M7G scoring annotation. The cells were divided into high-M⁷G group and low-M⁷G group.

WCGNA analysis filtered out modules and genes that are related to the m⁷G phenotype

The enrichment score (M7G score) of the M7G phenotype for each cell was first calculated by ssGSEA analysis. To further search for genes associated with the M7G phenotype, WGCNA analysis was performed. In Figure 3A, when setting a soft-zone threshold value reaches 9, R² > 0.8, suggesting that the data fit a power-law distribution and were suitable for WGCNA analysis. Moreover, the mean connectivity fluctuates little as the soft domain value increases. As shown in Figure 3B, all the genes were collectively clustered into six non-gray modules, and it was found that among these modules (Figure 3C), the turquoise module had the highest correlation with the M7G score, with R = 0.56 and p <0. 001. In Figure 3D, the correlation coefficient between module membership correlation and gene significance for body weight was 0.67, with p <0.001. Therefore, the 5,757 genes from the turquoise module were selected and included in the subsequent analysis.

Figure 3. Weighted co-expression network analysis. (A) When setting a soft-zone threshold value reaches 9, R² > 0.8, suggesting that the data fit a power-law distribution and were suitable for WGCNA analysis. Moreover, the mean connectivity fluctuates little as the soft domain value increases. (B, C) Clustering of modules. The turquoise module had the highest correlation with the M7G score. (D) The correlation coefficient between module membership correlation and gene significance for body weight was 0.67, with p <0.001.

The construction of a prognostic signature

After intersecting the 1153 genes obtained from the single-cell data analysis with those obtained from the WGCNA analysis described above, a total of 355 genes were obtained. Eventually, 315 genes were identified to be suitable for further analysis. We then conducted a univariate cox analysis to filter the prognosis-related genes. We found that 45 genes were of significant prognostic value (p<0.05) (Figure 4A).

Figure 4. Construction and evaluation of the prognostic signature. (A) Univariate Cox identified genes associated with prognosis. (B–D) Lasso regression was used to construct the prognostic model. (E–G) Prognosis assessment, ROC curve construction and principal component analysis in training cohort. (H–J) Prognosis assessment, ROC curve construction and principal component analysis in validation cohort.

Subsequently, in the TCGA cohort, we further screened for genes associated with prognosis using the LASSO regression method. (Figure 4B, 4C) When the optimum Lambda is 0.055, A total of eight genes significantly associated with UVM prognosis were obtained, genes CTSF, PAG1, TFAP2C, SAMD4A, DNAJA4, SRRM2, NDUFA13, and CHCHD10 were included in the prognostic signature. The signature was illustrated as follows:

$$\begin{array}{l}RiskScore=0.282\ast TP{M}_{DNAJA4}+0.196\ast TP{M}_{PAG1}+\\ 0.159\ast TP{M}_{CHCHD10}+0.048\ast TP{M}_{TFAP2C}+0.024\ast \\ TP{M}_{NDUFA13}+(-0.223)\ast TP{M}_{SRRM2}+(-1.058)\ast TP{M}_{CTSF}\end{array}$$

Figure 4D vividly depicted the genes and their corresponding coefficients. HR values of model genes were summarized in Supplementary Table 2.

All patients in the training cohort were divided into high risk-score and low risk-score groups according to the median value of the risk score. The high-risk group had significantly poorer outcomes (Figure 4E), and the model predicted that the area under the ROC curve (AUC) for patients’ 1,2,3,5-year mortality was 0.845, 0.819, 0.835, and 0.894, respectively (Figure 4F). Additionally, the model can precisely predict the outcomes of the patients with UVM in the validation cohort (downloaded from the GEO database GSE84976) (Figure 4H, 4I). Patients in the high-risk group had a poorer prognosis, and the model predicted the area under the ROC curve of 3,4,5 years was 0.823,0.742,0.882. PCA analysis further revealed the distribution of the patients in the PCA plot between the two groups varied both in the training cohort and the validation cohort. (Figure 4G, 4J).

Immune infiltration and immune checkpoint analysis showed the risk score had strong immune correlations

We further analyzed the tumor immune-microenvironment between the high-risk group and the low-risk group. Using multiple immune infiltration algorithms, we illustrated a case-based heatmap for the infiltration algorithms which scored significantly different between the two groups (Figure 5A). Furthermore, we compared the expression of 79 immune checkpoint genes acknowledged by previous researchers and performed Wilcoxon tests. (Figure 5B) [25]. As shown in Figure 5C, 5D, the expression of tumor necrosis factor and leukocyte antigen genes were all increased in the high-risk group when compared to the low-risk group.

Figure 5. Analysis of immune microenvironment. (A) Immune landscapes between high-risk and low-risk groups. (B) Expression levels of immune checkpoint related genes between high-risk and low-risk groups. (C) Tumor necrosis factor-related gene expression levels between high - and low-risk groups. (D) Expression levels of leukocyte antigen related genes between high - and low-risk groups.

Tumor mutation analysis

In Figure 6A, in the high-risk group, the top five genes mutated were, respectively, GNA11, GNAQ, BAP1, SF3B1, and C3. In the low-risk group, however, the top five genes were, GNA1, SF3B1, GNA11, EIF1AX, and BAP1 (Figure 6B). The commonly mutated genes are GNA11 and BAP1. Interestingly, gene GNAQ mutation, which is a classical alteration in the G_q pathways, is more commonly seen in the high-risk group, implying our signature’s underlying correlation with the G_q pathways.

Figure 6. Mutation analysis. (A) The mutant landscape of the high-risk group. (B) The mutant landscape of the low-risk group.

Construction of a nomogram

As shown in Figure 7A, the 1,3,5-year fatality of patient TCGA-VD-AA8N was 0.0439, 0.47, 0.939, respectively. Moreover, the ROC C-index was 0.849 95%CI:(0.771, 0.928), suggesting the high accuracy of this nomogram in assessing patient prognosis. Besides, the continuous prognosis ROC analysis found that the area under the ROC curve of the Nomogram evaluation for patient prognosis was above 0.8, which was higher than the other clinical characteristics (Figure 7B). Moreover, as shown in Figure 7C, the clinical interventions of the patients according to the Nomogram benefit the patients better than the other clinical characteristics, such as gender, age, tumor size, and Stages.

Figure 7. Construction and evaluation of the nomogram. (A) The nomogram. (B) The ROC curve. (C) The decision curve.

In vivo studies validated the oncogenic potential of gene PAG1

To further validate the robustness of our previously constructed risk model, we used in vivo experiments to validate the oncogenic role of a gene in the signature model in UVM patients, gene PAG1. Because the overexpression of PAG1 could most significantly elevate the risk score of our previously constructed model, genetic manipulation on the mRNA level was performed. Firstly, we found that the expression of the PAG1 gene was significantly suppressed after siRNA transfection in 3 variant cell lines (Figure 8A). Cell proliferation assays showed that after the knockdown of gene PAG1, MUM-2B, and OCM-1 cell lines showed a significantly lower rate of proliferation (Figure 8B, 8C). Scratch and wound healing assays showed that the migration capacity of MuM2B cells and OCM-1 cells were significantly compromised after the knockdown of gene PAG1 (Figure 8D). Furthermore, the transwell assay showed that the number of cells which migrated through the transwell plate significantly reduced after PAG1 knockdown (Figure 8E). Collectively, these results showed that the PAG1 gene plays a crucial role in the proliferation and migration of the UVM cell lines.

Figure 8. The function of PAG1 in UVM cell line was verified by cell assay. (A) Knockdown efficiency of PAG1 gene. The knockdown efficiency of si-PAG1-2 group was the highest. (B) Clone formation experiment. The proliferation ability of MUM-2B and OCM-1 cell lines decreased significantly after PAG-1 knockdown. (C) EdU assay. (D) Transwell experiment. The migration and invasion ability of MuM-2B cell lines and OCM-1 cell lines were significantly reduced after PAG-1 knockdown. (E) Wound healing experiment.