Pan-cancer γδ TCR analysis uncovers clonotype diversity and prognostic potential

Abstract

Gamma-delta T cells (γδ T cells) play a crucial role in both innate and adaptive immunity within tumors, yet their presence and prognostic value in cancer remain underexplored. This study presents a large-scale analysis of γδ T cell receptor (γδ TCR) reads from 11,000 tumor samples spanning 33 cancer types, utilizing the TRUST4 algorithm. Our findings reveal extensive diversity in γδ TCR clonality and gene expression, underscoring the potential of γδ T cells as prognostic biomarkers in various cancers. We further demonstrate the utility of TCR gamma (TRG) and delta (TRD) gene expression from standard RNA-sequencing (RNA-seq) data. This comprehensive dataset offers a valuable resource for advancing γδ T cell research, with implications for enhanced immunotherapy approaches or alternative therapeutic strategies. Additionally, our centralized database supports translational research into the therapeutic significance of γδ T cells.

Keywords: T cell receptor; TCR clonality; TRUST4; cancer immunotherapy; gamma-delta T cells; immune repertoire; prognostic biomarkers.

Graphical abstract

Figure 1

Comprehensive pan-cancer analysis of γδ TCRs 660 million RNA-seq reads from a total of 10,970 tumor samples in the TCGA dataset were reanalyzed for γδ TCRs using TRUST4. It led to the identification of 34,129 unique γδ TCR clones of 6,751 tumors. γδ-TCRs were further filtered to remove (1) samples without clinical information, (2) metastasis tumors for patients with multiple tumors available, (3) out-of-frame CDR3 sequences as well as CDR3 sequences carrying alpha-beta V genes. Similarity, approximately 40 million RNA-seq reads from 503 tumors downloaded from 8 immunotherapy studies were reanalyzed, generating 6,770 unique γδ TCR clones, of which 3,574 clones from 459 tumors were retained for downstream analysis. Comprehensive pan-cancer and cancer/study-specific γδ TCR repertoire analysis was conducted using the filtered γδ TCR sequences. These refined datasets, representing a significant resource in γδ TCR research, were made available to the research community, facilitating broader scientific exploration. See also Table S1.

Figure 2

TCR gamma chain expression and clonotype distribution across TCGA cancer types (A) Variability in TCR gamma gene enrichment and clonotype presence across 33 distinct TCGA cancer types. Top: gamma gene enrichment scores of tumors. Each dot represents a tumor, with y axis position indicating the normalized gamma gene enrichment score. Cancer types are ordered by their median enrichment scores, with the lowest scores on left and highest on right. The second panel shows the number of patients in each cancer types, with light green indicating patients without any γ TCR identified and dark green indicating patients with at least one γ TCR clone identified. The third panel shows number of unique γ TCR clones identified in each cancer type, with light purple representing clones with only 1 read and dark purple representing clones with >1 reads. Bottom: the proportion of unique γ TCR clones with specific clone sizes: = 1 read, = 2 reads, 3–10 reads, 11–100 reads, and >100 reads. (B) Correlation between TCR gamma gene enrichment scores and normalized gamma chain counts across cancer types. Each dot represents a cancer type with y axis indicating the mean normalized gamma chain counts and x axis indicating mean gamma gene enrichment scores across tumors. (C) Spearman correlation between gamma gene enrichment scores and normalized gamma chain counts within each cancer type, with correlation coefficients shown in different gradients of red. See also Figures S1–S3.

Figure 3

Prognostic implications of γδ TCRs in TCGA The forest plot displays the results of the Cox regression survival analysis on TCR γ and δ gene expression adjusted by covariates across various cancer types. It highlights genes with significant or marginally significant impacts on patient survival, indicated by FDR-adjusted Cox regression p values below 0.05. Each horizontal line in the forest plot represents a γ or δ gene in a cancer type, with the square and length of the line indicating the estimated hazard ratio (HR) and its 95% confidence interval, respectively. Horizontal lines are colored by cancer types. Within each cancer type, genes are ranked by HR from highest to lowest. A dashed vertical line at HR = 1 represents no effect.

Figure 4

γδ TCR-based prognostic signatures in HPV-positive vs. HPV-negative HNSC tumors in TCGA (A) Elevated TCR gamma and delta gene expression in patients with HPV-positive HNSC. Y axis represents the purity-adjusted gamma (left) and delta (right) gene enrichment scores. (B) Higher normalized entropy of TCR gamma signature in the HPV-positive HNSC. Forest plots show the estimated hazard ratio (square) and its 95% confidence interval (horizontal line) of gamma and delta gene enrichment by Cox regression adjusted for covariates in HPV+ and HPV− HNSC, respectively. A dashed vertical line at HR = 1 represents no effect. (C) Association of gamma and delta TCR signatures with improved survival in HPV-positive HNSC, absent in HPV-negative. p value was calculated by two-sided Wilcoxon rank-sum test. Violin plot shows the kernel probability density of the data with box in middle representing the median (central line), the 25% and 75% interquartile (IQR) (lower and upper hinges). See also Figures S4A–S4D.

Figure 5

γδ TCR-based prognostic signatures in TCGA COAD tumors (A) Significantly higher gamma and delta TCR enrichment scores in the MSI-high group. (B) Higher normalized entropy of the gamma TCR in the MSI-high group. (C) The gamma enrichment score associated with worse survival outcomes in the MSI-low group. Forest plots show the estimated hazard ratio (square) and its 95% confidence interval (horizontal line) of gamma gene enrichment by Cox regression adjusted for covariates in MSI-high and MSI-low COAD tumors, respectively. A dashed vertical line at HR = 1 represents no effect. For (A) and (B), p value was calculated by two-sided Wilcoxon rank-sum test. Violin plot shows the kernel probability density of the data with box in middle representing the median (central line), the 25% and 75% interquartile (IQR) (lower and upper hinges). See also Figures S4E–S4G.

Figure 6

γδ TCR-based prognostic signatures in immunotherapy studies (A) Prognostics γδ TCR signatures in the study by Gide. Violin plots on left compare gamma and delta gene enrichment between patients with complete response (CR) or partial response (PR) vs. patients with stable disease (SD) and progressive disease (PD), in pre-treatment and on-treatment samples, respectively. Number of samples in each category was labeled below x axis. Right: the (1) Kaplan-Meier overall survival plot by patients’ pre-treatment samples stratified based on optimal cutoff of gamma gene enrichment, with yellow representing high and blue representing low. (2) Forest plot of PFS (progression-free survival) HR (square) and its 95% confidence interval (horizontal line) of gamma gene enrichment by Cox regression adjusted for age and gender. A dashed vertical line at HR = 1 represents no effect. p values are shown on right. (B) Prognostics γδ TCR signatures in the study by Riaz. Left: violin plot with box shows elevated gamma and delta gene enrichment in patients with CP/PD vs. patients with SD/PD, in both pre- and on-treatment tumors. Right: violin plot with box shows elevated gamma and delta gene enrichment in CP/PD vs. SD/PD, in both pre- and on-treatment tumors within ipilimumab treatment-progressive (Ipi Prog) patients. (C) In the study by Riaz, Kaplan-Meier overall survival plot by patients’ pre-treatment samples shows that higher gamma gene enrichment is associated with better survival. Patients are stratified based on optimal cutoff of gamma gene enrichment. (D) Kaplan-Meier overall survival plot by patients’ pre-treatment samples stratified based on optimal cutoff of gamma true diversity in the study by Gide (left) and gamma normalized entropy in the study by Riaz (right), with orange representing high and purple representing low. For violin plots in A–C, p value was calculated by two-sided Wilcoxon rank-sum test. Violin plot shows the kernel probability density of the data with box in middle representing the median (central line), the 25% and 75% interquartile (IQR) (lower and upper hinges), and the ±1.5 IQR (Tukey whiskers). See also Figure S5.