Development of a Prognostic m7G Methylation-Related miRNA Signature
Overall, 792 m7G methylation-related miRNAs were isolated from the TargetScan database,47 of which were differently expressed miRNAs(28 upregulated and 19 downregulated). The volcano plot of miRNA differential expression and the heatmap of the top 20 miRNA differential expression were displayed by Fig. 2a-2b. To identify the survival-related miRNAs, miRNA differential expression was assessed via univariable COX analyses. Eventually, 7 miRNAs(hsa-miR-146b-3p, hsa-miR-181a-2-3p, hsa-miR-6860, hsa-miR-221-5p, hsa-miR-96-5p, hsa-miR-214-3p, hsa-miR-612) were highly correlated with OS (P < 0.05) (Fig. 3a).Then, these 7 miRNAs were involved in multivariable COX analyses,and the results showed that hsa−miR−214−3p and hsa−miR−612 were significant (P < 0.05) (Fig. 3b). These 7 miRNAs were utilized to construct a prognosis model through LASSO Cox regressive analyses and the model was constructed on the basis of the optimum score of λ (Fig. 3c-3d). Thus, the equation of our signature was: Risk Score = (-0.00001043 × expression hsa-miR-146b-3p) + (-0.000038×expressionhsa-miR-181a-2-3p) + (-0.07251 × expressionhsa-miR-6860)+(-0.004331 × expression hsa-miR-221-5p)+(0.004627 × expression hsa-miR-96-5p)+(0.0251 × expression hsa-miR-214-3p)+(0.3026 × expressionhsa-miR-612).The risk scoring of all sufferers was computed, and on the basis of the mid-value of risk scoring,we classfied sufferers into the riskhigh and risklow groups(Fig. 3e).The results of PCA and t-SNE analyses revealed that the cases in diverse groups had 2 diverse distribution(Fig. 4a-4b).Riskhigh sufferers exhibited a greater mortality vs risklow sufferers (Fig. 4c).Consistently, the Kaplan-Meier (K-M) curve analyses revealed that riskhigh sufferers displayed a remarkably inferior OS vs risklow sufferers(Fig. 4d, P < 0.001). The prediction ability of the OS risk scoring of THCA patients was assessed via time-reliant ROC curves,and the the region below the curve (AUC) was 0.743, 0.802, and 0.887 of 1-,3-, 5- years, separately(Fig. 4e). Meanwhile,the influence of the expression of 7 miRNAs on the prognosis of patients was analyzed by the Kaplan-Meier curve. The outcomes unveiled that the expressing levels of hsa−miR−214−3p and hsa−miR−612 considerably affected the OS(Fig. 5a-5b),which coincided with the outcomes of multivariable COX analyses.
Independent Prognostic Value of the Risk Signature
For the purpose of determining if the risk scoring can predict the OS independently, univariable and multivariable Cox regressive analyses were utilized.In univariable Cox regressive analyses, the risk scoring(HR= 23.5262, 95% CI =3.0887-179.1961, P=0.0023),age(HR= 1.1528, 95% CI =1.0951-1.2135, P＜0.0001), stage (HR= 7.2139, 95% CI =2.3137-22.4922, P=0.0007) were tightly correlated with OS ( Fig. 6a). Likewise, in multivariable Cox regressive analyses,the risk scoring (HR= 12.9671, 95% CI =1.6444-102.2522, P=0.015) and age(HR= 1.1387, 95% CI =1.0796-1.2010, P＜0.001) were tightly correlated with OS( Fig. 6b). The prognosis nomograph integrating independent factors of OS was displayed by Fig. 6c.The C-indexes for OS forecast were 0.959 (95%CI, 0.939‐0.979).The correction plot of the predictive ability at 1-, 3-, 5-years showed a good consistency bewteen the observed overall survival and nomogram−predicted overall survival( Fig. 6d).
Functional Analyses based on Risk Score and Immune Score
To unveil the biofunctions and paths related to the risk scoring and immunity scoring,our research group initially find out the DEGs between the riskhigh and risklow groups, and the DEGs between high immune score and low immune score were also analyzed. The DEGs existing in both the risk scoring and immunity scoring were selected to complete the GO and KEGG analyses.A total of 555 DEGs in the risk score and 1028 DEGs in the immune score were acquired, among which 270 DEGs were significant in both groups( Fig. 7a). Figure 7b displays the the top 10 BP terms, CC terms, and MF terms, the primary sponged Go terms were muscle filament sliding, actin-myosin filament sliding, muscle system process, contractile fiber, and hormone metabolic process. Unfortunately, no pathways were observed based on the KEGG analysis.
Comparison of Immune Activity
The enrichment scoring of immunocytes and immunity-associated pathways was caculated by ssGSEA.The heatmap showed the enrichment scoring of each immunocyte and immune pathway( Fig. 8a).Then, the relationship between immunocytes and immunity-associated pathways was studied.Among immunocytes, TIL and Th1 cells showed the strongest positive correlation, and TIL and Tfh also showed a strong positive correlation( Fig. 8b).Among immunity-associated pathways, CCR and Check-point showed the strongest positive correlation, and T_cell_co−inhibition also showed a strong positive correlation with check-point( Fig. 8c).Eventually, our team contrasted the immune status between the riskhigh and risklow groups.B cells, DCs, iDCs, Neutrophilic cells, Tfh, Th2_cells, and Treg were remarkably diverse between these 2 groups(P＜0.05, Fig. 9a). APC_co_stimulation and Type_II_IFN_reponse also demonstrated a notable diversity(P＜0.05, Fig. 9b).
Relationship among Prognostic DEGs and Immunocytes and Immunity-associated pathways
A total of 10 DEGs (FABP4, DPP6, PKHD1L1, TNNT3, EPYC, MYH2, MYL1, SLC25A15, TNNI2, ACTA1) were screened out throuth univariate cox analyses. EPYC, TNNT3, and TNNI2 displayed a strong positive relationship with immunocytes and immunity-associated pathways.Meanwhile, FABP4, DPP6, PKHD1L1 and SLC25A15 showed a negative correlation with immunocytes and immunity-associated pathways( Fig. 10).