Subtype recognition and identification of a prognosis model characterized by antibody-dependent cell phagocytosis-related genes in breast cancer

Data download and processing

Study cohorts A and B were based on data obtained from two publicly available datasets that derived from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Molecular Taxonomy of Breast Cancer International Consortium (METABRIC, https://www.cbioportal.org/study/summary?id=brca_metabric) databases, respectively [9, 10]. Cohort A consisted of 103 paracancerous samples and 903 breast cancer (BC) samples, which were used for gene expression profiling. On the other hand, cohort B included 154 paracancerous samples and 1826 BC samples for gene expression profiling. In order to evaluate the reliability of the prognostic model constructed in this study, cohort A was designated as the training set, while cohort B served as the validation set. This division allowed us to train the model on one dataset and validate its performance on an independent dataset.

Identification of the survival-related ADCP genes

To evaluate the prognostic significance of ADCP genes, Cox regression analysis was further adopted to evaluate the relationship of each gene with survival status in the TCGA cohort. To mitigate the risk of overlooking potentially important genes, a truncated p-value of 0.05 was set as a threshold, and a total of 130 survival-related genes were identified for further analysis.

Consensus clustering analysis of ADCP genes

The unsupervised clustering ‘Pam’ method was applied to the identification of different molecular subtypes based on survival-related ADCP gene expression [11]. The R package Consensus Cluster Plus performs this process and determines the number of clusters in Queue A, with 1000 times repeated to make sure classification stability [12]. Principal component analysis (PCA) was employed to illustrate the distribution differences. To evaluate the clinical significance of ADCP subtypes, their relationship with prognosis and other clinicopathological features (including age, stage and cancer subtypes) were examined. In different cohorts, the overall survival (OS) of different clusters was compared using the Kaplan-Meier (K-M) survival chart, with P ≤ 0.05 indicating statistical significance. Besides, the relationship between ADCP subtypes and other clinical variables was visualized using the Sankey plot and drawn using the ggalluvial R package [13]. The heat maps of different ADCP gene expression patterns were plotted by the pheatmap R package [14].

Identification and verification of the key ADCP genes

Based on the TCGA, the ‘limma’ R package was adopted to identify the differentially expressed genes (DEGs) between group1 and group2 [15]. The fold change cut off was 1, and the adjusted p value was less than 0.05. We used the glmnet R package in the TCGA cohort to use the LASSO Cox regression model to screen out the best prognostic biomarkers in 130 ADCP genes [16]. A total of 39 ADCP genes with non-zero coefficients were screened by 10-fold cross-validation. We overlapped the ADCP genes screened by DEGs and LASSO regression to screen out the key ADCP genes. Risk scoring formula:

$$\text{Risk}\text{score}={\sum }_{\text{i}=1}^{\text{n}}\text{coefi}\ast \text{id}$$ (1)

Where Coefi represents the coefficient and id represents the normalized count for each gene.

After constructing the risk score according to the LASSO regression model, the risk grouping and group were visualized through the Sankey diagram using the ggalluvial R package according to the survival situation.

Univariate Cox regression together with K-M plot were used to analyze the prognostic significance of key ADCP genes. P-values less than 0.05 indicated statistical significance. Subsequently, multivariate Cox regression model was utilized to explore whether ADCP genes could be combined with other clinicopathological features as independent prognostic factors.

Correlation of the ADCP genes with copy-number alterations and immune traits

Gene Set Cancer Analysis (GSCA) is a multi-omics online analysis tool based on TCGA data. It was used to analyze the expression and prognosis of pan-oncogenes in key ADCP genes (http://bioinfo.life.hust.edu.cn/web/GSCALite/) [17]. Results are displayed in the bubble plot by the ‘ggplot2’ package [18]. To analyze the expression characteristics of key ADCP genes in tumor and ADCP, t-test algorithm was used to compare their expression levels in different risk scores, groups and tumor groups in cohort A, respectively. The CAN (copy-number alterations) ratio as well as methylation level (β value) of key ADCP genes among TCGA-BC samples were verified by GSCA online tool. As a convenient online tool, TISIDB provides convenient access to correlation between genes and immune traits, such as lymphocytes, immunomodulators, and chemokines [19]. The scatter plots of ADCP genes with several tumor-infiltrating lymphocytes (TIL), immunosuppressants and immunostimulants were downloaded from TISIDB.

Construction and verification of risk model

For each patient, the risk score was predicted by the ‘predict’ function in the survival R package [20]. Patients were then divided into high- and low-risk groups according to the median value. To test the accuracy of the Cox regression model, C-indexes as well as receiver operating characteristic curves (ROC) were adopted. Prognostic factors, including C-indexes representing risk score, were visualized using histogram. ROC analysis was performed by the ‘pROC’ package [21], with the area under the curve (AUC) values calculated and compared to evaluate the performance of the model. In addition, the K-M survival analyses were conducted for risk scores. Based on training and test datasets, 1-, 3- and 5-year ROC curve were plotted using the R package timeROC [22].

Nomogram construction

To further investigate the prognostic significance of ADCP genes, the TCGA datasets were utilized to generate risk plots between high- and low-risk groups. To improve the predictive accuracy, a nomogram incorporating TNM stage, Cancer type and group was developed using survival and rms package [23]. To demonstrate the consistency between the predicted 1-, 3- and 5-year endpoint events and the corresponding actual outcomes, calibration plots were generated using rms package.

Evaluation of drug sensitivity

To investigate the molecular characteristics associated with drug sensitivity/resistance, the Genomics of Drug Sensitivity in Cancer (GSDC) database (https://www.cancerrxgene.org/) [24], a public resource for identifying cancer cell biomarkers, was utilized. The pRRophetic package was adopted to determine the sensitivity data of two BC groups towards various drugs [25]. Furthermore, the drug sensitivity of BC patients with different phenotypes was predicted from gene expression data.

Quantitative reverse transcriptase PCR (qRT-PCR)

The MCF-10A (human breast cancer paracancerous) and MDA-MB-453 (human breast cancer) cell lines were obtained (Shangcheng North Na Chuanglian Biotechnology Co., Ltd.). Control groups were created by combining MCF-10A with MDA-MB-453 cells to compare the expression of DEFB1, SIAH2, and SYT1 between normal and breast cancer cell lines. Total RNA was extracted using the Redzol kit (Beijing SBS Gene Technology Co., Ltd.), and qRT-PCR was performed using the SYBR® Premix Ex Taq™ II Kit. The relative mRNA expression levels were calculated using the 2^−ΔΔCt method, with β-actin as the internal reference gene. The forward primer sequences were as follows: (1) DEFB1: F-5′-CTCAGGTGGTAACTTTCTCA-3′; (2) SIAH2: F-5′-CACTTGACAGGCTGTTGCAC-3′; (3) SYT1: F-5′-CGCTTCGGCAGCACATATACTAAAATTGGAAC-3′. (4) GAPDH: F-5′-GAGTCAACGGATTTGGTCGT-3′. The reverse primer sequences were as follows: (1) DEFB1: R-5′-AAGCACTCCGGGTGATTCAG-3′; (2) SIAH2: R-5′-AAGCACTCCGGGTGATTCAG-3′; (3) SYT1: R-5′-TTGGTCAGCACAGATCATCG-3′. (4) GAPDH: R-5′-GATCTCGCTCCTGGAAGATG-3′.

Western blot (WB)

The WB analysis was performed following previously established protocols and the manufacturer’s recommendations. The concentrations of extracted proteins were determined using the Bradford method. Proteins were separated on a 10% SDS-polyacrylamide gel, transferred to PVDF membranes, and blocked with 5% skimmed milk. Primary antibodies against DEFB1 (1:500; catalog number ab115813; Abcam), SIAH2 (1:500; catalog number ab31234; Abcam), SYT1 (1:500; catalog number ab302627; Abcam), and GAPDH (1:500; catalog number ab181602; Abcam) were incubated with the membranes overnight at 4° C. Afterward, HRP-labeled secondary antibody was added and the membranes were washed with TBST. WB analysis was conducted using chemiluminescence, and the protein bands were visualized using a film.

Statistical analysis

Statistical analyses were conducted using R software (version 4.2.1) [26]. A two-sided P-value less than 0.05 was considered statistically significant. The K-M survival curve was constructed using the survival curve function ‘ggsurvplot’ in the R package survminer (0.4.2 version) [27], and the bilateral log-rank test was adopted to estimate the difference in OS between different groups. The cut-off value of ADCP genes, differentiating the high-risk group from the low-risk group, was determined as the median or calculated using the surv cutpoint function implemented in the R package survminer (0.4.2 version) [28]. The LASSO Cox regression model was used to identify the ADCP genes that exhibited a significant association with OS, and multivariate Cox regression analysis was performed using the identified ADCP genes to calculate hazard ratios (HR), 95 % confidence intervals (CI), and corresponding p-values to construct ADCP genes. Additionally, multivariate Cox regression analysis was used to evaluate the independent prognostic value of clinical indicators, including the ADCP genes, and a nomogram visual risk prediction map was formed after scoring each factor using the rms package of R software. To evaluate the consistency of nomogram, the calibration plot was generated using the rms package of R software [29].

Availability of data and material

Publicly available datasets were analyzed in this study. These data can be found here: https://portal.gdc.cancer.gov/, https://www.cbioportal.org/study/summary?id=brca_metabric. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Consent for publication

All authors approved the publication of the article.

Survival-related ADCP gene

To investigate the prognostic significance of ADCP genes, a gene signature was established using data from cohort A. The workflow for the construction of gene signature is depicted in Supplementary Table 1. Based on univariate Cox regression, 130 prognostic genes were screened out (p-value < 0.01, Supplementary Table 1).

Consensus clustering analysis of ADCP genes

Based on the expression levels of survival-related ADCP genes in the TCGA database, two distinct regulatory patterns were identified using an unsupervised clustering method. A total of 515 cases were classified into ADCP-related cluster 1 and 491 cases were classified into ADCP-related cluster 2 (Figure 1A). PCA analysis demonstrated that the patients can be classified into two distinct parts, providing further evidence of the existence of two significantly different subtypes (Figure 1B). To evaluate the survival difference between the two clusters, an analysis was conducted on a diverse range of patients from a variety of cohort A (TCGA) with complete survival information were analyzed. Subsequently, the relationship of these two clusters with various clinical features were analyzed (Figure 1D). The survival advantage of cluster 1 was higher compared to that of cluster 2. In cluster 1, the proportion of age less than 60, M0, N0-1, T1-2, Stagei-ii was higher. These suggest that ADCP genes may affect tumor development through some potential mechanisms. The heat map depicts the transcriptomic characteristics of different expressed ADCP genes in the two ADCP subtypes (Figure 1E).

Figure 1. Consensus clustering of ADCP-related genes in BC. (A) The consensus matrices and (B) PCA analysis were performed to assess the stability of clustering and explore coagulation subtypes. (C) Kaplan-Meier curves and (D) an alluvial diagram demonstrated the association between ADCP groups, tumor stage, age, grade, and survival state. (E) Unsupervised clustering of all coagulation-related genes in TCGA cohorts.

The immune landscape of ADCP subtypes

To explore the disparity in pathway enrichment analysis between the two clusters of ADCP in cohort A, GSEA was performed, with different immune infiltration patterns identified within the two subtypes. The enrichment histogram showed that cluster 1 significantly enriched hormone metabolism (including early estrogen response, late estrogen response) and inflammatory signal regulation pathways (allograft rejection, inflammatory response, kras signal and il6 jak stat3 signal, tnfa signal conduction through NFKB) (Figure 2A). At the same time, GSEA confirmed that there were differences in the immune pathways of ADCP clusters. The results showed that DEGs with higher expression levels in cluster 1 were significantly enriched in early estrogen response, allograft rejection pathway, inflammatory response, interferon response, kras signaling pathway and tnfa signal transduction through NFKB (Figure 2B). Given the strong correlation between ADCP subtypes and immune activity, the TME of the two clusters in cohort A was investigated (Figure 2C). Cluster 1 subtype is featured with high infiltration of Natural Kill cells (resting) and Macrophages, while cluster 2 subtype exhibits elevated infiltration levels of B cells, Plasma cells, T cells, Natural Kill cells, Dendritic cells, Mast cells and Neutrophils. Based on TCGA expression profile, the stromal score, immune score as well as ESTIMATE score of malignant tumor tissues were calculated through ESTIMATE algorithm. ESTIMATE produces a matrix score that measures the presence of tumor-associated matrices, as well as an immune score that reflects the level of immune cell infiltration. These scores are combined to produce an index called an ‘estimated score’ which provides a comprehensively estimation tumor purity. Compared with cluster 1, samples in cluster 2 also showed significantly higher estimated scores (Wilcoxon test, P < 0.05, Figure 2D). This trend was also observed for matrix scores as well as immune scores (Wilcoxon test, P < 0.05). Additionally, we studied the association between the two subtypes and major histocompatibility complex (MHC) and T cell stimulators. In addition to HLA-C, the expression levels of MHC as well as T cell stimulators exhibited a tendency to be higher in cluster 2 (Figure 3E, 3F).

Figure 2. Clinical significance and immune landscape of ADCP groups in the TCGA cohort. (A) GSVA and Gene Set Enrichment Analysis (GSEA) revealed the activation or inhibition of biological pathways and (B) significant enrichment in immune-associated processes. (C) Immune cell infiltration, (D) stromal and immune scores, and gene expression of HLA and MHC gene sets were analyzed between ADCP groups. (E) The expression box plot of HLA family genes between the two subtypes. (F) The expression box plot of inflammatory genes between the two subtypes. Statistical significance at the level of * <0.05, ** <0.01, and *** <0.001.

Figure 3. Identification of key ADCP targets. (A) Volcano plots and (B) LASSO regression model were used to identify differentially expressed genes (DEGs) and select key CRGs. (C) An alluvial diagram illustrated the changes in risk groups, ADCP groups, and survival state. (D) Univariate Cox regression analysis and (E) Kaplan-Meier curves evaluated the associations between infiltrate levels of ADCP targets and overall survival.

Recognition of key ADCP genes

Using the limma algorithm, a total of 432 degs (Differentially expressed gene) between group1 and group2 are identified under the filtering threshold of FDR q-value < 0.01 and the absolute value of logFC > 1 (Figure 3A and Supplementary Table 2). Prognostic models were constructed for 1006 BC patients with OS information in cohort A. LASSO Cox regression analysis was utilized to determine the best prognostic features based on 130 survival-related ADCP genes. After the variables were included in the LASSO Cox regression model with the smallest lambda, the genes of 39 ADCP-related features were chosen to construct the ADCP-related risk scoring model. Three key ADCP genes (DEFB1, SIAH2 and SYT1) by overlapped DEGs from cohort A and the model ADCP genes were identified (Figure 3B). The relationship of the expression levels with 3 ADCP-related signatures together with OS are also presented in the forest plot (Figure 3D). The expression of three key ADCP genes is of significance for survival (Figure 3E).

Multidimensional analysis of key ADCP genes

We studied the relationship of CNV with immune infiltration in BRCA. Additionally, we explored the association of gene methylation with immune infiltration. CNV and methylation of key ADCP genes were closely associated with the infiltration of key immune cells including T and B cells (Figure 4A, 4B). At the same time, based on cohort A, the gene expression levels of 3 key ADCP genes were analyzed. In the data set, only DEFB1 was down-regulated in group1, while SIAH2 and SYT1 were up-regulated in group2 (Wilcoxon test, P > 0.05). Only DEFB1 was up-regulated in the high-risk group, while in the low-risk group, SIAH2 together with SYT1 were down-regulated (Wilcoxon test, P > 0.05). Only DEFB1 was up-regulated in the tumor group; in contrast, SIAH2 and SYT1 were down-regulated in the low-risk group (Wilcoxon test, P > 0.05, Figure 4C). Among the three key ADCP genes, the frequency of copy number alterations (CNA) of SIAH2 was higher compared to that of the others. Specifically, CNA deletions were predominant among all types (Figure 4D). Methylation analysis showed that the beta value of SIAH2 was higher in the tumor group than in the normal group, and DEFB1 was the opposite (Figure 4E). It is well known that gene expression is negatively correlated with the level of methylation, while CNA has a positive effect on gene expression. In summary, CNA and methylation may be may be important in the up-regulation of SYT1 expression in BC.

Figure 4. Validation of key ADCP genes in expression level. (A, B) Correlation analyses were conducted between CNV/methylation and immune infiltrates in BC. (C) Expression levels and (D, E) CNA percentages of key ADCP genes were compared between ADCP groups, risk groups, and cancer types. Methylation levels of key ADCP genes in BC and normal samples were also examined. Statistical significance at the level of ns ≥ 0.05, *** <0.001 and **** <0.0001.

ADCP-related signatures for the prognostic prediction of BC

LASSO algorithm was used to establish a risk model. Finally, 39 genes related to prognosis were identified, and used to construct models based on risk scores using the training (n = 1006) and test (n = 1980) data sets of cohorts A and B patients, respectively. Survival analysis showed that a higher risk score of the training and the test sets was related with a lower survival rate (p < 0.0001) (Figure 5A, 5B). A time-varying ROC curve was generated for the assessment of the sensitivity of the model. The 3-year, 5-year, and 10-year AUCs were found to be 0.743,0.754, and 0.79 in the training set, respectively (Figure 5C); in contrast, they were 0.546,0.678, and 0.716 in the test set, respectively (Figure 5C).

Figure 5. Survival evaluation of ADCP-related risk scoring model in TCGA and METABRIC datasets. (A) Kaplan-Meier survival curves, (B) patient subgrouping based on risk scores, and (C) ROC curves assessed the performance of the ADCP-related risk scoring model in predicting overall survival at 1, 3, and 5 years.

The ADCP group served as an independent prognostic factor in BC

Since ADCP genes are significantly associated with high malignancy and advanced tumors of BC, univariate and multivariate Cox regression analysis were carried out to determine the prognostic significance of ADCP genes for BC patients. ADCP group, age, TNM stage, stage together with risk score were included as covariates. Results showed that ADCP group, age, TNM stage, stage as well as risk score were independent prognostic factors for BC patients (Figure 6A, 6B). we constructed a nomogram by combining independent prognostic factors, serving as a clinically relevant quantitative method tool for predicting the mortality of individual BC patient (Figure 6C). Based on the c-exponential curve of different variables over time in the TCGA cohort, nomogram performed best compared to other single factors (Figure 6D). Add up the scores of each prognostic parameter and assign a total score to each patient. The higher the total score, the worse the prognosis of patients. The modal diagram has similar performance to the ideal model (Figure 6E). In addition, the DCA curve also indicated that the nomogram had good stability and reliability net benefit curve in age compared with other clinical factors (Figure 6F).

Figure 6. A nomogram was developed to predict the probability of 1-, 3-, and 5-year overall survival (OS) in the training cohort. (A, B) Univariate and multivariate analyses were conducted, including the ADCP-related risk scoring model and clinical factors. (C) The comprehensive nomogram provided predictions for BC patients’ OS probabilities in the TCGA dataset. (D) The time-dependent c-index plot compared the performance of the nomogram with other clinical factors. (E) Calibration plots assessed the accuracy of the nomogram’s predictions for 1-, 3-, and 5-year OS in the TCGA cohort. (F) Decision curve analysis evaluated the clinical utility of the nomogram and other factors for 1-, 3-, and 5-year risk assessment.

ADCP groups for the prediction of the chemotherapeutic response

We evaluated the chemotherapy response and drug resistance of patients in the ADCP group. Figure 7 shows the sensitivity of two ADCP subtypes to six anticancer drugs (AZD8055, A.443654, AMG.706, AKT.inhibitor.VIII, ABT.888, ATRA). Results showed that the IC50 level of group2 was higher compared to group1 (Figure 7A–7F), and small molecule drugs with therapeutic effects on BC could be found according to the results of drug sensitivity. Three-dimensional structural tomography of AZD8055, A.443654, AMG.706, AKT.inhibitor.VIII, ABT.888 and ATRA was found in PubChem (Figure 7G).

Figure 7. The ADCP-related groups were evaluated for their ability to predict therapeutic benefit. (A–F) The chemotherapy response of two ADCP-related groups to six common chemotherapy drugs was analyzed. (G) The 3D structure tomographs of six candidate small-molecule drugs for BC were examined. Statistical significance at the level of *** <0.001.

Gene expression level verification via qRT-PCR and WB

We verified the mRNA and protein levels of DEFB1, SIAH2 and SYT1 in BC cell lines and adjacent cell lines by qRT-PCR and WB. Results of qRT-PCR is shown in Figure 8A. In contrast with the normal control, DEFB1 level in BC cell line MDA-MB-453 was lower, but the SIAH2 and SYT1 levels were higher. The protein expression levels of the three key ADCP-genes based on WB analysis were consistent with the results of qRT-PCR experiments (Figure 8B, 8C).

Figure 8. (A) RT-qPCR and (B, C) Western blotting analysis of DEFB1, SIAH2, and SYT1 expression in BC cell lines and adjacent cell lines. Statistical significance at different levels was reported for the results. Statistical significance at the level of * <0.05, ** <0.01, and *** <0.001.