We previously retrospectively analyzed two large prospectively collected datasets for early-stage patients presenting for primary tumor resection at Memorial Sloan Kettering Cancer Center: 740 patients with lung adenocarcinoma (MSK-LUAD) (11) and 1,157 patients with colon adenocarcinoma (MSK-COAD) (12). In the present work we revisited these datasets to further probe relationships between TMB, FGA, drugs, and outcomes.

The Cancer Genome Atlas (TCGA) (16) contains clinicopathologic and bulk RNA sequencing gene expression data for over 20,000 primary cancer samples from 33 cancer types. We focused on 18 TCGA solid tumor cohorts, corresponding to 6,488 patients (see Supplementary Table S1 for a breakdown of total samples by cancer type).

The MSK-LUAD cohort featured next-generation sequencing (MSK-IMPACT) (25) of tumor samples removed during surgery, which enabled calculation of TMB and FGA (among other genomic factors). The MSK-COAD cohort featured immunohistochemical determination of DNA mismatch repair (MMR) status [either deficient (dMMR) or proficient (pMMR)] (26). MMR status correlates with TMB, where dMMR patients have high TMB and pMMR patients have low TMB (27, 28). TCGA cohorts featured TMB and FGA values for patient tumor samples, which were obtained through cBioPortal (29–31). TMB in TCGA was determined by the total number of nonsynonymous coding mutations.

Distributions for TMB and FGA were calculated separately for the MSK and TCGA cohorts by plotting log TMB or unnormalized FGA vs. TMB/FGA rank. Since TMB distributions are known to differ by cancer type (32), we calculated separate distributions for the different TCGA cohorts. Consistent with existing literature and clinical algorithms for immunotherapy response prediction (33, 34), these distributions were not linear but rather featured breakpoints dividing the distributions into high (hypermutators) (35) and low TMB and high, mid, and low FGA. Three breakpoints were computed for TMB, while two were computed from FGA. Because TMB was binarized, the choice between the 2nd and 3rd breakpoints was manual; the 1st breakpoint was discarded as it represented samples with virtually no mutations and so was too low to be the hypermutation breakpoint. Breakpoints were calculated using the R package Segmented (36) with the mixed-linear segmented models.

The multivariable regression models developed previously (11–13) were employed to calculate predicted 5-year survival estimates for model patients for specific TMB and FGA values.

Bulk RNA sequencing data was analyzed for the same TCGA patients described above. HTSeq files containing raw RNA counts were downloaded from the National Institutes of Health/National Cancer Institute Genomic Data Commons public database (37). Standard DESeq2 workflow (38) was employed to normalize the data and calculate differential gene expression for the TMB and FGA categories described above (i.e., high vs. low TMB and high vs. low FGA). In brief, genes were prefiltered by cancer to only include genes where at least 50% of the samples had non-zero count data. Standard differential expression using DESeq2 was carried out and log₂ fold-change estimate shrinkage was carried out using the APEGLM method (39). P-values were automatically adjusted in DESeq2 for multiple testing correction. Count data was normalized using the variance-stabilized transform algorithm, which decreases the variance of counts as compared to the mean of counts for a given gene, along with outputting expression in log scale.

A list of 79 genes corresponding to target receptors (or subunits where applicable) for common anesthetic and analgesic drugs was compiled, listed in Supplementary Table S2.

Computational estimates for immune cell-type proportions for each patient sample in the TCGA cohorts were previously computed by Wang and colleagues (40) using the default LM22 immune composition expression matrices (41) in CIBERSORT (42).

Increased intraoperative opioid dose was previously found to be associated with worse overall survival (OS) in the MSK-LUAD cohort (11), while intraoperative ketamine (vs. either dexmedetomidine or no adjunct) (11) and intraoperative ketorolac were associated with higher recurrence-specific survival (RSS) (13). TMB was previously found to interact with opioid dose to modify the opioid-OS association (11), though not with ketorolac (13). Further analysis of this dataset in this work demonstrated that TMB interacted with ketamine to modify the ketamine-RSS association (see below). FGA was previously found to modify the opioid-OS association (11).

To probe these interactions further, we calculated predicted 5-year survival estimates for model patients with different TMB values, derived from the TMB distribution for MSK-LUAD cohort (Figure 1A) as explained in Figure 1B (TMB values: 1.8 (17th percentile), 3.9 (41st percentile), 10.8 (83rd percentile), and 16.7 (92nd percentile)). The 10.8 muts/mb cutoff separating hypermutators from non-hypermutators is consistent with the FDA-defined cutoff (10 muts/mb) above which patients are considered hypermutators and approved to receive pembrolizumab (20). TMB interacted with ketamine such that increased TMB [hypermutation (HM)] opposed ketamine's effect to reduce recurrence (Figure 1C, TMB = 1.8: predicted 5-year RSS = 0.94 with ketamine vs. 0.78 for no adjunct; 3.9: 0.93 vs. 0.78; 10.8: 0.87 vs. 0.78; 16.7: 0.80 vs. 0.77). TMB did not significantly interact with ketorolac, illustrated by the overlapping curves for TMB values (Figure 1D, TMB = 1.8: predicted 5-year RSS 0.81 with ketorolac vs. 0.71 without; 3.9: 0.80 vs. 0.71; 10.8: 0.79 vs. 0.70; 16.7: 0.78 vs. 0.69). TMB increase magnified the opioid-OS association such that OS further worsened at higher opioid dose (Figure 1E).

For FGA, the median of each group (high, mid, and low) was selected from the distribution (Figure 1F) for further analysis; these were 0.0135 (31st percentile), 0.1442 (78th percentile) and 0.3975 (96th percentile). Like TMB, FGA increase magnified the opioid-OS association such that OS further worsened at higher opioid dose (Figure 1G). Note that the error on these predicted values increases for higher opioid dose and FGA (Supplementary Figure S1) owing to the decreasing number of patients at both higher MME and high FGA, and so the observed near 100 percent mortality at high FGA and high opioid dose should be interpreted in that context, with a focus on the trend of increasing mortality for a given opioid dose as FGA increases.

In the MSK-COAD cohort, increased intraoperative opioid dose was associated with lower risk of recurrence (12) (highlighting that oncological effects of opioids may depend on cancer type/subtype). For this analysis, a correlate of TMB, MMR status, was used to determine hypermutation high-TMB status (Figure 1H). Hypermutation (dMMR) further decreased recurrence at higher opioid dose. (Figure 1I, modified from the original figure).

Distributions for TMB and FGA by cancer type are shown for the TCGA lung adenocarcinoma cohort (TCGA-LUAD) in Figures 2A,B and for the TCGA colon adenocarcinoma cohort (TCGA-COAD) in Figures 3A,B (and for all 18 TCGA cancer types in Supplementary Figure S2, S3 for TMB and FGA respectively with breakpoints provided in Supplementary Table S3, S4 for TMB and FGA respectively). For TCGA-LUAD, the breakpoint separating high and low TMB was 23.57 muts/mb. The breakpoints separating low to mid FGA and mid to high FGA were 0.318 and 0.494, respectively. For TCGA-COAD, the TMB breakpoint was 5.9 muts/mb, while the breakpoints separating low to mid FGA and mid to high FGA were 0.004 and 0.465, respectively. Note that these values are determined from whole-exome sequencing and do not correspond exactly to those calculated from a targeted oncology panel like MSK-IMPACT [the harmonization of different panel assays and whole exome sequencing is the subject of ongoing research (43, 44)].

Differential gene expression for TCGA-LUAD is illustrated in the volcano plot Figure 2C for high vs. low TMB and in Figure 2D for high vs. low FGA. For the TMB differential expression, 6,469 of a total 31,521 genes (20.5%) were significantly differentially expressed, with 18 of the 59 expressed anesthetic receptor genes (30.5%) being differentially expressed. For the FGA differential expression, 16,082 total genes were differentially expressed (51.0%) with 38 of anesthetic target receptors (64.4%) being differentially expressed. Differential gene expression for TCGA-COAD is shown similarly in Figures 3C,D. For TCGA-LUAD, 12,492 of 28,711 (43.5%) total genes were found to be differentially expressed for TMB with 32 of the 58 (55.2%) expressed anesthetic receptors being differentially expressed. In the FGA comparison, 7,479 (26.0%) of all genes were differentially expressed, and 12 of the 58 (20.7%) of the anesthetic genes had differential expression. Note that while this demonstrates overrepresentation of the receptor genes, more relevant to the current analysis is that there is differential expression of these genes (and the specific genes that are differentially expressed).

Anesthetic/analgesic target receptor genes are highlighted throughout. In TCGA-LUAD, an increase in TMB or FGA generally causes a subset of same genes to be upregulated (e.g., OPRD1, GRIN2D, GRIN1, GRIN2C) and downregulated (e.g., GRIN2A, GABRP, CHRNA2), though some may switch, such as CHRNA9. However, in TCGA-COAD, most of the genes switch up- or downregulation when comparing an increase in TMB vs. FGA (e.g., GABRB3, HTR7, CHRNA7, GABRE, CHRNG, GLRA2). It should also be noted that FGA for TCGA-LUAD and TMB for TCGA-COAD showed an increase in differential expression in both the background genes, as well as the target receptor genes, as compared to TMB and FGA, respectively. As a check, differential expression between FGA medium and high was carried out for TCGA-COAD and TCGA-LUAD (Supplementary Figure S4). This showed in general less differential expression than the FGA high vs. low comparisons, and no additional genes were found to be enriched.

Dot plots showing differential target receptor gene expression for all 18 TCGA cancer types for both TMB and FGA are displayed in Figure 4. For TMB (Figure 4A), the genes that are upregulated significantly most frequently are GRIN1 (5 cancer types), CHRNA5 (4), PTGS2, HTR3A, HTR1D, GRIN3A, GRIN2D, and GRIN2B (3). Similarly, the genes that are downregulated significantly the most frequently are HTR2B, GRIN2A (7), HTR3A, GABRP, GABRG3, GABRA2, and CHRNA1 (5). With the same determination for FGA (Figure 4B), the upregulated genes are GABRE (8), OPRK1, GRIN2D, GABRA3, CHRNA5 (7), HTR3A, GRIN1, ADRA2B, and GABRQ (6), and the downregulated genes are TLR4 (12), HTR7 (11), CNRIP1 (10), PTGS1, CNR2 (9), HT2RB (8), HT2RA, and GABRP (7). When comparing pan-cancer if genes were alternatively down- and upregulated by increasing TMB or FGA, 7 TCGA cancer types were found to have at least 5 significant genes that were common between the TMB differential expression experiment and the FGA differential expression experiment. The 7 types were BRCA, COAD, GBM, LUAD, LUSC, STAD, and UCEC. For these types, GBM had 0/10 genes (0%) swap, BRCA had 6/34 (17.6%), LUAD had 3/15 (20%), LUSC had 2/5 (40%), STAD had 12/21 (57.1%), UCEC had 17/21(81.0%), and COAD had 8/8 (100%) common genes swap whether they were up- or downregulated between high/low TMB/FGA. Differential expression results for all cancer types may be found in the Supplementary Data File.

Consistent with existing literature (45–47), variation in immune composition correlated with TMB and FGA. TCGA-LUAD immune composition is shown in Figures 5A,B for TMB and FGA variation, respectively, while TCGA-COAD immune composition is shown in Figures 5C,D, respectively. For hypermutated tumors (high TMB), most cancer types showed a significant increase in CD8+ T-cells and macrophages, while noting that many subtypes showed exception to this trend. For high FGA, CD8+ T-cells, macrophages (all subtypes but specifically M2 macrophages), and CD4+ memory resting T-cells are significantly decreased. One exception to this rule is TCGA-PRAD (figures and data for all cancer types can be found in Supplementary Data).

The anesthetic/analgesic target receptor genes found to be significantly correlated with TMB or FGA in TCGA-COAD and TCGA-LUAD were assessed for correlation with the 22 immune cell populations in each cancer type. Correlation was assessed using Spearman's correlation and p-values were computed using the permutation test from the R package coin (48). Comparison was done using Spearman's correlation so as to make no assumption of the linearity when analyzing the correlation of immune cell composition with the receptor genes. Data is provided in Supplementary Data File for the correlation between all significant genes and immune types with resulting p-values that were adjusted for multiple hypothesis testing within each immune cell type for each cancer type. Sample correlation plots are displayed in Figure 6. Data is shown for the expression of 3 example genes that were significantly differentially expressed for the TCGA-COAD TMB and TCGA-LUAD FGA categories (GRIN3A, OPRD1, and GRIN1). There is no clear functional significance at this time for selection of these genes; they are used to illustrate how differential expression may be affecting different immune cells in similar or different ways in various tumor types. These were compared to the presence of M2-type macrophages and CD8+ T-cells, also significantly different when comparing TCGA-COAD TMB and TCGA-LUAD FGA categories. As shown, GRIN3A is highly expressed in both TCGA-LUAD and TCGA-COAD across a wide variety of immune cells (Figures 6A–D, p < 1e-11 for all 4). However, OPRD1 shows a significant decreasing expression (p = 0.006) in TCGA-LUAD with an increase in M2 macrophages, but no change in expression in TCGA-COAD (p = 0.97). In contrast, GRIN1 shows no significant change in expression based on CD8+ T-cell level in TCGA-LUAD (p = 0.58), but an increase in expression is found with an increase in CD8+ T-cells in TCGA-COAD (p = 0.0004).

Heatmaps displaying Spearman's rho for TCGA-LUAD and TCGA-COAD are shown for the significant genes for TMB and FGA in Figure 7. These plots are both clustered by cell type and gene expression. As evidenced, specific genes that may be descriptive of a broader cell type such as GRIN3A or HTR7 for macrophages are generally upregulated in samples with estimated higher macrophages, while other genes that may be differentially expressed in a cancer type due to underlying genomic events in the tumor such as hypermutation or genomic alteration can cause significant differences to appear.