Among the 32 cancer types, SKCM had the highest percentage of high-TMB cancers (49.4%), consistent with previous studies (14, 15). Two lung cancers had the second and third highest percentage of high-TMB cancers (36.9% and 28.1% for LUAD and LUSC, respectively). The other cancer types encompassing more than 10% of high-TMB cancers included BLCA, UCEC, COAD, STAD, DLBC, HNSC, ESCA, PAAD, and CESC ( Figure 1A ). The cancer types having less than 1% of high-TMB cancers included THYM, PRAD, KIRC, KIRP, LAML, OV, PCPG, and THCA. Since the high TMB is a biomarker for cancer immunotherapy response, the cancer types with a high percentage of high-TMB cancers would be suitable for immunotherapy. Indeed, immunotherapies, including ICB and CAR T cell therapy, have successfully treated skin, lung, bladder, head and neck, kidney, lymphoma, and MSI/dMMR cancers [10], most of which embraced more than 10% of high-TMB cancers except kidney cancer. The percentage of high-TMB cancers in pan-cancer was 11.5%, suggesting that around 11.5% of cancer patients are appropriate to immunotherapy in terms of the TMB criterion.

We examined TMB in different subtypes of a single cancer type. As expected, MSI-high tumors often had high TMB. For example, within the MSI-high subtype, 98% of COAD, 100% of ESCA, 97% of STAD, and 56% of UCEC were high-TMB tumors ( Figure 1B ). In contrast, within the microsatellite instability (MSS)/MSI-low subtype, only 6% of COAD, 17% of ESCA, 2% of STAD, and 10% of UCEC had high TMB. In BRCA, triple-negative breast cancers (TNBCs) had significantly higher TMB than non-TNBCs (one-tailed Mann-Whitney U test, P = 4.36 × 10^-11) and harbored a higher proportion of high-TMB tumors (4.39% versus 2.30%) ( Figure 1C ). In addition, basal-like and HER2-enriched breast cancers had significantly higher TMB than luminal A&B (ER+) breast cancers (P < 0.001). In COAD, BRAF-mutated tumors had significantly higher TMB than BRAF-wildtype tumors (P = 1.37 × 10^-11) and harbored a higher proportion of high-TMB tumors (70.59% versus 12.64%; P = 2.67 × 10^-17, OR = 16.09) ( Figure 1D ). In LUAD, EGFR-mutated tumors had significantly lower TMB than EGFR-wildtype tumors (P = 2.13 × 10^-5) and harbored a higher proportion of high-TMB tumors (23.08% versus 38.28%; P = 0.03, OR = 0.49) ( Figure 1D ). However, KRAS-mutated and KRAS-wildtype LUADs showed no significant difference in TMB (high-TMB tumors: 35.98% in KRAS-mutated versus 36.15% in KRAS-wildtype). In addition, we analyzed a gastrointestinal (GI) pan-cancer (17), which included 79 esophageal, 383 gastric, 341 colon, and 118 rectal cancers belonging to five subtypes: Epstein-Barr virus (EBV), MSI, hypermutated-SNV (HM-SNV), chromosomal instability (CIN), and genome stable (GS). We found that TMB followed the pattern: MSI&HM-SNV > CIN > GS&EBV (one-tailed Mann-Whitney U test, P < 0.05). The proportions of high-TMB tumors in these subtypes were 99.04%, 80%, 3.70%, 2.94%, and 0% for MSI, HM-SNV, EBV, CIN, and GS, respectively ( Figure 1E ).

We identified 376 genes whose mutations significantly correlated with increased TMB in at least 10 of the 14 cancer types (one-tailed Mann-Whitney U test, false discovery rate (FDR) < 0.1) ( Supplementary Table S1A ). These genes included some critical tumor suppressor genes (such as BRIP1, SMARCA4, PALB2, and EP300) and oncogenes (such as AFF3, CNTRL, NCOA1, NIN, NOTCH2, NSD1, NUP214, PCM1, PDE4DIP, RNF213, ROS1, and TPR) ( Figure 2A ). Notably, many of these genes belong to the same families, including ABCA (ABCA4, 12&13), ARID1 (ARID1A&B), CHD (CHD3, 4, 6, 7&9), DNAH (DNAH1, 2, 3, 5, 7, 8, 10, 11&17), FAT (FAT1, 2, 3&4), MUC (MUC5B, 6, 16&17), NOTCH (NOTCH2&4), and ZEB (ZEB1&2). GSEA (18) revealed that these genes were significantly associated with the pathways regulating cellular activity, including focal adhesion, ECM-receptor interaction, calcium signaling, tight junction, gap junction, and ABC transporters. Also, we found 28 genes whose mutations correlated with reduced TMB in a single cancer type ( Supplementary Table S1B ). For example, GATA3 mutations connected with reduced TMB in BRCA and increased TMB in five other cancer types ( Figure 2B ). IDH1 mutations were associated with reduced TMB in GBM. They increased TMB in four different cancer types ( Figure 2C ). EGFR mutations were associated with reduced TMB in LUAD and increased TMB in nine other cancer types ( Figure 2D ). The negative correlation between mutations of these genes and TMB may explain why their mutations are associated with a better prognosis in these cancer types (19–21). Overall, these results indicate that the association between mutations in a single gene and TMB is likely positive in cancer. However, a negative association between them may occur in a few cases.

Since high TMB is associated with a favorable immunotherapy response in cancer patients (13), we anticipated that the patients harboring mutations in any of the 376 genes would have better survival than those without such mutations. As expected, among the 21 genes overlapping between the 376 genes and the 474 genes studied in the pan-cancer cohort [MSKCC-Samstein cohort (13)] receiving the ICB therapy, 12 genes exhibited a significant correlation of their mutations with better overall survival (OS) in the MSKCC-Samstein cohort (log-rank test, P < 0.05) ( Figure 3A ). The 12 genes included ARID1A, ARID1B, BRIP1, NOTCH2, NOTCH4, EPHA5, ROS1, FAT1, SPEN, NSD1, PTPRT, and ZFHX3. Among them, only the mutation of ZFHX3 was associated with better OS in another pan-cancer cohort [MSKCC-Zehir cohort (13)] not receiving the ICB therapy (P = 0.007) ( Supplementary Figure 1 ). It indicated that 11 of the 12 genes (except ZFHX3) had the characteristic of their mutations correlated with both high TMB and a favorable response to immunotherapy in cancer. Next, we used logistic models to explore the ability of mutations in the 11 genes to predict high- versus low-TMB in TCGA pan-cancer and the 14 individual cancer types. We first utilized the least absolute shrinkage and selection operator (Lasso) (22) to select variables (genes) and refitted logistic regression models with the variables selected by the Lasso. The Lasso selected all 11 genes in the pan-cancer analysis, and in analyzing individual cancer types, each gene was selected by the Lasso at least eight times ( Figure 3B ). The prediction accuracy was 92% in both training and test sets in the pan-cancer analysis and more than 80% in 10 individual cancer types. Furthermore, we analyzed the correlations between mutations of the 11 genes and OS in the MSKCC-Samstein cohort using the Cox proportional hazards model. We found three genes (ROS1, SPEN, and PTPRT) whose mutations had a significant association with better OS (ROS1: P = 0.019, HR = 0.67, 95% CI: [0.48, 0.94]; SPEN: P = 0.023, HR = 0.64, 95% CI: [0.44, 0.94]; PTPRT: P = 0.022, HR = 0.74, 95% CI: [0.57, 0.96]) ( Figure 3C ). We reanalyzed the correlations between mutations of the three genes and OS in the MSKCC-Samstein cohort using the Cox proportional hazards model and obtained their coefficients. Based on the mutation status of the three genes, we defined the TMB prognostic score (TMBPS) as: TMBPS = - (a × mutated (ROS1) + b × mutated (SPEN) + c × mutated (PTPRT)), where a = -0.5071, b = -0.5189, and c = -0.3485, which were the coefficients of the three variables in the Cox proportional hazards model; mutated(X) = 1 if the gene X is mutated in the tumor sample, otherwise mutated(X) = 0. We found that high-TMBPS (> median) tumors showed better OS than low-TMBPS tumors in the MSKCC-Samstein cohort (p = 1.85 × 10^-9) ( Figure 3D ), but high- and low-TMBPS tumors showed no significant different OS in the MSKCC-Zehir cohort without the ICB therapy ( Figure 3D ). We further confirmed that the TMBPS correlated positively with OS in two lung cancer cohorts [Rizvi cohort (10) and Hellmann cohort (23)] and two melanoma cohort [Hugo cohort (11) and Allen cohort (12)] receiving the ICB therapy ( Figure 3E ). In contrast, it showed no significant correlation with OS in the TCGA lung cancers ( Figure 3F ). These results demonstrate that the TMBPS can predict the response to ICB therapy.

We identified nine cancer-associated pathways whose activity exhibited a significant positive correlation with TMB in at least eight cancer types (Spearman’s correlation test, FDR < 0.1). These pathways included nucleotide excision repair, DNA replication, homologous recombination, base excision repair, mismatch repair, cell cycle, spliceosome, proteasome, and RNA degradation ( Figure 4 ). Notably, numerous DNA damage repair-associated pathways showed significant positive correlations with TMB ( Figure 4 ). This is consistent with previous studies showing that dMMR has a strong correlation with high TMB in cancer (15). Interestingly, the cell cycle pathway showed a marked positive correlation with TMB in nine cancer types, including (BLCA, BRCA, COAD, LUAD, LUSC, SARC, SKCM, STAD, UCEC) ( Figure 4 ), indicating that the proliferation of tumor cells increases TMB. The spliceosome pathway exhibited a significant correlation with TMB in eight cancer types ( Figure 4 ); a possible explanation is that the spliceosome is associated with DNA damage repair and cell cycle progression (24). Moreover, we found seven pathways showing significant inverse correlations with TMB in at least eight cancer types (Spearman’s correlation test, FDR<0.1) ( Figure 4 ). These pathways included Wnt, Hedgehog, PI3K-AKT, MAPK, neurotrophin, axon guidance, and pathways in cancer. Most of these pathways are overexpressed in cancer (25–29). Our data indicated that the hyperactivation of these pathways was associated with reduced TMB in cancer, suggesting that their activities inhibit the antitumor immune response. Interestingly, two neural development-associated pathways (neurotrophin and axon guidance) correlated significantly with TMB in diverse cancers. It suggests that the molecules involved in regulating neural development may play essential roles in modulating tumor immunity, as evidenced in a recent study (30).

We identified ten proteins having significant positive expression correlations with TMB in at least five cancer types (Spearman’s correlation test, FDR < 0.1) ( Figure 5 ). These proteins included ASNS, Cyclin_B1, Cyclin_E1, Cyclin_E2, Caspase-7, FoxM1, PCNA, Rad51, TFRC, and 4E-BP1. Many of these proteins were associated with cell cycle regulation, including ASNS, Cyclin_B1, Cyclin_E1, Cyclin_E2, and FoxM1. This is consistent with the fact that the cell cycle pathway correlates positively with TMB in multiple cancer types ( Figure 4 ). Besides, several proteins were involved in the regulation of DNA damage repair, including PCNA and RAD51, consistent with the significant correlation between dMMR and high TMB in cancer. Furthermore, we identified seven proteins with significant inverse expression correlations with TMB in at least five cancer types (Spearman’s correlation test, FDR < 0.1) ( Figure 5 ). These proteins included ACVRL1, Bcl-2, Caveolin-1, c-Kit, PKC-alpha, PR, and INPP4B. Among them, ACVRL1 is associated with the TGF-β pathway, which hampers antitumor immune response (31); Bcl-2 regulates the apoptosis pathway, which correlated inversely with TMB in four cancer types ( Figure 4 ).

We identified 356 genes from the immune response module (pink module) which enriched in higher-TMB cancers and 674 genes from the cell cycle/chromosome module (turquoise module) which enriched in lower-TMB cancers among pan-cancer (WGCNA-pink-turquoise genes) by WGCNA method ( Figure 6A ). Meanwhile, we identified 730 and 248 mRNAs upregulated and downregulated differentially expressed mRNA genes (DEmRNA) between high-TMB cancers and low-TMB cancers in at least five cancer types, respectively (Student’s t test, FDR < 0.05, |FC| > 1.5). WGCNA-DEmRNA were defined as the intersection of WGCNA-pink-turquoise genes and DEmRNA. ClueGO mainly enriched the 107 WGCNA-DEmRNA in the mitotic nuclear division, regulation of cell division, and positive regulation of cell cycle phase transition ( Figure 6B ). Furthermore, we identified 58 (or 18) miRNAs and 8 (or 16) lncRNAs upregulated and downregulated differentially expressed miRNAs (DEmiRNA) and lncRNAs (DElncRNA) between high-TMB cancers and low-TMB cancers in at least five cancer types, respectively (Student’s t test, FDR<0.1) ( Figures 6C, D ). Many of these miRNAs were involved in the KEGG pathway MicroRNAs in cancer pathway (https://www.genome.jp/kegg-bin/show_pathway?hsa05206), including hsa-mir-192, hsa-mir-7-2, hsa-mir-96, hsa-mir-25, hsa-mir-7-3, hsa-mir-107, hsa-mir-128-2, hsa-mir-128-1, hsa-mir-183, hsa-mir-200b, hsa-mir-615, and hsa-mir-335. In addition, hsa-mir-25, hsa-mir-3074, and hsa-mir-106b were involved in the DNA damage response pathways. To further explore the molecular features of TMB, a competing endogenous RNAs (ceRNA) network was constructed and visualized by Cytoscape based on the intersection of 24 DElncRNA–DEmiRNA pairs and 383 WGCNA-DEmRNA–DEmiRNA pairs that were identified from the three databases [miRecords (32), miRTarBase (33, 34), starBase (35)]. ( Figure 6E ). The ceRNA network contains 7 DElncRNAs, 20 DEmiRNAs, and 107 WGCNA-DEmRNA.

It has been recognized that high-TMB cancers are likely to yield more neoantigens to incite antitumor immune responses (7). We found that the enrichment levels of CD8+ T cells were higher in high-TMB than in low-TMB cancers in seven cancer types (COAD, CESC, UCEC, BLCA, LUSC, STAD, LUAD). The immune cytolytic activity was higher in high-TMB cancers in eight cancer types (COAD, CESC, STAD, BRCA, BLCA, UCEC, LUSC, LUAD) (one-tailed Mann-Whitney U test, FDR < 0.1) ( Figure 7A ). Meanwhile, in 10 cancer types (COAD, STAD, LUAD, BLCA, SKCM, CESC, BRCA, UCEC, SARC, PAAD), high-TMB cancers displayed significantly higher PD-L1 expression levels than low-TMB cancers (two-tailed Student’s t test, P < 0.1). It indicates that elevated TMB is associated with increases in both immunostimulatory and immunoinhibitory signatures. However, we found that the ratios of immunostimulatory/immunoinhibitory signatures (CD8+/CD4+ regulatory T cells, pro-/anti-inflammatory cytokines, and M1/M2 macrophages) were significantly higher in high-TMB than in low-TMB cancers (two-tailed Student’s t test, FDR < 0.1) in diverse cancer types ( Figure 7B ). The ratios were the average expression levels of immune-stimulatory signature marker genes over the average expression levels of immune-inhibitory signature marker genes. Furthermore, we found that high-TMB cancers encompassed more neoantigen loads than low-TMB cancers in all ten cancer types, in which the data of predicted neoantigen loads are available ( Figure 7C ) (36). These results confirmed that high TMB is associated with a more active antitumor immune response.

In seven of the 14 TCGA cancer types, high-TMB cancers displayed significantly higher expression levels of MKI67, a marker for tumor cell proliferation (37), and proliferation signature scores than low-TMB cancers (two-tailed Student’s t test, P < 0.1) ( Figure 8A ). It indicates that high TMB is associated with the increased proliferation potential of tumor cells. The association between TMB and survival prognosis is cancer type dependent. In LIHC and HNSC, high-TMB tumors had significantly worse survival (OS and/or disease-free survival (DFS)) prognosis than low-TMB tumors, while in BLCA, STAD, and UCEC, high-TMB tumors had significantly better survival prognosis (log-rank test, P < 0.1) ( Figure 8B and Supplementary Figure 2 ).

We analyzed associations of TMB with age, gender, height, weight, smoking, and race. We found that the high-TMB patients were significantly older than the low-TMB patients in pan-cancer and in six individual cancer types, including BRCA, CESC, ESCA, HNSC, SKCM, and STAD (one-tailed Mann-Whitney U test, P < 0.1) ( Figure 8C ). This is reasonable because tumor mutations accumulate with age. One exception was LUAD, in which the high-TMB patients were significantly younger than the low-TMB patients (P < 0.001) ( Figure 8C ). The main reason could be that smoking significantly contributes to gene mutations during the development of lung cancer. Indeed, the smoking load had a significant positive correlation n with TMB in LUAD (P < 0.001) ( Figure 8D ). In pan-cancer and three individual cancer types, including BLCA, SARC, and SKCM, males had significantly higher proportions of high-TMB tumors than females (Fisher’s exact test, P < 0.1) ( Figure 8E ). However, in STAD, females had significantly higher proportions of high-TMB tumors than males (Fisher’s exact test, P < 0.002). In pan-cancer and four individual cancer types, including BLCA, ESCA, LIHC, and SKCM, the high-TMB patients were significantly higher than the low-TMB patients (one-tailed Mann-Whitney U test, P < 0.05) ( Figure 8F ). In addition, in pan-cancer and three individual cancer types, including BLCA, ESCA, and SKCM, the high-TMB patients were significantly heavier than the low-TMB patients (one-tailed Mann-Whitney U test, P < 0.01) ( Figure 8F ). Interestingly, in certain cancer types, TMB showed a significant association with race. In SKCM, ESCA, and BLCA, the proportion of high-TMB tumors in the white was significantly higher than that in other races (Chi-square test, P < 0.1) ( Figure 8G ). In contrast, in LUAD, the white harbored a significantly lower proportion of high-TMB tumors than other races (P = 0.01). In LIHC, the American Indian had the highest proportion (50%) of high-TMB tumors compared to other races (less than 10%).

We downloaded gene somatic mutations (level 3), RNA-Seq gene expression profiles (level 3, RSEM-normalized), protein expression profiles (level 3), and clinical datasets for 32 cancer types from the TCGA data portal (https://portal.gdc.cancer.gov/). The cancer-associated pathways and genes involved in these pathways were obtained from KEGG (51). The somatic mutation profiles and clinical datasets for two other pan-cancer cohorts (MSKCC-Samstein and MSKCC-Zehir cohorts) were attained from their associated publication (13). In addition, we obtained the somatic mutation profiles and clinical datasets for two lung cancer cohorts (Rizvi cohort (10) and Hellmann cohort (23)) and two melanoma cohorts (Hugo cohort (11) and Allen cohort (12)) receiving ICB therapy from their associated publications.( Supplementary Table 2 )

For each tumor sample, its TMB was defined as the total count of non-synonymous somatic mutations/30Mb in the sample. We defined the high- and low-TMB according to the FDA standard. That is, if a tumor sample’s TMB is no less than 10, the tumor sample has a high-TMB. Otherwise, it has a low-TMB. The sample size for high- and low-TMB cancers is presented in Supplementary Table S3 .

For a pathway, immune signature, or biological process, we quantified its enrichment level in a tumor sample as the single-sample gene-set enrichment analysis (ssGSEA) (52) score of its marker genes. The marker gene sets representing different immune signatures were from several publications, including CD8+ T cells (36), CD4+ regulatory T cells (36), pro-/anti-inflammatory cytokines (36), and M1/M2 macrophages immune cytolytic activity (36), and proliferation signature score (53). These gene sets are listed in Supplementary Table S4 .

We explored the ability of gene mutations to predict high- versus low-TMB cancers in TCGA pan-cancer and in 14 individual cancer types. We separated all samples into a training set (70% of the samples) and a test set (30% of the samples). On the training set, we used the Lasso to select variables with 10-fold cross-validation (CV) and refitted a logistic regression model with the variables selected by the Lasso using a 10-fold CV. We then tested the model on the test set and calculated the area under the curve of the receiver-operating characteristic (ROC-AUC) as the prediction performance. We performed the Lasso and logistic regression analyses using the R package “glmnet” and “caret,” respectively, and calculated the standardized regression coefficients (β values) using the function “lm. beta” in the R package “QuantPsyc.”

The correlation between miRNA or lncRNA expression levels and TMB was assessed by comparing miRNA or lncRNA expression levels between high-TMB and low-TMB cancers in each caner type using Student’s t test. We identified miRNAs and lncRNAs whose expression levels were significantly higher (or lower) in high-TMB cancers than in low-TMB cancers in at least five cancer types as DEmiRNAs and DElncRNAs. We predicted DElncRNA–DEmiRNA pairs by the intersection of three databases [miRecords (32), miRTarBase (33, 34), starBase (35)] and using the R package “multiMiR”. WGCNA–DEmRNA–DEmiRNA pairs were identified by using starBase (http://starbase.sysu.edu.cn/) database. We establish and visualize the DElncRNA–DEmiRNA–WGCNA–DEmRNA ceRNA network based on the intersection of DElncRNA–DEmiRNA pairs and WGCNA–DEmRNA–DEmiRNA pairs (established by DElncRNAs, DEmiRNAs, and WGCNA–DEmRNA) by Cytoscape v3.7.2 WGCNA analysis was carried out using R package "WGCNA".

We using the first 70% mRNAs by variance comparison to build the gene modules that were differentially enriched between the higher TMB level and the lower TMB level tumors in pan-cancer. Then, we identified 14 co-expressed gene modules with softpower 7. WGCNA-pink-turquoise genes was defined as genes which enriched in immune response and cell cycle modules (pink and turquoise modules). The interaction of WGCNA-pink-turquoise genes and DEmRNA was defined as WGCNA – DEmRNA. Cytoscape plug-in (ClueGO) (54) provides functional annotation for WGCNA – DEmRNA genes.

In comparisons of two classes of data, we used the Mann-Whitney U test (for non-normally distributed data) or the Student’s t test (for normally distributed data). We used Spearman’s correlation test (for non-normally distributed data) or Pearson’s correlation test (for normally-distributed data) to evaluate the correlation between two groups of data. We used Fisher’s exact or Chi-square test to assess the association between two categorical variables. We utilized the FDR, evaluated by the Benjamini-Hochberg method, to adjust for multiple tests. We performed all statistical analyses in the R programming environment (version 4.0.2).

In the univariate survival analysis, we utilized the Kaplan-Meier curves to exhibit the survival time differences and the log-rank test to evaluate the significance of survival time differences. We performed the survival analyses using the function “survfit” in the R package “survival.” In the multivariate survival analysis, we used the cox proportional hazards model. We performed the multivariate survival analyses using the function “coxph” in R package “survival.”