17 lncRNAs were identified as autophagy-related
A total of 14142 lncRNAs were distinguished from the transcriptome matrix downloaded from the TCGA database. 258 autophagy-related mRNAs were obtained from the HADb database. The co-expression network was constructed to identify autophagy-related lncRNAs. Finally, 17 lncRNAs and 99 autophagy-related mRNA were screened out and considered as significant correlation (table S1, |R|>0.3, p value < 0.0001). We then used these RNAs as nodes to draw a network (Fig. 2a), the red ellipse indicating lncRNA, blue ellipse indicating autophagy-related mRNA, and the line between them indicating the co-expression relationship. Furthermore, we created this image (Fig. 2b)—a Sankey diagram—by plotting the connection within each node and connecting them to risk type which based on the HR value of multivariate Cox proportional hazards regression analyses (HR > 1 was considered as risk gene, HR < 1 was considered as protecting gene. table s2). Three of them were defined as protect genes (AC121338.2, EPB41L4A-DT, LINC01843) and others were all risk genes (AL391244.3, AC011462.4, AC103706.1, SNHG15, AL590094.1, AP003352.1, AC026356.2, SNHG17, LINC00460, HOTAIRM1, AC084876.1, AC027796.4, MELTF-AS1, AC010973.2).
Survival Analyses Of These 17 Autophagy-related Lncrnas
Given the 17 autophagy-related lncRNAs, we next explored their impact of expression level on OS (Fig. 3). For three protective autophagy-related lncRNAs, they showed a significant tendency that higher expression predicting better survival rate. In contrast, for 14 risky autophagy-related lncRNAs, higher expression may lead to poorer survival rates. The survival analysis results were consistent with the risk type classified by multivariate Cox proportional hazards regression analyses.
Risk Signature Assessment
After performing survival analyses of these autophagy-related lncRNAs, we next constructed the prognosis assessment of risk signature. First, we divided patients into high-risk group and low-risk group to compare their survival difference. As shown in Fig. 4a, the high-risk group displayed a lower survival rate than the low-risk group. Next, we plotted the expression patterns of the 17 autophagy-related lncRNAs, survival status dot plot, and risk score curve to visualize the association between these variables (Fig. 4b). The heatmap showed that 3 protective autophagy-related lncRNAs (LINC01843, EPB41L4A-DT, and AC121338.2) were down-regulated in high-risk group whereas other lncRNAs were overexpressed. Additionally, the survival status dot plot revealed that dead cases were almost concentrated in the high-risk group and more dead cases appeared with the increase of risk score.
In the ROC analysis (Fig. 4c), the prognostic model was evaluated using the area under the curve (AUC). The AUC for the risk score signature was 0.748, higher than age (0.580), gender (0.480), grade (0.743) and N stage (0.527, lymph node metastasis), but lower than stage (0.846), T stage (0.788, the extent of the primary tumor) as well as M stage (0.751, distant metastasis). The results suggested that the prognostic value of risk signature was better than age, gender, grade, and N stage.
For the independent prognosis analysis, age, gender, grade, stage, TNM system, and risk score were included in the independent prognostic factors. Univariate (Fig. 4d) and multivariate (Fig. 4e) analyses suggested that the risk signature was an independent predictive factor for KIRC patients (p < 0.001, HRs were 1.112 and 1.088, respectively). The details of analysis were shown in Table 1 and Table 2.
Table 1
Univariate Cox analysis of characteristics and risk score in KIRC.
id | B | SE | z | HR | HR.95L | HR.95H | pvalue |
age | 0.023401 | 0.009062 | 2.582371 | 1.023677 | 1.005656 | 1.042021 | 0.009812** |
gender | 0.017341 | 0.213742 | 0.081132 | 1.017492 | 0.669258 | 1.546923 | 0.935337 |
grade | 0.757602 | 0.143102 | 5.294148 | 2.133156 | 1.61144 | 2.823781 | 1.20E-07*** |
stage | 0.59752 | 0.097529 | 6.126595 | 1.817605 | 1.501353 | 2.200474 | 8.98E-10*** |
T | 0.614276 | 0.119009 | 5.161601 | 1.848318 | 1.463781 | 2.333873 | 2.45E-07*** |
M | 1.391904 | 0.222367 | 6.259497 | 4.022501 | 2.601461 | 6.219781 | 3.86E-10*** |
N | 1.110114 | 0.336707 | 3.296975 | 3.034706 | 1.568596 | 5.871135 | 0.000977*** |
riskScore | 0.106298 | 0.015117 | 7.031483 | 1.112153 | 1.079684 | 1.145598 | 2.04E-12*** |
**p < 0.01; ***p < 0.001. |
Table 2
Multivariate Cox analysis of characteristics and risk score in KIRC.
id | B | SE | z | HR | HR.95L | HR.95H | pvalue |
age | 0.037606 | 0.010293 | 3.653502 | 1.038322 | 1.017585 | 1.059482 | 0.000259*** |
gender | 0.262189 | 0.228932 | 1.145275 | 1.299773 | 0.829852 | 2.035796 | 0.252095 |
grade | 0.269428 | 0.172746 | 1.559684 | 1.309216 | 0.93319 | 1.836761 | 0.118835 |
stage | 0.416601 | 0.264258 | 1.57649 | 1.516797 | 0.903629 | 2.546036 | 0.114913 |
T | -0.14512 | 0.250034 | -0.5804 | 0.86492 | 0.529843 | 1.411903 | 0.561648 |
M | 0.509551 | 0.416789 | 1.222565 | 1.664544 | 0.735401 | 3.767613 | 0.221494 |
N | 0.428635 | 0.362005 | 1.184059 | 1.535161 | 0.755118 | 3.120992 | 0.23639 |
riskScore | 0.084718 | 0.018718 | 4.526048 | 1.088411 | 1.049204 | 1.129082 | 6.01E-06*** |
***p < 0.001. |
Clinical Value Of The Risk Signature For Kirc Patients
Subsequently, we evaluated the clinical value of the risk signature regarding age, gender, grade, stage, and TNM. As showed in Table 3, the risk score tended to elevate in higher grade (grade 3–4, p = 0.0001), advanced stage (stage Ⅲ-Ⅳ, p = 0.0001), higher T stage (p = 0.001) and N stage (p = 0.04). These results suggested that the risk signature might be associated with the progression of KIRC.
Table 3
Correlations between risk score signature and clinical features in the TCGA cohort.
| | | Risk score |
Clinical | Group | n | Mean | SD | t | P |
age | <=65 | 155 | 1.977 | 4.47 | 0.614461 | 0.54 |
| > 65 | 91 | 1.719 | 2.088 | | |
gender | Female | 97 | 1.803 | 2.73 | -0.28856 | 0.773 |
| Male | 149 | 1.933 | 4.316 | | |
grade | G1-2 | 109 | 0.973 | 0.87 | -3.84431 | 0.0001*** |
| G3-4 | 137 | 2.605 | 4.873 | | |
stage | Stage Ⅰ-Ⅱ | 133 | 1.036 | 0.978 | -3.64565 | 0.0001*** |
| stage Ⅲ-Ⅳ | 113 | 2.878 | 5.294 | | |
T | T1-2 | 145 | 1.098 | 1.045 | -3.40555 | 0.001** |
| T3-4 | 101 | 3.007 | 5.566 | | |
M | M0 | 205 | 1.463 | 1.899 | -2.0149 | 0.051 |
| M1 | 41 | 3.976 | 7.94 | | |
N | N0 | 232 | 1.659 | 3.442 | -2.27093 | 0.04* |
| N1-3 | 14 | 5.581 | 6.406 | | |
*p < 0.05; **p < 0.01; ***p < 0.001. |
Identified autophagy-related lncRNAs may function in tumor progression via Immunomodulation
To explore the potential mechanism of these 17 autophagy-related lncRNAs’ effect on tumor progression, we conducted Gene Ontology (GO) functional annotation with GSEA. Interestingly, immune-related pathways were markedly enriched in the high-risk group, including immune response mediated by circulating immunoglobulin, interferon-gamma production, regulatory T-cell differentiation, positive regulation of interferon-gamma production, immunoglobulin complex, B cell mediated immunity, regulation of humoral Immune response, and immunoglobulin receptor binding (Fig. 5a). The detailed enrichment results were summarized in Table 4. Given the pathway enriched in immune regulation, we next conducted analyses of the correlation between the expression of these lncRNAs and immune cells infiltration level in KIRC (Fig. 5b). The results revealed that the infiltration of T cell regulatory (Tregs), T cell NK and T cell CD4 + Th1 had a significantly positive correlation with lncRNAs expression. In contrast, the infiltration of neutrophil exhibited a negative correlation with the expression of lncRNAs.
Table 4
Gene set enrichment analysis (GSEA) results based on the risk signature of 14 autophagy related lncRNAs.
NAME | SIZE | ES | NES | P | FDR q-val | FWER p-val | RANK AT MAX | LEADING EDGE |
GO_HUMORAL_IMMUNE_RESPONSE_ MEDIATED_BY_CIRCULATING_ IMMUNOGLOBULIN | 139 | 0.600492 | 2.296511 | 0 | 0 | 0 | 18013 | tags = 77%, list = 33%, signal = 114% |
GO_INTERFERON_GAMMA_ PRODUCTION | 108 | 0.596715 | 2.199109 | 0 | 0 | 0 | 8083 | tags = 47%, list = 15%, signal = 55% |
GO_REGULATORY_T_CELL_ DIFFERENTIATION | 31 | 0.761541 | 2.293867 | 0 | 0 | 0 | 5855 | tags = 71%, list = 11%, signal = 79% |
GO_POSITIVE_REGULATION_OF_ INTERFERON_GAMMA_ PRODUCTION | 64 | 0.655313 | 2.235942 | 0 | 0 | 0 | 7431 | tags = 55%, list = 13%, signal = 63% |
GO_IMMUNOGLOBULIN_COMPLEX | 65 | 0.602932 | 2.09714 | 0 | 3.13E-04 | 0.006 | 14405 | tags = 65%, list = 26%, signal = 87% |
GO_B_CELL_MEDIATED_IMMUNITY | 208 | 0.511539 | 2.055525 | 0 | 5.55E-04 | 0.015 | 14596 | tags = 59%, list = 26%, signal = 80% |
GO_REGULATION_OF_HUMORAL_ IMMUNE_RESPONSE | 130 | 0.543822 | 2.056492 | 0 | 6.06E-04 | 0.015 | 18640 | tags = 77%, list = 34%, signal = 116% |
GO_IMMUNOGLOBULIN_RECEPTOR_ BINDING | 65 | 0.595271 | 2.021653 | 0 | 0.00122 | 0.04 | 17569 | tags = 71%, list = 32%, signal = 104% |
ES, enrichment score; NES, normalized enrichment score. FDR, false discovery rate; FWER, familywise-error rate; |