Identification of ARGs in HER2 + breast cancer cohort and healthy samples
To identify the differential expression of ARGs in HER2 + patients, we analysis 231 autophagy related genes in 53 samples from GSE29431, which included 41 HER2 + breast tissue samples and 12 normal breast tissue samples. The samples were categorized into the “Normal” group (n = 12) and “Her2” group (n = 41) based on the differential expression levels of autophagy-related genes. We performed PCA to evaluate the intra-group data, and the repeatability in GSE29431 was found to be satisfactory (Fig. 1A). Using the criteria of adjusted p < 0.05 and |log2FC| > 0.5, we identified 73 autophagy related genes. The volcano plots of the ARGs were presented in Fig. 1B (p < 0.05). Following the analysis of the GSE29431 with R software, 16 up-regulated genes and 57 down-regulated ARGs were presented in volcano plots (Fig. 1B and C). To validate the identified ARGs in the two cohorts, we conducted box plots to show up- or down-regulation of the most significant ARGs in terms of the basic characteristics of patients in the TCGA cohort, comparing HER2 + and normal samples. The top 10 up-regulated genes included BIRC5, ERBB2, CXCR4, VAMP7, ATF6, ATG16L1, TBK1, CDKN2A, CASP3, and PTK6, while top 10 down-regulated genes included FOS, CX3CL1, DLC1, BNIP3L, MYC, EGFR, FOXO1, HSPB8, GABARAPL1, and PPP1R15A (Fig. 2A and B).
Figure 1 Differentially expressed ARGs in HER2 + patient tissue and normal samples. (A) PCA for remarkable differences in ARGs between two cohorts. (B) Volcano plot of the differentially expressed ARGs. (C) Heatmap plot of the differentially expressed ARGs.
Figure 2 The boxplot of top 10 up or down-regulation ARGs in HER2 + patient tissue and normal samples. (A) The boxplot of up-regulated top 10 ARGs in HER2 + patient tissue and normal samples. (B) The boxplot of down-regulated top 10 ARGs in HER2 + patient tissue and normal samples (*p < 0.05, **p < 0.01, and ***p < 0.001. ns, not statistically significant).
Functional analysis of ARGs and validated hub genes
To further explore the functions of these ARGs, we conducted GO analysis with top 10 up or down-regulated genes. The top-ranked GO enrichment analysis revealed involvement in protein phosphatase binding, phosphatase binding, ubiquitin/ubiquitin-like protein ligase binding, autophagy, process utilizing autophagic mechanism, and macroautophagy. The pathway analysis showed that most ARGs were enriched in the process of autophagy and mitophagy, in terms of their biological process (BP), cellular component (CC), and molecular function (MF) (Fig. 3A, and B). To further analyze the potential biological relationships of these ARGs which played the most critical role, the KEGG analysis was performed (Fig. 3C and D). Additionally, the PPI interactions among differentially expressed ARGs were performed. The results showed these ARGs linker nodes CASP3, EGFR, ERBB2, FOS, FOXO1, MYC, CDKN2A, CXCR4, BIRC5, ATG16L1, GABARAPL1, PPP1R15A, TBK1, BNIP3L, PTK6, ATF6, CX3CL1, VAMP7 were interacted with each other (Fig. 3E and F). These top-ranked ARGs were mainly involved in various metabolic pathways, such as human cytomegalovirus infection, microRNAs in cancer, MAPK signaling pathway, endocrine resistance, indicating their potential as autophagy related hub genes.
Figure 3 GO, PPI, KEGG pathway analysis and identification of differentially expressed ARGs. (A-D) GO and KEGG pathway enrichments results of ARGs. (E-F) The PPI among differentially expression genes and interaction number of each ARGs.
Prognostic value of ARGs in TCGA
To assess the relationship between hub ARGs and prognosis of HER2 + patients, we conducted a prognostic correlation analysis. The prognostic value of the hub genes was evaluated by assessing the correlations between each gene and survival status in HER2 + TCGA cohort with employing Cox regression analysis. According to the survival data in HER2 + TCGA cohort, the Kaplan-Meier curves were used to compare the survival between low-risk and high-risk subgroups. The cohort offered DSS and OS information, allowing us to test the hub genes’ prognostic value. We found that most of ARGs had a significant impact on prognosis, and those could be potential risk factors. Among autophagy related hub genes, high expression of PPP1R15A (HR = 1.96, p = 0.041), VAMP7 (HR = 2.3, p = 0.035), PTK6 (HR = 2.42, p = 0.005), and CASP3 (HR = 2.1, p = 0.019) increased the risk of DSS, compared to the low expression group (Fig. 4A-D). Additionally, high expression of CASP3 (HR = 1.42, p = 0.035), PTK6 (HR = 1.58, p = 0.006), and VAMP7 (HR = 1.39, p = 0.046) increased the risk of OS compared to the low expression group. However, the expression of PPP1R15A did not show a significant difference between the low-risk and high-risk subgroups in terms of OS (Fig. 4E-H).
Figure 4 The prognostic value analysis of autophagy related hub genes in TCGA cohort. (A-D) Kaplan-Meier curves for comparison of the DSS between low- and high-risk groups. (E-H) Kaplan-Meier curves for comparison of the OS between low- and high-risk groups. (*p < 0.05, **p < 0.01, and ***p < 0.001. ns, not statistically significant).
Validation prognostic value of the hub genes in trastuzumab treatment cohort and its clinical characteristics
To validate the repeatability of the hub genes, we further assessed their prognostic robustness in a trastuzumab adjuvant treatment cohort of our clinical data. The expression levels of autophagy related genes PPP1R15A, VAMP7, PTK6, and CASP3 were further measured by qRT-PCR between pCR and non-pCR cohort, which received trastuzumab as neoadjuvant therapy for HER2 + breast cancer. As a result, expression levels of VAMP7 and PTK6 were significantly increased in HER2 + breast cancer samples in the non-pCR cohort (Fig. 5A). Additionally, VAMP7 and PTK6 staining were evaluated in a cohort of our breast cancer samples by immunohistochemistry. Consistent with the mRNA results, the positive staining of VAMP7 and PTK6 were mainly observed in the cytoplasm, and the expression levels of VAMP7, PTK6 were significantly increased in breast cancer patient in the non-pCR group (Fig. 5B-C). The clinicopathological variables of our clinical data were summarized in Table 1. In addition, we also made correlation analysis about VAMP7, PTK6 expression level in the tumor size, nodes status, ER status, PR status, metastasis status, vital and disease state. Results showed that high expression of VAMP7 and PTK6 were significantly associated with disease state to breast cancer (Table 1, p < 0.05). No associations were observed between VAMP7, PTK6 expression and other clinicopathological criteria.
Table 1
Clinical characteristics of the patients
Characteristic | levels | Numbers(N = 29) | VAMP7 | PTK6 |
+ | - | Pvalue | + | - | Pvalue |
Tumor size, n (%) | ≤ 20mm | 10 (34.4%) | 8 | 2 | 0.775 | 6 | 4 | 0.278 |
| > 20mm | 19 (65.6%) | 16 | 3 | | 15 | 4 | |
Nodes status, n (%) | 0 | 2 (6.9%) | 1 | 1 | 0.960 | 0 | 2 | 0.824 |
| 1 ~ 3 | 19 (65.5%) | 10 | 9 | | 6 | 13 | |
| 4 ~ 9 | 6 (20.7%) | 4 | 2 | | 2 | 4 | |
| ≥ 10 | 2 (6.9%) | 1 | 1 | | 1 | 1 | |
ER status, n (%) | Negative | 7 (24.1%) | 4 | 3 | 0.331 | 2 | 5 | 0.430 |
| Positive | 22 (75.9%) | 8 | 14 | | 10 | 12 | |
PR status, n (%) | Negative | 8 (27.6%) | 4 | 4 | 0.561 | 3 | 5 | 0.976 |
| Positive | 21 (72.4%) | 13 | 8 | | 8 | 13 | |
Metastasis, n (%) | No | 26 (89.7%) | 15 | 11 | 0.422 | 15 | 11 | 0.422 |
| Yes | 3 (10.3%) | 1 | 2 | | 1 | 2 | |
Disease state, n (%) | non-pCR | 15(51.7%) | 14 | 1 | 0.013* | 13 | 2 | 0.016* |
| pCR | 14(48.3%) | 6 | 8 | | 5 | 9 | |
Figure 5 Expression of PPP1R15A, VAMP7, PTK6, and CASP3 in HER2 + breast cancer patients treated with trastuzumab. (A) mRNA expression of PPP1R15A, VAMP7, PTK6, and CASP3 were measured using qRT-PCR. *p < 0.05, **p < 0.01, and ***p < 0.001. ns, not statistically significant. (B, C) Immunohistochemistry staining of VAMP7 and PTK6 in HER2 + breast cancer samples in the non-pCR group. Representative photomicrographs at 20× magnification and scale bar representing 100 µm.