www.nature.com/npjprecisiononcology ARTICLE OPEN Panels and models for accurate prediction of tumor mutation burden in tumor samples 1 2,3 ✉ 1,3 ✉ Elizabeth Martínez-Pérez , Miguel Angel Molina-Vila and Cristina Marino-Buslje Immune checkpoint blockade (ICB) is becoming standard-of-care in many types of human malignancies, but patient selection is still imperfect. Tumor mutation burden (TMB) is being evaluated as a biomarker for ICB in clinical trials, but most of the sequencing panels used to estimate it are inadequately designed. Here, we present a bioinformatics-based method to select panels and mathematical models for accurate TMB prediction. Our method is based on tumor-specific, forward-step selection of genes, generation of panels using a linear regression algorithm, and rigorous internal and external validation comparing predicted with experimental TMB. As a result, we propose cancer-specific panels for 14 malignancies which can offer reliable, clinically relevant estimates of TMBs. Our work facilitates a better prediction of TMB that can improve the selection of patients for ICB therapy. npj Precision Oncology (2021)5:31 ; https://doi.org/10.1038/s41698-021-00169-0 1234567890():,; INTRODUCTION Tumors with high TMB are likely to generate more neoantigens, Tumor cells escape immune surveillance through a variety of which can be displayed on HLA molecules on the surface of the mechanisms, including upregulation of three key immune cancer cells and be recognized by the immune system24–27. A checkpoint proteins, programmed cell death protein 1 (PD-1), significant number of retrospective studies, some of them PD-1 ligand (PD-L1), and cytotoxic T-lymphocyte associated performed in tumor samples from clinical trials, have found an protein 4 (CTLA-4)1. Treatment with antibodies targeting these association between TMB and outcome to ICIs23,28–30. However, some negative results have also been reported31,32 and recent proteins, known as immune checkpoint blockade (ICB) therapy, clinical trials have yielded contradictory results33. Some reasons can overcome escaping mechanisms, restoring T cell activity and can help to explain these inconsistencies. Different cut-off values allowing the adaptive component of the immune system to detect have been adopted in different trials for the definition of “high and kill the tumor cells. Clinical trials have shown that ICB therapy TMB”23, and the possibility that TMB might not have homo- improves overall survival (OS) of patients with advanced non- geneous predictive power across all malignancies has not been small-cell lung cancer (NSCLC), melanoma, renal cell carcinoma, considered34. More importantly, although some studies have used urothelial cancer, Hodgkin’s lymphoma, and other tumors2–7. whole-exome sequencing (WES) to calculate TMBs35,36, this However, not all metastatic patients respond to ICB therapy and technique cannot be easily applied in the clinical setting due to only 15–35% of cases derive a durable clinical benefit1. In its high cost, material needs, long turnaround time, and consequence, predictive biomarkers are needed8. complicated bioinformatics pipelines. In consequence, gene At present, evaluation of PD-L1 protein expression by immuno- panels targeting around 1 MB (3% of the total exome) are histochemistry (IHC) in tumor and inflammatory cells is the most currently being employed to estimate TMB in all types of solid widely used biomarker for the selection of cancer patients for ICB tumors. Two examples are the Foundation One (FO)37 and the therapy. It has received FDA-approval as a companion or MSK-Impact34 panels, which cover the exomes of 315 and 468 complementary test in certain malignancies, including NSCLC or cancer-related genes. Although these panels have demonstrated urothelial carcinoma8. However, IHC presents a number of draw- some predictive power for clinical benefit to ICB therapy24,28,34,38, backs. Several anti-PD-L1 antibodies and cut-off values are used for they were designed to detect driver mutations and other different immune checkpoint inhibitors (ICIs) and pathologies9, and biologically relevant genetic alterations, not to calculate TMB. To their concordance is far from being perfect10–12. In addition, the date, to the best of our knowledge, bioinformatic tools have been predictive power of PD-L1 IHC, although well established, is limited; rarely applied for the rational design of gene panels to accurately extrapolate TMB, and the results have never been clinically with only 35–45% of patients responding to ICB therapy in the validated. highest expression group9,11. Consequently, it is now generally Here, we present a bioinformatics-based approach for the assumed that PD-L1 testing alone is not sufficient and should be selection of gene panels and associated models to determine complemented with other markers8,13. To this end, transcriptomic TMBs. First, we studied the distribution of mutant genes and and epigenetic signatures14–16, the presence of tumor-infiltrating found a high variability across the different human tumors. lymphocytes17, gut microbiome18–20, oncogenic mutations16, mis- Consequently, a cancer-specific strategy was employed where match repair deficiency21,22, and tumor mutation burden (TMB)23 are relevant genes were selected and used to generate panels by a currently under investigation. Among them, TMB is probably the multiple linear regression algorithm. The TMB predicted by these biomarker that has received more attention, being the focus of panels was compared with experimental TMBs, the smaller panels multiple preclinical and clinical studies. with good accuracy were subjected to rigorous internal and 1 Bioinformatics Unit, Fundación Instituto Leloir, Buenos Aires, C1405BWE, Avda. Patricias Argentinas 435 C1405BWE, Ciudad Autonoma de Buenos Aires, Argentina. 2Laboratorio de Oncología/Pangaea Oncology, Hospital Universitario Quirón Dexeus, Barcelona, Spain. 3These authors contributed equally: Miguel Angel Molina-Vila, Cristina Marino-Buslje. ✉email: mamolina@panoncology.com; cmb@leloir.org.ar Published in partnership with The Hormel Institute, University of Minnesota E Martínez-Pérez et al. 2 a non-melanoma skin tumor melanoma 193 502 138 952 170 65 b 20000 49% 79% 46% 42% 68% 17% 32% 99% 51% 92% 93% 99% 97% 97% 86% 99% 100% 99% 99% 100% 100% colorectal carcinoma 2112 129 33 20% 12% 39% 29% 65% 18% 82% 56% 47% 95% 98% 95% 96% 99% 99% 98% 88% 93% 100% 100% 100% 697 98 7 10000 esophageal squamous carcinoma cervical carcinoma 185 97 10 5000 gastic carcinoma 637 86 22 2000 Number of mutations uninary tract other 273 85 7 urinary tract transitional cell 107 82 4 1000 endometrial carcinoma 303 79 14 esophageal other tumor 88 76 4 500 liver carcinoma 2065 76 30 ovarian other carcinoma 155 72 4 200 papillary cell renal carcinoma 205 64 2 100 head and neck tumor 666 61 8 clear cell renal carcinoma 800 47 2 50 pancreatic tumor 1090 45 7 290 40 3 20 esophageal adenocarcinoma small cell lung carcinoma 284 40 4 10 2301 40 8 Cancer type breast carcinoma 5 biliary tract carcinoma 355 39 6 ovarian serous carcinoma 543 38 2 2 non-small cell lung carcinoma 1324 34 6 glioma 1295 33 2 1 prostate carcinoma 1565 31 4 renal cell carcinoma 204 25 2 neuroblastoma wilms tumor small intestine tumor myeloid leukaemia pituitary tumor adrenal gland tumor germ cell tumor medulloblastoma meningioma soft tissue sarcoma salivary gland tumor parathyroid tumor lymphoid leukaemia pancreatic tumor papillary cell renal carcinoma pancreatic endocrine tumor chondro and osteosarcoma thyroid carcinoma mesothelioma renal cell carcinoma prostate carcinoma glioma non-small cell lung carcinoma ovarian serous carcinoma biliary tract carcinoma breast carcinoma small cell lung carcinoma esophageal adenocarcinoma clear cell renal carcinoma head and neck tumor ovarian other carcinoma liver carcinoma esophageal other tumor endometrial carcinoma urinary tract transitional cell uninary tract other gastic carcinoma cervical carcinoma esophageal squamous carcinoma colorectal carcinoma melanoma non-melanoma skin tumor mesothelioma 131 21 2 thyroid carcinoma 619 17 4 chondro and osteosarcoma 407 16 2 pancreatic endocrine tumor 150 16 2 lymphoid leukaemia 1803 15 2 parathyroid tumor 35 14 2 salivary gland tumor 78 14 2 soft tissue sarcoma 241 12 2 meningioma 65 8 1 medulloblastoma 488 7 1 germ cell tumor 86 7 1 adrenal gland tumor 316 7 1 pituitary tumor 57 6 1 myeloid leukaemia 972 6 2 small intestine tumor 50 4 2 Cancer type wilms tumor 183 4 1 neuroblastoma 356 3 1 1 2 5 10 20 50 100 200 500 1000 2000 5000 10000 20000 class DELETE KEEP Number of mutations Fig. 1 Mutations in the tumor samples of COSMIC v84 (training dataset). a Distribution of TMBs across 42 cancer types in the 24,726 samples of the training dataset. The “Y” axis corresponds to the frequency of mutant samples. Numbers on the right of the plot indicate: black, number of samples included in the analysis; red, median TMBs; blue, median of common mutant genes. Cancer types are 1234567890():,; ordered by median TMBs. b Number of mutations in the individual tumor samples of COSMIC v84. Each dot represents a sample. In blue, samples conserved for the development of panels for TMB prediction, in red samples discarded. Numbers on top of each cancer type are the percentage of conserved samples. external validation, and those showing the best performance were these “gene-sets” were used for further analysis (Fig. 1b). The finally selected. As a result of this analysis, we propose cancer- “gene-sets” were submitted to a multiple linear regression specific panels for 14 malignancies and we demonstrate that algorithm in order to generate panels and associated mathema- these panels, which range between 0.23 and 1.49 Mb (33 and 231 tical models. The TMB predicted using these panels and models genes, respectively), can reliably estimate TMBs. was compared with experimental TMBs and panels with the minimum number of genes and good accuracy were used to design cancer-specific consensus panels. Consensus panels were RESULTS then subjected to rigorous internal and external in silico Distribution of TMB across cancer types validations and those showing the highest accuracy were finally The aim of our study was to rationally select panels, with a limited selected (see the “Methods” section). As a result of this process, we number of genes or exons, which could be used to predict TMB obtained 1446 panels and associated models for 14 cancer types using appropriate mathematical models. To this aim, we first that allow for rigorous prediction of TMB. The 14 cancer types analyzed the distribution of experimental TMBs in our training include the most common solid malignancies, such as non-small- dataset of 24,726 tumor samples from 42 cancer types with WES cell lung carcinoma (NSCLC), colorectal carcinoma (CRC), mela- or WGS data available in COSMIC_v84 (Fig. 1a). COSMIC is the noma, prostate, and breast carcinoma. The total number of panels most extensive mutation database in cancer and, since 2010, per cancer type, as well as the number of Mb and genes that integrates the samples from The Cancer Genome Atlas Project should be sequenced in order to estimate TMBs, are shown in (TCGA), the International Cancer Genome Consortium (ICGC), and Table 1 and Fig. 2a; while Fig. 2b shows the correlation (R2) with the Cancer Genome Project (CGP)39,40. Remarkable differences experimental TMBs and the error (expressed as RSE) of the TMBs were observed in median TMBs of COSMIC_v84 samples, which predicted according to our panels. showed a >100-fold variation, ranging from 3 in autonomic Next, we applied the same strategy to “our exon-sets” and to ganglia to 502 in skin carcinoma. In addition, some tumor types the genes and exons of the FO panel and the Cancer Gene Census presented an apparently unimodal distribution of TMBs, while (CGC). The FO and CGC sets allowed us to generate valid bimodal distributions were observed in others, such as stomach, consensus panels only for 3 and 8 cancer types, respectively. endometrium, or lung small cell and non-small cell carcinomas. Models based on “our gene-sets” systematically predicted TMB The distribution of common mutant genes across samples was with more accuracy than models based on the genes or exons of also analyzed. Although some malignancies such as skin the FO-panel or CGC (p < 0.05 in Kruskal–Wallis and p > 0.05 in carcinoma, melanoma, large intestine or liver carcinoma showed Dunn tests; Supplementary Data 1). Also, models based on genes a median of more than 30 mutant genes in common between were systematically superior to models based on exons (Supple- individual tumors, in most cancer types this percentage ranged mentary Figs. 2 and 3). from 0 to 10 (blue numbers in Fig. 1a). Due to different reasons, we could not generate panels for TMB prediction in 28 out of the 42 initial cancer types. Some tumors (n = Development of panels for TMB prediction 12) had an insufficient number of samples after data curation. In In view of the significant variability among different types of other cases (n = 5), the models and panels obtained did not qualify tumor, a cancer-specific strategy was used to design panels for to generate consensus models. Finally, consensus models for some TMB prediction (Supplementary Fig. 1). After excluding the largest cancer types did not pass the internal (n = 5) or external (n = 6) genes, we selected for each cancer type the genes with mutations validations due to high error (RSE) or low correlation (R2) with in ≥1% of the samples. Samples with mutations in at least one of experimental TMB. npj Precision Oncology (2021) 31 Published in partnership with The Hormel Institute, University of Minnesota E Martínez-Pérez et al. 3 Table 1. Characteristics of the total (left) and suggested panels (right, bold) for accurate prediction of TMB in 14 cancer types. Type of tumor Number of panels Genes (range) Mb (range) Suggested panel Genes Suggested panel Mb Breast carcinoma 161 130 (48–210) 1.05 (0.48–1.56) 199 1.49 Glioma 68 63 (21–101) 0.34 (0.11–0.52) 98 0.52 Endometrial carcinoma 25 28 (16–40) 0.21 (0.14–0.28) 33 0.23 Colorectal carcinoma 280 173 (33–323) 1.06 (0.26–1.81) 185 1.10 Liver carcinoma 250 150 (26–278) 1.30 (0.23–1.74) 231 1.48 Non-small cell lung carcinoma 107 130 (77–185) 0.79 (0.51–1.02) 175 0.97 Ovarian serous carcinoma 26 52 (34–67) 0.35 (0.23–0.43) 61 0.40 Pancreatic tumor 131 19 (17–31) 0.66 (0.25–1.11) 111 0.76 Prostate carcinoma 108 103 (49–170) 0.76 (0.37–1.17) 124 0.90 Non-melanoma skin carcinoma 31 24 (9–40) 0.14 (0.07–0.25) 40 0.25 melanoma 110 70 (15–126) 0.51 (0.17–0.79) 126 0.79 Gastric carcinoma 68 63 (30–99) 0.41 (0.24–0.62) 80 0.49 Head and neck tumor 58 67 (36–96) 0.44 (0.23–0.59) 89 0.56 Urinary tract other 23 30 (19–41) 0.20 (0.13–0.28) 40 0.28 Bold values indicate the number of genes and size of the only panel for TMB estimation we suggest for each malignancy. a 1.50 b breast carcinoma colorectal carcinoma melanoma non-small cell lung carcinoma 1.00 1.00 0.50 0.99 Mb R2 0.25 0.98 0.10 0.97 160 breast carcinoma colorectal carcinoma endometrial carcinoma gastic carcinoma glioma head and neck tumor liver carcinoma melanoma non-melanoma skin tumor non-small cell lung carcinoma ovarian serous carcinoma pancreatic tumor prostate carcinoma uninary tract other 120 RSE 80 40 0.10 0.25 0.50 1.00 1.50 0.10 0.25 0.50 1.00 1.50 0.10 0.25 0.50 1.00 1.50 0.10 0.25 0.50 1.00 1.50 Number of Mb Cancer type strategy CGC - exon CGC - gene FO-panel - gene our - exon our - gene Fig. 2 Characteristics of the panels for TMB prediction. a Size (in Mb). The red line indicates the size of the FO-panel. The center line in the boxes corresponds to the median, the bounds of the boxes to the first and third quartiles, and the whiskers extend to the largest and lowest values no further than 1.5 times the inter-quartile range. Data beyond the end of the whiskers are called “outlying” points and plotted individually. b R2 and RSE of the panels vs. size (in Mb). Only 4 representative cancer types are shown; breast carcinoma, colorectal carcinoma, melanoma, and non-small cell lung carcinoma (for other tumors, see Supplementary Fig. 4). The red line indicates the size of the FO-panel. Each dot represents a model, colors indicate initial datasets. Overlap between panels NSCLC (n = 34). First, we predicted the TMB for each sample using Next, we analyzed the overlap of the panels for TMB prediction our suggested panels of 126 genes (0.79 Mb) for melanoma and generated from “our gene-sets”. The percentage of genes shared 175 genes (0.97 Mb) for NSCLC. As expected, predicted TMBs between the panels for different tumors ranged from 0 to 40% and, showed an excellent correlation with the published, experimental in some cases such as urinary tract or non-melanoma skin tumor, TMBs (R2 = 0.96 and 0.84 for melanoma and lung, respectively; never reached 10% (Fig. 3a). Then, we compared the panels for TMB Fig. 4a–b). Regarding the response to ICB therapy, we evaluated prediction with the 314 genes of the FO-panel and the 719 genes of three different cut-off values for predicted TMB; 100, 150, and 200 the CGC. Most of the genes in our panels were not contained in the mutations. The 150 threshold showed the best correlation with FO-panel or in the CGC (Fig. 3b). When discriminating by cancer clinical outcomes (Supplementary Figs. 4–6). Response rate of type, the overlap was always <20% (Fig. 3c). melanoma patients with >150 predicted mutations was 51%, compared with 18% in those with <150 predicted mutations (p < 0.05 in a z-score test) (Supplementary Fig. 7). In addition, patients Suggested panels and models and their correlation with with high predicted TMB had a median OS of 24.3 months, response to immunotherapy significantly longer than the 10.1 months of patients with low In order to facilitate implementation in the clinical setting, we predicted TMB (HR = 0.66; CI95% = 0.44–0.93; P = 0.02) (Fig. 4c). In made a final selection of one panel and associated model for each the case of the NSCLC cohort, the 150 cut-off value also had the of the 14 cancer types mentioned above. These best panels best performance (Supplementary Fig. 8), with 56% of patients with showed optimal performances in terms of correlation with high predicted TMB showing partial responses to ICB therapy, experimental TMB and low error, having sizes between 0.23 and compared to 17% in patients with low TMB (p < 0.05 in a z-score 1.49 Mb (Table 1 and Supplementary Data 2). test) (Supplementary Fig. 7). Only progression-free survival (PFS) but To evaluate the association of the TMBs estimated with our not OS data was available for this cohort. Despite the relatively low suggested panels and outcome to ICB therapy, we used two number of patients, the median PFS of those patients with >150 published cohorts of metastatic melanoma (n = 174) and advanced predicted mutations was significantly longer, 14.5 months vs. Published in partnership with The Hormel Institute, University of Minnesota npj Precision Oncology (2021) 31 E Martínez-Pérez et al. 4 a uninary tract other b our-gene CGC c uninary tract other prostate carcinoma prostate carcinoma pancreatic tumor pancreatic tumor percent 69 ovarian serous carcinoma 1095 417 ovarian serous carcinoma percent <1 non-small cell lung carcinoma non-small cell lung carcinoma 64 <1 1-5 11 169 Cancer type Cancer type non-melanoma skin tumor non-melanoma skin tumor 1-5 melanoma 5-10 melanoma 70 5-10 liver carcinoma 10-15 liver carcinoma 10-15 head and neck tumor 15-20 head and neck tumor 15-20 glioma 20-40 glioma Fo-panel >20 gastic carcinoma gastic carcinoma >40 endometrial carcinoma endometrial carcinoma colorectal carcinoma colorectal carcinoma breast carcinoma breast carcinoma breast carcinoma colorectal carcinoma endometrial carcinoma gastic carcinoma glioma head and neck tumor liver carcinoma melanoma non-melanoma skin tumor non-small cell lung carcinoma ovarian serous carcinoma pancreatic tumor prostate carcinoma uninary tract other CGC FO-panel Percent of genes in Cancer type Fig. 3 Gene overlap between panels. a Percentage of shared genes between the panels for TMB prediction based on “our gene-set”. b Euler- Venn diagram of genes in the panels for TMB prediction based on “our set”, CGC genes and genes in the FO-panel. c Percentage of shared genes between panels for TMB prediction, CGC and FO-panel, per cancer type. a 4000 Lineal regression c 1.00 d 1.00 R2 = 0.96 3500 RSE = 220.98 ++ experimental TMB 3000 Spearman rank correlation ρ = 0.94 Log-rank p = 0.02 Log-rank p = 0.0019 0.75 0.75 Probability of PFS 2500 P.value = 6.9e-85 Probability of OS 2000 ++ + ++++ 1500 0.50 ++++ 0.50 ++ 1000 +++ 500 + ++ + + +++++ +++ + 0 0.25 ++++++++ + + 0.25 0 500 1000 1500 predicted TMB 2000 2500 3000 + + 0.00 0.00 b 1750 Lineal regression 0 10 20 30 40 50 60 70 0 5 10 15 20 25 R2 = 0.84 1500 Months Months RSE = 136.95 experimental TMB 1250 Spearman rank correlation 1000 ρ = 0.92 P.value = 1.6e-14 TMB (150 mutations) + High + Low TMB (150 mutations) + High + Low TMB (150 mutations) TMB (150 mutations) 750 Patients at risk Patients at risk 500 High 100 64 54 32 18 14 7 4 High 16 10 5 1 0 0 250 Low 74 37 24 18 9 4 3 3 Low 18 7 2 1 1 1 0 0 100 200 300 400 500 600 700 0 10 20 30 40 50 60 70 0 5 10 15 20 25 predicted TMB Months Months Fig. 4 Predicted TMB and response to immunotherapy. a, b Correlation between experimental TMB and the TMB predicted according to our suggested panels for two reported cohorts of (a) melanoma (n = 174) and (b) NSCLC (n = 34) patients treated with ICB therapy (two melanoma samples with >6000 mutations are omitted for clarity). c Kaplan–Meier plot of OS in the 174 metastatic melanoma patients, stratified according to predicted TMB (cut-off, 150 mutations). d Kaplan–Meier plot of PFS in the 34 advanced NSCLC patients, stratified according to predicted TMB (cut-off, 150 mutations). 3.5 months for patients with low predicted TMB (HR = 0.27; been designed to extrapolate TMB. Instead, they target genes CI95% = 0.10–0.57; P = 0.002) (Fig. 4d). carrying driver mutations or other clinically relevant alterations, The performance of predicted TMB to evaluate clinical benefit which are not necessarily representative of the total number of was also analyzed using ROC curves (Supplementary Figs. 9–10). mutations in the entire exome. Furthermore, the same NGS panels The areas under the curve (AUC) (0.67 for the melanoma and 0.70 are used in all solid tumors, ignoring the fact that mutant genes for the NSCLC cohort) were similar to those obtained when using are not the same and do not present the same frequency of experimental TMBs (0.67 and 0.76, respectively). mutation across different malignancies (Fig. 1). Finally, TMB estimates based on NGS panels are subjected to biases that can derive in the misclassification of up to 19–21% of samples when DISCUSSION using large panels (1.1–1.4 Mb), or up to 33–36% for smaller TMBs should ideally be determined by WGS or WES but the panels41. In consequence, the development and validation of gene widespread implementation of these techniques in the clinical panels for accurate estimation of TMB is of particular interest. setting is difficult due to long turnaround time, complexity, cost, However, to the best of our knowledge; rational, bioinformatics- and tissue consumption. In consequence, NGS panels34,37 ranging based strategies have rarely been attempted to achieve this aim. from 1.1 to 3 Mb are increasingly used to estimate TMB in tumor Here, we present a data-driven strategy that allowed us to samples, and two of them have recently received FDA approval develop cancer-specific panels and associated mathematical after demonstrating some predictive power in several malignan- models for accurate prediction of TMB in 14 types of tumors; cies23. However, the majority of NGS panels in clinical use have not including the most common malignancies such as NSCLC, CRC, npj Precision Oncology (2021) 31 Published in partnership with The Hormel Institute, University of Minnesota E Martínez-Pérez et al. 5 breast carcinoma, melanoma or prostate carcinoma. These panels, vs. 4–25% response rates obtained when the FO or MSK-Impact which were submitted to rigorous internal and external in silico panels were employed23,28,49. Regarding survival, patients with validations, range from 0.23 to 1.49 Mb, a size comparable with the high predicted TMB according to our models showed 14.5 months panels currently used in the clinical setting (0.9 Mb for FO and PFS in the lung and 24.3 months OS in the melanoma cohorts; 1.22 Mb MSK-Impact). A recent report42 concluded that 1.5–3 Mb compared to 3.5 and 10.1 months respectively for patients with panels are better suited to calculate TMB, whereas smaller panels low TMB. Again, these values are in the range of the figures can yield considerably imprecise estimates. However, we have reported in the literature for patients stratified according to WES demonstrated that the <1.5 Mb panels selected using our data; 7.1–17.1 months for high TMB lung and 24–25 months for bioinformatics-based strategy can accurately predict TMB (Fig. 2a). high TMB melanoma vs. 1.4–3.7 and 9–10 months, respectively, for Several studies have reported the development of methods for patients with low TMB24,46–48. TMB calculation based on commercial NGS panels43,44. However, to In summary, we have developed a new bioinformatics-based the best of our knowledge, there is only one study where “in silico” method to select NGS panels and associated mathematical models tools were used to select genes for TMB prediction45. As a result, for accurate prediction of TMB, which we have validated using the authors proposed a 525 gene pan-cancer panel, which is clinical data. The panels we suggest are tumor-specific, can be easily significantly larger than the 33 to 231 gene cancer-specific panels implemented in the experimental setting and outperform the driver- we suggest (Supplementary Data 2). In some cancer types, this 525 oriented NGS panels currently used to estimate TMB in tumor gene panel showed an average correlation with experimental TMB samples. In addition, the TMB predicted by our panels correlates (R2) comparable to our panels. However, the variability in the error with clinical outcome to ICIs. Our work can facilitate a more accurate of the samples (RSE), which can lead to sample misclassification, prediction of TMB which, in turn, can improve the discrimination was not reported by the authors. In addition, the 525 panel was not between patients responding and not responding to ICB therapy. validated with clinical data to determine its performance in terms of prediction of clinical outcome to ICB therapies. Finally, the common genes among our suggested tumor-specific panels range METHODS from 1 to 40%, indicating that a single selection of genes is not Training dataset optimal to estimate the TMB due to the heterogeneity of mutant We downloaded the WGS and WES datasets of the Cosmic version 84 genes across the spectrum of human malignancies. (Cosmic_v84)39 for human grch38 assembly, which contained 25,533 tumor During our analysis we also found that, although commercial samples classified into cancer types according to the primary site, primary NGS panels can offer acceptable estimations of TMB in certain histology, and histology subtype. Then, we excluded (i) benign tumors, (ii) malignancies, they are not optimal for this purpose. The cancer- samples with mutations without annotation of genomic coordinates, (iii) specific panels we suggest outperform the FO panel and the CGC samples with non somatic mutations (labeled as “Variant of unknown origin” or “Not specified”), (iv) samples with mutations exclusively in non gene set for TMB estimation (Fig. 2, Supplementary Fig. 2). For coding regions. Finally, cancer types with <10 samples were also excluded. some cancer types, only our set of genes generated acceptable Our training dataset contained 24,726 samples of 42 cancer types. For models for TMB prediction while for other malignancies, although detailed information of samples and filtering see Supplementary Data 3. sometimes FO or CGC-based estimates were acceptable, our panels showed a significantly better R2 and lower RSE. We also found most of the best genes for TMB calculation are not included External validation dataset in the FO panel or CGC set, being the overlap <20% (Fig. 3b–c). The external validation dataset was based on the new WGS and WES These findings further support the idea that panels composed of samples added to COSMIC_v90 (n = 3144), together with samples with WGS and WES data retrieved from 133 articles publicly available in the oncogenes and tumor suppressor genes systematically introduce literature (n = 4773) (Supplementary Data 4). All samples were mapped to biases and have a higher variance than those based on a random the grch38 genomic coordinate and filtered according to the procedure selection of genes42. described for the training dataset (Supplementary Data 1). Overall, the In our study, TMB was considered as the total number of exonic external validation dataset contained 7917 samples from 40 cancer types. somatic mutations; a definition aligned with the scientific literature Parathyroid and pituitary cancers were not represented. and published clinical trials. Although the analysis of variants in non-coding regions might be of interest, we excluded them for Study population in the clinical validation set several reasons. The non-coding segments of the genome are We used 174 metastatic melanoma35,36 and 35 NSCLC24 patients much longer than the exons, the number of possible mutations is (Supplementary Data 3) with anonymized WES data and clinical outcome enormous and a vast majority of them are probably silent, unlikely to ICB therapy publicly available. All patients had been enrolled in clinical to generate neoantigens. Also, the number of samples with whole- trials outside of our institution. The patients had signed informed consent genome sequencing (WGS) and clinical data available is relatively under institutional review board (IRB) approved protocols. In addition, the small, probably insufficient to obtain statistically significant results. results of the genetic testing and outcome to ICIs had been previously Therefore, testing if a “WGS TMB” can be superior to the commonly published by the corresponding investigators. No additional human clinical used “WES TMB” would require a comprehensive study, where data was used. For these reasons, the IRB of the Quirón Dexeus University WGS should be performed in a significant number of patients Hospital granted a waiver for this study. receiving ICB therapy with response data available. Clinical validation of our results in a prospective study would Outline of the development of panels and models for TMB add considerable value to our work. However, in the 14 cancer prediction types where we could select panels and models, our predicted TMB was defined as the total number of exonic mutations in a given TMB highly correlated with the experimental TMB. Also, the TMBs sample. Consequently, only mutations in exons were considered. Our aim calculated using our suggested panels closely correlated with the was to select cancer-specific panels, with a limited number of genes or TMB experimentally obtained by WES in two small cohorts of exons, which could be used to predict TMB using appropriate mathema- melanoma and lung cancer patients with publicly available data tical models (Supplementary Fig. 1). First, we excluded from our training database the largest genes (11 genes >18928 bp; 0.05% of the genome), (Fig. 4). In addition, when patients were classified according to our some of which had been also removed or penalized by length in the estimates, response rates were 51–56% for the high TMB vs. previous studies50. Next, for each cancer type, we selected all genes 17–18% for the low TMB group. Both values are in excellent harboring mutations in at least 1% of the samples and obtained what we agreement with the 46–59% vs. 12–23% reported in several called “our gene-sets” and “our exon-sets”. These sets were the input for clinical trials of ICB therapies in lung and melanoma tumors where developing gene and exon models, respectively. Regarding samples, those TMB was determined by WES24,46–48 and outperform the 20–45% with mutations in at least one of “our genes” were selected. We found that Published in partnership with The Hormel Institute, University of Minnesota npj Precision Oncology (2021) 31 E Martínez-Pérez et al. 6 some cancer types had <80% of samples with mutations in at least one of Generation of initial panels and associated linear models for the genes of “our” sets (Fig. 1b and Supplementary Fig. 2). Therefore, they TMB prediction were excluded from further analysis. For the rest of the cancer types, we A forward step algorithm was used for generating the first set of panels generated panels to predict TMB using multiple linear regression and associated linear regression models for TMB prediction. The accuracy algorithms that assigned appropriate coefficients to each gene. Predicted of the panels and associated models was measured with four evaluators. TMB accuracy was assessed by the coefficients of correlation (R2) and Three of them were metrics; namely R2, adjusted R2, and RSE; the fourth relative standard error (RSE) when compared with experimental TMB. was the sum of the samples selected for the development of the model. Panels with the minimum number of genes, R2 > 0.6 and RSE lower than The R2 or coefficient of determination was defined as the proportion of the the median TMB by cancer type were subsequently used to generate variance in the TMB (dependent variable) that could be predicted from the cancer-specific consensus panels and to develop associated models. Then, independent variables (genes/exons in the models). The R2 was adjusted consensus panels and models were subjected to internal and external by the number of variables in the models to obtain the adjusted R2. Finally, validation and those showing the highest R2 and lowest RSE were finally RSE was defined as the square root of the estimated variance of the selected (Supplementary Data 5). For comparison purposes, we applied the random error. same procedure to the genes of the FO-panel and the CGC51. A detailed In each step of the iteration, we used the four evaluators to select the description of the development of cancer-specific panels and models for 100 best panels and associated models. Then, we deleted the so-called TMB prediction is presented in the sections below. “non-feasible” panels, defined as those having ≥1 gen/exon not statistically relevant. The remaining panels and models were subsequently employed Selection of genes and samples as seeds for the next step. The forward step algorithm ended when it could Within each cancer type, genes were ordered according to the number of not add more genes/exons to the panel. All final panels and associated mutant samples and the top 10 genes were selected. A total of 157 models generated in each step were considered for further analysis; the different genes appeared in this analysis, including, TP53, TTN, MUC16, total number was 731,958 for 30 cancer types. Then, we selected those MUC4, KRAS, PIK3CA, SYNE1, APC or OBSCN, (Supplementary Fig. 11a). Some models with R2 > 0.6 and RSE < median TMB by cancer type; ending up of them were long genes codifying for large proteins, such as TTN and with 354,644 panels. MUC16, which have coding regions of 109 and 44 Kb, respectively. In consequence, similarly to the previous studies50, we penalized genes by Generation of consensus panels and associated models length (see below). After the penalization, 57 genes no longer appeared Due to a large number of panels and associated models obtained as a among the top 10 genes most frequently mutated by cancer type, result of our analysis, we decided to generate cancer-specific consensus including TTN or MUC16; while and 100 genes remained, such as TP53, panels. Some genes/exons, which only appeared in a few panels, were KRAS, PIK3CA, BRAF, or APC (Supplementary Fig. 11b). eliminated to select a core of genes per cancer type according to the The human genome assembly Grch38 comprises 20,291 coding genes, following strategy. First, the genes/exons present in the initial panels were which were obtained using the “gtf data” from the Ensembl version ordered by the number of occurrences. Then, they were added in a grch38.91. For each gene, the number of coding bases was calculated; the forward step manner to obtain a series of consensus panels. Those with ≥1 median length was 1341 bp; Q1 and Q3 were 828 and 2139 bp, gene/exon with no statistical significance were deleted, leaving a total of respectively. Among the 20,291 coding genes of the human genome, 7,242 consensus models in 25 cancer types. These consensus models were 499 (2.46%) showed a disproportionate large size of the coding region submitted to 1000 repetitions of bootstrap internal validation and then to (Supplementary Fig. 12), higher than three times the interquartile distance an external validation using samples obtained from a newer version of (6072 bp), and were considered extreme outliers. To avoid the bias due to COSMIC and the literature (see External Validation Dataset, in the the gene length, the number of mutant samples per gene (G) and cancer “Methods” section). The rules used to select the acceptable consensus type (C) was corrected as shown in Eq. 1. For example, the FLG gene has panels and the final numbers of consensus panels after each step of the 12,186 coding bp and 379 mutant samples among large intestine tumors. validation are shown in Supplementary Data 3 and Supplementary Table 1. After applying the equation, they counted as 188 mutant samples [379/ A total of 2301 panels and associated models for 14 cancer types passed (12,186/6072)]. The equation did not penalize genes with <6072 coding bp. the internal and external validation process. Panels developed from “our samples in C gene-set” showed the best results in terms of quality parameters (R2 and samples penalized ðG; CÞ ¼   (1) RSE in the internal and external validations), as demonstrated with the max 1; coding base 6072 pairs of G Kruskal–Wallis and posthoc Dunn statistical tests (Supplementary Data 1). In order to develop the cancer-specific panels for TMB prediction, we To facilitate the implementation in the clinical setting, we made a final selected the genes with mutations in >5 samples representing >1% of the selection of the best panels and associated models for each of the 14 total number of samples in a particular cancer type; generating what we cancer types. We also selected 25 panels with almost optimal performances called “our gene-sets”. Regarding samples, those with mutations in at least of different sizes (in Mb). The genes contained in the 14 + 25 panels and one of “our genes” were selected. Some cancer types had <80% of samples their associated models are presented in Supplementary Data 2. with mutations ≥1 of the genes in “our set” and were excluded from further analysis (Fig. 1b). The final number of samples and cancer types Clinical validation. Correlation between response to selected for further analysis is presented in Supplementary Data 3. immunotherapy and predicted TMB Using a similar strategy, we generated “our exon-sets”. First, to avoid the bias due to the gene length, the number of mutant samples per gene (G), First, we evaluated the correlation between experimental and predicted exon (E) and cancer type (C) was corrected according to Eq. 2. Then, we TMB by Spearman’s rank correlation and linear regression in the cohorts of selected the exons with mutations in >2 samples representing >1% of the patients mentioned above. Then, the association of predicted TMBs with total number of samples in a particular cancer type; generating “our exon- OS and PFS to ICB therapy was analyzed at cut-off values of 100, 150, and sets”. Finally, only samples with mutations in at least one of “our exons” 200 mutations. The correlation between predicted TMB and type of clinical were selected (Supplementary Fig. 13 and Supplementary Data 3) response was also evaluated. SPSS and GraphPad Prism v6.0 were used to generate Kaplan–Meier plots and perform log-rank tests in order to samples penalizedðG; CÞ compare OS and PFS data; while ROC curves were made with plotROC52 samples penalized 2 ðE; G; CÞ ¼   (2) coding base pairs of G and Venn diagrams with VennDiagram53. coding base pairs of E For comparison purposes, we applied the same strategies to the genes Reporting summary and exons of the Cancer Gene Census (CGC)51 and the Foundation One panel (FO-panel)37. The CGC was downloaded from the Cosmic_v84 version Further information on research design is available in the Nature Research and contained 719 genes. Six of them (IGH, IGK, IGL, TRA, TRB, and TRD) were Reporting Summary linked to this article. “macrogenes” and we replaced them with all possible Ensembl gene names within the same genomic coordinates, ending with a dataset of 1080 Ensembl genes. Regarding the 315 genes of the FO-panel, one of them had DATA AVAILABILITY STATEMENT no coding sequence and was excluded from the analysis. Using the same The data generated and analysed during this study are described in the following strategies discussed above, we selected the samples for further analysis data record: https://doi.org/10.6084/m9.figshare.1407445154. The datasets and code (Supplementary Fig. 13, Supplementary Data 3, Supplementary Table 1). are available as part of this data record, and also via https://gitlab.com/ npj Precision Oncology (2021) 31 Published in partnership with The Hormel Institute, University of Minnesota E Martínez-Pérez et al. 7 bioinformatics-fil/predict_tmb. A list of the files underlying the figures, tables and 23. Chan, T. A. et al. Development of tumor mutation burden as an immunotherapy supplementary tables of the related manuscript is available as part of the data record biomarker: utility for the oncology clinic. Ann. Oncol. 30, 44–56 (2019). in the file ‘Martinez-Perez_et_al_underlying_data_list.xlsx’. For details about the 24. Rizvi, N. A. et al. Cancer immunology. Mutational landscape determines sensitivity software used, see the methods of the related publication. to PD-1 blockade in non-small cell lung cancer. Science 348, 124–128 (2015). 25. Alspach, E. et al. MHC-II neoantigens shape tumour immunity and response to immunotherapy. Nature 574, 696–701 (2019). CODE AVAILABILITY 26. Gubin, M. M., Artyomov, M. N., Mardis, E. R. & Schreiber, R. D. Tumor neoantigens: The codes and scripts generated and/or analyzed during the current study are also building a framework for personalized cancer immunotherapy. J. Clin. Invest. 125, available at https://gitlab.com/bioinformatics-fil/predict_tmb. All the programming 3413–3421 (2015). and related statistics were performed with R. The following R libraries were used; 27. Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G. & Hacohen, N. Molecular and Genomic Features55, biomaRt156, GenomicRanges, and rtracklayer57. Translation genetic properties of tumors associated with local immune cytolytic activity. Cell databases hg18ToHg38.over.chain and hg19ToHg38.over.chain from UCSC web page 160, 48–61 (2015). (http://genome.ucsc.edu/)58 were employed for the filtering datasets process. All 28. Hellmann, M. D. et al. Nivolumab plus ipilimumab in lung cancer with a high graphics were made with ggplot259 Library caret was used for internal and external tumor mutational burden. N. Engl. J. Med. 378, 2093–2104 (2018). validation. Wilcoxon test was employed to compare gene and exon models; the 29. Wang, F. et al. Safety, efficacy and tumor mutational burden as a biomarker of overall Kruskal–Wallis test to compare models obtained using different strategies (our-gene, survival benefit in chemo-refractory gastric cancer treated with toripalimab, a PD-1 our-exon, CGC-gene, CGC-exon, FO-panel-gene, and FO-panel-exon) and Dunn Test antibody in phase Ib/II clinical trial NCT02915432. Ann. Oncol. 30, 1479–1486 (2019). to identify differences between strategies when the Kruskal–Wallis test was 30. Ready, N. et al. First-line Nivolumab plus ipilimumab in advanced non-small-cell significant. lung cancer (checkmate 568): outcomes by programmed death ligand 1 and tumor mutational burden as biomarkers. J. Clin. Oncol. 37, 992–1000 (2019). 31. Mishima, S. et al. Clinicopathological and molecular features of responders to Received: 16 November 2020; Accepted: 10 March 2021; Nivolumab for patients with advanced gastric cancer. J. Immunother. Cancer 7, 24 (2019). 32. Morrison, C. et al. Predicting response to checkpoint inhibitors in melanoma beyond PD-L1 and mutational burden. J. Immunother. Cancer 6, 32 (2018). 33. Buttner, R. et al. Implementing TMB measurement in clinical practice: con- REFERENCES siderations on assay requirements. ESMO Open 4, e000442 (2019). 1. Ribas, A. & Wolchok, J. D. Cancer immunotherapy using checkpoint blockade. 34. Samstein, R. M. et al. Tumor mutational load predicts survival after immu- Science 359, 1350–1355 (2018). notherapy across multiple cancer types. Nat. Genet. 51, 202–206 (2019). 2. Hodi, F. S. et al. Improved survival with ipilimumab in patients with metastatic 35. Snyder, A. et al. Genetic basis for clinical response to CTLA-4 blockade in mela- melanoma. N. Engl. J. Med. 363, 711–723 (2010). noma. N. Engl. J. Med. 371, 2189–2199 (2014). 3. Brahmer, J. et al. Nivolumab versus docetaxel in advanced squamous-cell non- 36. Van Allen, E. M. et al. Genomic correlates of response to CTLA-4 blockade in small-cell lung cancer. N. Engl. J. Med. 373, 123–135 (2015). metastatic melanoma. Science 350, 207–211 (2015). 4. Borghaei, H. et al. Nivolumab versus docetaxel in advanced nonsquamous non- 37. Singal, G. et al. Association of patient characteristics and tumor genomics with small-cell lung cancer. N. Engl. J. Med. 373, 1627–1639 (2015). clinical outcomes among patients with non-small cell lung cancer using a clin- 5. Bellmunt, J. & Bajorin, D. F. Pembrolizumab for advanced urothelial carcinoma. N. icogenomic database. JAMA 321, 1391–1399 (2019). Engl. J. Med 376, 2304 (2017). 38. Balar, A. V. et al. Atezolizumab as first-line treatment in cisplatin-ineligible 6. Motzer, R. J. et al. Nivolumab versus everolimus in advanced renal-cell carcinoma. patients with locally advanced and metastatic urothelial carcinoma: a single-arm, N. Engl. J. Med. 373, 1803–1813 (2015). multicentre, phase 2 trial. Lancet 389, 67–76 (2017). 7. Ansell, S. M. et al. PD-1 blockade with Nivolumab in relapsed or refractory 39. Tate, J. G. et al. COSMIC: the Catalogue of Somatic Mutations in Cancer. Nucleic Hodgkin’s lymphoma. N. Engl. J. Med. 372, 311–319 (2015). Acids Res. 47, D941–D947 (2019). 8. Havel, J. J., Chowell, D. & Chan, T. A. The evolving landscape of biomarkers for 40. Forbes, S. A. et al. COSMIC: mining complete cancer genomes in the Catalogue of checkpoint inhibitor immunotherapy. Nat. Rev. Cancer 19, 133–150 (2019). Somatic Mutations in Cancer. Nucleic Acids Res. 39, D945–D950 (2011). 9. Mino-Kenudson, M. Programmed cell death ligand-1 (PD-L1) expression by 41. Budczies, J. et al. Optimizing panel-based tumor mutational burden (TMB) immunohistochemistry: could it be predictive and/or prognostic in non-small cell measurement. Ann. Oncol. 30, 1496–1506 (2019). lung cancer? Cancer Biol. Med. 13, 157–170 (2016). 42. Buchhalter, I. et al. Size matters: dissecting key parameters for panel-based tumor 10. Tsao, M. S. et al. PD-L1 immunohistochemistry comparability study in real-life mutational burden analysis. Int. J. Cancer 144, 848–858 (2019). clinical samples: results of blueprint phase 2 project. J. Thorac. Oncol. 13, 43. Fancello, L., Gandini, S., Pelicci, P. G. & Mazzarella, L. Tumor mutational burden 1302–1311 (2018). quantification from targeted gene panels: major advancements and challenges. J. 11. Kerr, K. M. et al. Programmed death-ligand 1 immunohistochemistry in lung Immunother. Cancer 7, 183 (2019). cancer: in what state is this art?. J. Thorac. Oncol. 10, 985–989 (2015). 44. Heydt, C. et al. Analysis of tumor mutational burden: correlation of five large gene 12. Hirsch, F. R. et al. PD-L1 immunohistochemistry assays for lung cancer: results panels with whole exome sequencing. Sci. Rep. 10, 11387 (2020). from phase 1 of the blueprint PD-L1 IHC assay comparison project. J. Thorac. 45. Xu, Z. et al. Assessment of tumor mutation burden calculation from gene panel Oncol. 12, 208–222 (2017). sequencing data. Onco Targets Ther. 12, 3401–3409 (2019). 13. Gibney, G. T., Weiner, L. M. & Atkins, M. B. Predictive biomarkers for checkpoint 46. Riaz, N. et al. Tumor and microenvironment evolution during immunotherapy inhibitor-based immunotherapy. Lancet Oncol. 17, e542–e551 (2016). with Nivolumab. Cell 171, 934–949 (2017). 14. Cristescu, R. et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade- 47. Carbone, D. P. et al. First-line Nivolumab in stage IV or recurrent non-small-cell based immunotherapy. Science 362, eaar3593 (2018). lung cancer. N. Engl. J. Med. 376, 2415–2426 (2017). 15. Auslander, N. et al. Robust prediction of response to immune checkpoint 48. Hellmann, M. D. et al. Genomic features of response to combination immunotherapy blockade therapy in metastatic melanoma. Nat. Med. 24, 1545–1549 (2018). in patients with advanced non-small-cell lung cancer. Cancer Cell 33, 843–852 (2018). 16. Jiang, P. et al. Signatures of T cell dysfunction and exclusion predict cancer 49. Rizvi, H. et al. Molecular determinants of response to anti-programmed cell death immunotherapy response. Nat. Med. 24, 1550–1558 (2018). (PD)-1 and anti-programmed death-ligand 1 (PD-L1) blockade in patients with 17. Tumeh, P. C. et al. PD-1 blockade induces responses by inhibiting adaptive non-small-cell lung cancer profiled with targeted next-generation sequencing. J. immune resistance. Nature 515, 568–571 (2014). Clin. Oncol. 36, 633–641 (2018). 18. Roy, S. & Trinchieri, G. Microbiota: a key orchestrator of cancer therapy. Nat. Rev. 50. Tan, H., Bao, J. & Zhou, X. Genome-wide mutational spectra analysis reveals Cancer 17, 271–285 (2017). significant cancer-specific heterogeneity. Sci. Rep. 5, 12566 (2015). 19. Routy, B. et al. The gut microbiota influences anticancer immunosurveillance and 51. Sondka, Z. et al. The COSMIC Cancer Gene Census: describing genetic dysfunc- general health. Nat. Rev. Clin. Oncol. 15, 382–396 (2018). tion across all human cancers. Nat. Rev. Cancer 18, 696–705 (2018). 20. Routy, B. et al. Gut microbiome influences efficacy of PD-1-based immunotherapy 52. Sachs, M. C. plotROC: a tool for plotting ROC curves. J. Stat. Softw. 79, 2 (2017). against epithelial tumors. Science 359, 91–97 (2018). 53. Chen, H. & Boutros, P. C. VennDiagram: a package for the generation of highly- 21. Le, D. T. et al. PD-1 blockade in tumors with mismatch-repair deficiency. N. Engl. J. customizable Venn and Euler diagrams in R. BMC Bioinforma. 12, 35 (2011). Med. 372, 2509–2520 (2015). 54. Martínez-Pérez, E., Molina-Vila, M. A. & Marino-Buslje, C. Metadata record for the 22. Le, D. T. et al. Mismatch repair deficiency predicts response of solid tumors to PD- manuscript: Panels and models for accurate prediction of tumor mutation burden 1 blockade. Science 357, 409–413 (2017). in tumor samples. figshare https://doi.org/10.6084/m9figshare14074451 (2021). Published in partnership with The Hormel Institute, University of Minnesota npj Precision Oncology (2021) 31 E Martínez-Pérez et al. 8 55. Lawrence, M. et al. Software for computing and annotating genomic ranges. PLoS ADDITIONAL INFORMATION Comput Biol. 9, e1003118 (2013). Supplementary information The online version contains supplementary material 56. Durinck, S., Spellman, P. T., Birney, E. & Huber, W. Mapping identifiers for the available at https://doi.org/10.1038/s41698-021-00169-0. integration of genomic datasets with the R/Bioconductor package biomaRt. Nat. Protoc. 4, 1184–1191 (2009). Correspondence and requests for materials should be addressed to M.A.M.-V. or 57. Lawrence, M., Gentleman, R. & Carey, V. rtracklayer: an R package for interfacing C.M-B. with genome browsers. Bioinformatics 25, 1841–1842 (2009). 58. Kent, W. J. et al. The human genome browser at UCSC. Genome Res. 12, 996–1006 Reprints and permission information is available at http://www.nature.com/ (2002). reprints 59. Ito, K. & Murphy, D. Application of ggplot2 to pharmacometric graphics. CPT Pharmacomet. Syst. Pharm. 2, e79 (2013). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. ACKNOWLEDGEMENTS E.M.P. and C.M.B. are CONICET researchers. This work was partially funded by The Agencia Nacional de Promoción Científica y Tecnológica. Argentina, Grant/Award Number: Préstamo BID PICT 1787 2014–2017. All the icons included in Supplemen- Open Access This article is licensed under a Creative Commons tary Fig. 1 were taken from https://melbournechapter.net. The picture of Attribution 4.0 International License, which permits use, sharing, immunologic treatment was taken from https://commons.wikimedia.org/wiki/File: adaptation, distribution and reproduction in any medium or format, as long as you give CAR_T-cell_Therapy.svg. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless AUTHOR CONTRIBUTIONS indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory E.M.P., M.A.M.V., and C.M.B. conceived and designed the study; E.M.P. conducted the regulation or exceeds the permitted use, you will need to obtain permission directly experiments; E.M.P., M.A.M.V., and C.M.B. analyzed data and wrote the manuscript. from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. COMPETING INTERESTS The authors declare no competing interests. © The Author(s) 2021 npj Precision Oncology (2021) 31 Published in partnership with The Hormel Institute, University of Minnesota