{'fraction_non_alphanumeric': 0.04231320306069567, 'fraction_numerical': 0.028442520501173674, 'mean_word_length': 5.050165990409443, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 3, 'https://': 7, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A gene expression signature in HER2+ breast cancer patients related to neoadjuvant chemotherapy resistance, overall survival, and disease-free survival

Breast cancer ranks first in terms of mortality and incidence rates worldwide among women. The HER2+ molecular subtype is one of the most aggressive subtypes; its treatment includes neoadjuvant chemotherapy and the use of a HER2 antibody. Some patients develop resistance despite positive results obtained using this therapeutic strategy. Objective. To identify prognostic markers for treatment and survival in HER2+ patients. Methods. Patients treated with neoadjuvant chemotherapy were assigned to sensitive and resistant groups based on their treatment response. Differentially expressed genes (DEGs) were identified using RNA-seq analysis. KEGG pathway, gene ontology, and interactome analyses were performed for all DEGs. An enrichment analysis Gene set enrichment analysis was performed. All DEGs were analyzed for overall (OS) and disease-free survival (DFS). Results. A total of 94 DEGs were related to treatment resistance. Survival analysis showed that 12 genes (ATF6B, DHRS13, DIRAS1, ERAL1, GRIN2B, L1CAM, IRX3, PRTFDC1, PBX2, S100B, SLC9A3R2, and TNXB) were good predictors of disease-free survival, and eight genes (GNG4, IL22RA2, MICA, S100B, SERPINF2, HLA-A, DIRAS1, and TNXB) were good predictors of overall survival (OS). Conclusion: We highlighted a molecular expression signature that can differentiate the treatment response, overall survival, and DFS of patients with HER2+ breast cancer.


Introduction
Breast cancer is a heterogeneous disease characterized by abnormal and uncontrolled growth of malignant breast cells. Among all types of cancer, this disease ranks first in mortality and incidence rates in women over 25 years of age worldwide (Sung et al., 2021). In 2000, Perou et al. reported different molecular expression patterns in patients with breast cancer, and these patterns were subsequently used to classify breast cancer into distinct molecular subtypes (Perou et al., 2000;Sorlie et al., 2003). According to this classification, cancer cells that express human epidermal growth factor 2 (ERBB, formerly HER2) and not estrogen receptors (ER) are identified as the HER2+ molecular subtype, which represents 15%-30% of breast cancer patients, is an aggressive phenotype, and a predictor of poor outcome (Ban et al., 2020).
The treatment of HER2+ breast cancer includes the administration of chemotherapy and trastuzumab, a monoclonal antibody against the HER2 receptor (Abal et al., 2003;Harbeck and Gnant, 2017;Waks and Winer, 2019). Conventional neoadjuvant chemotherapy involves anthracyclines followed by taxane application. Anthracyclines work by joining DNA and suppressing the binding of DNA polymerase, thereby preventing DNA replication (McGowan et al., 2017). Taxanes affect mitotic spindle formation by binding to tubulin dimers, thereby preventing the division of tumor cells (Yardley, 2013;Harbeck and Gnant, 2017). Furthermore, adding trastuzumab in conventional chemotherapy helps block HER2 receptor-induced cell growth signaling (Maximiano et al., 2016). Despite the positive results obtained with this therapeutic strategy, some patients develop resistance. The molecular mechanisms underlying resistance are not fully understood; therefore, there is a lack of predictive biomarkers that are helpful in the prognosis and prediction of chemotherapy response (Iwamoto et al., 2020).
This study aimed to evaluate the transcriptome of HER2+ breast cancer patients and, according to their response to chemotherapy (sensitivity or resistance), to identify differentially expressed genes (DEGs) that could be useful in predicting patient outcomes after neoadjuvant chemotherapy treatment.

Materials and methods
Sample selection, chemotherapy treatment, and study design Patients aged 18 years and older with a diagnosis of breast cancer, HER2+/PR-/ER-, tumor size >2 cm, and positive nodes, candidates to receive neoadjuvant chemotherapy, and without previous therapy against cancer were recruited for this study. Patients with metastatic cancer, those with insufficient breast cancer biopsy tissue for pathological analysis, or those with RNA extraction were excluded. All participants provided written informed consent prior to enrolment. The data were deposited in the Gene Expression Omnibus (GEO) repository under the number GSE162187. Samples were separated by pathologic response into two groups: pathological complete response (pCR) was considered the sensitive group, and those in the non-pCR group were considered the resistant group.
Additionally, we used and analyzed data from the GSE163882 study, and HER2+/PR-/ER-samples were selected. The results obtained from both databases were compared.
Finally, from the TCGA breast ductal carcinoma database, HER+/PR-/ER-breast cancer samples were selected for analysis of overall survival.

Ethics and informed consent statements
The study was conducted in accordance with the guidelines of the Declaration of Helsinki and the ethical standards of the institutional and/or national research committee. This study was approved by the Ethical and Research Committee of the Instituto Mexicano del Seguro Social (IMSS) (number R-2013-785-061). Informed consent was obtained from all subjects involved in the study.

Quality control, alignment, and differential expression
The FASTQ files were analyzed with the Flexbar software tool version 3.5.0 (https://github.com/seqan/flexbar/releases/tag/v3.5. 0) (Dodt et al., 2012;Roehr et al., 2017) to remove Illumina adapters and to filter reads by a Phred score >30. To quantify the RNA-seq data, a pseudo-alignment was performed using Kallisto software version 0.46.1 (https://pachterlab.github.io/kallisto/ download.html) (Bray et al., 2016) with the default parameters and the GRCh38 human genome reference (GRCh38. p12). The DESeq2 package version 1.28.1 (https://bioconductor.riken.jp/ packages/3.0/bioc/html/DESeq2.html) (Love et al., 2014) was used for the analysis of abundance tables and the identification of differentially expressed genes (DEGs) for comparing resistant and sensitive samples (set as the reference group). The Ensembl database was used for the annotation of genes. To decrease the false discovery rate, the Benjamini-Hochberg correction test was applied to obtain adjusted p-values.
The Panther database v.16.0 (Mi et al., 2021) web tool was used for Gene Ontology enrichment analysis using Fisher's exact test and false discovery rate (FDR), with a threshold of p < 0.05, which was considered to be significant for each of the three categories, that is, molecular function, cellular component, and biological process.
To perform an interactome analysis, the DEGs were filtered by adjusted p-value <0.05 and analyzed using STRING-DB v11.0 software (https://string-db.org/) (Szklarczyk et al., 2019), the confidence score was set up at 0.7 to represent protein-protein associations.

Principal component analysis and heatmap representation
The geometric mean of the counts for each gene was used as a normalization factor. Once normalized, the principal component analysis (PCA) and heatmap representation were performed with the prcomp package and heatmap functions, respectively, with the default parameters in R v.4.0.2 ("Taking off Again") using as variables the normalized counts of the DEGs with an adjusted p-value of <0.05.

Survival analysis
Furthermore, a database of 109 patients obtained from the TCGA breast ductal carcinoma study with HER+/ER-/PR-was analyzed (TCGA Research Network: https://www.cancer.gov/ tcga) at 60 and 120 months to analyze overall survival (OS) and disease-free survival (DFS). Gene expression levels were determined according to normalized Log2-read counts for each gene. Median and quartiles were used for determining high-and low-expression groups. DEGs with an adjusted p-value < 0.05 were analyzed. Statistical significance was set at p < 0.05.

Results
A previous study was conducted to determine biomarkers of response to neoadjuvant chemotherapy in patients with breast cancer (GSE162187) (Barron-Gallardo et al., 2022). In this study, HER2+ samples were taken and included in the RNA-seq analysis; five samples were from patients categorized as resistant to treatment, and three samples were from patients sensitive to treatment. This small subset was used as training data. The results were validated using the GSE163882 dataset, which included information from 222 patients with breast cancer. Patients were over 33 years old; the mean ages for the resistant and sensitive groups were 52.2 (±12.15) and 62 (±6.24) years, respectively, with no statistical differences. The diagnostic status of all the patients was invasive ductal carcinoma breast cancer. The histological grades for tumor biopsies based on SBR (Scarff-Bloom-Richardson) parameters were five SBRII, two SBRIII, and one with non-available information (Supplementary Table S1).
The transcriptomic pattern was studied to determine the variables that specifically discriminated HER2+ patients according to their neoadjuvant treatment response. Despite having two groups defined by their pathological response, principal component analysis (PCA) with all genes detected by RNA-seq showed that the samples did not form specific clusters. Moreover, the distribution of the samples followed a heterogeneous pattern, indicating that HER2+ breast cancer patients may have high variability in gene expression (Supplementary Figure S1A).
Determination of DEGs related to treatment resistance in the training data HER2+ patients were grouped into sensitive and resistant to neoadjuvant chemotherapy groups to obtain DEGs related to treatment response. The transcriptomes of both groups were compared, and a total of 1383 genes were observed to be differentially expressed (p < 0.05), of which 719 were subexpressed and 664 were overexpressed in the resistant group as compared with the sensitive group ( Figure 1A). To diminish the inclusion of false-positive DEGs, the Benjamini-Hochberg post hoc test was applied; among the 1383 DEGs, only 94 maintained statistical significance with an adjusted p-value (p-adj) <0.05 (45 subexpressed and 49 overexpressed genes) ( Figure 1B).
Thereafter, we investigated whether this set of 94 DEGs could adequately classify the samples as sensitive and resistant to treatment. Therefore, principal component analysis (PCA) with only 94 genes was performed again. The results showed two clusters defined by principal component 1 (PC1) with 55.35% and principal component 2 (PC2) with 15.79% of the data variance (Supplementary Figure S1B). Therefore, the selection of the 94 DEGs included genes capable of clustering HER2+ patients into resistant and sensitive groups. As shown in Supplementary Table S2,  To visualize the expression pattern (color density) and distribution (clustering) of the 94 DEGs, heatmap analysis was performed ( Figure 2). The column dendrogram results showed two clusters belonging to the sensitive and resistant groups. The row dendrogram shows four clusters of genes with similar expression patterns.

Pathways and enrichment analysis of DEGs related to chemotherapy resistance
The 94 DEGs were analyzed using the KEGG Mapper search pathway tool and DAVID v6.8. Among the 318 KEGG pathways, seven were statistically modulated (FDR <0.05) ( Figure 3A), including graft-versus-host disease, allograft rejection, type I diabetes mellitus, autoimmune thyroid disease, viral myocarditis, antigen processing and presentation, and cell adhesion molecules. In addition, GO analysis results showed that the biological processes enriched by DEGs were related to the interferon-gamma-mediated signaling pathway. The cellular components in which the DEGs were included were associated with the MHC class II protein complex, luminal side of the endoplasmic reticulum membrane, endoplasmic reticulum (ER)to-Golgi transport vesicle membrane, and extracellular space. Finally, the modulated molecular functions were MHC class II receptor activity and peptide-antigen binding ( Figure 3B).
Enrichment analysis showed a total of 40 gene sets enriched with a p-value < 0.05 (35 positively and five negatively), which belongs to C1 (1 enriched set), C2 (6 enriched sets), C3 (2 enriched sets), C5 (19 enriched sets), C7 (7 enriched sets), and C8 collections (5 enriched sets). The C4, C6, and H collections did not contain enriched sets. From the enriched sets, we found two related to therapy resistance (Massarweh tamoxifen resistance and Creighton endocrine therapy resistance gene sets) and three related to the immune system, such as GOBP immune response, Goldrath antigen response, and CHR6P21, which is a location for genes related to the immune system (HLA-DQA1, HLA-DRB1, HLA-B, and MICA) (Supplementary Table S2).

Determination of interactions clusters between DEGs
An analysis of 94 DEGs was performed to determine the molecular interactions between them. The STRING-DB tool was used to set an interaction score with high confidence (0.7). The results showed 18 edges (genes) distributed among seven clusters: one with five genes, one with three genes, and five with two genes. Among the seven clusters, the cluster with five edges was related to the interferon-gamma-mediated signaling pathway, MHC class I/II-like antigen recognition protein, and cell adhesion molecules, including the HLA-A, HLA-DQA1, TRIM26, HLA-B, and HLA-DRB1 genes ( Figure 3C).

Evaluation of DEGs for survival prediction
To evaluate whether the expression of the 94 DEGs was related to survival prediction, measured as DFS or OS, we analyzed the 94 DEGs individually using a database with Frontiers in Genetics frontiersin.org

FIGURE 2
Expression patterns of DEGs in resistance and sensitivity. The row Z-score of the normalized read counts of DEGs with p-adj < 0.05 are plotted in the heatmap. The red color indicates a row Z-score >0, and the blue indicates a row Z-score <0. Columns represent each patient, and each row represents a gene. The dendrogram at the top of the heatmap clusters the patients according to their gene expression pattern, while the dendrogram at the left side of the heatmap groups the genes with similar expression patterns. Columns 1 to 3 represent sensitive patients, and columns 4 to 8 represent resistant patients.
Frontiers in Genetics frontiersin.org 05 expression information of 109 patients with HER+/ER-/PR-, and data were obtained from TCGA breast ductal carcinoma study. A total of 12 DEGs predicted the DFS. The high expression of ATF6B, DHRS13, DIRAS1, ERAL1, GRING2B, IRX3, PRTFDC1, and PBX2 was found to be an excellent prognostic of DFS at 5 years; on the other hand, a high expression of L1CAM was associated with lower DFS at 5 years ( Figure 4). We found that low expression of TNXB and SLC9A3R2 and high expression of S100B were associated with better DFS in the long term (10 years) ( Figure 5).
According to the OS analysis, groups with high expression of GNG4, IL22RA2, S100B, and SERPINF2 were associated with better OS at 5 years; the same was true for HLA-A and DIRAS1 at 10 years. In contrast, high expression of MICA and TNXB was related to lower OS times at 5 and 10 years, respectively ( Figure 6).

FIGURE 3
KEGG, GO, DEGs interaction related to chemotherapy resistance. The 94 DEGs with a p-adj < 0.05 were analyzed to know their contribution to KEGG pathways, gene ontology, and the interaction clusters (A) KEGG enrichment analysis (B) GO enrichment analysis. Each bar represents the fold enrichment value for KEGG and GO. The x-axis plots the fold enrichment values, and the y-axis shows the pathway's name or GO terms. In GO enrichment analysis, the plot is divided into three categories: biological process (red bars), cellular component (green bars), and molecular function (blue bars). (C) Interactome analysis. Only DEGs that interact with each other were plotted in the graph. Network nodes represent proteins encoded by DEGs; colors represent the category to which encoded proteins belong; edges represent protein-protein interactions. Line colors indicate the type of interaction reported.
Frontiers in Genetics frontiersin.org

FIGURE 5
DEGs related to the prediction of DFS at 10 years. DEGs that meet the criteria of p-adj<0.05 were analyzed using the TCGA ductal breast cancer database (https://www.cancer.gov/tcga) to determine their association with DFS at 10 years. The red line indicates the high expression group, the blue line represents the low expression group. Y-axis shows the disease-free survival percentage; X-axis shows the time in years. (A) S100B, (B) SLC9A3R2, (C) TNXB.
Frontiers in Genetics frontiersin.org 07 Odds ratio analysis was performed to determine if there is any difference at the end point of five of 10 years in OS or DFS. The results of the odds ratio analysis showed a similar prognostic pattern for each gene compared with the results obtained from the log-rank analysis, except for GRIN2B, PRTFDC1, SLC9A3R2 in DFS and HLA-A, IL22RA2 in OS whose p values were greater than 0.05 (Supplementary Table S3).
Furthermore, univariate and multivariate Cox analyses were performed. The results show some genes in which the expression can be considered a predictor variable associated with survival time. In univariate cox analysis for DFS, the coefficients were negative for DHRS13, GRIN2B and positive for L1CAM with p < 0.05. When applying the univariate cox analysis for OS, DIRAS had a negative coefficient, and MICA had a positive coefficient p < 0.05 (Table 1). We performed multivariate cox analysis using age, pathologic stage, radiation therapy, and the expression level as variables. The results show that as higher the pathologic stage, the hazard to disease recurrence increases for ATF6B, DHRS13, DIRAS1, ERAL1, GRIN2B, L1CAM, PBX2, PRTFDC1, SLC9A3R2, and TNXB, furthermore, increase the risk of death when analyzed DIRAS1, IL22RA, MICA, and S100B. Alongside, neoadjuvant radiation therapy was correlated with decrease recurrence risk in ATF6B, DHRS13, DIRAS1, LICAM, and SLC9A3R2, and decreases dead risk when analyzed S100B, MICA, and DIRAS1 (Table 2).

FIGURE 6
DEGs related to overall survival. DEGs with p-adj<0.05 were contrasted with overall survival data at 5 and 10 years in the TCGA ductal breast cancer database (https://www.cancer.gov/tcga). The red line indicates the high expression group, and the blue line represents the low expression group. Y-axis shows the overall survival percentage; X-axis shows the time in years. Panels for (A) GNG4, (B) IL22RA2, (C) MICA, (D) S100B, and (E) SERPINF2 represent OS at 5 years, while panels for (F) HLA-A, (G) DIRAS1, and (H) TNXB showed OS data for 10 years.
Frontiers in Genetics frontiersin.org

Analysis of the 94 DEGs in other studies highlight similar DEGs as possible biomarkers
Finally, to evaluate whether the data from GSE162187 (94 DEGs between resistant and sensitive patients) has a consistent expression with other studies, we analyzed the data of the study GSE163882, which aimed to predict pCR to neoadjuvant therapy in breast cancer patients. Data from GSE162187 were used as training data, and data from the GSE163882 dataset were used as corroboration data. This analysis discard 84 DEGs and identified 10 DEGs in common (ATF6B, ERAL1, CRYM, MUC16, SOX10, MICA, PDE2A, TMEM97, SDF2, and BICDL2) that could discriminate patient  outcomes. Therefore, these 10 DEGs were considered possible biomarkers of pCR and neoadjuvant chemotherapy response. Moreover, three DEGs (ATF6B, ERAL1, and MICA) have a strong correlation with DFS and OS.

Discussion
One fundamental aspect of treating breast cancer patients is the knowledge of their molecular subtypes. This information has per se a prognostic value for predicting patient treatment response (von Minckwitz et al., 2012), which can be evaluated according to the criteria for the diagnosis of pCR. Achieving pCR has been associated with better overall survival (Broglio et al., 2016;Spring et al., 2020); however, the intrinsic factors involved in pCR have not been clarified. There is still controversy on whether standard adjuvant therapy increases pCR (Mauri et al., 2005). The percentage of patients who achieve pCR ranges between 27%-47% (Muller et al., 2021;Xin et al., 2021). In this study, 37.5% of patients achieved pCR. Therefore, HER2+ breast cancer patients were categorized as sensitive (pCR achieved) or resistant (pCR did not achieve) to neoadjuvant chemotherapy and were used as an RNA-seq strategy to identify predictors of pCR. It should be considered that the diagnosis of HER2+ breast cancer was because more than 10% of the tumor cells present detectable HER2 expression; therefore, there is a large percentage of cells that do not express HER2. This highlights the heterogeneity of this breast cancer subtype (Ng et al., 2015;Chen et al., 2020). When the transcriptional profiles of all patients were compared to determine clusters, this heterogeneity was emphasized (Supplementary Figure S1A). The implications of molecular differences in HER2+ breast cancer patients are not fully understood and may be relevant to prognosis and treatment response.
The DEGs found in HER2+ breast cancer patients sensitive and resistant to neoadjuvant chemotherapy were mainly related to plasma membranes, vesicles, and extracellular space and were involved in different biological processes, such as cellular response to chemical stimulus, cell adhesion, and signal regulation. Variations in the protein components of the extracellular matrix have been reported in breast tumors of different origins (Borghesi et al., 2021). In addition, extracellular components such as the extracellular matrix, vesicles, and plasma membranes can be modified by cancerassociated fibroblasts, leading to a tumor microenvironment involved in cancer development and drug resistance (Mashouri et al., 2019;Helal-Neto et al., 2020;Lugo-Cintron et al., 2020).
One of the most enriched pathways is related to cell adhesion molecules involved in tight junctions of epithelial and endothelial cells, such as claudins, which participate in epithelialmesenchymal transition (EMT) and chemoresistance (Hewitt et al., 2006;Agarwal et al., 2009;Gowrikumar et al., 2019).
According to KEGG enrichment analysis, the AMPK signaling pathway is involved in the resistance process; this pathway is considered a double-edged sword that protects and promotes cancer progression (Jeon and Hay, 2015). Sensitization of breast cancer cells to chemotherapy by activating AMPK signaling by CTAB has been observed (Pan et al., 2019). Similarly, histological evaluations have reported altered AMPK signaling in breast cancer samples (Hadad et al., 2009), and this pathway is considered a therapeutic target for breast cancer treatment (Hadad et al., 2008). However, it has been hypothesized that once cancer has developed, AMPK promotes the survival of cancer cells by protecting them against DNA damage, nutritional stress, and hypoxia (Russell and Hardie, 2020). Further studies are needed to delineate the role of the AMPK pathway in breast cancer and the development of chemotherapy resistance.
Another enriched pathway was cyclic guanosine 3,5monophosphate (cGMP) and protein kinase G (PKG). The cGMP-PKG pathway has been associated with the modulation of apoptosis and growth inhibition in MCF-7 and MDA-MB-468 breast cancer cell lines (Fallahian et al., 2011). An essential component of this pathway is the protein kinase cGMPdependent 2 (PRKG2), which was found to be downregulated in the resistant group in this study. Our results correlate with those of Karami-Tehrani et al. (Karami-Tehrani et al., 2012), who observed lower expression of PRKG2 protein in breast tumor samples. In addition, it has been reported that PRKG2 inhibits EGF-induced MAPK/c-Jun N-terminal kinase (JNK) signal transduction in human breast cancer cells (Lan et al., 2012) and also inhibits the activation of EGFR and HER2 in gastric cancer cells (Zhu et al., 2016;Lan et al., 2019). PRKG2 inhibits the migration, invasion, and proliferation of cancer cells and activates CREB, which modulates anti-apoptotic genes, such as BCL2 (Shankar et al., 2010), which are overexpressed in the resistant group, thereby contributing to the survival of cancer cells in the resistant group.
In this study, many DEGs related to resistance were identified. With the dimensional reduction, samples clustered better, highlighting the possibility of using these genes to predict the response to treatment. An interesting finding in our results was a group of DEGs that interacted with each other, including HLA-A, HLA-DQA1, HLA-DRB1, HLA-B, and TRIM26, which are components of the MHC protein complex, except TRIM26. These DEGs were found to be overexpressed in the resistant group. The upregulation of classical and non-classical HLA-I molecules has been reported to acquire a "protective" phenotype in melanoma cells (Balsamo et al., 2012). HLA molecules play a role in self-recognition by immune cells, which is essential for hematopoietic and healthy cells to avoid their destruction, and the loss, alteration, or absence of HLA molecules can cause susceptibility to NK cell attack (Ljunggren and Karre, 1990;Moretta et al., 2004). HLA molecules interact with inhibitory receptors such as killer cell immunoglobulin-like receptors (KIRs), leukocyte immunoglobulin-like receptors (LIRs), and Frontiers in Genetics frontiersin.org natural killer group 2A (NKG2A) on the NK surface, avoiding its activation (Khan et al., 2020). Overexpression of HLA by cancer cells has been reported as a mechanism for evading the immune response of NK cells and is termed immune checkpoint inhibition (Bi and Tian, 2019). From this group of genes, a variant of HLA-A (ENSEMBL ID ENSG00000235657) was observed subexpressed in the resistance group. In addition, low expression of this gene was associated with a worse prognosis for OS. This gene has already been reported to predict treatment response and OS (Sinn et al., 2021;Barron-Gallardo et al., 2022). ATF6B has two ensemble IDs (ENSG00000228628 and ENSG00000213676).
ENSG00000228628 ID was overexpressed. Nevertheless, ENSG00000213676, which corresponds to the primary assembly of this gene, was found to be sub-expressed, and low expression was related to lower OS and worse DFS. Variants of this gene have been associated with an increased risk of breast cancer development (Dierssen-Sotos et al., 2018).
DIRAS1 was found to be sub-expressed in the resistance group. Subexpression of this gene was correlated with lower OS. This gene has tumor-suppressive activity by binding to SmgGDS, which blocks the interactions of small GTPases, such as Rho and K-Ras4B. The expression of DIRAS1 is downregulated in most types of breast cancer (Bergom et al., 2016).
Other sub-expressed genes in resistant treatment and lower DFS were GRIN2B, GNG4, and IRX3. GRIN2B is involved in breast cancer progression and acts as a promoter of CpG islands (Park et al., 2011;Park et al., 2012). GNG4 is hypermethylated in breast cancer; however, when comparing all molecular subtypes, the HER2 subtype shows the highest expression levels for this gene (Fernandez-Nogueira et al., 2016;Mao et al., 2021). IRX3 plays an important role in obesity and type 2 diabetes; however, it plays an important role in the adaptability of tumor cells to metabolic challenges, a process that has a parallelism with the development of chemotherapeutic resistance (Singh et al., 2016).
A set of genes that showed high expression in the resistant group, which were related to lower OS and worse DFS, were L1CAM, MICA, and TNXB. The expression of L1CAM is increased in luminal B breast cancer, and its expression is related to disease recurrence and higher levels of Ki-67 expression (Moisini et al., 2021). A soluble form of L1CAM has been found in HER2-enriched primary breast cancer patients (Wu et al., 2018). There are reports that inhibition of L1CAM reverses cisplatin resistance in triple-negative breast cancer cells (Zhang et al., 2022). MICA is overexpressed in breast cancer when compared to normal tissue and is considered an indicator of poor prognosis (Madjd et al., 2007). It is an activation ligand of NK cells, which induces the lysis of cells that express it. However, there is a soluble form of MICA (sMICA) that decreases the expression and presentation of NKG2D, a natural cytotoxic receptor in natural killer cells, thus sMICA helps cancer cells to evade immune cell attack (Pan et al., 2017) and contributing to a worse prognosis in cancer (Roshani et al., 2016). In this study, high expression of MICA was observed in the resistant and lower OS groups; however, further studies are needed to determine the role of MICA or sMICA in chemotherapy resistance. In the case of TNXB, the expression of this gene has been analyzed in breast cancer, and a correlation between high TNXB expression and good survival prognosis has been found (Liot et al., 2020). Its expression decreases at late stages, major tumor grade, and node status of the disease (Liot et al., 2020), however, its expression in the HER2 molecular subtype and in relation to chemotherapy resistance has not been evaluated.
In contrast, genes with high expression but related to better OS and DFS were IL22RA2, PRTFDC1, PBX2, S100B, SERPINF2, DHRS13, ERAL1, and SLC9A3R2. IL22RA2 expression decreases in luminal A, B, and triplenegative breast cancers (Fu et al., 2015); however, but HER2+ breast cancer has not been reported. PRTFDC1 has been associated with the triple-negative basal-like immunesuppressed breast cancer subtype (TNBC-BLIS), which is considered one of the worst prognoses (Yin et al., 2020). The most highly expressed gene is PBX2. This gene was found to be upregulated in breast lesions and has been proposed along with other genes as a candidate biomarker for distinguishing breast cancer lesions (Hou et al., 2020). It has been showed that the overexpression of PBX2 increases the tumorigenic properties of SkBr3 breast cancer cell line when transfected with HoxB7 (Fernandez et al., 2008). S100B expression has been negatively correlated with lymph node metastasis (Wang et al., 2021), inhibition of cell migration, better overall survival in luminal B breast cancer patients, and being a good distant metastases-free survival biomarker (Yen et al., 2018). SERPINF2 is differentially expressed in breast cancer tissues compared with normal tissues (Malvia et al., 2019). The protein product of SERPINF2 has been found in the serum of breast cancer patients when evaluating treatment response; however, this protein appeared in both resistant and sensitive groups (Chantada-Vazquez et al., 2022).
Finally, DHRS13, ERAL1, and SLC9A3R2 could predict treatment response and survival; however, there are no reports related to breast cancer and its possible function in this disease.

Conclusion
This study underlines a molecular expression pattern related to the response of patients with HER2-positive breast cancer to neoadjuvant chemotherapy. Differentially expressed genes highlight the involvement of pathways, such as extracellular components, adhesion molecules, and immune responses, in the process of resistance to chemotherapy. Some differentially expressed genes can be Frontiers in Genetics frontiersin.org used as biomarkers of overall survival and disease-free survival in breast cancers.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

Ethics statement
The studies involving human participants were reviewed and approved by Comite Nacional de Investigación Científica del Instituto Mexicano del Seguro Social. The patients/participants provided their written informed consent to participate in this study.

Author contributions
CB-G and MG-C. contributed to sample recruitment and processing. CB-G performed the bioinformatic analysis, interpretation of data, and drafting of the manuscript. MG-C. contributed to data analysis and funding acquisition. RD-C, AM-M, MM-S, and MV-G were involved in patient recruitment, obtaining clinical information, and data interpretation. AA-L. and LJ-S. conceived the study, advised, analyzed the results, contributed to funding acquisition, and wrote and revised the manuscript. All authors substantively revised the manuscript, suggested modifications, and approved the final version.

Funding
This work was funded by the Fondo de Investigación en Salud, Instituto Mexicano del Seguro Social with grant number FIS/IMSS/PROT/PRIO/14/030 to LJ-S.