Gathering and processing public data

We obtained target tRNAs associated with the prognosis of BC from the OncotRF database (http://bioinformatics.zju.edu.cn/OncotRF). We downloaded the GSE20711 and GSE20685 datasets, which include the expression profiles of 92 and 327 breast cancer patients’ tumor and normal tissues, from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo//). For the purpose of developing and validating our risk model, we accessed data on BC patients from TCGA-BRCA, including 64 patients with metastases and 19 patients with original tumors, from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). In order to investigate potential relationships between t-mRPM (t-mRNA related prognostic model) and its model genes across various cancers, we also obtained the processed RNA-Seq data of 32 solid malignant tumor types from the UCSC Xena database (https://xena.ucsc.edu/). For an extensive pan-cancer analysis, the 32 processed scRNA-Seq datasets and annotations of 32 solid malignant tumor types were additionally acquired from the TISCH2 database (http://tisch.comp-genomics.org). The TISCH2 database has already undergone quality control, normalization, unsupervised clustering, and cell type annotation for these datasets.

Screening and functional analysis of t-mRNA related to BC prognosis

Using the OncotRF database, we screened tRNAs differentiating breast cancer from normal tissues that influence prognosis, selecting from 3ʹ-tRF and 5ʹ-tRF sequences. We calculated Spearman correlation coefficients between these tRNAs and mRNA expression profiles, defining t-mRNAs as those with a correlation coefficient > 0.4 and p < 0.01, which were then visualized using Cytoscape [23]. We used Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to comprehend the biological roles of these mRNAs [24].

Identification of differentially expressed t-mRNA related to BC prognosis

We first screened prognosis-related genes from t-mRNAs using the Boruta algorithm [25], a random forest-based method for feature ranking and selection, to screen prognosis-related genes from t-mRNAs. Next, we examined the mRNA expression matrix between BC samples and normal breast samples using the “limma” package [26]. A false discovery rate (FDR) < 0.05 and an absolute log2 (fold change) > 1 were the requirements assigned for finding differentially expressed mRNAs (DEmRNAs). From these results, we obtained BC-related genes and the selected t-mRNA genes. By intersecting these genes, we acquired DEt-mRNAs for further analysis.

Construction of t-mRPM using machine learning-based integrative approaches

Ten machine learning algorithms and their 101 combinations were applied to guarantee great accuracy and stability of our t-mRPM. Among them were the following: supervised principal components (SuperPC), generalized boosted regression modeling (GBM), random survival forest (RSF), elastic net (Enet), Lasso, Ridge, stepwise Cox, CoxBoost, and survival support vector machine (survival-SVM). Based on the highest average concordance index (C-index), the best model was chosen. We used three datasets (TCGA-BRCA, GSE20711, and GSE20685) for Kaplan-Meier analysis in order to verify the stability and reproducibility of the model. The results were displayed using survival curves. Additionally, we plotted the model’s 1-year, 3-year, 5-year, 7-year, and 9-year ROC curves to evaluate the predictive performance.

Association of clinical-pathological features and construction and validation of the nomogram

We examined the relationship between the risk model and pathological features in order to assess the clinical value of t-mRNA. Additionally, we analyzed the gene expression patterns across several pathological features in the risk model, which were shown in a heatmap. Furthermore, the TCGA BC metastasis database was used to evaluate the predictive potential for metastasis of BC patients. The model was then integrated with clinical-pathological variables to produce the nomogram, and calibration and decision curves were used to assess the nomogram’s clinical efficacy.

Gene set variation analysis (GSVA) and single-sample gene set enrichment analysis (ssGSEA)

The R package “GSVA” [27] was used in our investigation to calculate scores for fifty HALLMARK pathways. The “limma” program was then used to examine pathways that showed notable variations between high- and low-risk groups. In order to determine the contribution of tumor-related molecular mechanisms in different risk groups, we used the R package “clusterProfiler” [28] to perform a GSEA of HALLMARK gene sets with FDR < 0.25 and |NES| > 1.

Immune feature correlation analysis

Using ssGSEA, we evaluated the variations in immune-related pathways between the groups [29]. The CIBERSORT and xCell algorithms [30] were used to quantify the relative abundance of tumor-infiltrating immune cells (TIICs) within tumor samples. Based on recent advances in immunotherapy for BC [31], we used the IMgor210 dataset to predict the effects of immunotherapy in BC patients using our model.

Performance of the risk model and its genes in pan-cancer

We sought to further assess the performance of the model and its genes regarding expression, mutations, copy number variations, and methylation across various cancers. Additionally, we examined the impact of the t-mRPM on the prognosis of different cancers, affirming the model’s broad applicability across various tumor types. We also assessed the relationship between the model and features of tumor malignancy such as the cell cycle, angiogenesis, and epithelial–mesenchymal transition (EMT).

Expression of model genes at the pan-cancer single-cell level

We obtained single-cell gene expression information from the TISCH2 database, which was used to investigate the relationship between the tumor microenvironment and t-mRNA by examining the single-cell expression of model genes in 32 distinct solid malignant cancer types.

Quantitative real-time PCR (RT-qPCR)

Human mammary epithelial cells (HS578BST) and breast cancer cell lines (MDA-MB-231 and MCF-7) were procured from Fenghui Biotechnology Co., Ltd. (Hunan, China). These cells were maintained in DMEM medium supplemented with 10% fetal bovine serum (Gibco) and incubated at 37°C with 5% CO₂. Total RNA was isolated from cultured cells using TRIzol reagent (Tiangen, Beijing, China) according to the manufacturer’s protocol. Then the total RNA was reverse-transcribed into cDNA using RevertAid Reverse Transcriptase (Thermo Fisher Scientific, Waltham, MA, USA). The RT-qPCR analysis was subsequently carried out on the ABI QuantStudio 1 Plus instrument using PerfectStart Green qPCR SuperMix. The reaction conditions were as follows: initial denaturation at 94°C for 35 s, annealing at 60°C for 15 s, and extension at 72°C for 10 s. Relative gene expression levels were quantified using the 2^–∆∆Ct method. The primers used in this study were designed by Pulateze Biotech (Hunan, China), as listed in Table I.

Statistical analysis

R software (version 4.1.1) was used for the analysis. Pearson’s χ² test was utilized to assess categorical data, and one-way ANOVA was employed to investigate continuous variables. A p-value of less than 0.05 was considered statistically significant. The OS of individuals with breast cancer was examined using the Kaplan-Meier method. Cox regression models, both univariate and multivariate, were used to find independent predictive markers and clinical features that varied significantly between patient groups. Using the Mann-Whitney U test, differences in immune cell infiltration were evaluated. Statistical significance was defined as p-values < 0.05 (*p < 0.05, **p < 0.01, ***p < 0.001, ns: not significant).

Identification of DEt-mRNAs

Using the OncotRF database, we first chose a prognostically significant tRNA (5ʹ-M-tRNA-Tyr-GTA-7-1_L23) in BC. We next calculated the Spearman correlation coefficient to find 92 mRNAs that were significantly associated with this tRNA (Figure 1 A). Meanwhile, we also found that these t-mRNAs are mainly involved in pathways that include the cell cycle, infection of the human T-cell leukemia virus, the p53 signaling pathway, maturation of oocytes through progesterone, meiosis of oocytes, cellular senescence, leukocyte transendothelial migration, and the FoxO signaling pathway (Figure 1 B). We identified 62 genes associated with BC prognosis using the Boruta method (Figures 1 C, D). In TCGA-BRCA cohorts, we found 296 DE-mRNAs, which were visualized in a volcano plot (Figure 1 E). Notably upregulated genes in this plot include MMP11, COL10A1, COMP, COL11A1, and S100P. Finally, these 62 prognostic genes were intersected with the DE-mRNAs using a Venn diagram (Figure 1 F), yielding 22 key genes that were designated as DEt-mRNAs.

Figure 1

Identification of DEt-mRNA. A – t-mRNA related network diagram. B – Biological process involvement by t-mRNA. C – t-mRNAs related to BC prognosis: identified by the Boruta algorithm. D – Feature selection detail diagram: shows t-mRNAs selected by the Boruta algorithm across different iterations E – Volcano plot: displays DE-mRNAs between normal and tumor tissues in BC patients. F – Venn diagram: used to identify common genes between differentially expressed genes (DEGs) and t-mRNAs in BC

Construction and validation of t-mRPM

An integrative technique based on machine learning was applied to the 22 DEt-mRNAs in order to generate a consensus tRNA-related mRNA signature. Even though Ridge had the greatest average C-index (0.616), we decided that the Lasso and CoxBoost combo was the best algorithm pair. This combination achieved the greatest average C-index value, 0.611, among all combinations that incorporated Lasso, which affected the decision. Additionally, the combination was required to obtain coefficient values from Lasso regression for a subsequent pan-cancer study (Figure 2 A). According to the LOOCV framework (Figure 2 B), the optimal λ for the Lasso regression was found at the point where the partial likelihood deviance was lowest. A definite collection of 12 t-mRNAs was identified by applying CoxBoost proportional hazards regression analysis to t-mRNAs with nonzero Lasso coefficients. These t-mRNAs are MELK, CENPF, TSPAN7, BIRC5, NEK2, TOP2A, GPIHBP1, COX7A1, SCN4B, ANLN, UBE2T, and KIF4A (Figure 2 C). Next, the expression levels of these 12 t-mRNAs, weighted by their regression coefficients in the CoxBoost model, were used to calculate the risk score for each BC patient, referred to as the t-mRPM score. The ‘survminer’ R package was implemented to determine the optimal cutoff value, classifying BC patients into high-risk and low-risk groups. As demonstrated in the TCGA training dataset and the GSE20585 validation set, overall survival was significantly lower in the high-risk group compared to the low-risk group (p < 0.05) (Figures 2 D–F). However, in the GSE20711 validation set, there was no significant difference in survival rates between the high- and low-risk groups (p = 0.21), likely due to the smaller sample size of this validation cohort. ROC analysis assessed the discrimination capability of t-mRPM, with 1-, 3-, 5-, 7-, and 9-year AUC values of 0.72, 0.66, 0.66, 0.66, and 0.63, respectively, in TCGA-BRCA, indicating the model’s effective predictive power for the OS of BC patients across these time points (Figure 2 G).

Figure 2

Construction and validation of t-mRPM. A – 101 different kinds of prediction models were created using the LOOCV framework. B – Logarithmic sequence graph (λ) illustrating the t-mRNA by LASSO regression analysis. C – CoxBoost regression coefficients for the last set of 12 t-mRNAs D–F – Kaplan-Meier curves for OS based on t-mRPM in several datasets. G – The optimal model’s ROC curves

Correlation between t-mRPM score and clinical-pathological features

The TCGA dataset’s analysis of the relationship between t-mRPM risk scores and clinical-pathological features showed significant differences (p < 0.05) in TNM staging, pathological grading, survival status (fustat), and gender between the high- and low-risk groups (Figures 3 A, B). As seen in Figure 3 C, Kaplan-Meier curve analysis indicated that the high-risk group had a worse prognosis than the low-risk group (log-rank test, p < 0.01). Since early BC (stages I and II) still accounts for a significant portion of diagnoses, we first looked at how this risk model varied across different T stages. Patients in later stages (T3-4) showed greater risk ratings than those in early stages (T1-2) (Figure 3 D; p < 0.01, Wilcox test). Furthermore, ROC curve analysis showed that the AUC value for t-mRNA predicting M staging was 0.668 (Figure 3 E), suggesting that the t-mRPM score might also predict the development of distant organ metastasis in BC patients. t-mRPM was also found to be predictive of M staging (M0 and M1). We investigated the predictive power of model genes for BC metastasis using the TCGA-BRCA dataset. While individual genes were not very predictive, the risk model that consists of these genes showed better predictive power in predicting metastases of breast cancer (p = 0.032) (Figure 4). According to this model, metastasis is more likely to occur in patients with higher scores. All things considered, the t-mRPM score enhances the prediction for BC survival and is a useful predictor of cancer metastasis.

Figure 3

Analysis of correlations between various clinicopathological features and the prognostic model in the TCGA cohort. A – Pie chart showing the relationship between clinical-pathological traits and high- and low-risk groups. B – Heatmap showing how six distinct clinicopathological traits are distributed along with each patient’s risk score determined by their signature C – Kaplan-Meier OS curves for breast cancer in its early stages based on t-mRPM (log-rank test: p = 0.0067). D – Risk scores were considerably greater in patients in stages T3–4 than in those in stages T1-2 (p = 0.0018). E – ROC curves for BC patients showing M stage prediction (M0 and M1, AUC = 0.668)

Figure 4

Expression levels of individual model genes in the BC metastasis cohort. Green represents the non-metastasis group (primary) and orange represents the metastasis group (Metastasis)

Development and evaluation of nomogram

To confirm the independent predictive capability of t-mRPM, both univariate and multivariate Cox regression analyses were conducted (Figures 5 A, B). In the univariate analysis, all clinical-pathological characteristics, except for gender (HR = 0.815, 95% CI = 0.114–5.836, p = 0.839), i.e. age, TNM staging, pathological staging, and risk score, were significantly correlated with prognosis, each indicating a poor prognosis (all p < 0.05). These included age (HR = 1.034, 95% CI = 1.02–1.047, p < 0.001), T stage (HR = 1.763, 95% CI = 1.225-2.538, p = 0.02), M stage (HR = 5.649, 95% CI = 3.387–9.423, p < 0.001), N stage (HR = 2.189, 95% CI = 1.533–3.127, p < 0.001), overall stage (HR = 2.628, 95% CI = 1.878–3.678, p < 0.001), and risk score (HR = 2.576, 95% CI = 1.887–3.517, p < 0.001) (Figure 5 A). In the multivariate analysis, the risk score (HR = 2.531, 95% CI = 1.802–3.556, p < 0.001) and M stage (HR = 3.541, 95% CI = 1.883–6.659, p < 0.001) emerged as independent prognostic factors for predicting adverse OS in BC patients (Figure 5 B). A predictive nomogram including clinical-pathological variables was developed to predict an individual’s overall survival at 1, 3, and 5 years in order to enhance clinical decision-making. This nomogram closely matched the optimal prediction model (Figures 5 C, D), demonstrating high accuracy in predicting 1-year, 3-year, and 5-year OS. The decision curve analysis (DCA) (Figure 5 E) and calibration chart (Figure 5 D) demonstrated that t-mRNA is a reliable prognostic indicator for BC patients, confirming the nomogram’s accuracy in predicting survival probability with real outcomes. This enhances predictive models and facilitates clinical judgment.

Figure 5

Analysis of the clinical pathological characteristics of t-mRPM and its independent prognostic potential. A, B – Outcomes of the univariate and multivariate Cox regression analyses. C – A nomogram that can be used to forecast a patient’s overall survival. D – The nomogram’s decision curve values. E – The nomogram model calibration curve

To explore the difference of biological functions in two risk groups and the research value of model genes in multiple cancer types, we conducted pathway enrichment analysis, immune correlation analysis, and pan-cancer analysis, as detailed in the supplementary materials.

Validation of the 12 model genes with RT-qPCR assay

To delineate differential expression patterns of model genes between malignant and normal mammary epithelial cells, we conducted quantitative real-time PCR (RT-qPCR) assays. As shown in Figure 6, compared with normal breast cells (HS578BST), the expression of CENPF, TSPAN7, NEK2, TOP2A, COX7A1, SCN4B, and ANLN was downregulated in tumor cell lines (MCF-7 and MDA-MB-231). Moreover, distinct expression patterns of the signature genes were observed across diverse breast cancer cell subtypes. MELK, BIRC5, UBE2T, and KIF4A showed high expression levels in MDA-MB-231, but with low expression in MCF-7. In contrast, GPIHBP1 exhibited an inverse expression pattern.

Figure 6

Validation of model genes in malignant and normal mammary epithelial cells