Identification of a Risk Predictive Signature Based on Genes Associated with Tumor Size and Lymph Node Involvement in Breast Cancer

Abstract

Background:

Breast cancer is a heterogeneous disease. Small tumors with extensive lymph node involvement (STEL) in breast cancer often reflect a biologically aggressive phenotype and poor prognosis. The aim of this study was to identify key genes associated with STEL and investigate their prognostic values in breast cancer.

Methods:

RNA sequence data from breast cancer specimens were acquired from The Cancer Genome Atlas (TCGA) database for differential analysis. Weighted gene correlation network analyses (WGCNA) were performed to identify coexpressed gene modules associated with tumor size and lymph node metastases. Gene set enrichment analysis (GSEA) was employed to investigate the biological functions of the identified genes. A combination of LASSO and Cox regression analyses was conducted to establish a risk predictive signature, and time-dependent receiver operating characteristic (tdROC) and Kaplan—Meier analyses were used to evaluate its prediction precision. Quantitative RT-PCR was employed to validate the expression levels of the key genes from the signature set.

Results:

A total of 2777 genes from three coexpressed gene modules were identified by WGCNA, and 880 differentially expressed genes were identified by transcriptome analyses. The 63 overlapping genes identified by both methods were considered STEL-associated genes, and a 9-gene risk-predictive signature was established based on them, with AUCs at 3, 5, and 7 years reaching 0.810, 0.811, and 0.753, respectively.

Conclusion:

This study demonstrated the transcriptomic profile of STEL breast cancer and successfully established a risk predictive signature with satisfactory accuracy. These findings may provide insights in to the genetic etiology of breast cancer.

Introduction

Breast cancer is one of most common malignancies and the most serious threat to women's health worldwide (Siegel et al., 2021). Although the mortality of breast cancer is decreasing due to improvements in diagnostic and therapeutic methods, its morbidity continues to increase each year; thus, breast cancer remains a leading cause of cancer-associated death in women (Allemani et al., 2018; Cardoso et al., 2018). Breast cancer is widely recognized as a heterogeneous disease, and the prognosis varies greatly among patients (Brenton et al., 2005); therefore, metastasis and outcomes of breast cancer are difficult to predict.

To date, many efforts have been devoted to the identification of prognostic factors for breast cancer. Several indicators have been well established to be associated with breast cancer-specific mortality (BCSM), including local tumor size, degree of axillary lymph node (LN) involvement, histological grade, lymphovascular invasion (LVI), and hormone receptor and HER2 status (Perou et al., 2000; Schoppmann et al., 2004; Goldhirsch et al., 2005). Of these, local tumor size and regional LN status are two of the most important prognostic determinants and constitute the American Joint Committee on Cancer (AJCC) staging system (Singletary et al., 2002), along with the presence of distant metastases.

In most cases, increasing tumor size is closely associated with a high risk of regional LN involvement (Gajdos et al., 1999); however, the relationships between these two determinants are not always linear. Accumulating evidence has revealed that the acquisition of metastatic potential and distant spread may occur during the early stage of cancer development (Husemann et al., 2008; Wo et al., 2011). In the context of extensive LN involvement, a small tumor size is paradoxically a surrogate for aggressive biological disease and poor prognosis (Wo et al., 2011; Ryu et al., 2016). Conversely, large tumors (e.g., ≥5 cm) that fail to metastasize to regional LNs may reflect a biologically indolent phenotype and thus can be associated with a relatively lower risk of distant spread and BCSM (Yu et al., 2012).

It could reasonably be inferred that the gene expression profiles of the two distinct types of breast cancer may differ considerably. The differentially expressed genes (DEGs) may play pivotal roles in the acquisition of biologically aggressive phenotypes and ultimately contribute to tumor development and metastasis in breast cancer. However, few studies have focused on transcriptome analysis between small tumor with extensive LN involvement (STEL), and large tumor with negative LN count (LTNL) breast cancers. Therefore, we carried out this study and performed a series of bioinformatics analyses to identify STEL-associated genes and investigate their prognostic values in breast cancer.

Materials and Methods

Data acquisition, processing, and differential expression analysis

The RNA-seq data and corresponding clinical information of STEL and LTNL breast cancer samples from The Cancer Genome Atlas (TCGA) database were downloaded from the GDC (Genomic Data Commons) Data Portal (URL: https://portal.gdc.cancer.gov/). A total of 24 STEL and 47 LTNL breast cancer samples from TCGA database were included in this study, and the clinicopathological characteristics of the patients are shown in Table 1. All the patients involved in this study were primarily diagnosed by surgery or biopsy, with no neoadjuvant treatment. STEL and LTNL were defined by the patients' T and N stage.

Table 1.

Characteristics of Patients from The Cancer Genome Atlas Database

Characteristics	STEL (n = 24; Mean ± SD/%)	LTNL (n = 47; Mean ± SD/%)	p
Age	59.88 ± 13.17	60.79 ± 14.68	0.799
Menopause status
Peri	0 (0%)	2 (4.3%)	0.172
Post	20 (83.3%)	30 (63.8%)
Pre	4 (16.7%)	9 (19.1%)
Unknown	0 (0%)	6 (12.8%)
Vital status
Alive	21 (87.5%)	44 (93.6%)	0.004
Dead	3 (12.5%)	3 (6.4%)
Survival days	523 ± 551.31	767.34 ± 1188.91	0.344
Histological type
IDC	16 (66.7%)	21 (44.7%)	0.106
ILC	6 (25%)	24 (51.1%)
Other	2 (8.3%)	2 (4.3%)
ER status
Negative	6 (25%)	7 (14.9%)	0.208
Positive	15 (62.5%)	38 (80.9%)
Unknown	3 (12.5%)	2 (4.3%)
PR status
Negative	7 (29.2%)	13 (27.7%)	0.183
Positive	13 (54.2%)	32 (68.1%)
Unknown	4 (16.7%)	2 (4.3%)
HER2 status
Equivocal	2 (8.3%)	4 (8.5%)	0.876
Negative	15 (62.5%)	28 (59.6%)
Positive	2 (8.3%)	7 (14.9%)
Unknown	5 (20.8%)	8 (17%)

ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; IDC, infiltrating ductal carcinoma; ILC, infiltrating lobular carcinoma; LNM, lymph node metastasis; PR, progesterone receptor; SD, standard deviation; STEL, small tumor with extensive lymph node involvement.

Specifically, patients with pathological T0 ∼ 1, N2 ∼ 3, and M0 were considered to have STELs, while patients with pT3N0M0 or pT4N0M0 were considered to have LTNLs. Inflammatory breast cancer samples were excluded. The raw count values of RNA-seq data were converted to transcripts per million (TPM) values for normalization, and the missing values were imputed by the K-nearest neighbor (KNN) imputation procedure. Differential expression analysis was conducted by the R package DESeq2 (Love et al., 2014). We used absolute log₂-fold change (FC) and the false discovery rate (FDR) as evaluation indicators for DEGs, and |log2FC| > 1 and FDR <0.05 were set as cutoff criteria for DEGs.

Weighted gene correlation network analysis

To explore the roles of STEL-associated genes in breast cancer, we performed a multistep bioinformatic analysis, and the design of this study is briefly shown as a flowchart in Figure 1. First, a gene coexpression network was constructed by conducting weighted gene correlation network analysis (WGCNA) on 24 STEL and 47 LTNL breast cancer samples. The whole set of coding and noncoding genes with at least one transcript was involved in this analysis. A scale-free network was built by using the R package WGCNA (Langfelder and Horvath, 2008) for coexpression gene module identification. A hierarchical clustering analysis was performed to show the adjacency between different module eigengenes (MEs).

FIG. 1.

Flowchart of the study design and analysis.

Moreover, the correlations between LN metastasis, tumor size, and MEs were calculated with a Pearson correlation analysis, and MEs with absolute coefficient r > 0.3 and p-value <0.05 were considered key modules associated with STEL. Finally, the overlapping genes from both key modules and DEGs were considered STEL-associated genes.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was used to further investigate the biofunctions and signaling pathways associated with STEL. The R package clusterProfiler (Wu et al., 2021) was used to conduct this analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set information was acquired from the Molecular Signatures Database (MSigDB). All genes involved in GSEA were reordered by their log₂FC according to the differential expression analysis. Thresholds of absolute normalized enrichment score (NES) >1 and FDR <0.05 were used to identify significantly enriched GO terms and KEGG pathways.

Training and testing of the risk predictive model

To further investigate the prognostic values of the identified STEL-associated genes in breast cancer, we performed a combination of least absolute shrinkage and selection operator (LASSO) and Cox regression analyses. The vital status and survival time information of the patients were acquired from TCGA database, and samples with incomplete clinical information were removed from the cohort. Subsequently, the whole set of breast cancer patients (n = 681) was randomly divided into two subsets for risk model training (n = 341) and validation (n = 340).

To narrow the range of the target genes, LASSO regression was conducted by using the R package glmnet, and the penalty parameter lambda (λ) was automatically determined by machine learning. Thereafter, multivariate Cox regression analysis was performed for risk predictive signature identification. A formula was established for risk score calculation: Risk score = $\sum C o e f i \times x i$ . Coef and x represent the coefficients obtained from multivariate Cox regression and the expression levels of the selected genes, respectively. The cutoff value between high and low risk was the median risk score of the training set. Moreover, time-dependent receiver-operating characteristic (tdROC) and Kaplan‒Meier analyses were used to evaluate the precision of the risk model in both the training and validation sets. Finally, univariate and multivariate Cox regression analyses with risk scores and other clinical features were carried out to determine independent risk factors.

Clinical sample collection and quantitative real-time PCR

Clinical samples of 35 patients, including 16 STELs and 19 LTNLs, were collected from patients who were primarily diagnosed between 2019 and 2021 at the Affiliated Hospital of Jiangnan University. Specimens were kept in the gas phase of liquid nitrogen immediately after surgery or biopsy. Total RNA was extracted from the specimens by using the E.Z.N.A. Total RNA Kit I (Omega) following the manufacturer's protocol. Reverse transcription PCR was performed to synthesize complementary DNA by using the Prime Script RT Reagent Kit (Takara Bio, Japan). Real-time PCR was performed on a CFX96 Touch Real-Time PCR Detection System (Bio-Rad) to measure the relative expression levels of the target genes using the comparative Ct method. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as a normalization control.

The sequences of the primers are shown in Supplementary Table S1. Each patient enrolled in this study signed an informed consent form, and this study was approved by the Ethics Committee of the hospital.

Statistics

Kaplan‒Meier analysis was used to evaluate the overall survival (OS) of different groups of patients. Univariate and multivariate Cox regression analyses were performed to identify prognostic factors among the clinicopathological features. The significance of differences between gene expression levels obtained from qRT-PCR was assessed by one-way analysis of variance. All statistical analyses were conducted in R software (R version 4.1.2), and a p-value (or adjusted p-value) <0.05 was considered to be statistically significant.

Results

Identification of STEL-associated genes

WGCNA and differential expression analysis were used to identify genes associated with STEL. The whole sets of 19213 protein-coding and 20730 noncoding genes were included in the analyses. To construct a scale-free gene network, the soft threshold was selected as 6 according to Figure 2A. The result of hierarchical clustering analysis is shown as a dendrogram (Fig. 2B). A total of 47 coexpression gene modules were identified, and the adjacency between different modules is shown in Figure 2C. Furthermore, Pearson correlation analysis was conducted to evaluate the associations between gene modules and clinical traits, including tumor size and lymph node metastasis (LNM). Three gene modules, MEdarkgrey (n = 568), MEdarkgreen (n = 600), and MEbrown (n = 1609), were ultimately identified, with a significant positive correlation with LNM and a negative correlation with tumor size (Fig. 2D). Thus, genes from these three modules were considered to be associated with STEL.

FIG. 2.

Identification of STEL-associated genes by WGCNA and differential expression analysis. (A) Screening of the soft threshold for WGCNA. (B) Identification of coexpression gene modules and visualized by a cluster dendrogram. (C) Eigengene adjacent heatmap shows the association between different modules. (D) Correlation analysis between clinical traits and gene modules; MEdarkgrey, MEdarkgreen, and MEbrown are significantly correlated to tumor size and LNM. (E) Volcano plot of the differential expression analysis between STEL and LTNL breast cancer samples; red dots represent for upregulated genes, and blue dots represent for downregulated genes. (F) Upset plot shows the intersection of different gene sets. LNM, lymph node metastasis; LTNL, large-sized tumor with negative lymph node count; ME, module eigengene; STEL, small tumor with extensive lymph node involvement; WGCNA, weighted gene correlation network analysis.

Differential expression analysis was subsequently carried out to identify DEGs between STEL and LTNL breast cancer samples. According to the criteria mentioned above, 885 DEGs were identified, including 367 upregulated and 518 downregulated genes (Fig. 2E). A total of 63 overlapping genes between the identified coexpression gene modules and DEGs were considered STEL-associated genes (Fig. 2F) and were further used for functional enrichment and prognostic analyses. In addition, the expression levels of the 63 identified genes in different clinicopathological groups are presented in a heatmap (Fig. 3).

FIG. 3.

Heatmap of the STEL-associated gene expression in TCGA breast cancer samples. TCGA, The Cancer Genome Atlas.

Functional and pathway enrichment analyses

To characterize the functional roles of the STEL-associated genes, GSEA was conducted for enrichment analysis. According to the results, multiple terms and pathways associated with cell adhesion, immune response, and tumorigenesis were significantly enriched (Fig. 4A, B). The running enrichment score curves of representative GO terms and KEGG pathways are shown in Figure 4C and D, respectively. For GO terms under the biological process (BP) category, regulation of cell adhesion (NES = −1.669), mammary gland development (NES = −1.799), and activation of immune response (NES = −1.500) were included in the top significant terms, and the P53 signaling pathway of KEGG (NES = −1.660) was also significantly enriched. These terms and pathways were all negatively correlated with STEL-associated genes, indicating the pivotal roles of these genes in breast cancer tumorigenesis and metastasis.

FIG. 4.

Gene set enrichment analysis for STEL-associated genes. (A) Dot plot of top 15 enriched GO terms. (B) Dot plot of top 15 enriched KEGG pathways. (C, D) Representative RES curves of GO terms and KEGG pathways significantly enriched, including regulation of cell adhesion, mammary gland development, activation of immune response, and P53 signaling pathway. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RES, running enrichment score.

Screening of the risk predictive signature in STEL-associated genes

To explore the prognostic values of the STEL-associated genes in breast cancer, a combination of LASSO and Cox analyses was conducted. LASSO regression was performed to narrow the range of the 63 STEL-associated genes (Fig. 5A), and multivariate Cox regression was subsequently used to identify a gene signature that could predict the prognosis of breast cancer. Thus, a risk predictive model based on nine STEL-associated genes was successfully constructed, and the formula for risk score calculation was as follows: risk score = (−0.00167*KRT6A)+(0.14177*LINC00839)+(0.03247*MMADHC)+(1.09877*NR0B2)+(−6.75219*PHB2)+(0.00947*SOD1)+(0.60041*SUN3)+(0.03683*VPS35)+(0.02745*YKT6). The cutoff value between high and low risk was set to the median risk score of samples from the training set (RS cutoff = 4.05820, Table 2).

FIG. 5.

Training of the risk predictive signature based on STEL-associated genes. (A) The process of STEL-associated gene selection using LASSO regression; the penalty parameter lambda (λ) was automatically decided by machine learning. (B) Risk score and vital status distribution in low- and high-risk groups from the training set according to the risk model. (C) Kaplan-Meier analysis of the low- and high-risk groups. (D) Time-dependent ROC curves for 3-, 5-, and 7-year survival prediction. LASSO, least absolute shrinkage and selection operator; ROC, receiver-operating characteristic.

Table 2.

Coefficients of the Five Identified Genes from the Risk Model

Genes	Coefficients
KRT6A	−0.00167
LINC00839	0.14177
MMADHC	0.03247
NR0B2	1.09877
PHB2	−6.75219
SOD1	0.00947
SUN3	0.60041
VPS35	0.03683
YKT6	0.02745

The cutoff value of risk score is 4.05820.

The distribution of risk score and survival status is shown in Figure 5B. Kaplan‒Meier analysis showed that the OS of high-risk patients was significantly worse (p < 0.001, Fig. 5C), and tdROC analysis verified the predictive accuracy of the risk model, with area under curves (AUCs) for 3-, 5-, and 7-year survival reaching 0.810, 0.811, and 0.753, respectively (Fig. 5D).

Validation of the prognostic signature

Kaplan‒Meier and tdROC analyses were also performed in the validation set to confirm the prognostic value of the identified prognostic signature. The distribution of the risk score and vital status of patients in the validation set is shown in Figure 6A. For Kaplan‒Meier analysis, the OS of the high-risk group was significantly worse than that of the low-risk group (p = 0.0034, Fig. 6B), and the precision of the signature in predicting the prognosis of breast cancer was satisfactory according to the ROC curves. The AUCs for 3-, 5-, and 7-year survival were 0.828, 0.674, and 0.684, respectively (Fig. 6C).

FIG. 6.

The validation of the risk predictive signature. (A) Risk score and vital status distribution in the validation set. (B) Kaplan-Meier survival analysis of the low- and high-risk groups from the validation set. (C) ROC curves for 3-, 5-, and 7-year survival prediction in the validation set. (D) Univariate and Multivariate Cox regression involving clinicopathological features and the risk score.

Moreover, to investigate whether the risk predictive signature was an independent prognostic factor of breast cancer, univariate and multivariate Cox regressions were performed, and multiple clinicopathological features were involved (Fig. 6D). The results of univariate (hazard ratio [HR] = 5.166, 95% confidence interval [CI]: 2.346-11.158; p < 0.001) and multivariate (HR = 4.839, 95% CI: 1.944-12.044; p < 0.001) analyses confirmed the prognostic value of the signature.

Furthermore, Kaplan‒Meier analysis was performed to evaluate the OS of the high- and low-expression groups of each gene from the risk predictive signature (Fig. 7A-F). The results showed that the OS of the KRT6A and PHB2 high-expression groups was more beneficial than that of the low-expression groups, while the high-expression groups of MMADHC, NR0B2, VPS35, and YKT6 had a worse prognosis than the low-expression groups. In addition, the expression profiles of the nine genes from the signature were validated by qRT‒PCR in clinical specimens. The relative expression levels of LINC00839, MMADHC, NR0B2, SOD1, SUN3, VPS35, and YKT6 were significantly higher in STEL breast cancer samples, and the expression levels of KRT6A and PHB2 were significantly lower (Fig. 7G). These results were consistent with the differential expression analysis.

FIG. 7.

Survival analysis and qRT-PCR of the nine genes from the signature. (A-F) Kaplan-Meier survival analysis of KRT6A, MMADHC, NR0B2, PHB2, VPS35, and YKT6. (G) Relative expression of the nine genes from the signature in clinical samples. qRT-PCR, quantitative reverse transcriptional polymerase chain reaction (*p < 0.05; **p < 0.01; ***p < 0.001).

Discussion

Breast cancer is the most prevalent tumor type worldwide. Despite distinct progress in treatment, breast cancer remains a leading cause of female cancer-related death (Sung et al., 2021). Breast cancer is known as a heterogeneous disease comprising discrepant subtypes that vary evidently in accordance with molecular biology and clinical features (Pellacani et al., 2019). Currently, breast cancer is divided into four subtypes, luminal A, luminal B, Her-2 positive, and triple-negative, based on immunohistochemical results of specific molecules, including estrogen receptor, progesterone receptor, HER-2, and Ki-67 (Waks and Winer, 2019). Such biological heterogeneity among different subtypes has brought about distinct differences in clinical outcome, response to therapeutic approaches, and patient survival (Wilcock and Webster, 2021).

The cancer staging system was formulated based on tumor size (T), lymph node metastasis (N), and distant metastasis (M) by the AJCC and the Union for International Cancer Control (UICC). Tumor size and LN metastasis have been the most commonly used methods and the most valid prognostic indicators in assessing the prognosis of breast cancer (Giuliano et al., 2018). Due to the prognostic value, an increasing number of researchers have evaluated the correlation between tumor size and nodal metastasis. The correlation between tumor size and LN status varies significantly among different breast cancer subtypes (Min et al., 2021). The tumor areas of 106 patients with breast invasive ductal cancer measured with contrast-enhanced ultrasound could predict regional LNM and N stage (Wang et al., 2014).

A retrospective study of 819,647 patients with invasive breast tumors showed that the relationship between tumor size and LN status is not linear (Sopik and Narod, 2018). However, a larger tumor size might not always reflect an increased risk of axillary LN metastasis, specifically in triple-negative breast cancer subtypes (Yaman et al., 2012; Altundag, 2019). Patients with small tumors and extensive LN status or with large tumors and negative LN status have been commonly found in clinical practice. The traditional TNM stages might not always be suitable for all patients due to individual differences. Therefore, finding genes driving these unusual clinical features could provide a basis for individualized diagnosis and treatment.

In this study, to identify genes driving different pathological phenotypes of STEL and LTNL breast cancer, we performed a series of bioinformatics analyses of breast cancer datasets from TCGA. First, three gene modules positively correlated with LNM and negatively correlated with tumor size were identified by WGCNA. Then, in combination with DEG analysis, 63 overlapping genes between the identified coexpression gene modules and DEGs were determined to be STEL-associated genes. GSEA was used to portray the functional roles of the STEL-associated genes. The results showed that the STEL-associated genes were significantly enriched in the regulation of cell adhesion, mammary gland development, activation of the immune response, and P53 signaling pathway BP, which indicated that all the above signaling pathway dysregulation may cause abnormal pathological phenotypes, such as enhanced LNM.

To evaluate the value of STEL-associated genes in predicting the clinical outcome of breast cancer, a risk predictive model was constructed through LASSO regression and multivariate Cox regression and included the following nine genes: KRT6A, LINC00839, MMADHC, NR0B2, PHB2, SOD1, SUN3, VPS35, and YKT6.

Keratin 6A (KRT6A), a type II cytokeratin, was recently found to be hyperexpressed in various types of cancer. However, a study reported that Krt6a expression in mammary epithelial progenitors did not favor tumorigenesis initiated by ErbB2 (Yang et al., 2020) in mice. In our study, we found that KRT6A expression was lower in STEL patients, and high KRT6A expression predicted good clinical outcome in breast cancer. NR0B2 (nuclear receptor subfamily 0 group B member 2), an orphan nuclear receptor, is a favorable prognostic factor in patients with liver, kidney, lung, and urinary bladder cancers but is a negative factor in patients with colon, thyroid, uterine, and head and neck cancers (Zhu et al., 2021).

NR0B2 functions as a STEL phenotype-associated gene and is negatively related to patient prognosis in breast cancer. Prohibitin 2 (PHB2) has been studied as a tumor suppressor in breast cancer (Kim et al., 2015), consistent with our results. Superoxide dismutase 1 (SOD1) might assist in eliminating excessive ROS produced during the process of metabolic reprogramming of cancer cells, which is conducive to the malignant transformation of tumors (Papa et al., 2014). Upregulated SOD1 in STEL breast cancer patients might also be involved in the occurrence of LNM, a possibility that needs further experimental confirmation. Vacuolar protein sorting-associated protein 35 (VPS35) was identified as an autophagy-related encoding gene in breast cancer.

VPS35 upregulation indicated worse prognosis in both progression-free and OS and promoted the progression of breast cancer (Li et al., 2021). YKT6 (YKT6 V-SNARE Homolog) was reported to be a tumor promoter in nonsmall cell lung cancer (Lee et al., 2021), breast cancer (Ooe et al., 2007), pancreatic cancer (Sun et al., 2020), hepatocellular carcinoma (Xu et al., 2021) and oral squamous cell carcinoma (Yang et al., 2021). In accordance with previous reports, we also found that the expression levels of VPS35 and YKT6 were increased in STEL breast cancer patients and correlated with poor clinical outcome in breast cancer, indicating their involvement in the development of this tumorigenic phenotype. LINC00839, the only noncoding RNA identified, was reported to play protumor roles in different kinds of cancer (Yu et al., 2021; Zhou et al., 2021; Zhang et al., 2022).

LINC00839 was also upregulated in chemoresistant breast cancer cells and tissues and was associated with a poor prognosis (Chen et al., 2020). There has been no report regarding the roles of MMADHC and SUN3 in breast cancer, which need to be explored in the future. Combining the results of the present study, all nine genes were strongly associated with the STEL phenotype, and the risk model based on these genes provided high precision in predicting the prognosis of breast cancer. This evidence suggests that genes aberrantly expressed in the STEL phenotype play pivotal roles in breast cancer development and lead to poor outcomes in breast cancer. Further exploration is still needed on the aspects of function and mechanism at the experimental level.

To the best of our knowledge, our study is the first to evaluate the association between extensive LN involvement, breast cancer, and candidate genes. Moreover, a risk model based on these genes was constructed with high precision and accuracy to predict prognosis for breast cancer. These findings may provide a new theoretical basis for the individual diagnosis and treatment of breast cancer. However, our study had several limitations. First, since the proportions of STEL and LTNL breast cancer patients were small, the number of clinical samples analyzed by qRT‒PCR was limited. Multicenter clinical studies are still required to further validate the effects of these genes on STEL breast cancer. Second, the STEL-associated genes were primarily identified using bioinformatics methods in the present study. Wet-lab experiments should be performed to further investigate the functional roles of these genes in vitro and in vivo.

Footnotes

Ethical Approval

All procedures performed in studies involving human participants were in accordance with the Declaration of Helsinki 1964 and its later amendments or comparable ethical standards. This article does not contain any studies with animals performed by any of the authors.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

This work was supported by grants from the National Natural Science Foundation of China (81902932).

Supplementary Material

References

Allemani

, Matsuda

, Di Carlo

, et al. (2018) Global surveillance of trends in cancer survival 2000-14 (CONCORD-3): analysis of individual records for 37 513 025 patients diagnosed with one of 18 cancers from 322 population-based registries in 71 countries. Lancet, 391:1023-1075.

Altundag

(2019) Larger tumor size detected by sonography might not always reflect increased risk of axillary lymph node metastasis in patients with breast cancer. J Ultrasound Med, 38:2521.

Brenton

, Carey

, Ahmed

, et al. (2005) Molecular classification and molecular forecasting of breast cancer: ready for clinical application?. J Clin Oncol, 23:7350-7360.

Cardoso

, Senkus

, Costa

, et al. (2018) 4th ESO-ESMO International Consensus Guidelines for Advanced Breast Cancer (ABC 4) dagger. Ann Oncol, 29:1634-1657.

Chen

, Shen

, Zhu

, et al. (2020) A nuclear lncRNA Linc00839 as a Myc target to promote breast cancer chemoresistance via PI3K/AKT signaling pathway. Cancer Sci, 111:3279-3291.

Gajdos

, Tartter

, Bleiweiss

(1999) Lymphatic invasion, tumor size, and age are independent predictors of axillary lymph node metastases in women with T1 breast cancers. Ann Surg, 230:692-696.

Giuliano

, Edge

, Hortobagyi

(2018) Eighth Edition of the AJCC cancer staging manual: breast cancer. Ann Surg Oncol, 25:1783-1785.

Goldhirsch

, Glick

, Gelber

, et al. (2005) Meeting highlights: international expert consensus on the primary therapy of early breast cancer 2005. Ann Oncol, 16:1569-1583.

Husemann

, Geigl

, Schubert

, et al. (2008) Systemic spread is an early step in breast cancer. Cancer Cell, 13:58-68.

10.

Kim

, Yoshimaru

, Chen

, et al. (2015) BIG3 inhibits the estrogen-dependent nuclear translocation of PHB2 via multiple karyopherin-alpha proteins in breast cancer cells. PLoS One, 10:e0127707.

11.

Langfelder

, Horvath

(2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 9:559.

12.

Lee

, Park

, Sung

, et al. (2021) Tumor suppressor miR-584-5p inhibits migration and invasion in smoking related non-small cell lung cancer cells by targeting YKT6. Cancers (Basel), 13:1159.

13.

, Cao

, Yu

, et al. (2021) A novel autophagy-related genes prognostic risk model and validation of autophagy-related oncogene VPS35 in breast cancer. Cancer Cell Int, 21:265.

14.

Love

, Huber

, Anders

(2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15:550.

15.

Min

, Lee

, Woo

, et al. (2021) Relation between tumor size and lymph node metastasis according to subtypes of breast cancer. J Breast Cancer, 24:75-84.

16.

Ooe

, Kato

, Noguchi

(2007) Possible involvement of CCT5, RGS3, and YKT6 genes up-regulated in p53-mutated tumors in resistance to docetaxel in human breast cancers. Breast Cancer Res Treat, 101:305-315.

17.

Papa

, Hahn

, Marsh

, et al. (2014) SOD2 to SOD1 switch in breast cancer. J Biol Chem, 289:5412-5416.

18.

Pellacani

, Tan

, Lefort

, Eaves

(2019) Transcriptional regulation of normal human mammary cell heterogeneity and its perturbation in breast cancer. EMBO J, 38:e100330.

19.

Perou

, Sorlie

, Eisen

, et al. (2000) Molecular portraits of human breast tumours. Nature, 406:747-52.

20.

Ryu

, Lee

, Yoon

, et al. (2016) Breast cancer-specific mortality in small-sized tumor with node-positive breast cancer: a nation-wide study in Korean breast cancer society. Breast Cancer Res Treat, 159:489-498.

21.

Schoppmann

, Bayer

, Aumayr

, et al. (2004) Prognostic value of lymphangiogenesis and lymphovascular invasion in invasive breast cancer. Ann Surg, 240:306-312.

22.

Siegel

, Miller

, Fuchs

, et al. (2021) Cancer statistics, 2021. CA Cancer J Clin, 71:7-33.

23.

Singletary

, Allred

, Ashley

, et al. (2002) Revision of the American Joint Committee on Cancer staging system for breast cancer. J Clin Oncol, 20:3628-3636.

24.

Sopik

, Narod

(2018) The relationship between tumour size, nodal status and distant metastases: on the origins of breast cancer. Breast Cancer Res Treat, 170:647-656.

25.

Sun

, Wang

, Dong

, et al. (2020) LncRNA PVT1 promotes exosome secretion through YKT6, RAB7, and VAMP3 in pancreatic cancer. Aging (Albany NY), 12:10427-10440.

26.

Sung

, Ferlay

, Siegel

, et al. (2021) Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin, 71:209-249.

27.

Waks

, Winer

(2019) Breast cancer treatment. JAMA, 321:316.

28.

Wang

, Zhou

, Liu

, et al. (2014) Tumor size of breast invasive ductal cancer measured with contrast-enhanced ultrasound predicts regional lymph node metastasis and N stage. Int J Clin Exp Pathol, 7:6985-6991.

29.

Wilcock

, Webster

(2021) The breast cancer drug market. Nat Rev Drug Discov, 20:339-340.

30.

, Chen

, Neville

, et al. (2011) Effect of very small tumor size on cancer-specific mortality in node-positive breast cancer. J Clin Oncol, 29:2619-2627.

31.

, Hu

, Xu

, et al. (2021) clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb), 2:100141.

32.

, Jiang

, Xu

, et al. (2021) High expression of YKT6 associated with progression and poor prognosis of hepatocellular carcinoma. Scand J Gastroenterol, 56:1349-1354.

33.

Yaman

, Gumuskaya

, Ozkan

, et al. (2012) Lymphatic and capillary invasion patterns in triple negative breast cancer. Am Surg, 78:1238-1242.

34.

Yang

, Zhang

, et al. (2020) KRT6A Promotes EMT and cancer stem cell transformation in lung adenocarcinoma. Technol Cancer Res Treat, 19:1533033820921248.

35.

Yang

, Yan

, Zheng

, et al. (2021) YKT6, as a potential predictor of prognosis and immunotherapy response for oral squamous cell carcinoma, is related to cell invasion, metastasis, and CD8+ T cell infiltration. Oncoimmunology, 10:1938890.

36.

, Jiang

, Chen

, et al. (2012) Effect of large tumor size on cancer-specific mortality in node-negative breast cancer. Mayo Clin Proc, 87:1171-1180.

37.

, Jiang

, Hu

, et al. (2021) LINC00839/miR-519d-3p/JMJD6 axis modulated cell viability, apoptosis, migration and invasiveness of lung cancer cells. Folia Histochem Cytobiol, 59:271-281.

38.

Zhang

, Wei

, Li

, el al. (2022) LINC00839 promotes neuroblastoma progression by sponging miR-454-3p to up-regulate NEUROD1. Neurochem Res, 47:2278-2293.

39.

Zhou

, Chang

, Zhu

, et al. (2021) LINC00839/miR-144-3p/WTAP (WT1 Associated protein) axis is involved in regulating hepatocellular carcinoma progression. Bioengineered, 12:10849-10861.

40.

Zhu

, Tu

, Chang

, et al. (2021) The orphan nuclear receptor gene NR0B2 is a favorite prognosis factor modulated by multiple cellular signal pathways in human liver cancers. Front Oncol, 11:691199.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.01 MB