The RNA-seq data and clinical information of patients were obtained from the TCGA and GTEx databases (http://xena.ucsc.edu/). The BC cohort contains 1097 cancer samples and 120 paracancer samples. The pan-cancer dataset contains 10,534 samples covering 33 different types of cancers. All the TCGA files are stored in FPKM standard format. The validation dataset was downloaded from the Metabric database (https://www.cBioPortal.org), which provides RNA-seq data and related clinical information of 1918 BC patients. In addition, we downloaded the scRNA-seq profiles of BC from the GEO database (https://www.ncbi.nlm.nih.gov/GSE114727) and selected 11 normal breast single-cell samples and 30 tumor single-cell samples from it for subsequent analysis.

The Glycosylphosphatidylinositol (GPI) - anchor biosynthesis related genes were obtained from the relevant literature as indicated in Supplementary Table 1 (19). The scoring genes for cytotoxic T lymphocytes (CTL) are sourced from the String database (https://cn.string-db.org/) ( Supplementary Table 2 ). The chemokine score consists of 41 known chemokines ( Supplementary Table 3 ). All the scores are determined by the expression of scoring genes with the ssGSEA Algorithm (20). The cytolytic activity score (CYT-score), which quantifies the cell-killing function, is based on the geometric mean of gene expressions of GZMA and PRF1 (21).

Package ‘GSVA’, was employed to calculate GPI-anchor biosynthesis score (GPI-score) for patients in the TCGA pan-cancer dataset. Differences in GPI-scores between various tumors and normal samples were assessed using the “ggpubr” software package. The association between GPI-scores and prognostic indicators of patients was investigated by univariate Cox regression analysis (uniCox).Prognostic indicators include: overall survival (OS); disease-free interval (DFI); progression-free interval (PFI); and disease-specific survival (DSS) (22). And the Kaplan-Meier method was employed to assess differences in OS across distinct groups (23).

The differential gene analysis among various groups was conducted using the “limma” package, and | log2 fold change (FC) |≥ 0.5 and p< 0.05 were the standard for the differential genes. Finally, KEGG and GSEA enrichment analysis of the differential genes was performed on SangerBox portal (https://www.sangerbox.com) (24).

Apply the “ssGSEA” algorithm to assess the enrichment scores of 22 immune cell types in every individual, while employing the deconvolution method with LM22 immunological feature matrix in CIBERSORT for precise calculation of infiltrating immune cell components (25). In addition, Stromal score and immune score were calculated by “ESTIMATE” algorithm (26).

Previous studies have demonstrated that interleukin-2 (IL-2), tumor necrosis factor (TNF), interferon-gamma (IFN-γ), and the cytotoxic T lymphocyte (CTL) can serve as quantitative markers for assessing the abundance of T cell exhaustion (27, 28). Consensus clustering analysis was performed on the TCGA BC cohort based on the aforementioned characteristics of T cell exhaustion with K-Means clustering algorithm. The R package “ConsensusClusterPlus” was used to determine the number of clusters with the parameters were set as reps=50, maxK=9, pFeature=1, pItem=0.8. And 1000 repetitions were carried out to ensure the stability of the subtype (29).

The “Seurat” package was used for quality control of data, which involved excluding cells exhibiting ≥5% expression of mitochondria-related genes and those with fewer than 50 detected genes. Further normalization of the gene expression matrix is conducted based on the top 3000 variable genes. The significant principal components (PCs) were computed by principal component analysis (PCA), and the dimensionality was subsequently reduced for the first 20 PCs. Then, a cluster classification analysis was conducted on all units. With the help of the CellMarker database and FindAllMarkers function, we annotated different cell clusters. There are significant differences in the chromosomal copy number variations (CNV)between malignant cells and normal cells in solid tumors. The InferCNV algorithm is widely used for exploring scRNA-seq data to provide evidence of large-scale chromosomal copy number alterations, thereby effectively identifying malignant cells (30). We utilized the InferCNV algorithm to identify malignant cells within clusters of epithelial cells. The GPI-score in different cell clusters was calculated by the AddModuleScore algorithm implemented in Seurat. The Monocle 2 algorithm was employed for pseudo-time trace analysis, and the plot_pseudotime_heatmap function was utilized to construct heatmaps illustrating the dynamic expression of genes (31).

We performed ssGSEA analysis on the characteristics of T-cell exhaustion (including IL2, TNF, IFNG, TGFB1, IL10, TCF7, TBX21, TOX expression abundance, CTL score, Chemokines score, and CYT score) to derived a T-cell exhaustion score (TEX-score). Based on the median TE-score, we divided BC patients into two types (Severe exhaustion type and mild exhaustion type). Different exhaustion states are included as the outcomes of the model. The XGBoost algorithm and the Logistic regression algorithm were used for model training, and the TEX-score of patients was used as the gold standard to evaluate the predictive value of these two models for T-cell exhaustion. The ROC analysis for outcome prediction was conducted using pROC’s ROC function, and the final AUC result was determined by utilizing pROC’s ci function to obtain both AUC and confidence intervals.

Integrating GPI-anchor biosynthesis related genes and T cell exhaustion related genes (IL2, TNF, IFNG) as input features, we arranged 101 combinations of 10 algorithms in the training dataset for variable selection and model construction based on a ten-fold cross-validation framework (32). Optimal modeling method selection is based on c-din values, followed by feature gene screening, calculation of the proportion of feature genes in the prognostic model, and subsequent construction of a risk score.

25 BC patients were recruited from the First Affiliated Hospital of Wenzhou Medical University ( Supplementary Table 4 ). Breast tumor tissue was processed into 4μm thick paraffin sections and subsequently dewaxed. The sections were then incubated in 0.01M sodium citrate at room temperature for 100 minutes to eliminate endogenous enzyme activity. Subsequently, a solution containing 5% goat serum was applied and incubated at room temperature for 60 minutes, followed by overnight incubation with DAPI (1:100), CD8 (1:100), PIGU (1:100), and GPAA1 (1:100) antibodies, respectively. On the following day, after washing with 0.02% PBST, a secondary antibody solution (concentration of 1:200) was added under light-avoiding conditions and incubated at room temperature for 2 hours. Finally, the sections were washed with 0.02% PBST, mounted, and observed using confocal optical microscopy to visualize the staining results.

Small interfering RNA (siRNA) was designed by Tongyong (Shanghai)Company to knock down PIGU mRNA in BRCA cell lines. Total RNA was extracted by Trizol reagent (RNAiso Plus-9109, Takara), and cDNA was synthesized through PCR thermal cycler (Eppendorf, Germany) using reverse transcription reagent (RR036A, Takara) according to the instructions. An appropriate amount of cDNA was used as a template to verify the efficiency of knock down for qRT-PCR through TB Green dye (RR820A, Takara) on Mastercycler real-time PCR system (LightCycler, Roche). After targeted gene knock down, siNC and siPIGU cells were seeded in 96-well plates and cultivated for 48h. CCK8 reagent (CK04, Dojindo) was added and incubated for 2h. The absorbance was measured at 450nm by microplate reader (Tecan, Switzerland).

The t-test was applied to the normally distributed data in this study, while the Wilcoxon rank sum test was used for analyzing the non-normally distributed data. The chi-square test was employed to compare categorical variables. In order to compare the disparities among small sample sizes, a non-parametric test is selected. Pearson correlation analysis is used to measure the strength and direction of the linear relationship between two variables. All the analysis were conducted on the R software (version 4.1.2).And adhering to standard practice, the results are statistically significant with p< 0.05.

We quantified the intensity of GPI-anchored biosynthesis in patients with GPI-scores and evaluated the differences in GPI-anchored biosynthesis between cancer and paracancer tissues in the TCGA pan-cancer dataset. Figure 1A shows that there are 20 cancer types exhibit significant changes in GPI-anchored biosynthesis. The levels of GPI-anchored biosynthesis in the tumor tissues of LGG, TGCT, LUAD, READ, BRCA, THCA, LIHC, LUSC, SKCM, PRAD and GBM were higher than those in normal tissues. Conversely, the levels of GPI-anchored biosynthesis in the tumor tissues of DLBC, OV, KIRC, KICH, STAD, ESCA, KIRP, PAAD and ACC were lower than normal tissues. As the result of unicox analysis shown in the forest diagram ( Figure 1B ), GPI-anchored biosynthesis is a risky factor of prognosis in LGG, HNSC, GBM, CESC and BRCA, while it is a protective factor in READ, THYM, OV, MESO, LUSC and KIRC. The results from the KM plotter demonstrated that the enhancements to the GPI-anchored biosynthesis was associated with poor prognosis in HNSC, BRCA, LGG and LUAD ( Figure 1C ). However, in COAD, LIHC, LUSC, MESO, READ and KIRC, the stronger the GPI-anchored biosynthesis, the more favorable the prognosis for patients. Combined the results of the cox regression, K-M survival curve and the expression of GPI-scores in tumor and normal samples, we finally determined that the GPI-anchored biosynthesis may exert a pivotal role in the development of breast cancer.

The TCGA BC patients were divided into two cohorts with high and low GPI-anchored biosynthesis based on the best cutoff values of the GPI-score. The heat map illustrates the distinct gene patterns within various patient subgroups ( Figure 2A ). KEGG functional enrichment analysis of different clusters of patients showed that the patients with high GPI-anchored biosynthesis had a worse prognosis accompanied by a lower immune response, the specific manifestations include the down-regulation of Cytokine-cytokine receptor interaction, Helper T lymphocytes (Th17, Th1 and Th2) cell differentiation, Primary immunodeficiency, and Antigen processing and presentation ( Figure 2B ).Interestingly, the patients in the high GPI-anchored biosynthesis group exhibited augmented metabolic pathways, characterized by the Estrogen signaling pathway and AMPK signaling pathway, as well as the enhanced Fatty acid metabolism and Cortisol synthesis secretion ( Figure 2B ). In addition, GSEA analysis showed that as the GPI-score increased, protein secretion and estrogen response were elevated, peroxisome metabolism and oxidative phosphorylation pathway were enhanced. However, IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, TNFA signaling, Interferon Gamma response and Inflammatory response were downregulated ( Figure 2C ). We also examined the abundance of 22 immune cells in different groups. Surprisingly, the group of patients with high GPI-score showed a decrease in adaptive immune function, with a decrease in T cells CD8 ,T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs),and T cells gamma delta. It was also accompanied by an increasing of M2 macrophages (tumor-promoting) and monocytes ( Figure 2D ) (33). Simultaneously, the analysis of ESTIMATE scores and immune checkpoints revealed that patients with higher GPI-anchored biosynthesis exhibited a diminished immune status and response ( Figure 2E ). Additionally, these patients exhibited decreased stromal cell infiltration in their tumor tissues but elevated levels of tumor purity ( Figures 2F, G ).

GPI-score is negatively correlated with T cell CD8, T cell CD4 memory activation, Treg, B cell memory, and T-cell follicular helical (P<0.001), while showing a positive correlation with M2 macrophages (P<0.001) ( Figure 3A ). As it is depicted in Figure 3B , the expression of the majority of genes involved in GPI-anchored biosynthesis exhibited a significant negative correlation with diverse immune cell populations. Then, consensus clustering analysis was employed to categorize BC patients based on the specific TEX pathway which includes TNF, IL-2, IFN-γ, and CTL. We identified 4 clusters of patients with different T-cell exhaustion status in the TCGA-BC cohort (C1-C4) ( Supplementary Figure 1 ). The heat map also shows the change of other molecular pathways of TEX in different groups, including TGFB, IL10, chemokines, transcription factors (TCF7, TBX21, Tox), CYT scores, and the abundance of infiltrating lymphocytes (Activated B cell, Activated CD4 T cell, Activated CD8 T cell, Effector memeory CD4 T cell, Effector memeory CD8 T cell, Th1, Th2, and Th17) ( Figure 3C ). Consistently, these characterized pathways of TEX including CYT and CD4 T cell and CD 8 T cell showed a decreasing trend in C1-C4 (Figures 3D–F), suggesting that the four TEX subgroups identified in this study accurately represent the biology of the TEX grading stage. Interestingly, the intensity of GPI-anchored biosynthesis increases with the severity of the TEX ( Figure 3G ). Based on the KEGG enrichment analysis, we compared the differences in biological function between the groups with the greatest differences in the degree of TEX (C1 and C4), and surprisingly, alterations in GPI-anchored biosynthesis were the most significant change between the two groups, topping the list in terms of magnitude of change ( Figures 3H, I ).

We identified 14 major cell types using annotations from the Cell Marker database, including CD4 T cells, CD8 T cells, Exhausted CD8 T cells, B cells, Epithelial cells, Fibroblasts, Myeloid cells, Innate Lymphocytes (ILCs), Dendritic cells, Monocytes, Neutrophils, Mast cells, Endothelial cells and Macrophages ( Figure 4A ). Projecting the clinical information of the samples onto various cell subpopulations to trace the cellular origin, the exhausted CD8 T cells were predominantly present in tumor tissues ( Figure 4B ). Compared with the normal mammary microenvironment, the proportion of various types of cells in TME was significantly altered, in which the proportion of B cells, endothelial cells, macrophages, and monocytes was elevated while the CD4 T cells and CD8 T cells were decreased and replaced by a substantial increase in exhausted CD8 T cells ( Figure 4C ). This observation suggests that the T-cell exhaustion is a rare occurrence in normal breast tissue, but frequently manifests in tumor tissue. Then, we defined 127 malignant cells in the tumor epithelial cell population using the InferCNV algorithm ( Figures 4D, E ; Supplementary Table 5 ) and compared the differences of GPI-anchored biosynthesis between the malignant epithelial cells and normal breast epithelial cells. Interestingly, the GPI-anchored biosynthesis was instead reduced in malignant epithelial cells ( Figure 4F ). We compared the alterations of GPI-anchored biosynthesis in other tumor-derived cell populations, we found that remarkably differences could be only observed between CD8 T cells and exhausted CD8 T cells ( Figure 4G ; Supplementary Figure 2 ). The GPI-anchored biosynthesis was much stronger in exhausted CD8 T cells compared to normal CD8 T cells. Similarly, GSVA analysis of KEGG biological pathways showed there is no difference in GPI-anchored biosynthesis between malignant and non-malignant epithelial cells ( Figure 4H ), whereas this difference was evident in CD8 T cells and exhausted CD8 T cells ( Figure 4I ).

The differentiation trajectories of CD8 T cells and epithelial cells were constructed to observe the process of CD8 T cell exhaustion and metabolic remodeling of GPI-anchored biosynthesis during tumorigenesis. In the pseudotiming analysis, CD8 T cells were initially located along the trajectory path and gradually transitioned into exhausted CD8 T cells ( Figures 5A, B ). We calculated the dynamic expression levels of GPI-anchored biosynthesis related genes in CD8 T cells and exhausted CD8 T cells along the pseudotemporal transition, and the gene heat map suggested that PIGP, PIGC, PIGU, PIGW, PIGV, and GPAA1 were highly expressed at the end of the trajectory path ( Figure 5C ), which means that these genes are highly expressed in exhausted CD8 T cells. Two-dimensional map shows the dynamic expression of genes in CD8 T cells (green) and exhausted CD8 T cells (blue) ( Figure 5D ). These results suggest that the expression of GPI-anchored biosynthesis related genes is heterologous during the CD8 T cell state transition. Figures 5E, F represent the pseudotemporal developmental trajectories of epithelial cells. We calculated the GPI-scores in the cells and present them in a two-dimensional plot, where we can see that the overall level of GPI-anchored biosynthesis in exhausted CD8 T cells is always higher than that in non-exhausted CD8 T cell, and that metabolic intensity increases as CD8 T cell become more exhausted ( Figure 5G ). However, the intensity of GPI-anchored biosynthesis is lower in tumor epithelial cells than in normal epithelial cells ( Figure 5H ). Thus, it is worth considering whether GPI-anchored biosynthesis related genes could serve as a potential biomarker for predicting T cell exhaustion in BC patients. Interestingly, we did not observe this phenomenon in the metabolic landscape of exhausted CD4 T cells. It is worth noting that GPI-anchored biosynthesis levels are lower in exhausted CD4 T cells compared to non-exhausted CD4 T cells ( Supplementary Figure 3 ).

Correlation analysis was performed between the GPI-score and TEX-score in Metabric and TCGA BC cohort, with correlation coefficients (r-values) of -0.49 and -0.51 respectively ( Figures 6A, B ). Then, a variational analysis was performed in Metabric BC cohorts with the model_profile function in the XGBoost algorithm, assessing the association between the expression of GPI-anchored biosynthesis related genes and the states of T cell exhaustion in BC patients. The Partial Dependence profile demonstrates the impact of individual gene expression on the predictive ability of the XGboost model ( Figure 6C ). Figure 6D illustrates the magnitude of the impact of each GPI-anchored biosynthesis related gene on the states of T cell exhaustion in Metabric BC patients, it can be observed that PIGQ exerts the most significant influence on severe T cell exhaustion states in patients, followed by GPAA1. Finally, XGBoost and Logistic Regression were used to construct a T-cell exhaustion diagnostic model associated with GPI-anchored biosynthesis related genes in two BC cohorts. We evaluated the accuracy of these two models in predicting the occurrence of severe T-cell exhaustion in patients. The area under the ROC curve represents the accuracy of the model. Compared to the XGBoost algorithm, the Logistic regression model has greater AUC values in both the Metabric and TCGA cohorts, with values of 0.688 and 0.684 respectively ( Figures 6G, H ), while the AUC values for the XGBoost model are 0.637 and 0.623 ( Figures 6E, F ). Our results demonstrated that the GPI-anchored biosynthesis related genes are valuable as a tool to predict whether a patient has developed severe T cell exhaustion.

To validate the prognostic value of GPI-anchored biosynthesis related genes in BC, we develop a gene signature (GPIS) based on 21 GPI-anchored biosynthesis related genes and the characteristic factors of T cell exhaustion (TNF, IL-2, IFN-γ).The Metabric dataset, with its high number of patients, served as the training set, while the TCGA cohort, with its low number of samples, served as the validation set. We employed a ten-fold cross-validation framework to train and validate 101 predictive models, subsequently determining the C-index for each set. Models with C-index less than 0.5 indicate obvious overfitting problems in the construction process and should be excluded. Interestingly, the optimal model was Random Forest (RSF) with the highest average C-index (0.767) Figure 7A . Then, the RandomForestSRC algorithm was utilized to optimize the results of the RSF. Figure 7B represents the weight ratio of the characterized genes in GPIS, which shows that PIGV, PIGU, GPAA1, and PGAP1 are the core genes of the model. By weighting the 13 gene expression values and their regression coefficients in RandomForestSRC, a risk score was obtained for each patient ( Figure 7C ; Supplementary Table 6 ). With the optimal thresholds for risk scores determined by the “surviminer” package, patients were classified into high-risk and low-risk groups. K-M curves indicate that patients categorized as high-risk exhibited notably reduced OS duration compared to individuals in the low-risk cohort ( Figures 7D, E ). ROC analysis measured the discrimination of GPIS, and the 3-year, 4-year and 5-year AUC of training dataset were 0.62, 0.64 and 0.65 ( Figure 7F ). The 1-year, 3-year and 5-year AUC in validation dataset were 0.68, 0.60 and 0.59 ( Figure 7G ).The Accuracy of Random Forest model is 0.685, KS score is 0.510, F1 score is 0.645, and Precision is 0.672 ( Supplementary Figure 5 ).

The multivariable Cox regression analysis revealed that GPAA1, PIGU, PIGS, PIGQ, and PIGW independently serve as risk factors for OS ( Figure 7H ). In the feature weight analysis of the T cell exhaustion diagnostic model and prognosis model, further comprehensive investigation is necessary for GPAA1 and PIGU. We employed the Timer method from the IOBR software package to conduct immune cell infiltration scoring on 9,406 tumor samples covering 38 cancer types within the TCGA pan-cancer dataset. The heat map reveals a strong correlation between GPAA1 and immune infiltration in 32 cancer types, while PIGU exhibits a significant association with immune infiltration in 33 cancer types ( Figure 8A ). The infiltration of CD8+T cells in breast cancer exhibits a negative correlation with the expression levels of GPAA1 and PIGU. In addition, the K-M curve indicates that patients with high expression of GPAA1 and PIGU have a worse OS ( Figures 8B, C ). Compared to normal breast cells (MCF10A), the expression of GPAA1 and PIGU is relatively high in MCF7, MDA-MB-231 and BT-474 cells ( Figure 8D ). SiRNA was employed to downregulate the expression of PIGU in MCF7 and MDA-MB-231 cells, and the ability of PIGU to effectively retard proliferation of MCF7 and MDA-MB-231 cells was confirmed by CCK8 assay ( Figure 8E ).We also employed a multiplex immunofluorescence staining technique to detect the expression levels of GPAA1 and PIGU molecules, as well as the degree of CD8 T cell infiltration in 25 local breast cancer tissues. Surprisingly, In multiple microscopic views, the expression of GPAA1 and PIGU in breast cancer cells was negatively correlated with the expression of CD8 molecules in neighboring lymphocytes ( Figures 8F, G ). When tumor cells high expressed GPAA1 and PIGU, the fluorescence intensity of CD8 molecules in peripheral infiltrating lymphocytes noticeably decreased, indicating a reduction in the number of CD8+ T cells.