Correlation between S100A7 and immune characteristics, methylation, tumor stemness and tumor heterogeneity in pan-cancer and its role in chemotherapy resistance in breast cancer

Data collection

The Cancer Genome Atlas-breast cancer cohort (TCGA-BRCA) and pan-cancer RNA-seq data were collected from the University of California Santa Cruz database (https://xenabrowser.net/) and converted into the transcripts per million formats. Data on the clinical follow-up survival and clinicopathological characteristics were simultaneously downloaded. To reduce bias in the statistical analysis, BRCA and pan-cancer patients with missing overall survival data or a follow-up time <30 days and male BRCA patients were excluded. The tumor mutational burden (TMB) and mutant-allele tumor heterogeneity (MATH) scores were calculated.

Differential and survival analysis

To ensure the accuracy of the following analysis, we deleted duplicated and missing RNA-seq data, and converted the remaining data into the log₂(transcripts per million + 1) format. Differential S100A7 expression between normal and tumor tissues was identified using Wilcoxon rank-sum test, and samples were classified into two groups based on the median S100A7 expression.

The relationship between S100A7 expression and pan-cancer prognosis was examined using log-rank analysis. Cox analysis was performed to ensure the precision of the results, which were depicted as heatmaps and forest plots. The “Survminer” and “Survival” packages were used to perform survival analysis and plot survival curves.

Enrichment analysis

To clarify the impact of S100A7 expression on some potential pan-cancer pathways, Gene Set Enrichment Analysis (GSEA) was performed, and the results were displayed in the form of a heatmap.

We also explored the impact of S100A7 expression specifically in breast cancer. Gene Ontology analysis was performed to identify the molecular functions, biological processes, and cellular components that were enriched in breast cancer. Kyoto Encyclopedia of Genes and Genomes enrichment analysis and GSEA could predict the pathways related to S100A7 alteration-associated genes in breast cancer. The protein–protein interaction network was also elucidated to understand the interactions between proteins associated with S100A7 alterations in breast cancer. Finally, we compared the scores of the common pathways between the S100A7-high-expression and the S100A7-low-expression groups.

Immunomodulatory, immune cell infiltration, immune checkpoint, and immune activity score analysis

We assessed the correlations between S100A7 expression and immunomodulators, immune checkpoints, and immune cell infiltration. Lists of genes encoding immune checkpoints and immunomodulators were downloaded from TISIDB (http://cis.hku.hk/TISIDB/download.php). We also compared immunomodulatory scores between the S100A7-high-expression and the S100A7-low-expression groups to further define the relationship between S100A7 expression and immunomodulators.

The immune cell infiltration score was calculated by multiple methods (TIMER, EPIC, MCP-COUNTER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, XCELL, TIDE) and downloaded from TIMER2.0 (http://timer.cistrome.org/). Heatmaps were used to illustrate the relationship between immune cell infiltration and S100A7 expression pan-cancer. The ESTIMATE method was also used to perform the same analysis pan-cancer as well as in breast cancer.

Tracking Tumor Immunophenotype (http://biocc.hrbmu.edu.cn/TIP/) was used to calculate the immune activity score. We compared the immune activity scores between the S100A7-high and -low expression groups to establish the relationship between S100A7 expression and immune activity.

Finally, we performed a pan-cancer correlation analysis between S100A7 expression and immune checkpoints. We also compared immune checkpoint scores between the S100A7-high and -low expression groups.

Association with m⁶A, m⁵C, m¹A, and DNA methylation

A correlation analysis was performed between S100A7 expression and m⁶A, m⁵C, m¹A, and DNA methylation. DNA methylation data were downloaded from the Gene Set Cancer Analysis database (http://bioinfo.life.hust.edu.cn/GSCA/). The m⁶A, m⁵C, and m¹A gene sets were downloaded from the literature (https://doi.org/10.3389/fimmu.2022.918140).

Gene mutation, tumor stemness, and tumor heterogeneity analysis

TMB and MATH were used to predict the efficacy of tumor immunotherapy and tumor heterogeneity. The “maftools” package was used for calculating TMB and MATH scores in breast cancer [18]. We also downloaded microsatellite instability, neoantigen, and purity scores from the literature to examine the relationship between S100A7 expression and tumor heterogeneity [19, 20].

Six tumor stemness scores—RNA expression-based stemness score (RNAss), epigenetically regulated RNA expression-based stemness score (EREG.EXPss), DNA methylation-based stemness score (DNAss), epigenetically regulated DNA methylation-based stemness score (EREG-METHss), differentially methylated probes-based stemness score (DMPss), and enhancer elements/DNA methylation-based stemness score (ENHss)—were downloaded from the literature [21]. The relationship between S100A7 expression and tumor stemness was analyzed using Spearman’s method.

We downloaded the level 4 single nucleotide variation dataset for TCGA samples processed using Mutect2 [22] and the level 4 copy number variation dataset for all TCGA samples processed using GISTIC [23]. Samples with synonymous mutations were filtered, and log₂(× + 0.001) transformation was applied to each expression value. Finally, tumors with less than three samples in a single cancer were excluded. The data so processed were used to study the relationship among single nucleotide variation, copy number variation, and S100A7 expression.

We also evaluated the relationship between gene mutation and S100A7 expression. The “maftools” package was used to analyze mutation-related data and to visualize this relationship in breast cancer. All these results were visualized using the “ggplot2” and “ggpubr” packages.

Experimental material

Cell Counting Kit-8 (CCK8) was purchased from MedChemExpress (Monmouth Junction, NJ, USA). CellTiter 96^® Aqueous One Solution Reagent Kit for measuring the cell viability was obtained from Promega Co. (Madison, WI, USA). Lipofectamine 3000 was bought from Invitrogen Co. (Carlsbad, CA, USA). The S100A7 expression plasmid with the green fluorescent protein tag (pGensil-1-S100A7) was designed by Sangon Biotech Co. (Shanghai, China).

Cell culture

MCF7 and MDA-MB-231 cells were purchased from the Cell Resource Center, Institute of Basic Medicine, Chinese Academy of Medical Sciences and Peking Union Medical College. They were cultured in DMEM (containing 4.5 g of glucose; Gibco, Waltham, MA, USA) containing 10% fetal bovine serum (FBS; Gibco) and 1% penicillin/streptomycin (Gibco) at 37°C in a humidified incubator having 5% CO₂.

Cell transfection

MCF7 and MDA-MB-231 cells were seeded in 6-well plates and transfected at 90% confluence with the indicated plasmid using Lipofectamine 3000, according to the manufacturer’s instructions. Briefly, 2.5 μg of the plasmid was mixed with 5 μL of P3000 in 125 μL of serum-free DMEM. Separately, 5 μL of Lipofectamine 3000 was mixed in 125 μL of serum-free DMEM. These two mixtures were complexed for 15 min at room temperature before they were added to the cells. The transfection efficiency was evaluated 48 h later under an inverted fluorescence microscope.

Reverse transcription-quantitative PCR

Total RNA was extracted from MCF7 and MDA-MB-231 cells using TRIzol reagent (Takara, Dalian, China), and cDNA was synthesized according to the instructions of the Reverse Transcription Kit (Takara). Based on the Ct value, the expression of the target gene was calculated by the 2^−ΔΔCt method.

CCK8 assay

MCF7 or MDA-MB-231 cells in their growth phase were seeded in 96-well cell culture plates at a density of 5000 cells/100 μL of medium and cultured for 24 h at 37°C. The drug concentrations tested on MDA-MB-231 cells were 0, 1, 2, 4, 5, and 8 mg/L, whereas those tested on MCF7 cells were 0, 0.1, 0.2, 0.25, 0.4, and 0.5 mg/L. After 48 h of culture, the cells were incubated with 10 μL of CCK8 reagent for 2 h, following which the absorbance was measured at 450 nm. The inhibition rate was calculated according to the following formula: inhibition rate = (1-experimental group/control group). Finally, the half-maximal inhibitory concentration (IC₅₀) was calculated.

We also used the CCK8 assay to evaluate cell proliferation. MCF7 and MDA-MB-231 cells were seeded in 96-well cell culture plates as described above. The absorbance was measured at 450 nm at 0 h, 24 h, 48 h, and 72 h to evaluate the susceptibility of S100A7 expression to breast cancer cell proliferation and chemosensitivity. We calculated the percentage of cell viability. Cell viability was calculated as follows:

$$Cellviability\left(\%\right)=\frac{OD\text{450(}treatmentgroup\text{)}-OD\text{450(}blank\text{)}}{OD\text{450(}controlgroup\text{)}-OD\text{450(}blank\text{)}}\times \mathrm{100\%}$$

Treatment group: Containing cells and different concentrations of drugs (include 0) in knockdown or control groups.

Control group: Only containing cells in knockdown or control groups.

Blank: Both cells and drugs were not contained in knockdown or control groups.

Transwell assay

To test the migration ability of breast cancer cells, suspensions of MCF7 and MDA-MB-231 cells were prepared in DMEM without FBS at a density of 2 × 10⁴ cells/mL. The upper chamber of the Transwell was filled with 100 μL of this suspension, while the lower chamber was filled with 500 μL of DMEM containing FBS. The excess cells were retrieved with a sterile cotton swab, fixed in 4% formaldehyde, and stained with 0.1% crystal violet to count the number of migrating cells.

To investigate the invasiveness of MCF7 and MDA-MB-231 cells, the stroma was diluted in DMEM without FBS at a ratio of 1:7 (Beyotime, Guangzhou, China). An aliquot of 50 μL was aspirated into the upper chamber and allowed to rest for 3–4 h. The rest of the procedure was the same as that described for the migration assay.

Wound healing assay

For this experiment, 1 × 10⁶ cells per well were seeded in 6-well plates. After cell transfection and overnight incubation, the cells were scratched using a 10-μL pipette tip and photographed using a powered microscope with a 4× magnification at 0 h and 24 h.

Statistical analysis

T-test was used to compare two groups of measurement data. The Wilcoxon rank-sum and Kruskal tests were employed to identify the differences in three clusters. Univariate Cox regression was used to assess the differences between counting data. Log-rank test was used to evaluate the prognosis of cancer. P < 0.05 was considered statistically significant.

Availability of data and materials

The data can be downloaded from TCGA (https://portal.gdc.cancer.gov/). The codes used during the current study are available from the corresponding author on reasonable request.

S100A7 was highly expressed and associated with poor prognosis in pan-cancer

We found that S100A7 was expressed at a significantly higher level pan-cancer (in glioblastoma multiforme, glioma, brain lower grade glioma, uterine corpus endometrial carcinoma, breast invasive carcinoma, esophageal carcinoma, stomach and esophageal carcinoma, colon adenocarcinoma, prostate adenocarcinoma, stomach adenocarcinoma, lung squamous cell carcinoma, thyroid carcinoma, rectum adenocarcinoma, ovarian serous cystadenocarcinoma, testicular germ cell tumors, all, kidney chromophobe; Figure 1A). S100A7 was highly expressed in breast cancer as well (Figure 1B–1D).

Figure 1. The pan-cancer expression of S100A7. (A) S100A7 expression from TCGA + GTEx microarray data. (B) S100A7 expression from GSE15852 data. (C) S100A7 expression in breast cancer and normal breast tissues from GSE10797 data. (D) S100A7 expression in breast cancer and normal breast tissues from TCGA-BRCA data. Abbreviations: TCGA: The Cancer Genome Atlas; GTEx: Genotype-Tissue Expression.

We also explored the relationship between S100A7 expression and pan-cancer prognosis. Log-rank test and univariate Cox regression were used to conduct survival analysis. We found that S100A7 was related to the prognosis of multiple types of cancer (Figure 2A). S100A7 expression was closely related to the overall survival of patients with skin cutaneous melanoma (Figure 2B), kidney renal clear cell carcinoma (Figure 2C), liver hepatocellular carcinoma (Figure 2D), and the disease-specific survival of patients with kidney renal clear cell carcinoma (Figure 2E) and skin cutaneous melanoma (Figure 2F). What’s more, S100A7 expression was closely related to the overall survival of patients with bladder urothelial carcinoma (Figure 2G).

Figure 2. Relationship between S100A7 expression and the prognostic value of S100A7 pan-cancer. (A) The pan-cancer prognostic value of S100A7 (univariate Cox and log-rank analysis). Prognostic value of S100A7 expression in terms of OS in (B) SKCM, (C) KIRC, (D) LIHC, and DSS in (E) SKCM and (F) KIRC. It also reflected the prognostic value of S100A7 expression in terms of OS in (G) BLCA. Abbreviations: OS: overall survival; DSS: disease-specific survival; SKCM: skin cutaneous melanoma; KIRC: kidney renal clear cell carcinoma; LIHC: liver hepatocellular carcinoma; BLCA: bladder urothelial carcinoma.

Enrichment analysis revealed that S100A7 was related to multiple cancer development pathways

We performed enrichment analyses to investigate the potential pathways affected by S100A7 pan-cancer and in breast cancer. GSEA showed that S100A7 was related to multiple cancer development pathways, such as Notch signaling (Figure 3A). The same analysis was performed in breast cancer, and all potential pathways as well as those with the most significant changes were displayed (Figure 3B, 3C). We also compared the differences in potential pathways between the S100A7-high and -low expression groups. Overall, we discovered that S100A7 was closely related to various biological activities, such as the cell cycle, metabolism of various substances, and cancer development (Figure 3D).

Figure 3. GSEA of DEGs between S100A7-high and -low expression groups pan-cancer. (A) Pan-cancer GSEA. (B) The GSEA results in breast cancer were shown using a word cloud map. (C) Top three potential pathways according to GSEA results in breast cancer. (D) Pathway score in the S100A7-high and -low expression groups in breast cancer. Abbreviations: GSEA: Gene Set Enrichment Analysis; DEGs: differentially expressed genes.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses showed that the upregulated genes were closely related to system development, the extracellular region, regulation of molecular function, and the interleukin-17 pathway (Figure 4A, 4B), whereas the downregulated genes were closely related to RNA processing, the nucleolus, signaling receptor activity, and the nicotine addiction pathway (Figure 4C, 4D). The protein–protein interaction network of all differentially expressed genes was charted to understand the interaction among them (Figure 4E).

Figure 4. GO and KEGG pathway enrichment analyses of upregulated and downregulated genes in the two groups based on the median S100A7 expression in breast cancer. (A) GO functional enrichment analysis of upregulated genes based on the S100A7 median expression grouping. (B) KEGG pathway enrichment analysis of upregulated genes based on the S100A7 median expression grouping. (C) GO functional enrichment analysis of downregulated genes based on the S100A7 median expression grouping. (D) KEGG pathway enrichment analysis of downregulated genes based on the S100A7 median expression grouping. (E) PPI network of the DEGs between S100A7-high and -low expression groups pan-cancer. Abbreviations: GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein–protein interaction.

S100A7 affects immunomodulators, immune cell infiltration, immune activity score, and immune checkpoints

We examined the relationship between S100A7 and tumor immunity. Spearman’s correlation analysis showed that S100A7 was closely related to immunomodulators pan-cancer (Figure 5A). In breast cancer, S100A7 was associated with immunomodulators (Figure 5B), the major histocompatibility complex (Figure 5C), chemokines (Figure 5D), and receptors (Figure 5E).

Figure 5. Correlation between S100A7 expression and immunomodulators in pan-cancer. (A) Correlation between S100A7 expression and immunomodulators in pan-cancer. The expression of genes related to immunomodulators (B), MHC (C), chemokines (D), and receptors (E) in the two groups based on the median S100A7 expression in breast cancer. “ns” represents not significant. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001. Abbreviation: MHC: major histocompatibility complex.

Multiple methods were used to assess the correlation between immune cell infiltration and S100A7 expression pan-cancer. S100A7 was associated with B cells, CD4⁺ T cells, CD8⁺ T cells, cancer-associated fibroblasts, macrophages, neutrophils, and mast cells (Figure 6). The immune, stromal, and ESTIMATE scores were also related to S100A7 expression pan-cancer (Figure 7A) as well as in breast cancer (Figure 7B–7D).

Figure 6. Correlation between S100A7 expression and immune cell infiltration pan-cancer. “×” represents P > 0.05, which was not statistically significant.

Figure 7. Correlation between S100A7 expression and immune infiltration as well as immune checkpoints pan-cancer, and the relationship between S100A7 expression and immune activity scores in breast cancer. (A) Correlation between S100A7 expression and immune infiltration pan-cancer. Correlation between S100A7 expression and the ESTIMATE score (B), immune score (C), and stromal score (D) in breast cancer. (E) Relationship between S100A7 expression and immune activity score in breast cancer. STEP1: Release of cancer cell antigens; STEP2: Cancer antigen presentation; STEP3: Priming and activation; STEP4: Trafficking of cells to tumors; STEP5: Infiltration of immune cells into tumors; STEP6: Recognition of cancer cells by T cells; STEP7: Killing of cancer cells. (F) Correlation between S100A7 expression and immune checkpoints pan-cancer. The expression of genes related to inhibitors (G) and stimulators (H) of immune checkpoints in the two groups based on the median S100A7 expression in breast cancer. “ns” represents not significant. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001.

Further, the immune activity score was higher in the S100A7-high than in the S100A7-low expression group (Figure 7E).

Finally, correlation analysis showed that S100A7 was related to immune checkpoints pan-cancer (Figure 7F). A comparison of the expression of immune checkpoints between the S100A7-high and -low expression groups further verified the close association between S100A7 expression and immune checkpoints (Figure 7G, 7H).

S100A7 is associated with methylation, gene mutation, tumor heterogeneity, and stemness

We also studied the relationship between methylation and S100A7 expression in cancer. S100A7 expression was closely correlated with the m⁶A, m⁵C, and m¹A modifications (Figure 8A). It was also associated with the DNA methylation sites cg00325910, cg02892624, and cg17421062 pan-cancer (Figure 8B). Moreover, the expression of m⁶A-, m⁵C-, and m¹A-related genes differed between the S100A7-high and -low expression groups (Figure 8C).

Figure 8. Correlation between methylation and S100A7 expression pan-cancer. (A) Correlation between S100A7 expression and m⁶A, m⁵C, and m¹A pan-cancer. (B) Correlation between DNA methylation and S100A7 expression pan-cancer. (C) The expression of genes related to m⁶A, m⁵C, and m¹A in the two groups based on the median S100A7 expression in breast cancer. “ns” represents not significant. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001. Abbreviations: m6A: N⁶-methyladenosine; m1A: N¹-methyladenosine; m5C: 5-methylcytosine.

Next, we analyzed the differences in the gene mutation landscape of S100A7 between the S100A7-high and -low expression groups (Figure 9A). Copy number variation was similar between the two groups pan-cancer (Figure 9B), but single nucleotide variation was associated with S100A7 expression (Figure 9C). In the S100A7-high expression group, the genes most susceptible to mutation were TP53, PIL3CA, TTN, GATA3, and CDH1 (Figure 9D).

Figure 9. Relationship between gene mutation and S100A7 expression. (A) The pan-cancer gene mutation landscape of S100A7. The relationship between S100A7 expression and CNV (B) and SNV (C) pan-cancer. (D) The gene mutation landscape in the two groups based on the median S100A7 expression in breast cancer. “ns” represents not significant. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001. Abbreviations: CNV: copy number variation; SNV: single nucleotide variation.

Finally, we explored the relationship between S100A7 expression and tumor stemness and heterogeneity. S100A7 expression was associated with tumor heterogeneity (TMB, MATH, microsatellite instability, neoantigen, purity, ploidy, homologous recombination deficiency (HRD), and loss of heterozygosity (LOH) scores; Figure 10A). In breast cancer, the tumor heterogeneity indicators that most correlated with S100A7 expression were HRD, TMB, and LOH (Figure 10A). The tumor stemness indicators RNAss, ENHss, DMPss, DNAss, EREG-METHss, and EREG-EXPss were also associated with S100A7 expression (Figure 10B). In breast cancer, DMPss and EREG-EXPss were associated with S100A7 expression (Figure 10C).

Figure 10. Correlation between S100A7 expression and the tumor stemness and heterogeneity pan-cancer. (A) Correlation between S100A7 expression and tumor heterogeneity pan-cancer and in breast cancer. (B) Correlation between S100A7 expression and the tumor stemness score pan-cancer. (C) Tumor stemness score in the two groups based on the median S100A7 expression in breast cancer. “×” represents P > 0.05, which was not statistically significant. “ns” represents not significant. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001.

S100A7 promotes breast cancer cell proliferation, invasion, migration, and chemoresistance in vitro

We performed in vitro experiments to verify the functions of S100A7 in cancer. The CCK8 proliferation assay showed that S100A7 knockdown inhibited the proliferation of MCF7 and MDA-MB-231 cells (Figure 11A), which was accompanied by a downregulation of the proliferation biomarkers Ki67 and proliferating cell nuclear antigen (Figure 11B).

Figure 11. S100A7 promotes the proliferation, migration, and invasion of breast cancer cells in vitro. (A) The CCK-8 assay showed that knockdown S100A7 reduced the proliferation of MCF7 and MDA-MB-231 cells. (B) RT-qPCR showed that knockdown of S100A7 reduced the expression of the proliferation biomarkers Ki67 and PCNA. (C) The Transwell assay showed that knockdown of S100A7 attenuated the migration and invasion of MCF7 and MDA-MB-231 cells. (D) RT-qPCR showed that knockdown of S100A7 reduced the expression of the migration biomarkers MMP2 and MMP9. (E) The wound healing assay showed that knockdown of S100A7 decreased the migration of MCF7 and MDA-MB-231 cells. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001. All experiments were repeated thrice. Abbreviations: CCK-8: Cell Counting Kit-8; RT-qPCR: reverse transcription-quantitative PCR; PCNA: proliferating cell nuclear antigen; MMP: matrix metalloproteinase.

The Transwell assay demonstrated that S100A7 was associated with the invasion and migration of breast cancer cells (Figure 11C). The migration biomarkers matrix metalloproteinases-2 and -9 were also downregulated after S100A7 was knocked down (Figure 11D). The wound healing assay showed that S100A7 knockdown inhibited the migration of MCF7 and MDA-MB-231 cells (Figure 11E).

We also used Taxol to explore the relationship between S100A7 expression and sensitivity to chemotherapy. The CCK8 assay showed that the IC₅₀ of breast cancer cells to Taxol diminished when S100A7 was knocked down (Figure 12A). The wound healing assay further showed that S100A7 knockdown inhibited the migration of breast cancer cells after Taxol treatment, indicating that S100A7 contributed to chemotherapy resistance in breast cancer (Figure 12B, 12C).

Figure 12. S100A7 reduced the chemosensitivity of breast cancer cells in vitro. (A) The CCK-8 assay showed that knockdown of S100A7 reduced the IC₅₀ of MCF7 and MDA-MB-231 cells towards Taxol. The wound healing assay showed that knockdown of S100A7 inhibited the migration of (B) MCF7 and (C) MDA-MB-231 cells after Taxol treatment (MCF7: 0.05 mg/L, MDA-MB-231: 1 mg/L). ^*P < 0.05; ^**P < 0.01; ^***P < 0.001. All experiments were repeated thrice. IC₅₀, half-maximal inhibitory concentration.