3.1 Genetic alterations of granzyme genes in pan-carcinoma
To investigate the pattern of familial dysregulation of the GZM gene in cancer, we analyzed genomic data from tumors and normal tissues for 20 different cancer types. This analysis included genetic variation, somatic copy number changes (SCNA), mRNA expression, and DNA methylation data. Our examination of nonsynonymous mutations in the GZM gene across the 20 cancer types revealed consistently high mutation frequencies in nearly all cancers, with particularly noteworthy frequencies in SKCM and UCEC (Fig. 1A). Notably, family mutations in the GZM gene have shown strong associations with the prognosis of various cancers. Additionally, we observed relatively high mutation frequencies for GZMA and GZMK in multiple cancer types (Fig. 1A). To further explore genetic aberrations in the GZM gene family, we analyzed the percentage of SCNA. Our results demonstrated that SCNA occurred more frequently (in over 5% of all samples) in most cancer types (Fig. 1B). However, all GZM gene families in THCA displayed a lower frequency of SCNA, which may be specific to THCA itself. Moreover, the SCNA signatures of the GZM gene family varied across different cancers. For instance, GZMA and GZMK exhibited similar SCNA frequencies in major cancers, while GZMM showed an opposite pattern to GZMA in GBM and LUAD. To further investigate copy number variations, we examined the correlation between mRNA expression of the GZM gene family and CNV data from the TCGA database. Our findings indicated a positive correlation between mRNA expression and CNV, particularly for GZMA, GZMK, and GZMM (P < 0.05, Fig. 1C) in BLCA, CESC, HNSC, STAD, and UCEC. However, a negative correlation was observed between GZMH and GZMB in BRCA, CESC, PAAD, PRAD, STAD, and THCA. We then proceeded to explore the methylation levels of the GZM gene family to understand its epigenetic regulation. Our analysis revealed highly heterogeneous methylation patterns for the GZM gene families across different tumors. In BRCA, CHOL, KIRP, LUAD, LUSC, PRAD, and THCA, a higher proportion of genes displayed high methylation levels compared to low methylation levels. Conversely, in BLCA and LICH, there were more genes with hypomethylation than hypermethylation. In most cancers, GZMM and GZMH exhibited hypermethylation (P < 0.05, Fig. 1D), while GZMK displayed hypomethylation (P < 0.05, Fig. 1D). The correlation analysis between methylation levels and mRNA expression revealed a negative correlation for most genes, except for GZMK, which showed a positive correlation with methylation (P < 0.05, Fig. 1E). Overall, these results suggest that mutations, CNV, and methylation of the GZM gene family may play a crucial role in cancer progression.
3.2 Differential expression and prognostic value of GZM gene families
First, we describe the genomic location of the granzyme genes. GZMA and GZMK are situated on chromosome 5, while GZMB and GZMH are located on chromosome 14, and GZMM is situated on chromosome 19 (Fig. 2A). We conducted Pearson correlation analysis and discovered significant correlations between PD-1 and GZMK (r = 0.724, P < 0.01), GZMH (r = 0.636, P < 0.01), GZMM (r = 0.591, P < 0.01), GZMA (r = 0.447, P < 0.01), and GZMB (r = 0.410, P < 0.01) (Fig. 2B). To investigate the central protein in the granzyme interactome, we performed a functional similarity analysis, which revealed that GZMK, GZMB, and GZMA may play pivotal roles, GZMK exhibits the highest critical value, which stands at 0.90(Fig. 2C). Next, we provide a detailed description of granzyme expression levels in various tumors. Generally, granzyme expression is downregulated in most cancers, except for renal chromocytosis (KICH) and brain glioma GBM, where it is highly expressed. Additionally, certain granzymes exhibit unique expression patterns in multiple cancer types that have not been well-described previously. For instance, GZMK is downregulated in most cancer types but significantly upregulated in KICH and GBM. Furthermore, we note that GZMB displays contrasting expression profiles in different tumor subtypes, suggesting diverse roles for granzymes in different cancers (Fig. 2D; Fig. 2E; Figure S1A-E). Survival analysis demonstrated that the GZM gene family primarily serves as a risk factor in LGG, GBM, THYM, and UVM, and is associated with poorer survival. On the other hand, in HNSC, CESC, and BRCA, the GZM gene family mainly acts as a protective factor and is associated with better survival (P < 0.05, Fig. 2F; Figure S1A-E). Moreover, correlation analysis of GZMA and GZMK with overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) across all carcinoma types yielded similar results, indicating that GZMA and GZMK are strongly associated with the prognosis of BRCA, LGG, SKCM, THYM, UCEC, and UVM cancers (Figure S2). Using the Cancer Genome Atlas (TCGA) data, we compared the differences in GZM scores between tumor and normal tissues (Fig. 2E). Importantly, GZM scores were significantly higher in GBM, KIRP, BRCA, and KIRC tumors, and significantly lower in LUSC, LIHC, and LUAD.
3.3 Prognosis and pathway of pan-carcinogenesis of GZMscores
GZMscores exhibit significant organ specificity with lower scores observed in brain tumors and higher scores in tumors of secretory glands such as thymus, testicular, cervical squamous cell carcinoma, and adenocarcinoma (Fig. 3A). Patients were categorized into high-risk and low-risk groups based on whether their GZMscores in the TCGA training cohort were above or below the population median. All-cancer patients with higher GZMscores exhibited poorer survival outcomes compared to the low-risk group, including disease-free survival (DFS), overall survival (OS), and progression-free survival (PFS) (Fig. 3A). To explore the cellular signaling pathways associated with GZMscores, gene set enrichment analysis (GSEA) was performed on transcriptome data from tumor groups with the top 30% and bottom 30% of GZMscores. Notably, the GZM gene family was found to play a significant role in signaling pathways pertaining to cancer, such as TNFA_SIGNALING_VIA_NFKB, INTERFERON_GAMMA_RESPONSE, INTERFERON_ALPHA_RESPONSE, INFLAMMATORY_RESPONSE, IL6_JAK_STAT3_SIGNALING, IL2_STAT5_SIGNALING, TNFA_SIGNALING_VIA_NFKB, and ALLOGRAFT_REJECTION (Fig. 3B). These results suggest that the GZM gene family is involved in regulating pathways related to the cancer immune response. Furthermore, a separate analysis revealed a positive correlation between the GZM gene family and INFLAMMATORY_RESPONSE in the entire cancer cohort (Fig. 3B). Malignant cell transformation is characterized by rapid proliferation, active epithelial-mesenchymal transformation (EMT), and angiogenesis. To investigate the relationship between GZMscores and these malignant features, the z-score algorithm was employed to quantify the tumor's contribution to angiogenesis, EMT, and cell cycle in GZMscores. The findings revealed a significant positive association between GZMscores and angiogenesis (R = 0.4, P < 0.001) as well as EMT (R = 0.45, P < 0.001) in the overall TCGA cohort or most tumor types (Figure S3). In essence, tumors with higher potential for GZM scores often exhibit angiogenesis and more aggressive tumor cells in a highly active microenvironment.
3.4 Expression distribution of pan-carcinoca single-cell atlas of GZM gene family
The expression analysis of GZM gene family in myeloid cell subsets of pan-cancer tissues showed that the expression of GZM gene family in pDC_LILRA4 cells was generally higher in BRCA, ESCA, CRC, UCEC, LYM, MEL, MYE, NPC, HCC, LUNG and THCA (Fig. 4; Figure S4; Figure S5). Similarly, in LYM, MYE, NPC, HCC, and THCA, the GZM gene family is also expressed higher in cDC1_CLEC9A cells (Figure S4; Figure S5). In addition, we found that the expression of GZMA and GZMK in myeloid cell subsets was also similar. However, there are also special cases in individual tumors, such as in BRCA, where the GZM gene family is more expressed in cDC3_LAMP3 and Marco_SPP1 cells (Fig. 4A). The expression analysis of GZM gene family in T cell subsets of pan-cancer tissues showed that GZM gene family was expressed in Tm.CD52, Tem.CXCR5, Tem.GZMK, Temra.CX3CR1, Tk.TYROBP, Tm.ZNF683, Tex.PDCD1 and other cell subsets in CD8 + T cells. In CD4 + T cells, the GZM gene family was also expressed in cell subsets such as Tm.TNF, Tm.ANXA1, Tm.ANXA2, Tm.CREM, Tm.CCLS, Tm.CAPG, Tm.GZMA, Tem.GZMK, Temra.CX3CR1, Th17.SLC4A10, Th17.IL23R, TfhTh1.CXCL13 and other cell subsets (Figure S6).
3.5 Immunoassay of GZM gene family and GZMscores in breast cancer
In this study, we assessed annotations for various factors, including immune cell infiltration, immune checkpoints, cytolytic activity, IFN response, and Siglec. To capture the true extent of immune cell infiltration in both the high and low score groups, we employed eight algorithms—CIBERSORT, MCPcounter, xCell, EPIC, Estimate, TIMER, quantiseq, and IPS—to analyze immune cell infiltration in BRCA. The results from the analysis of these eight algorithms, as depicted in Fig. 5A, demonstrated that nearly all immune cells exhibited significantly higher infiltration in the high GZMscores group. This finding implies a strong association between this score and the immune cell activity of patients. Additionally, we utilized epigenetics and mRNA mechanisms to investigate the expression of IM genes and SCNA, along with their respective expression controls. The expression of the IM genes, as presented in Fig. 5B, experiences mutation due to variants in immunoisotypes. This mutation is primarily utilized by immunoisotypes to differentiate tumors, which suggests the potential role of the GZM gene family in both the formation of the tumor microenvironment (TME) and immunotherapy. Furthermore, we analyzed the association between immune checkpoints, cytolytic activity, IFN response, Siglec markers, and scores, and discovered that all of these markers were significantly upregulated in the high GZMscores group (Fig. 5C). The correlation analysis between the GZM gene family and PD-1 in BRCA revealed a strong positive association between PD-1 and GZMA (r = 0.838, P < 0.01), GZMM (r = 0.822, P < 0.01), GZMK (r = 0.811, P < 0.01), GZMB (r = 0.800, P < 0.01), and GZMH (r = 0.761, P < 0.01) (Fig. 5D).
3.6 GZMscores-related clinical prediction models in breast cancer
The survival curves for breast cancer patients revealed a 5-year survival rate of 76%. Significant changes and fluctuations were observed during the 2–3 year, 4–5 year, and 6–8 year survival periods. In the first 5 years, the chance of survival decreased each year compared to the overall survival time. However, by year 6, the survival rate had recovered to 96% before declining again until year 9. Thus, it can be concluded that the 5-year and 8-year marks are important milestones for breast cancer patients (Fig. 6A). The survival factor association plot indicated a close correlation between the GZM gene family, GZMscores, and breast cancer prognosis in BRCA (Fig. 6B). To predict the probability of adverse events in a readable and quantifiable manner, a comprehensive network scoring curve was constructed by combining GZMscores with various clinicopathological features, such as patient age and cancer stage (Fig. 6C). The calibration curves for 3-year, 5-year, and 8-year survival closely aligned with the standard curve, which is a straight line with a slope of 1 through the origin of the coordinate axis. This alignment indicates that the network scoring curve is highly accurate in predicting the actual survival probability (Fig. 6D). Moreover, the results of the time-dependent C-index suggest that the model performs well when GZMscores are combined with clinical prognostic factors (Fig. 6E). In conclusion, incorporating GZMscores with other clinical features is critical for accurate prognostic prediction.
3.7 Single-cell distribution of GZMA, GZMK and GZMscores-related candidate genes and transcription factor regulatory network activity in breast cancer
Analysis results based on The Human Protein Atlas database reveal that in the context of BRCA, GZMA and GZMK are predominantly expressed in T cells (Figure S7). Through DepMap screening of candidate genes associated with GZMscores, we have identified a set of potential CRISPR candidate genes, namely POLR2F, MYC, SRP9, NVL, POLR3K, SNRNP25, ATP6V0B, LSM4, CCT3, GINS4, MTBP, TONSL, RCC1, CHAF1B, GINS2, CDT1, MCM2, POLE2, CENPK, GINS1, PCNA, DONSON, TIMELESS, SMC4, WDHD1, TOP2A, CDC6, KIF23, KIF11, ESPL1, CDK1, NDC80, NCAPG, CCNA2, NUF2, TPX2, PLK1, RAD51, CDC45, SPC24, CDC20, AURKB, RRM2, BIRC5, PKMYT1, DTL, SPC25, MAD2L1, AURKA, PRC1, RACGAP1, CDC7, CHEK1, ORC6, and TICRR (Figure S8A). Furthermore, the results of the transcriptional regulatory network analysis suggest that FOXM1, HIF1A, EGFR, GATA6, ESR2, RARB, AR, PGR, RXRA, STAT3, FGFR1, PPARG, TP63, KLF4, RARG, RXRB, ERBB2, ERBB3, FGFR3, ESR1, GATA3, FOXA1, and RARA might serve as potential regulators (Figure S8B).
3.8 Expression distribution of GZM gene family in TNBC single-cell atlas
The expression levels of the GZM gene family in different single-cell subsets of triple-negative breast cancer (TNBC) tissues were further analyzed. In the T cell subsets of TNBC tumor tissue, the expression of the GZM gene family was consistently higher in specific cell types, namely CD8_MAIT-KLRB1, CD8_Teff-GNLY, CD8_Tem-GZMK, CD8_Trm-ZNF683, CD8-CXCL13, and Tprf-MKI16 (Fig. 7A). Likewise, in B cell subsets, the GZM gene family showed higher expression in Bfoc-MKI67, Bmem-MKI67, and pB-IGHG1 cells (Fig. 7B). In myeloid cell subsets, the GZM gene family was generally expressed in cDC2-CD207 and pDC-LILRA4 cells (Fig. 7C). Among NK cell subsets, higher expression of GZM gene families was observed in ILC1-CD160, ILC1-CX3CR1, ILC1-FGFBP2, ILC1-GZMK, ILC1-IFNG, ILC1-IL32, ILC1-SELL, ILC1-SPON2, ILC1-VCAM1, and ILC1-ZNF683 cells (Fig. 7D). Overall, the GZM gene family was found to be commonly expressed in both T cells and NK cells in TNBC (Fig. 7).