4.1 Mutation and copy number analysis of genes related to m6A modification
Genetic changes were analyzed by mutation and copy number. In terms of mutation, 39 (8.58%) of 452 samples in TCGA-PCA data set had RNA modification mutation, among which ZC3H13 had the highest mutation rate, reaching 21% (Fig.1A). In order to explore whether the genetic change of m6A affects the progression of prostate cancer, all samples were divided into two groups according to whether the mutation occurred. The survival analysis of these two groups showed that the survival difference between the two groups was significant (P = 0.021), and the prognosis of the mutation group was poor, which indicated that the mutation of m6A gene could affect prostate cancer (Fig.1B). GSVA showed that the mutation group of m6A gene was more concentrated in the cancer hallmark gene set, such as RTK and EASMAPK signaling pathways. This mutation of m6A gene may lead to functional changes, which may affect the prognosis of prostate cancer (Fig.1C).
Copy number amplification generally corresponds to increased expression, while copy number deletion generally corresponds to decreased expression. To further explore whether somatic copy number changes affect the prognosis and progression of prostate cancer, our study showed that the change of copy number of each kind of m6A(Fig.2A), and the expression of each kind of m6A in tumor and paired normal(Fig.2B). It was found that the high frequency of CNV gain of RBM15B, HNRNPC and YTHDF3 in PCA corresponded to the increased expression, which suggested that CNV might be a regulator of the expression of m6A mRNA. However, the CNV loss of some m6A modified genes corresponded to increased expression. The samples were divided into three groups (CNV gain, CNV loss, none CNV) according to the CNV of m6A modified gene. The expression of FTO,YTHDF3, YTHDC2, ZC3H13 and IGF2BP2 in CNV loss group was significantly lower than that in normal group (Fig.2C).
4.2. Different RNA modification patterns related to hallmark in cancer
The original count and clinical information of RNA sequencing data of m6A modification related genes were obtained from the Cancer Genome Atlas (TCGA) data set. In the present study, a total number of 452 PCA patients were included from TCGA. The clinicopathologic characteristics of the patients are listed in Table 1. In order to explore different RNA modification patterns, correlation analysis was conducted on the expression data of 20 m6A modified genes, and it was found that there was interaction between different types of m6A modified genes (Fig.3A). Next, all samples were clustered according to the expression of these 20 m6A modified genes, and two clusters were obtained (Cluster_1 and cluster_2)(Fig 3B, 3C). Survival analysis showed that there was a significant difference in the prognosis between the two groups (Fig.3D). So what do these two types affect PCA? Two clustering gene expression heat maps and volcano maps (Fig 3E, 3F) were drawn. Go and KEGG methods were used to analyze the two types of differential genes. The GO of up-regulated gene was mainly enriched in stone modification and covalent chromotin modification, and the KEGG of up-regulated gene was mainly enriched in proteoglycans in cancer and focal adhesion. The GO of down_regulated genes was mainly enriched in oxidative phosphorylation, while KEGG of down-regulated genes was mainly enriched in thermogenesis and oxidative phosphorylation (Fig.3G).
4.3. Construct the signature of m6A RNA modification
Firstly, PCA transcription samples were divided into two groups by using 20 m6A modifying genes and clinical information(Cluster_1 and cluster_2). Using these two groups of samples, 3022 differentially expressed genes (DEGs) were screened, and then the samples were divided into two groups by unsupervised clustering (gene.cluster_A and gene.Cluster_B) , and the heat maps and volcano maps(Fig.4A, 4B) were drawn. Finally, it was found that the samples in gene.cluster_A and Cluster_1 were consistent, and the samples in gene.cluster_B and Cluster_2 were consistent, and the prognosis of gene.cluster_A group was worse than that of gene.cluster_B group (Fig.4C), which also verified that the 3022 DEGs could represent the m6A RNA modification pattern. Univariate analysis of prognosis (Fig.4D) was performed for 3022 differentially expressed genes, and a scoring model (Risk_ Score) was constructed using genes with prognostic value, which was used to quantify RNA modification patterns in prostate cancer patients. The risk score of cluster_1 was significantly higher than that of cluster_2 (Fig.4E), and the risk score of high-grade (T and N) tumors was higher than that of low-grade (Fig.4F, 4G). The PFS survival of risk_score-low group was higher than that in risk_score-high group (Fig.4H). In order to test whether the risk_score can be used as an independent prognostic indicator, multivariate Cox risk regression was conducted and it was found that the risk_score group is a poor and independent prognostic indicator (Fig.4I). The cancer related pathways of GSVA showed that risk_score-high group was negatively correlated with Cellcycle, EMT and hormone Er, and positively correlated with RTK, EMT and RasMAPK. On the contrary, risk_score-low was negatively correlated with RTK, RasMAPK and EMT, and positively correlated with hormone ER and dnadamage (Fig.4J). Immunocorrelation of GSVA analysis showed that high risk scores were positively correlated with itreg cells, but positively correlated with B cell, monocyte and CD8+T were negatively correlated; On the contrary, low risk scores were associated with B cell, macrophage, monocyte and CD8+, and negatively correlated with Treg cells (Fig.4K).
4.4. Prediction of immune microenvironment and immunotherapy response by risk_score group
A large number of articles have reported that immune checkpoint genes and infiltrating immune cells are related to the modification of m6A RNA. Heat map of immune checkpoint gene expression showed that immune checkpoint gene was highly expressed in high risk_score group (Fig.5A). Among them, SSIGLEC15,TIGIT, CD274, HAVCR2, LAG3 and PDCD1LG2 were highly expressed in high risk_score group (Fig. 5B). CIBERSORT method was used to get the proportion of various kinds of cells in each sample, and we found that m6A RNA modification may be strongly related to the invasion of TME cells. According to the high and low score classification, differential expression analysis of the obtained cell data showed that the infiltration of B cells, T cells CD4+, endothelial cells and T cells CD8+ was higher in risk_score-high group (Fig. 5C), which may be the reason for its better prognosis. The result of patients' response to ICB treatment predicted by risk_score showed that risk score was positively correlated with TIDE score (Fig.5D), and the TIDE score of high risk group was higher than that of low risk group (Fig.5E). In addition, MSI was closely related to tumor occurrence and tumor immunotherapy, and the MSI score of high risk group was significantly higher than that of low risk group (Fig.5F). Combined with the lower TME of B cells, T cells CD4+ and T cells CD8+ in the high risk group, we can conclude that the risk_score model can well reflect the tumor immune infiltration pattern and guide the immunotherapy. In order to evaluate the role of risk score in drug response, the relationship between risk score and drug response of cancer cell lines was explored. By using Spearman correlation coefficient, the combination of drug sensitivity between risk score and CTRP database was identified, and 30 kinds of drugs related to risk score were screened out. Among them, the high score group was mainly related to drug resistance, and the low score group was mainly related to drug sensitivity (Fig.5G).