Characterization of gene expression profiles of esophageal cancer patients with different nonsynonymous tumor mutation burden.

BACKGROUND: Nonsynonymous tumor mutation burden (NSTMB) could affect the prognosis of esophageal cancer (EC) patients, but differentially expressed genes between EC patients with different NSTMB have not been explored. Our study aimed to compare differentially expressed genes between EC patients with different NSTMB (high vs. low).
METHODS: RNA-seq data for EC patients were downloaded from The Cancer Genome Atlas (TCGA). The edgeR package was used to identify differentially expressed genes between patients with different NSTMB. Cell type identification by estimating relative subsets of known RNA transcripts (CIBERSORT) software was employed to underscore immune cell differences between patients with different NSTMB.
RESULTS: In total, we discovered 2215 differentially expressed genes between patients with different NSTMB, among which 842 genes were upregulated and 1373 downregulated in patients with high NSTMB. The differentially expressed genes were enriched in pathways such as heme binding and structural molecule activity. We built a logistic model that may be used to predict patients' NSTMB. We found that tumors with high NSTMB had a significantly higher percentage of resting natural killer (NK) cells than those with low NSTMB (P = 0.028). The percentages of regulatory T (Treg) and CD8+ T cells were also higher in those with high NSTMB, although it was not statistically significant (P = 0.064 for Treg cells and P = 0.12 for CD8+ T cells).
CONCLUSIONS: NSTMB may cause changes in gene expression and immune cell infiltration in EC patients, and affect the overall survival of EC patients. KEY POINTS: Significant findings of the study This study found differentially expressed genes and differences in infiltration of immune cells between esophageal cancer (EC) with different NSTMB. What this study adds This study highlights differences between EC patients with different NSTMB.

Introduction

Esophageal cancer (EC) was the sixth leading cause of death among malignant tumors worldwide in 2018. Current treatment mainly focuses on surgery, chemotherapy and radiotherapy. Nonetheless, the five‐year overall survival rate remains low (less than 10% for overall survival and 15%–20% for post esophagectomy survival). In recent years, immunotherapy has shown some promising effects for these patients. Some studies have reported that many cancer patients with high tumor mutation burden (TMB) can benefit from immunotherapy, , , suggesting the existence of differences in immune landscape associated with TMB. However, the immune landscape of EC patients between tumors with different nonsynonymous tumor mutation burden (NSTMB) remains to be determined.

In this study, we explored the differentially expressed genes between tumors with different NSTMB distinct from the traditional perspective, which mainly focuses on PD‐1/PD‐L1. We downloaded RNA‐seq data from The Cancer Genome Atlas (TCGA) and sorted patients into two groups: one group with high NSTMB and another with low NSTMB. We analyzed differentially expressed genes between these two groups and performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to uncover the underlying pathways. We built a logistic model using five differentially expressed long intergenic noncoding RNAs (lincRNAs) that may be used to predict NSTMB in EC patients. Finally, we explored immune infiltration differences between these two groups. Overall, our study characterized the gene expression profiles of EC patients with different NSTMB.

Methods

NSTMB calculation method

TCGA data were obtained from the Genomic Data Commons(GDC)application program interface(API)with the Bioconductor R package “TCGAbiolinks”. Aggregated whole‐exome sequencing somatic mutation information from the project “TCGA‐ESCA” was downloaded as a Mutation Annotation Format (MAF) file with the VarScan pipeline used to perform somatic mutation calling. Mutations with a coverage <10× or a variant allele fraction <5% were filtered. NSTMB was calculated as the total number of nonsynonymous somatic mutations divided by the length of the whole exome (35 M). Variant classifications with high/moderate variant consequences (“Frame_Shift_Del”, “Frame_Shift_Ins”, “Splice_Site”, “Translation_Start_Site”, “Nonsense_Mutation”, “Nonstop_Mutation”, “In_Frame_Del”, “In_Frame_Ins” and “Missense_Mutation”) were considered nonsynonymous mutations.

Overall survival

Overall survival (OS) was defined as the time from diagnosis to death or the last follow‐up.

Original RNA‐seq dataset

We downloaded EC patients' RNA‐seq count data from TCGA with the following filtering criteria: (i) the project was “TCGA‐ESCA”; (ii) the workflow type was “HTSeq‐Counts”; (iii) the data category was “Transcriptome Profiling”; (iv) the data type was “Gene Expression Quantification”; and (v) the sample type was “Primary Tumor”. We obtained information on 161 patients. Among these patients, the original site of two patients was the stomach, and one patient did not have NSTMB information, so these patients were excluded. In total, we obtained data on 158 patients. After the data were downloaded, we extracted the mRNA expression information from the files and merged them with the clinical information.

Data analysis

Differentially expressed genes (DEGs) were analyzed using the edgeR (version 3.22.5) package. Genes with a false discovery rate (FDR) <0.05 and |fold change| >2 were considered DEGs.

GO and KEGG analyses

After the DEGs were obtained, we assessed these genes using Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) and selected pathways that were significant (FDR < 0.05).

Logistic model

Since the model should be used to predict NSTMB so that it could be predicted when sequencing cannot be performed, we selected genes in a stricter way. We considered FDR < 0.01 and |fold change| > 5 to be significant. After that, we selected lincRNAs that are logistically significant with NSTMB (Fig 3[a]). We used these lincRNAs to build a logistic model.

Figure 3

Logistic model used to predict nonsynonymous tumor mutation burden (NSTMB). (a) Workflow of selecting logistically significant lincRNAs. (b) Receiver operating curves of the logistic model for diagnosing patients with high NSTMB.

Cell infiltration prediction

Cell type identification by estimating relative subsets of known RNA transcripts (CIBERSORT) software (https://cibersort.stanford.edu/) was used to predict the infiltration rate of immune cells in tumors, with 100 permutations. Data used to run CIBERSORT were downloaded from TCGA. Before running CIBERSORT, we used variance modeling at the observation level (voom) to render our RNA‐seq data more similar to microarray data. All filtering criteria were the same as those used to find differentially expressed genes except that the workflow type was “HTSeq‐FPKM”. The LM22 signature matrix was used as a reference.

Statistical analysis

R (version 3.5.3) and SPSS software (version 24.0) were used to perform statistical analyses. Kaplan‐Meier analysis was used to compare the prognosis of patients. The differences in immune cell infiltration between patients with different NSTMB were compared using the Wilcoxon test. The edgeR package in R was used to find DEGs. Except for genes selected for the logistical model (with P < 0.01), P < 0.05 or FDR < 0.05 was considered significant.

Results

Baseline characteristics of patients

In total, we enrolled 158 patients into our study. To determine the correlation between NSTMB and overall survival, we divided the patients into two groups according to the cutoff value determined by the maximally selected log‐rank statistic with the R package “survminer”. Patients with NSTMB higher than 3.228 571 were considered the NSTMB high group; patients with a tumor NSTMB value lower than 3.228 571 were considered the NSTMB low group. Finally, we obtained 96 EC patients with low NSTMB and 62 with high NSTMB. Among the enrolled patients, all of the tumor primary sites were the esophagus. The detailed baseline characteristics of the patients are listed in Table 1.

Table 1

Baseline characteristics of patients selected

	NSTMB low	NSTMB high
Gender
Male	82	53
Female	14	9
Primary site
Esophagus	96	62
Other	0	0
Vital status
Alive	64	32
Dead	32	30
Age at initial pathological diagnosis	59.87 ± 11.85	65.97 ± 11.55
(mean ± SD)
Pathological stage
I	11	5
II	42	26
III	27	20
IV	6	2
Histological type
Adenocarcinoma	36	42
Squamous cell carcinoma	60	20
History of neoadjuvant treatment
Yes	0	0
No	96	62

Patients with low NSTMB have better prognosis than those with high NSTMB

Since NSTMB can affect the immunogenicity of a tumor, which may influence prognosis, we sought to determine whether differences in prognosis existed between the two groups. Patients with low NSTMB had a significantly better prognosis than those with high NSTMB (χ2 = 6.764, log‐rank P = 0.009: see Fig 1[a], and Table 2, 3). Then, we asked whether the prognosis of patients with different NSTMB was different in different pathological stages. In stages II and III, patients with high NSTMB had a worse overall survival than patients with low NSTMB (χ2 = 4.134, log‐rank P = 0.042 for stage II, and χ2 = 4.914, log‐rank P = 0.027 for stage III: see Fig 1[c] and [d]). However, in stages I and IV, the overall survival of patients with different NSTMB was not different statistically (χ2 = 2.200, log‐rank P = 0.138 for stage I, and χ2 = 0.102, log‐rank P = 0.749 for stage IV: see Fig 1[b] and [e]).

Figure 1

The overall survival of patients with different nonsynonymous tumor mutation burden (NSTMB) in different groups. (a) The overall survival of all patients with different NSTMB. (b) The overall survival of patients with different NSTMB in stage I. (c) The overall survival of patients with different NSTMB in stage II. (d) The overall survival of patients with different NSTMB in stage III. (e) The overall survival of patients with different NSTMB in stage IV. : low NSTMB, : high NSTMB.

Table 2

Mean survival time (months) of patients with different nonsynonymous tumor mutation burden (NSTMB)

			95% Confidence interval
NSTMB	Estimate	Std. error	Lower bound	Upper bound
Low	43.969	4.756	34.647	53.292
High	28.208	4.022	20.325	36.092
Overall	38.669	3.618	31.578	45.761

Table 3

Median survival time (months) of patients with different nonsynonymous tumor mutation burden (NSTMB)

			95% Confidence interval
NSTMB	Estimate	Std. error	Lower bound	Upper bound
Low	45.367	11.642	22.548	68.185
High	20.000	2.684	14.739	25.261
Overall	28.500	3.614	21.417	35.583

Comparison of gene expression between patients with different NSTMB

Since there is a difference in the survival time between patients with different NSTMB, we assessed differences in gene expression between them. We downloaded RNA‐seq data from TCGA and compared gene expression patterns using the edgeR (version 3.22.5) package in R. Genes with an FDR < 0.05 and |fold change| >2 between the two groups were considered differentially expressed. We found that compared with patients with low NSTMB, there were 842 genes that were upregulated and 1373 genes that were downregulated in patients with high NSTMB (Fig 2[a]).

Figure 2

RNA‐seq data for patients with different nonsynonymous tumor mutation burden (NSTMB). (a) Gene expression of tumors with different NSTMB. Red, significantly upregulated gene expression in tumors with high NSTMB; green, significantly downregulated gene expression in tumors with high NSTMB. (b) Significantly enriched GO terms in tumors with low NSTMB. (c) Significantly enriched KEGG pathways in tumors with low NSTMB.

We next examined the functions of these genes by submitting the genes to DAVID (https://david.ncifcrf.gov/) for GO and KEGG analyses. As shown in Fig 2, genes downregulated in tumors with high NSTMB were enriched in 12 GO and 4 KEGG terms; genes upregulated in tumors with high NSTMB were enriched in only one GO term (GO: extracellular space, FDR = 0.04837) and no KEGG terms.

Differentially expressed lincRNAs between tumors with different NSTMB

LincRNAs are involved in many biological processes in tumors. Since lincRNAs can be used to predict patient prognosis, , , we asked whether there are differentially expressed lincRNAs that could be used to predict NSTMB. In total, we selected five lincRNAs that are significantly associated with NSTMB. These lincRNAs are LINC00200, LINC01206, LINC01043, LINC01019 and LINC01580. Interestingly, we found that these lincRNAs are specifically expressed in the testis (Fig S1). Detailed information on these lincRNAs is listed in Table S1. [Correction added on 3 July 2020, after first online publication: the variable ‘LINC01026’ has been corrected to ‘LINC01206’ throughout the article.]

We then calculated the predicted probability and used it to predict NSTMB. The P‐value was calculated as:

Correct index = sensitivity + specificity – 1.

We set the cutoff value by which the correct index calculated as above was the largest. After calculation, we chose the cutoff value as 0.466. A P‐value greater than 0.466 could be considered high NSTMB in patients. As shown in Fig 3(b) and Table 4, the area under the curve (AUC) is 0.772, suggesting that our model may be used to predict NSTMB.

Table 4

Statistical analysis of variables used to predict nonsynonymous tumor mutation burden (NSTMB)

				95% CI
Test variable	AUC	SE	P‐value	Lower bound	Upper bound
LINC00200	0.375	0.044	0.008	0.289	0.462
LINC01206	0.338	0.043	0.001	0.254	0.423
LINC01043	0.387	0.044	0.016	0.300	0.474
LINC01019	0.417	0.045	0.078	0.328	0.505
LINC01580	0.490	0.048	0.835	0.396	0.584
Pre‐1	0.772	0.037	0.000	0.699	0.845

AUC, area under the curve; CI, confidence interval; SE, standard error.

Immune landscape of tumors with different NSTMB

NSTMB may affect the prognosis of EC patients, and we accordingly evaluated immune infiltrate differences between patients with different NSTMB. We downloaded LM22 from CIBERSORT (https://cibersort.stanford.edu/) and obtained the expression profile of 22 immune cells (Fig 4[a]). We used CIBERSORT to calculate the immune infiltration percentage of each cell type based on the LM22 signature matrix. We excluded patients with CIBERSORT P‐values greater than 0.05 and ultimately obtained 76 patients, among whom 51 had low NSTMB and 25 had high NSTMB. To our surprise, we did not observe many differences between the two groups, except for resting natural killer (NK) cells (Fig 4[b] and [c]). The percentages of regulatory T (Treg) and CD8+ T cells were also higher in the NSTMB high group than in the NSTMB low group, although the differences were not statistically significant (P = 0.064 for Treg cells, P = 0.12 for CD8+ T cells: see Fig 4([b]). Correlation analysis showed that resting NK cells had a negative correlation with activated NK cells, and CD8+ T cells had a strong correlation with activated memory CD4+ T cells (Fig 4[d]).

Figure 4

High nonsynonymous tumor mutation burden (NSTMB) is associated with a higher percentage of resting NK cell infiltration. (a) Gene expression of 22 immune cells based on the LM22 file. The rows are genes, and the columns are cells. (b) Violin plot showing the CIBERSORT estimation of immune cell infiltration based on the LM22 file. Red, high NSTMB; blue, low NSTMB. Data were analyzed using the Wilcoxon test. (c) Heatmap showing the CIBERSORT estimation of immune cell infiltration based on the LM22 file. Low, patients with low NSTMB. High, patients with high NSTMB () low, () high. (d) Pearson correlation analysis showing the correlation of the 22 immune cell types in these patients.

Discussion

In this study we first compared the overall survival of patients with different NSTMB and found that patients with high NSTMB had worse overall survival than patients with low NSTMB. We then compared the overall survival in different pathological stages and found that in stages II and III, patients with high NSTMB had a worse overall survival than that in patients with low NSTMB. However, the prognosis was not different between patients with different NSTMB in stages I and IV. This may be caused by the small number of patients in stages I and IV (Table 1). Therefore, more patients are needed in further studies to verify this hypothesis. We then compared the differences between EC patients with different NSTMB (high vs. low) tumor from a new perspective. We found 2215 genes that were differentially expressed in tumors with different NSTMB. We performed GO, KEGG, and CIBERSORT analyses based on these genes. We selected genes to make a logistic model that could be used to predict patients' NSTMB in a stricter way (selection criteria: FDR < 0.01, |fold change| >5) since this model may be used to predict NSTMB. Interestingly, lincRNAs for the prediction in the model are specifically expressed in the testis, suggesting that these lincRNAs are cancer‐testis genes. The differentially expressed genes are enriched in terms such as oxygen binding, integral component of membrane. We also found that the proportion of resting NK cells was significantly higher in patients with high NSTMB than those with low NSTMB. Moreover, the proportion of Treg cells was higher in tumors with high NSTMB, although the difference was not statistically significant (P = 0.064).

As shown in Fig 2, 12 GO and four KEGG terms were enriched in genes downregulated in tumors with high NSTMB, suggesting that NSTMB may affect the prognosis of EC patients in these ways. Nevertheless, elucidation of the exact mechanisms requires further investigation. As shown in Fig 4(b), the proportion of resting NK cells was significantly higher in tumors with high NSTMB, suggesting that higher NSTMB could inhibit NK cell activation, thus causing the poorer prognosis of patients with high NSTMB. The proportion of Treg cells in tumors with high NSTMB was also higher (P = 0.064). Since resting NK and Treg cells are related to immune inactivation to a certain degree, we speculated that the antitumor immunity in tumors with high NSTMB was weaker than that in tumors with low NSTMB, causing patients with high NSTMB to have a relatively poorer prognosis.

Unlike our observation, Li et al. found that the proportion of many immune cells, including CD8+ T cells and some kinds of macrophages, was different between tumors with different NSTMB in non‐small cell lung cancer (NSCLC) patients. This may be due to the fact that our discovery was based on EC, yet our study also showed that the proportion of CD8+ T cells was higher in patients with high NSTMB, even though this difference was not statistically significant (P = 0.12, Fig 4[b]). In addition, we found that the proportion of resting NK cells was higher in tumors with high NSTMB, suggesting that high NSTMB may cause NK cell inactivation to be common in malignant tumors, but more data for other kinds of tumors are needed. Since NSTMB is correlated with the amount of tumor neoantigens, which is correlated with the immunogenicity of tumors, we speculated that the relatively higher proportion of CD8+ T cells in patients with high NSTMB is caused by the higher number of tumor neoantigens, which in turn activated the body's immune response more strongly. Interestingly, we found that the proportion of CD8+ T cells had a positive correlation with that of activated memory CD4+ T cells. CD4+ Tcells might help the differentiation of tumor antigen‐associated CD8+ T cells, verifying our speculations that the antitumor immune response of patients with high NSTMB is more strongly activated. However, the proportion of Treg cells in the high NSTMB group was also higher (P = 0.064), suggesting that the immune suppressive response of patients with high NSTMB was also higher. Taken together, our results showed that both the antitumor and immune suppressive responses might be higher in patients with high NSTMB.

To our knowledge, this is the first study to uncover expression differences between tumors with different NSTMB in EC patients. We also built a logistic model that may be used to predict the EC patients' NSTMB, and this may be used to guide further treatments, predict patients' prognosis, or guide proper follow‐up intervals for each patient; however, more research is needed when used in clinical settings. The proportion of resting NK cells was higher in tumors with high NSTMB, suggesting that we may use NK cell activation strategy to obtain a better prognosis. Regardless, further research is needed.

Our study has some limitations. First, our sample size was small. We obtained data for only 158 patients. Increasing the number of patients in this analysis may reveal different results. Second, as all of our studies were based on bioinformatic analyses without experimental verification, more experimental results are needed in further research. Third, part of our discussions was based on results with P‐values greater than but close to 0.05, including those for CD8+ T cells and Treg cells. However, these P‐values were calculated based on current data and methods mentioned in our study, and they were very close to 0.05, indicating that although their P‐values were more than 0.05, it does not mean they have no difference in real situations. Finally, the result of the immune cell infiltration rate was obtained based on a CIBERSORT estimation, which may vary slightly from real situations.

In summary, we found that NSTMB is associated with different gene expression and immune cell infiltration in EC patients, and we created a model that may be used to predict EC patients' NSTMB. These findings highlight the differences between EC patients with different NSTMB, indicating that further research is needed to uncover the mechanism by which NSTMB functions in EC patients.

Disclosure

The authors declare no conflicts of interests.

NSTMB, nonsynonymous tumor mutation burden

Figure S1 Five lincRNAs are specifically expressed in testis. (a) LINC01206. (b) LINC00200. (c) LINC01043. (d) LINC01019. (e) LINC01580.These data were analyzed in NCBI (https://www.ncbi.nlm.nih.gov/).

Click here for additional data file.

Table S1 Position of logistically significant lincRNAs

Click here for additional data file.