Exploration of genes and tumor infiltrating lymphocytes in female lung adenocarcinoma microenvironment that predicted prognosis.

ABSTRACT: The tumor microenvironment has an important impact on tumor growth, invasion, metastasis, anti-tumor immune tolerance, and prognosis. The present study aimed to explore female lung adenocarcinoma microenvironment-associated tumor infiltrating lymphocytes (TILs) and genes that predict prognosis in The Cancer Genome Atlas (TCGA) database. Gene expression profiles of female patients with lung adenocarcinoma were downloaded from TCGA. Base on the CIBERSORT algorithm, we determined the fractions of TILs. By applying the ESTIMATE algorithm, immune scores and stromal scores were derived. According to the immune and stromal scores, we categorized the female patients with lung adenocarcinoma into high and low score groups. We also identified the fractions of TILs and differentially expressed genes (DEGs) that were significantly related with prognosis. The proportion of M1 macrophages was significantly negatively related to overall survival in female patients with lung adenocarcinoma. There were 269 upregulated genes and 35 downregulated genes both in immune scores and stromal scores. PTPRC (protein tyrosine phosphatase receptor type C) and GIMAP6 (GTPase, IMAP family member 6) were not only hub genes, but also were significantly related to overall survival in the Kaplan-Meier Plotter online and TCGA databases. In summary, our study provided new insight into the tumor microenvironment-related cellular and molecular mechanisms of women with lung adenocarcinoma. The results will be useful for future clinical studies.

Introduction

Lung cancer is the leading cause of cancer death in worldwide.[ Based on its biology, therapy, and prognosis, the World Health Organization (WHO) divides lung cancer into 2 major classes: Non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC).[ Adenocarcinoma is the most common subtype of NSCLC.[ The 5-year survival rate of lung adenocarcinoma is only 18%.[ Good prognostic factors in lung adenocarcinoma include good performance status, diagnosis at early stage, no significant weight loss (<5%), and female sex.[ Different clinical characteristics and outcomes of lung adenocarcinoma are related to sex.[ Moreover, women with lung adenocarcinoma tend to be younger and have more advanced stage disease.[ Prognostic factors and pathogenesis of women with lung adenocarcinoma are subject to intensive research.

Yu et al showed that microRNAs (miRNAs) miR-144-5p and miR-218-3p are tumor suppressor miRNAs that are significantly related with prognosis in female patients with lung adenocarcinoma.[ NF1 (Neurofibromin 1) is a tumor suppressor gene, and mutation of NF1 in female patients with lung adenocarcinoma resulted in significantly short disease-free survival and overall survival.[ Previous studies have shown that tumor infiltrating lymphocytes (TILs) and the extent of immune scores or stromal scores contribute significantly to prognosis in several tumors.[ However, the utility on TILs, immune scores, and stromal scores of female patients with lung adenocarcinoma have not been investigated in detail.

In the present study, based on CIBERSORT (Cell Type Identification By Estimating Relative Subsets Of RNA tissues using Expression data) and ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor) algorithms, we compiled a list of microenvironment-associated TILs and genes that predicted prognosis from The Cancer Genome Atlas (TCGA) database. Moreover, we validated these genes in an independent database: Kaplan–Meier Plotter online.

Material and methods

Data collection

To analyze the TILs and genes associated with immune or stromal scores in women with lung adenocarcinoma, tumor tissue and adjacent normal tissue data were downloaded from TCGA database (https://portal.gdc.cancer.gov/), on March 2020. The workflow type was set as Fragments Per Kilobase Million (FPKM) to meet the needs of subsequent analysis. The corresponding patient follow-up data was also acquired from TCGA.

Evaluation of tumor-infiltrating immune cells

CIBERSORT (http://cibersort.stanford.edu/) can accurately estimate TILs in tumor samples profiled by microarrays or RNA-Seq.[ In this study, by applying the CIBERSORT algorithm, the normalized gene expression data downloaded from TCGA were used to infer the relative proportions of 22 types of TILs. A P value less than .05 was used a filter in the follow-up analysis. The proportions of the TILs and clinical follow-up were evaluated.

Evaluation of immune scores and stromal scores

The ESTIMATE algorithm is used to predict tumor purity, and calculate immune and stromal scores.[ In this study, we applied the normalized gene expression data downloaded from TCGA to infer the immune and stromal scores in women with lung adenocarcinoma. The correlation of immune or stromal scores and overall survival were analyzed using a Kaplan–Meier plot.

Comparison of the gene expression profile with immune scores and stromal scores

We divided the cases into the top and bottom halves by their immune and stromal scores. We compared the RNA-seq data between the top and bottom immune or stromal scores, using the R software (Limma package). The intersection of differentially expressed genes (DEGs) between the top immune scores and top stromal scores were used in subsequent analysis. The intersection of DEGs between the bottom immune scores and bottom stromal scores were also used in subsequent analysis.

Enrichment analyses of DEGs

To carry out a functional enrichment analysis of the DEGs, we applied the online tool DAVID (https://david.ncifcrf.gov/). Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes pathways that were enriched for the DEGs were identified using the WEB-based Gene Set Analysis Toolkit, P < .05.

Analysis of protein-protein interaction networks

The STRING database (https://string-db.org/cgi/input.pl) provides experimental and predicted interaction information between gene-encoded proteins. We uploaded the DEGs to the STRING database, and then downloaded identified the PPI network. The downloaded PPI network was manipulated using Cytoscape software. Hub genes are those that encode proteins with the most connections with the other proteins in the PPI network. A plug in of Cytoscape, CytoHubba, was used to find the hub genes in the PPI network.

Correlation of DEGs expression with overall survival

We divide the expression of DEGs in tumor patients into high-expression group and low-expression group based on the average value. The association between DEGs and corresponding clinical follow-up were analyzed using Kaplan–Meier curves and evaluated using the log-rank test. P < .05 was considered to be statistical significantly.

Validation in the Kaplan–Meier plotter database

Kaplan–Meier plotter (http://kmplot.com) is an independent online database that includes both clinical and expression data.[ We identified the DEGs with prognostic significance using Kaplan–Meier plotter.

Statistical analysis

All statistical analysis was performed using the R software. TILs in normal and cancer tissues were compared using the Wilcoxon signed rank test. The DEGs screening criteria were: adjusted P values <.05 and absolute value of Log2Fold Change (FC) >2. Kaplan–Meier analysis and the log-rank test were used to evaluate the relationship between TILs or DEGs and overall survival (OS). A P value <.05 was considered statistically significant.

Results

Composition of immune cells in tumor and normal tissues

Total of 254 female patients diagnosed with lung adenocarcinoma were included in this study. The median age at diagnosis was 63 years old (range from 33 to 87 years old). The clinical features grouped by age of the 254 female lung adenocarcinoma patients were shown in Table 1. Using the CIBERSORT algorithm, there were 242 patients that were selected for research into the composition of their immune cells.

Table 1

Clinical features grouped by age.

Group	<65 Y	>=65 Y	χ2	P
T1&2	90	124	6.620	.008
T3&4	24	13
N1&2	95	115	0.003	.550
N3&4	16	19
M0	69	94	0.627	.314
M1	6	5
StageI&II	90	109	0.014	.514
StageIII&IV	24	28

We first investigated the differences among 22 subpopulations of immune cells between tumor and normal lung tissues. The proportions of naïve B cells, memory B cells, plasma cells, CD4+ memory T cells, follicular helper T cells, regulatory T cells, gamma delta T cells, monocytes, M1 macrophages, and resting dendritic cells were significantly higher in lung tumor tissues than in normal tissues. The proportions of CD4- memory T cells, resting natural killer (NK) cells, M0 macrophages, M2 macrophages, resting mast cells, and neutrophils were significantly lower in lung tumor tissues than in normal tissue, as shown in Figure 1A. The correlations between the proportions of the 22 types TIL were weak (Fig. 1B).

Figure 1

The differences among 22 subpopulations of immune cells between tumor and normal lung tissue. (A) Violin plots illustrating whether there was a significantly different fraction between tumor and normal lung tissue. (B) The correlations between the proportions of the 22 TIL subsets were weak.

The proportion of M1 macrophages correlates negatively with overall survival

The correlations between the proportions of TILs and clinicopathological characteristics or overall survival were also analyzed. The proportions of M0 macrophages, M2 macrophages, neutrophils, and resting NK cells were significantly lower in N0&N1 lymph nodes than in N2&N3 lymph nodes. The proportions of plasma cells and CD4+ memory T cells were significantly higher in N0&N1 lymph nodes than in N2&N3 lymph nodes, as shown in Figure 2B–2G. The proportion of gamma delta T cells was significantly higher in patients aged ≥ 65 years old than in patients aged <65 years old, as shown in Figure 2H. More importantly, the proportion of M1 macrophages was significantly negatively correlated with overall survival, as shown in Figure 2A.

Figure 2

The correlation between the proportions of tumor infiltrating lymphocytes (TILs) and the patients’ clinicopathological characteristics or overall survival. (A) The proportion of M1 macrophages was significantly negatively related with overall survival. (B) The fractions of Plasma cells were significantly higher in N0&N1 lymph nodes than in N2&N3 lymph nodes. (C) The fractions of M2 macrophages were significantly lower in N0&N1 lymph nodes than in N2&N3 lymph nodes. (D) The fractions of neutrophils were significantly lower in N0&N1 lymph nodes than in N2&N3 lymph nodes. (E) The fractions of Resting natural killer (NK) cells were significantly lower in N0&N1 lymph nodes than in N2&N3 lymph nodes. (F) The fractions of plasma cells were significantly higher in N0&N1 lymph nodes than in N2&N3 lymph nodes. (G) The fractions of CD4+ memory T cells were significantly higher in N0&N1 lymph nodes than in N2&N3 lymph nodes. (H) The fraction of T gamma delta cells was significantly higher in patients ≥65 years old than in patients <65 years old.

Comparison of gene expression profiles with immune and stromal scores

We analyzed the correlation of the immune and stromal scores with overall survival. Based on the immune and stromal scores, we divided the female patients with lung cancer into top and bottom halves (high scores vs low scores). Median overall survival was longer in cases with low immune scores than that in high immune scores (775.7 d vs 750.1 d, P = .35 in log-rank test), as shown in Figure 3A. Median overall survival was longer in cases with low stromal scores than that in high stromal scores (814.9 d vs 708.3 d, P = .81 in log-rank test), as shown in Figure 3B. However, neither the immune scores nor the stromal scores were statistically correlated with overall survival.

Figure 3

The correlation of immune scores or stromal scores with overall survival. (A) Median overall survival was longer in cases with low immune scores than in cases with high immune scores, P = .347. (B) Median overall survival was longer in cases with low stromal scores than in cases with high stromal scores P = .81.

To explore the DEGs between high immune/stromal scores and low immune/stromal scores, we constructed a heatmap of DEGs associated with immune scores and stromal scores, as shown in Figure 4A and 4B. There were 573 upregulated DEGs and 105 downregulated DEGs between the high and low immune scores. There were 655 upregulated DEGs and 82 downregulated DEGs between high and low stromal scores, as shown in Figure 4C and D. We selected 269 upregulated genes (intersection of immune scores and stromal scores) for follow-up studies. We also choose 35 downregulated genes (Intersection of immune scores and stromal scores) for follow-up studies.

Figure 4

The differentially expressed genes (DEGs) between high immune/stromal scores and low immune/stromal scores. (A) A heatmap of gene expression according to the immune scores. (B) A heatmap gene expression according to the stromal scores. (C) There were 269 upregulated genes that were associated with both the immune scores and stromal scores. (D) There were 35 downregulated DEGs that were associated with both the immune scores and stromal scores.

The common up and downregulated genes in immune scores and stromal scores were used as DEGs for enrichment analyses. According to the GO hierarchy, we grouped the 304 DEGs into BP (biological process), MF (molecular function), and CC (cellular component) categories. The most enriched GO terms under BP were “immune response,” “inflammatory response,” “adaptive immune response,” and “innate immune response.” The most significant GO terms under MF were “chemokine binding” and “IgG binging.” The most important GO terms under CC were “collagen trimer,” “membrane raft,” and “extracellular space.” According to the Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis, among the DEGs the most enriched in pathways were “phagosome” and “cell adhesion molecules,” Figure 5.

Figure 5

Enrichment analyses of differentially expressed genes (DEGs). Gene ontology (GO) categories:(A) BP, biological process; (B) CC, cellular component. (C) MF, molecular function; (D) KEGG: Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis.

Protein-Protein Interaction network among DEGs

The common up and down-regulated genes in the immune scores and stromal scores were used as DEGs to construct a protein-protein interaction (PPI) network. CytoHubba (a plug-in of Cytoscape) was used to find the major genes that have most interactions. The upregulated Hub-genes were PTPRC, IL10, ITGAM, TLR8, CCR5, CSF1R, CD163, C3AR1, FCER1G, FPR1, FCGR3A, IGLL5, CD86, CCR2, LILRB2, CD180, CXCL11, CCL1, PIK3AP1, GIMAP6, ADORA3, P2RY12, CXCL10, CCR1, CXCL13, P2RY13, FPR3, CNR2, CCR8, LY86, FCGR1A, CCL19, CCL13, CXCL9, CXCR5, CCR4. The downregulated hub-genes were NPW, F2, and GNRH2, as shown in Figure 6.

Figure 6

The protein-protein interaction network and hub genes. The small circle represents the hub genes. The large circle plus the small circle are the differentially expressed genes. Red represents upregulation and green represents downregulation.

Prognosis-related DEGs

Kaplan–Meier survival analysis was applied to compare overall survival in data from TCGA. Among the 304 DEGs, 30 genes were significantly related to overall survival, as shown in Table 2.

Table 2

Thirty DEGs significantly related with overall survival in TCGA database.

Group	Differentially expressed genes
Upregulated genes associated with good prognosis	LILRB5, B4GALNT4, GNRH2, CR1, CCR2, CD33, MPEG1, CCR6, GIMAP8, CEBPE, PLEK, CRHBP, NLRP3, ARHGAP15, FCRL3, CLEC10A, SIGLEC1, TREM2, HPR, CD200R1, FYB1, PTPRC, SLCO2B1, NAIP, RCSD1, IKZF1, GIMAP6, EOMES, APOC4-APOC2
Downregulated genes associated with good prognosis	F2

To validate the association of the 30 genes with OS, we analyzed them using an independent online database, Kaplan–Meier Plotter. There were data for 286 women with lung adenocarcinoma in the Kaplan–Meier Plotter online database. Among the 30 genes, 8 DEGs were validated significantly related with OS, both in TCGA database (Fig. 7) and Kaplan–Meier Plotter online database (Fig. 8), but this difference is barely noticeable. PTPRC (protein tyrosine phosphatase receptor type C), and GIMAP6 (GTPase, IMAP family member 6) were not only identified as hub genes, but also were significantly related to OS in the Kaplan–Meier Plotter online and TCGA databases.

Figure 7

Eight differentially expressed genes were significantly related to overall survival, in TCGA database.

Figure 8

Eight differentially expressed genes were significantly related to overall survival, in Kaplan–Meier Plotter online database.

Discussion

Understanding the tumor microenvironment-related TILs and genes remains fundamental for the elaboration of the pathogenesis, prognostic factors, and novel therapeutic methods for women with lung adenocarcinoma. In this study, we explored tumor microenvironment-related TILs and genes that contribute to OS in women with lung adenocarcinoma in data from TCGA. Using the CIBERSORT algorithm, we determined the proportions 22 immune cell populations that acted as TILs in women with lung adenocarcinoma and evaluated their clinical follow-up. Applying the ESTIMATE algorithm allowed us to calculate the immune and stromal scores for women with lung adenocarcinoma. We then categorized the patients into high and low immune or stroma score groups. Furthermore, we explored the DEGs whose expression was significantly associated with prognosis between the high and low immune/stromal score groups. Finally, we validated the prognosis-related DEGs using an independent online database: Kaplan–Meier plotter.

CIBERSORT showed that the proportion of M1 macrophages was significantly higher in tumors than in normal tissues. Furthermore, M1 macrophage levels were closely associated with poor prognosis. Macrophages, as key regulators of host immunity, play an important role in the tumor microenvironment.[ Macrophages can differentiate into 2 phenotypes, namely the classically activated macrophages (M1) and alternatively activated macrophages (M2).[ A study of NSCLC[ showed that the number of M1 macrophages was significantly higher in tumor tissue than in healthy lung tissue from the control group. Consistent with this study, our results revealed that M1 macrophage levels were significantly higher in females with lung adenocarcinoma than in normal tissue. One possible reason might be that monocyte chemotactic protein-2 (CCL2, also known as C-C motif chemokine ligand 2), which is produced by tumor cells, can recruit circulating M1 macrophages to the tumor site.[ In many tumors, the proportion of M1 macrophages is related to prognosis. In colorectal cancer, tumor infiltrating M1 macrophages were associated with poor prognosis.[ In breast cancer, the tumor-associated macrophage density correlated with poor prognosis.[ The role of M1 macrophages in NSCLC patient’ survival was a matter of controversy in previous studies. In a study by Liu et al, neither M1 nor M2 macrophages influenced the prognosis of patients with lung adenocarcinoma.[ In contrast, Jurgita et al showed that M1 macrophages were associated with improved NSCLC patient survival.[ The controversy might reflect the fact that during tumor progression, the macrophage phenotype changes from the classically to the alternatively activated form.[ However, there have been few studies about the relationship between M1 macrophages and the OS of women with lung adenocarcinoma. Our results showed that M1 macrophage levels were related to poor prognosis in women with lung adenocarcinoma. This result might provide new target for immunotherapy of women with lung adenocarcinoma.

Based on the ESTIMATE algorithm, we identified 304 DEGs from comparison of high versus low immune or stromal scores, and found that most of the DEGs were involved in the tumor microenvironment. PTPRC and GIMAP6 were not only identified as hub genes, but also were significantly related to OS in the Kaplan–Meier Plotter online and TCGA databases.

Protein tyrosine phosphatase receptor type C (PTPRC) is a biomarker of T cells.[ In gastric cancer, PTPRC was highly expressed, and the overexpression of PTPRC could promote the development of gastric cancer.[ PTPRC was also overexpressed in renal cell carcinoma, and participate in the progress of cell adhesion.[ In triple-negative breast cancer, PTPRC acted as an immune gene, was associated with overall survival of patients.[ The expression of PTPRC is related to the prognosis of lung adenocarcinoma,[ as confirmed by the results of the present study. GIMAP6 is expressed at high levels in cells of the immune system, and has been associated with immunological functions, such as thymocyte development, apoptosis of peripheral lymphocytes, and T helper cell differentiation.[ Studies on GIMAP6 in lung adenocarcinoma are rare; however, both PTPRC and GIMAP6 might be potential therapeutic targets.

The present study had limitations. The fraction of the TILs, the expression of DEGs, the identities of the hub genes, and the association between the hub genes and OS, all require experimental validation.

In conclusion, based on the CIBERSORT algorithm, this study demonstrated that low levels of tumor infiltrating M1 Macrophages were associated with better prognosis in women with lung adenocarcinoma. Applying the ESTIMATE algorithm, we identified that PTPRC and GIMAP6 were hub-genes that were significantly related with overall survival. We believe that this study will provide new insights into the cellular and molecular mechanisms of lung adenocarcinoma in women.

Acknowledgments

We thank the authors who provided the TCGA public data set and Kaplan–Meier Plotter online database.

Author contributions

Funding acquisition: Hao Wang.

Investigation: Song Wei.

Resources: Song Wei.

Software: Lijun Wang, Hao Wang, Song Wei.

Supervision: Zhihong Zhang.

Validation: Zhihong Zhang.

Visualization: Zhihong Zhang.

Writing – original draft: Lijun Wang.

Writing – review & editing: Hao Wang, Zhihong Zhang.