In July 2021, the TCGA-BRCA dataset (Version 3; The Cancer Imaging Archive. http://doi.org/10.7937/K9/TCIA.2016.AB2NAZRP) (16) was retrieved via the University of California Santa Cruz Xena platform (https://xenabrowser.net), which includes RNA sequencing data, survival statistics and clinical details for both cancerous and normal breast tissue samples. It enables the segmentation of cancer cases into four molecular categories based on the status of ER, PR, and HER2: Luminal A (positive for ER or PR, and negative for HER2), luminal B (positive for ER or PR, and HER2), HER2⁺ (negative for ER and PR, and positive for HER2) and triple-negative (negative for ER, PR and HER2) (2). For the present analysis, the Gene Annotation by the ENCODE Consortium (GENCODE) human gene annotation (release 22) (https://www.gencodegenes.org/human/release_22.html) was used, selecting ‘protein-coding’ tags for mRNA research, and classifying several non-coding RNA types, such as ‘antisense’, ‘sense-intronic’, ‘lincRNA’, ‘ncRNA’, ‘sense-overlapping’ or ‘processed-transcript’, as lncRNAs. The data was organized into separate matrices for mRNA and lncRNA evaluation. Additionally, immune gene data were compiled from the ImmPort database (https://immport.niaid.nih.gov) to create a dedicated immune gene matrix, setting the stage for the subsequent analysis.

In the analysis, differential expression of immune-related genes across luminal A, luminal B, triple negative breast cancer (TNBC) and normal samples compared with HER2⁺ samples was assessed using the ‘limma’ package in R (version 3.5.1; http://bioconductor.org/packages/release/bioc/html/limma.html) (17). The screening criteria of P<0.05 and log₂fold change (FC) >0.5 were selected to identify significantly differentially expressed genes. A Venn diagram analysis was performed using the online tool Venn (http://bioinformatics.psb.ugent.be/webtools/Venn/) to extract the overlapping set of differentially expressed immune genes across the different conditions. This identical approach was also applied to identify differentially expressed lncRNAs, adhering to the same screening criteria of P<0.05 and log₂FC >0.5, to ascertain the common subset of lncRNAs for further investigation.

The intersecting differentially expressed immune genes and lncRNAs identified from the previous steps were analyzed using Pearson correlation analysis to ascertain significant relationships. Significantly correlated pairs exhibiting an absolute value of r>0.3 and P<0.05 were earmarked for further assessment. The outcomes of this screening were graphically represented utilizing the ‘ggplot2’ package in R (version 3.3.2; http://ggplot2.tidyverse.org/), facilitating a visual interpretation of the data's underlying patterns and associations.

The ‘pROC’ package within R software (version 1.12.1; http://www.rdocumentation.org/packages/pROC/versions/1.18.5) was used to create receiver operating characteristic (ROC) curves (18), comparing the performance of immune-related lncRNAs in HER2⁺ breast cancer samples and normal controls. The diagnostic capabilities of these lncRNAs were quantified by measuring the area under the curve (AUC). lncRNAs with AUC>0.7 were deemed as significant, highlighting their potential as effective markers for differentiating between healthy and cancerous conditions.

The optimal cut-off value of each diagnostic marker was determined by the surv_cutpoint function, enabling the division of HER2⁺ breast cancer samples into groups based on high or low expression levels of specific lncRNAs. To assess how the expression levels of these lncRNAs correlated with OS among HER2⁺ patients, the ‘survival’ package in R (version 3.5–3; http://CRAN.R-project.org/package=survival) was used to create Kaplan-Meier survival curves. lncRNAs with P<0.05 were considered significant and were therefore identified as potential prognostic biomarkers for further analysis in the present research.

The correlation between biomarkers and clinicopathological characteristics was analyzed by creating box plots to visualize the differences in biomarker expression across several clinicopathological categories. This was accomplished using the ‘ggplot2’ R package, allowing for a clear and informative representation of the data trends and distributions relative to clinical attributes.

The subcellular locations of prognostic biomarkers were identified using the LncLocator (http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/index.html#) and iLoc-lncRNA database (http://lin-group.cn/server/iLoc-LncRNA/predictor.php), respectively.

The Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) approach (19) was used to perform immune cell infiltration analysis in patients with HER2⁺ breast cancer. Following this, the R package ‘ggstatsplot’ (version 3.6.1; http://indrajeetpatil.github.io/ggstatsplot/) was used to calculate the Spearman's rank correlation coefficient between the identified biomarkers and 22 types of infiltrating immune cells. This analysis aimed to assess the relationship between biomarkers and the presence of immune cells in the tumor environment. Immune cells that demonstrated a correlation coefficient of r>0.3 and P<0.05 were recognized as significantly correlated with the biomarkers in question.

A Pearson correlation analysis was performed to identify the relationships between nuclear and cytoplasmic lncRNAs and immune-related genes, selecting only lncRNA-mRNA pairs with r>0.3 and P<0.05. The lncbaseV2 database (https://dianalab.e-ce.uth.gr/html/diana/web/index.php?r=lncbasev2) was used to predict lncRNA-miRNA relationships, with significance set at a prediction score >0.7. Furthermore, the miRWalk 2.0 platform (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) was used to identify miRNA-mRNA interactions. Only interactions recognized by at least 4/6 reference databases (miRWalk, miRanda, miRDB, PITA, RNA22 and Targetscan) were considered valid. A lncRNA-miRNA-mRNA regulatory network was constructed based on cytoplasmic co-expression relationships. Cytoscape software (version 3.8.0; http://cytoscape.org/) facilitated the visualization of the lncRNA-mRNA co-expression and ceRNA networks for nuclear and cytoplasmic lncRNAs, respectively (20). Finally, functional characterization of mRNAs present in both networks was achieved using Gene Ontology (GO) (https://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) pathway enrichment analyses, using the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.ncifcrf.gov/tools.jsp).

Two tissue microarrays (cat. nos. HBreD080CS01 and TFBrec-01) were purchased from Shanghai Outdo Biotech Co., Ltd., containing 217 breast cancer tissues and 24 adjacent non-tumor tissues. The expression of lncRNA CTC-537E7.2 was evaluated utilizing the ISH Test Kit (Boster Biological Biotechnology; cat. no. MK11221) according to the manufacturer's protocols. In brief, the paraffin-embedded tissue sections were deparaffinized and rehydrated. Endogenous peroxidase activity was blocked by 3% hydrogen peroxide for 10 min and then digested with protease K. The digoxigenin-labeled oligonucleotide probe was added at a concentration of 1.5–2 µg/ml and incubated overnight at 37°C. Next, the sections were treated with a blocking reagent at 37°C for 30 min, followed by incubation with biotinylated rat anti-Digoxigenin at 37°C for 60 min. Subsequently, the sections were treated with Sterptavidin-Biotin-Amplified Complex at 37°C for 20 min, followed by the application of biotinylated peroxidase at 37°C for 20 min. Finally, sections were stained using 3,3′-diaminobenzidine. Oligonucleotide probes marked with DIG-dUTP at the 3′ end were purchased from Boster Biological Biotechnology (cat. no. MK1122). The sequences for these probes were as follows: 5′-AGGTACAGTCATGTGCCGCATAATGACATTCAGTCAATGA-3′; 5′-ATATATGAGGTGGTCCTGTAAGGTGATAATGGAGCTGAAA-3′; and 5′-TGTAACATCATAGCACAACCAATCACCTTTTCTATATTTA-3′.

ISH results were captured using a tissue microarray scanner (Aperio VERSA 8; Leica Biosystems). The captured images were analyzed using ImageScope software (v.11.2.0.780; Leica Biosystems), which automatically detects regions of interest on the tissue sections based on staining intensity: Dark brown signifies strong positivity; light brown indicates moderate positivity; light yellow shows weak positivity; and blue marks negative staining. Quantitative analyses were performed based on the intensity levels, encompassing both the area (in pixels) and the % of positively stained regions. To quantify the expression of lncRNA CTC-537E7.2, a semi-quantitative method known as the histoscore (H-score) was used (21). This scoring system combines the staining intensity (with scores ranging from 0 for negative, 1 for weak, 2 for moderate and 3 for strong) and the % of positive cells at each intensity level. Specifically, the H-score was calculated as follows: H-SCORE=∑(I × Pi)=(% of cells of weak intensity × 1) + (% of cells of moderate intensity × 2) + (% of cells of strong intensity × 3), where I=intensity of staining and Pi=% of stained tumor cells (22,23). This method produces a score ranging from 0–300, serving as a comprehensive quantitative indicator of lncRNA expression.

Statistical analyses were performed using R (version 3.5.1; http://www.r-project.org/foundation/) and incorporated several packages such as ‘limma’, ‘pROC’, ‘ClusterProfiler’, ‘ggplot2’, ‘Cytoscape’, ‘ggstatsplot’ and ‘survival’. The differential expression of mRNAs and lncRNAs was determined using the ‘limma’ package. Correlation analysis was performed using either Pearson or Spearman's correlation, depending on the data distribution and the analysis requirements. ROC curve analysis was used to evaluate the predictive power of each lncRNA regarding clinical outcomes. The Wilcoxon rank-sum test was used to assess the statistical differences in the expression of biomarkers across several clinicopathological categories between two groups. For comparisons involving >2 groups, the one-way ANOVA test with Bonferroni's correction was used. Kaplan-Meier survival curves were generated to compare survival times among groups based on lncRNA expression levels, identifying those with statistically significant impacts. The Kruskal-Wallis test was used to assess differences between multiple groups in experimental data, followed by Dunn's post hoc test for pairwise comparisons. The results are presented as median (interquartile range). P<0.05 was considered to indicate a statistically significant difference, unless stated otherwise.

A total of 803 samples were extracted from the TCGA database, which included 430 luminal A samples, 124 luminal B samples, 37 HER2⁺ samples, 113 TNBC samples and 99 normal samples. Out of these, 689 samples had comprehensive survival and clinical information available. A list of immune-related genes was also compiled from the ImmPort database. After removing duplicates, a list of 1,811 unique genes was finalized for subsequent analysis.

Differential gene expression analysis was performed comparing the HER2⁺ breast cancer group with the luminal A, luminal B, TNBC and normal sample groups. A total of 485 immune genes were identified as differentially expressed in the HER2⁺ breast cancer group compared with the normal sample group, consisting of 185 upregulated and 300 downregulated genes. In comparison with the luminal A breast cancer group, 193 immune genes exhibited differential expression (76 upregulated and 117 downregulated). Similarly, when compared with the luminal B breast cancer group, 135 immune genes showed differential expression (79 upregulated and 135 downregulated). Compared with the TNBC group, 152 immune genes demonstrated differential expression, with 43 upregulated and 109 downregulated genes (Table I). The distribution of these differential genes among each sample group was visually represented using volcano plots (Fig. 1A-D). A Venn diagram analysis using online tool Venn (http://bioinformatics.psb.ugent.be/webtools/Venn/) revealed 23 shared differential immune genes among the four groups (Fig. 1E). Similarly, after analyzing the differential lncRNAs, 662 lncRNAs were found to be differentially expressed in the HER2⁺ breast cancer group compared with the normal sample group. This included 158 upregulated and 504 downregulated lncRNAs. In comparison with the luminal A breast cancer group, 233 lncRNAs exhibited differential expression (55 upregulated and 178 downregulated). Similarly, compared with the luminal B breast cancer group, 138 lncRNAs showed differential expression (36 upregulated and 102 downregulated). In contrast with the TNBC group, 184 lncRNAs demonstrated differential expression, with 68 upregulated and 116 downregulated lncRNAs (Table II). The distribution of these differential lncRNAs was illustrated using volcano plots (Fig. 2A-D). By comparing the differential lncRNAs across the four groups, 22 commonly differentially expressed lncRNAs were identified (Fig. 2E).

In total, 74 significantly correlated pairs (differentially expressed immune-related genes and differentially expressed lncRNAs exhibiting a correlated relationship meeting the standard criteria of r>0.3 and P<0.05) were identified. These pairs included 20 differential immune genes, encompassing C-X-C motif chemokine ligand 17, Dermcidin, Dickkopf WNT signaling pathway inhibitor 1, Fibroblast Growth Factor Receptor 4, Galanin, Growth Differentiation Factor 15, Indoleamine 2,3-Dioxygenase 1, Insulin-Like Growth Factor 1 Receptor, Lipopolysaccharide-Binding Protein, Mesenchymal-Epithelial Transition Factor, Phosphatidylinositol 3, Prolactin Receptor, Proteasome 26S Subunit, Non-ATPase 3, S100 Calcium-Binding Protein A7, S100 Calcium-Binding Protein A7A, S100 Calcium-Binding Protein A8, S100 Calcium-Binding Protein A9, Syndecan 1, Semaphorin 3C, Solute Carrier Family 40 Member 1 and Thymosin Beta 15A, and 19 differential lncRNAs (TMEM92-AS1, ST8SIA6-AS1, RP11-95M15.1, RP11-783K16.5, RP11-612B6.2, RP11-510J16.5, RP11-45215.2, RP11-428L9.2, RP11-369C8.1, RP11-287D1.4, RP11-20F24.2, RP11-206M11.7, RP11-19E11.1, MNX1-AS1, LINC01213, LINCO0993, LA16C-380H5.4, HOTAIR, CTC-537E7.2, CTA-384D8.35, AC087491.2 and AC008268.1). As a result, a set of 19 immune-related differential lncRNAs was obtained. The aforementioned findings were visualized using the R package ‘ggplot2’, as depicted in Fig. 3.

Diagnostic lncRNAs were identified by constructing ROC curves for 19 immune-related differential lncRNAs between normal control and HER2⁺ samples. Out of these, 13 immune-related differential lncRNAs exhibited an AUC >0.7, pointing to a promising predictive value for breast cancer survival (Table III and Fig. S1). Therefore, all 13 immune-related differential lncRNAs with diagnostic potential were included in the analysis, namely AC008268.1 (AUC=0.709), CTA-384D8.35 (AUC=0.909), CTC-537E7.2 (AUC=0.788), HOTAIR (AUC=0.942), LA16c-380H5.4 (AUC=0.879), LINC00993 (AUC=0.717), RP11-95M15.1 (AUC=0.826), RP11-287D1.4 (AUC=0.915), RP11-510J16.5 (AUC=0.854), RP11-612B6.2 (AUC=0.883), RP11-783K16.5 (AUC=0.963), ST8SIA6-AS1 (AUC=0.913) and TMEM92-AS1 (AUC=0.844).

To elucidate prognostic biomarkers for patients with HER2⁺ breast cancer, the optimal threshold value of each diagnostic lncRNA was used to stratify the HER2⁺ samples into high and low expression groups. Subsequently, Kaplan-Meier survival curves were generated using the R package ‘survival’. Notably, among these lncRNAs, CTC-537E7.2 was significant associated with prognosis in HER2⁺ patients (P=0.04), establishing itself as a promising prognostic biomarker (Fig. 4A). Patients with elevated levels of CTC-537E7.2 exhibited notably extended overall survival compared with those with lower expression.

To further assess the association between prognostic biomarkers and clinicopathological factors, a correlation analysis was performed. Using the R package ‘ggplot2’, box plots were generated to visualize the expression values of the biomarker across different clinicopathological factors (Fig. 4B). The results revealed significant differences (P<0.05) in the expression of lncRNA CTC-537E7.2 across several sample groups concerning six clinicopathological factors, namely ER, PR, HER2, age, group and ethnicity.

According to the subcellular localization prediction of lncRNA CTC-537E7.2, the LncLocator database assigned the highest score to the cytoplasm (0.598853; Fig. 4C). Conversely, the iLoc-lncRNA database identified the nucleus as having the highest score (0.905448) for this lncRNA (data not shown). These results indicate that lncRNA CTC-537E7.2 is expressed in both the cytoplasm and nucleus.

To evaluate the potential association between lncRNA CTC-537E7.2 and immune cells in HER2⁺ breast cancer samples, the CIBERSORT algorithm was initially used to assess the proportions of immune cell populations within the samples. This analysis provided insights into 22 distinct immune cell types across all 37 HER2⁺ breast cancer samples. Following this, correlation analyses were performed to evaluate the relationship between the expression levels of lncRNA CTC-537E7.2 and each specific immune cell type. Noteworthy results revealed significant associations between lncRNA CTC-537E7.2 and four immune cell types: M0 Macrophages (r=−0.386; P=0.018), monocytes (r=0.340; P=0.039), neutrophils (r=−0.332; P=0.045) and M2 macrophages (r=0.331; P=0.046; Table IV). Neutrophils and macrophages serve dual roles in both promoting and suppressing cancer development (24,25). Macrophages, derived from monocytes (26), can transition from inactive M0 macrophages to polarized phenotypes (27), primarily classified as classically activated (M1) and alternatively activated (M2) macrophages (28). M2 macrophages, predominant among tumor-associated macrophages, markedly contribute to tumor progression and metastasis (24). Notably, the present study found lncRNA CTC-537E7.2 shows a positive correlation with M2 macrophages, which are known to be activated by cytokines such as interleukin (IL)-4, IL-13, IL-10 and transforming growth factor-β (TGF-β) (28). To further assess the relationship between lncRNA CTC-537E7.2 and macrophages, the correlation between lncRNA CTC-537E7.2 and phenotypic factors of M2 macrophages (IL-10, TGF-β, IL-4, and IL-13) were analyzed. The correlation plot revealed a positive correlation between lncRNA CTC-537E7.2 and IL-4, whilst demonstrating negative correlations with other phenotypic factors, including IL-10, TGF-β and IL-13 (Fig. 4D).

Based on the subcellular localization findings for lncRNA CTC-537E7.2, lncRNA-mRNA and ceRNA networks were constructed. Nuclear lncRNAs directly modulate target gene expression through binding, whilst cytoplasmic lncRNAs act as ceRNAs, influencing gene expression via miRNA interactions (29). Pearson correlation analysis was performed on lncRNA CTC-537E7.2 and immune-related genes, whereby, based on the parameters of r>0.3 and P<0.05, a total of 51 lncRNA-mRNA pairs were identified, consisting of 1 lncRNA and 51 mRNAs. The LNCbaseV2 database was used for predictive analysis for selected lncRNA CTC-537E7.2 and miRNAs relationship pairs with a score >0.7, resulting in 13 pairs comprising 1 lncRNA and 13 miRNAs. The miRWalk2.0 database, with data integrated from six databases: miRWalk, miRanda, miRDB, PITA, RNA22 and Targetscan, was used to predict miRNA-mRNA relationship pairs for the 51 mRNAs associated with lncRNA CTC-537E7.2. This approach identified 6,456 miRNA-mRNA relationship pairs, involving 48 mRNAs and 1,284 miRNAs. Finally, by integrating the co-expression relationships, a ceRNA network was formulated with 24 pairs, consisting of one lncRNA, 8 unique miRNAs and 11 distinct mRNAs. Visualization was achieved using Cytoscape software (Fig. 5).

The results of GO and KEGG enrichment analyses were considered significant when P<0.05 and counts ≥2. Within the lncRNA-mRNA co-expression network, mRNAs were enriched across 40 GO Biological Processes (BP), 8 Cellular Components (CC), 7 Molecular Functions (MF) and 18 KEGG pathways (Fig. S2). In the ceRNA network, mRNAs showed significant enrichment in 12 GO BP, 2 CC, 5 MF and 5 KEGG pathways, including pathways such as Janus kinase-signal transducer and activator of transcription (Jak-STAT) signaling related to macrophage polarization (Fig. S3). The top 10 GO and KEGG terms, selected by count value, are presented due to the extensive number of enriched pathways.

ISH assays were performed to detect the expression of lncRNA CTC-537E7.2 in tissue microarrays. Cases with severe detachment during immunohistochemistry and those lost to follow-up were excluded from the statistical analysis due to tissue chip detachment issues. After the removal of invalid cases, 160 breast cancer tissues and 24 paired cancer-adjacent tissue samples were included in the statistical analysis. The analysis revealed distinctive color reactions for lncRNA CTC-537E7.2 in both breast cancer and adjacent tissues, primarily localized within the nuclei of tumor cells and cancer-adjacent tissues, and occasionally observed in the cytoplasm. The expression of lncRNA CTC-537E7.2 across different groups is summarized in Table V. Notably, compared with other subtypes of breast cancer and adjacent normal tissues, HER2⁺ breast cancer tissues exhibited significantly higher staining signal intensity for lncRNA CTC-537E7.2 (P<0.05; Fig. 6A and B).