Assessment of tumor mutation burden calculation from gene panel sequencing data

Introduction

Immunotherapy has profoundly changed cancer treatment, providing durable benefits across a wide range of cancer types once considered untreatable.^1–7 To date, more than 1,600 clinical trials on PD-1/PD-L1 blockade have been registered in clinicaltrials.gov. Patients with advanced non-small cell lung cancer (NSCLC) treated with nivolumab (Opdivo), an anti-PD-1 monoclonal antibody checkpoint inhibitor, achieved a five-year overall survival (OS) rate of 16%, comparing with an OS of 4% with standard-of-care chemotherapy.⁸ However, patient response rates remain relatively low at 15–40%,⁹ indicating a need for effective biomarkers to identify patients who are likely to respond to PD-1/PD-L1 blockade.

Many studies have attempted to identify reliable biomarkers in patients most likely to benefit from immune checkpoint blockade. Immunohistochemistry (IHC) for PD-L1 protein expression was the first diagnostic biomarker approved by the FDA for PD-1/PD-L1 treatment. However, the high false-negative rate and lack of standardized procedures mean that this assay is still far from being an ideal biomarker for immunotherapy response prediction.^10,11 In recent years, some progress has been made in the discovery of multiple genome stability-related biomarkers that correlate highly with the response of immunotherapy in multiple cancer types.^{1,6,7,9,12–14} These include tumor mutation burden (TMB) or tumor mutation load, neoantigen burden, DNA mismatch repair deficiency, and high microsatellite instability. Of these, TMB is most robust, effective and clinical verifiable biomarker in multiple cancer types.^9,12,15,16 Previous studies in small-cell lung cancer and NSCLC found a remarkable association between TMB values and the production of neo-antigen peptides, thought to be highly correlated with tumor immune response.^17,18 Moreover, the clinical trial CheckMate-032 study showed that TMB level has strong distinguishing power in both Nivolumab single agent treatment and Nivolumab plus Ipilimumab combination therapy.^19,20

TMB of a tumor sample is calculated by the number of non-synonymous somatic mutations (single nucleotide variants and small insertions/deletions) per mega-base in coding regions.²¹ TMB represents both the stability level of the tumor genome and heterogeneity of the tumor micro-environment. However, methods to calculate TMB differ greatly in scientific research and clinical practice. Usually, whole genome sequencing, WES and panel sequencing are used to calculate TMB, with various approaches to filter somatic mutations. In research, the most common method is WES, which sequences only the coding regions of the genome, with average coverage of 120X to 150X. However, in clinical application, some studies use panel sequencing to estimate TMB.²² Although only a small number of genes (about 300–500) are covered, gene panel sequences can identify rare somatic mutations because of the higher sequencing depth compared to WES. For example, the FDA recently approved two cancer-related genes panels, FoundationOne CDx (F1CDx) and MSK-IMPACT.^11,21 F1CDx adopts a tumor-only sequencing strategy while MSK-IMPACT requires sequencing of both the tumor and its matched normal sample. In F1CDx, synonymous mutations are counted while hotspot driver mutations are excluded when calculating TMB. With MSK-IMPACT, TMB is calculated with similar filtering criteria to those used in WES. In addition to these two panels, there are also other large gene panels from different institutes that approximate TMB, with the size of the panel another important factor influencing TMB approximation.

Further efforts are still needed for more accurate TMB approximation from targeted sequencing data of cancer-related gene panels. In this work, we attempted to investigate the TMB distributions in WES level across different cancer types using WES data from the TCGA data-set. We studied frequently mutated genes and hotspot mutations and their associations with TMB. We also examined the correlation between TMB approximation from F1CDX, MSK-IMPACT, their union panel (F1CDX+MSK) and WES. Furthermore, based on hotspot-mutated gene sets in different cancers, we proposed an optimized pan-cancer panel for better TMB estimation. This panel could serve as a reference data-set for gene selection in TMB-oriented panels designed for immunotherapy and other clinical applications, such as targeted therapy for solid tumors.

Material and methods

Test data-set

VCF files of somatic mutations from WES sequencing data called by Mutect2 for 15 different cancer types were downloaded from the TCGA database (https://portal.gdc.cancer.gov/projects). The data-set contained a total of 6,625 samples, of which 567 samples are lung adenocarcinoma (LUAD, 567), lung squamous cell carcinoma (LUSC, 492), bladder urothelial carcinoma (BLCA, 412), colon adenocarcinoma (COAD, 399), breast invasive carcinoma (BRCA, 986), skin cutaneous melanoma (SKCM, 467), stomach adenocarcinoma (STAD, 438), kidney renal papillary cell carcinoma (KIRC, 281), liver hepatocellular carcinoma (LIHC, 364), glioblastoma multiforme (GBM, 393), uterine corpus endometrial carcinoma (UCEC, 530), ovarian serous cystadenocarcinoma (OV, 530), prostate adenocarcinoma (PRAD, 498), esophageal carcinoma (ESCA, 184) and pancreatic adenocarcinoma (PAAD, 178). For this analysis, we used these somatic mutations as raw data to evaluate the correlation between TMB estimated by cancer-related genes panels FoundationOne CDx (F1CDx, 315 genes), MSK-IMPACT (468 genes), union set of F1CDx and MSK-IMPACT (F1CDx+MSK, 525 genes) and WES. The list of oncogenes and suppressor genes was obtained from the Oncokb database (http://oncokb.org/) and Cosmic (https://cancer.sanger.ac.uk/cosmic), which are manually curated.

TMB calculation

TMB was defined as the number of somatic mutations in the coding region per megabase, including single nucleotide variants (SNVs) and small INDELs (insertions and deletions, usually less than 20 bases). However, the means to determine reliable somatic coding mutations for TMB calculation was not trivial. In the F1CDx approach, synonymous mutations were included. However, stop-gain mutations in tumor suppressor genes and hotspot driver mutations were not included in order to reduce bias due to enrichment of cancer-related genes in the F1CDx panel. Here, we followed the F1CDx approach to calculate the TMB for our panel, defining the cutoff values as TMB-high (≥20 mutations/Mb), TMB-medium (<20 mutations/Mb ≥10 mutations/Mb) and TMB-low (<10 mutations/Mb). Others filtering parameters of somatic mutations identified by WES for TMB calculation were a mutated allele frequency greater than 5% and a sequence depth greater than 20X in tumor samples greater or 10X in normal samples.

Pan-cancer TMB panel

Genes with a high mutation frequency in cancer patients were selected for a pan-cancer panel to estimate TMB. We calculated the proportion of patient samples from 15 different cancers with mutations in each gene. For each cancer, genes were sorted in descending order of the previous proportion and the cumulative value of the proportion was calculated. Based on a cumulative value of 30%, we selected genes from all 15 cancers and added genes appearing in more than 10 cancers to the pan-cancer TMB panel. In total, 524 genes were identified and included in the panel.

Improved correlation analysis

As TMB calculation from panels and WES both use somatic mutations from TCGA samples, a simple correlation analysis might cause “over-fitting”. To avoid this problem, we further trained a linear regression model on a set (approximately 80%) of randomly selected samples from TCGA for every cancer. The remaining samples were used as the test data-set. During the training process, we fit linear regression models as WES TMB = a + b * Panel TMB for 15 cancers. The purpose of these models was to calculate the predicted TMB for samples in the test data-set. In this way, we obtained the correlation between predicted TMB and real TMB. By comparing the correlation between panels and WES, performance of different panels and the pan-cancer TMB panel was evaluated.

Results

TMB distribution in WES level

We first examined the distribution of TMB at the WES level for each cancer type as the somatic mutations were called from WES data in the TCGA data-set. Based on the WES TMB distribution results shown in Figure 1A and B, we found that TMB distributions were very different in different cancer types. For example, SKCM had the highest TMB level, with the third quartile reaching 20 mutation/Mb (the threshold for TMB-high samples in the F1CDx study²³). Using this threshold, 530 of 6,625 (7.2%) samples were identified as TMB-high. In BLCA, COAD, LUAD, LUSC, SKCM, STAD and UCEC, the percentage of TMB-high samples ranged from 5–30%. However, the percentage of TMB-high samples in the other eight cancer types (BRCA, KIRC, LIHC, GBM, OV, PRAD, USCA and PAAD) were very low (0–4%). Though KIRC and LIHC had low TMB distributions, a relatively high objective response rate to immunotherapy for these two cancers were still observed.⁹ These results indicated that different cancer types have different TMB distribution level, and that a simple cutoff value for TMB-high samples (defined here as 20 mutations/Mb) is not suitable in predict efficacy of immunotherapy in every solid tumor type. For cancers with low TMB distributions, more clinical trials are needed to confirm TMB as a good marker, and a suitable TMB cutoff must be determined using a large number of patient samples with known corresponding clinical response results to immunotherapy.

Figure 1 Tumor mutation burden (TMB) distribution in whole exome sequencing (WES) level. (A) The histogram plots show the frequency of TMB high, medium and low in WES level (B) Violin plots show the distribution of the TMB for 15 tumor types in WES level. Abbreviations: LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma ; BLCA, bladder urothelial carcinoma; COAD, colon adenocarcinoma; BRCA, breast invasive carcer; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; KIRC, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; GBM, glioblastoma multiforme; UCEC, uterine corpus endometrial carcinoma; OV, ovarian serous cystadenocarcinoma; PRAD, prostate adenocarcinoma; ESCA, esophageal carcinoma; PAAD, pancreatic adenocarcinoma.

Top mutated genes in WES level

The top 20 most commonly mutated genes in each cancer were calculated and shown in Figure S1. The gene USH2A, TTN, TP53, SYNE1, RYR2, MUC16, CSMD3 showed highly mutation frequencies in at least ten cancer types. The APC gene had the highest mutation frequencies (270/399) in COAD and BRAF was most frequently mutated (228/467) in the SKCM than the other cancer type. Some well-established cancer-associated genes, such as PI3KCA and TP53, were frequently mutated in most cancer types. For example, the mutation ratio for PI3KCA in seven cancer types ranged from 2.4–42%; for TP53, 9.2–64.0%. For PTEN, the mutation ratio in the same seven cancer types was relatively low for most cancer types, starting from 1.25%, but reaching 50.9% in UCEC. Moreover, the top mutated genes in different TMB level data-set (TMB-high, TMB-medium and TMB-low) were interesting as well. However, the number of samples with high TMB in BRCA(29/986), KIRC(1/281), LIHC(10/364), GBM(28/393), OV(20/530), PRAD(4/498), ESCA(7/184) and PAAD(2/178) were relatively low (Figure 1). Thus, we only calculated the top 20 mutated genes in the other seven cancer types (BLCA, COAD, LUAD, LUSC, SKCM, STAD and UCEC) from three data-sets: TMB-high, TMB-medium and TMB-low samples. Comparing the results in these three data-sets ( Figure S1), we found that the top genes in these three data-sets were different for every cancer and the impact of these genes on TMB approximation was different. Figure 2 shows the number of samples with mutations in nine representative genes in LUAD (Figure 2A) and SKCM (Figure 2B). TP53, EGFR, PIC3CA and KRAS were often mutated in TMB-low samples in LUAD. Conversely, these genes were mutated in TMB-high samples in SKCM. Interestingly, TTN and MUC16 also appeared in the top ten mutated genes for six cancers with various TMB values (Figure S1). TTN, which encodes for titin, is a tumorigenic gene with a coding length of 80781 bp, likely making it susceptible high mutation rates because of random DNA repair errors.²⁴ Therefore, the TTN gene is frequently mutated both in TMB-high samples and TMB-low samples. MUC16 (coding length, 43524 bp) encodes a transmembrane glycoprotein with a molecular weight of 2000 kDa, meaning that it also has a high risk of mutation.²⁵

Figure 2 The mutated frequency of 9 genes in WES level for data-sets of all samples, TMB-high samples, TMB-medium samples and TMB-low samples. (A) LUAD. (B) SKCM.Abbreviations: LUAD, lung adenocarcinoma; SKCM, skin cutaneous melanoma; TMB, tumor mutation burden; WES, whole exome sequencing.

Association of hotspot mutations with TMB

TMB is defined as the number of somatic coding mutations per million bases. As the coding length of WES data is generally about 33 million bases, many mutations from different genes contribute to the calculation of TMB from WES. As shown above, some genes were frequently mutated in particular cancer types. However, the relationship of those hotspot mutations with TMB are still not well-understood. Here, we investigated whether hotspot somatic mutations might be associated with TMB in specific cancers.

Firstly, we defined hotspot mutations as highly frequent mutations found in at least ten samples in the complete set of cancer WES data (Figure S2). Next, we focused on those hotspot mutations which occurred in at least three samples in one cancer type. Finally, a total of 150 unique mutations were investigated for their association to TMB level (327 cancer type-specific mutations). For these mutations, we calculated the P-value using a Mann-Whitney U-test and the fold-change by the median TMB value in mutation-positive samples (samples containing the mutation) and mutation-negative samples (samples without the mutation), as listed in Figure S3. Figure 3A shows that most of these hotspot mutations were associated with high TMB (101/327) and only a small number of mutations were associated with low TMB (5/327). We found that some hotspot mutations occurred in only one cancer type. For example, the mutations P.S71L in OR2A5, P.G362E in CNTNAP2, P.F17F in BCL2L12 and P.T554I in SLC27A5 only occurred in SKCM and were associated with high TMB, with an adjusted-log10 (P-value)>1.3 and log2 (fold-change)>1. However, other hotspot mutations occurred in at least two cancer types. For example, the mutation KRAS:P.G12V appeared in BRCA, COAD, LUAD, OV, PAAD, STAD and UCEC, but only in OV was this mutation was positively related to low TMB, with a -log10 (P-value) of 2.29 and log2 (fold-change) of −2.26. As shown in Figure 3A–B, only in LUAD were the mutations P.E746_A750del and P.L858R in EGFR significantly associated with low TMB (adjusted P-value <0.00006 and log2 (fold-change) of 4.28). In other words, samples with these two hotspot mutations were usually TMB-low in LUAD and would be targetable by first-generation tyrosine kinase inhibitors. Most colorectal patients with the BRAF P.V600E mutation were TMB-high, with a P-value of 8.21E-19 (fold-change: 31.9 vs 3.4). However, the BRAF V600E mutation was associated with TMB-low in LUAD patients with P-value 0.015 (fold change: 1.9 vs 5.9). In summary, our results indicated that some hotspot mutations were strongly associated with TMB level, either low or high.

Figure 3 (A) The volcano figure of the -log10 (P-value) and the log2 (fold-change) in the mutation-positive samples (samples with this mutation) and mutation-negative samples (samples without this mutation) for 150 hotspot mutations. (B) The number of samples with 15 typical hotspot mutations in different cancer types. Abbreviations: LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma ; BLCA, bladder urothelial carcinoma; COAD, colon adenocarcinoma; BRCA, breast invasive carcer; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; KIRC, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; GBM, glioblastoma multiforme; UCEC, uterine corpus endometrial carcinoma; OV, ovarian serous cystadenocarcinoma; PRAD, prostate adenocarcinoma; ESCA, esophageal carcinoma; PAAD, pancreatic adenocarcinoma.

The impact of gene set on TMB estimation

Though WES is the best way to calculate TMB, gene panel sequencing has become increasingly popular and practical in routine molecular diagnostics. To demonstrate the capability of TMB estimation from gene panels, we carried out in silico TMB estimation for various panels using WES mutations from the public TCGA data-set. A total of 6,625 WES samples from the TCGA were filtered and analyzed in silico using the cancer-related gene panels F1CDx, MSK-IMPACT and F1CDX+MSK. TMB estimated from simulated F1CDx and MSK-IMPACT panels showed a high correlation to TMB estimated from WES, with R² correlation values of 0.95 and 0.94, respectively, for total mutations (Figure 4A and B). However, for most cancer types, there were obvious differences observed between WES and the two panels (Figure 5). F1CDx, MSK-IMPACT and F1CDX+MSK had differences in R² when correlating the data to WES, with F1CDX+MSK having the best R², followed by MSK-IMPACT and F1CDx. Despite this, the TMB calculated from these three panels correlated with TMB calculated from WES in most cancer types. However, they did not perform well in some cancer types (ESCA, KIRC and PAAD.), particularly KIRC, with low R² correlation values below 0.4. The results indicated that these three panels were only suitable for TMB calculations for a subset of cancer types. As the above three panels were not specially designed to calculate TMB values, we attempted to propose a pan-cancer TMB panel as described in Methods. Finally, 524 high-frequency mutated genes were selected (Figure S4). Comparing this pan-cancer TMB panel to the three panels above improved the R² correlation values in all 15 cancer types (Figure 5).

Figure 4 The correlation of estimated TMB by two FDA-approved panels with WES. (A) F1CDx with WES (R2=0.95). (B) MSK-IMPACT with WES (R2=0.94). Abbreviations: TMB, tumor mutation burden; WES, whole exome sequencing.

Figure 5 Correlation analysis between panels (F1CDx, MSK-IMPACT, F1CDx-MSK and our optimized pan-cancer panel) and WES. Abbreviation: WES, whole exome sequencing.

To avoid the problem of “over-fitting”, we further fitted linear regression models using approximately 80% of the samples to predict TMB values for the remaining 20% of samples for each cancer type. We assessed the correlation between the predicted and real TMB for the samples in the test data-set. Furthermore, 300, 400 or 500 genes were randomly selected five times and their average R² values were compared (Figure 6). Similar to previous results, F1CDx, MSK-IMPACT and F1CDX+MSK performed well in some cancer types, like BLCA, COAD, SKCM, STAD and UCEC, while results were poorer in ESCA, KIRC and PAAD. Excitingly, using the pan-cancer TMB panel, the correlations improved tremendously in these cancer types. Besides, we found that randomly selected gene sets had little difference between F1CDx, MSK-IMPACT and F1CDX+MSK in almost all cancer types. In general, the pan-cancer TMB panel we proposed in this study had a higher correlation with WES for all cancer types.

Figure 6 Correlation analysis between the predicted TMB and the real TMB in different panels. Abbreviations: TMB, tumor mutation burden; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma ; BLCA, bladder urothelial carcinoma; COAD, colon adenocarcinoma; BRCA, breast invasive carcer; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; KIRC, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; GBM, glioblastoma multiforme; UCEC, uterine corpus endometrial carcinoma; OV, ovarian serous cystadenocarcinoma; PRAD, prostate adenocarcinoma; ESCA, esophageal carcinoma; PAAD, pancreatic adenocarcinoma.

Discussion

Here, we used somatic mutations identified from WES data from the TCGA data-set to investigate the TMB distribution level across 15 cancer types. Frequently mutated genes and hotspot mutations and their associations with TMB level were then systematically analyzed. However, determining TMB by WES is not feasible for larger patient cohorts in current clinical practice due to its high costs, substantial turn-around time and limited availability of fresh cancer tissue samples.²⁶ Recent research efforts have demonstrated that targeted sequencing panels used routinely in the clinic may provide a reasonable estimate of TMB from only a few hundred genes of interest.^11,21,27 Our study on the correlation between TMB calculated from WES and that of different panels systematically addresses the question of whether TMB estimated from panels is equivalent to TMB estimated using public WES data from TCGA.

Recently, many organizations have provided panel sequencing for cancer patients with solid tumors for the purpose of diagnosis, targeted therapy and prognosis. Generally, such panels consist of a few genes to a few hundred genes, depending on the clinical purpose. As new target drugs are approved by FDA and new guidelines from the National Comprehensive Cancer Network are released, panels should be updated with new target genes with actionable mutations. In the meantime, the development of an increasing number of immune checkpoint inhibitors requires a more accurate pan-cancer panel to estimate TMB and determine patient suitability for treatment. The pan-cancer TMB panel proposed here provides a reference data-set for gene selection in any TMB-oriented panel design. However, we still need more clinical data with immunotherapy response results to refine this pan-cancer TMB panel. We anticipate that these genes can be added into existing panels and are capable of estimating TMB to predict patient response to cancer immunotherapy.