Sage Journals: Discover world-class research

Abstract

BACKGROUND:

Leukocyte infiltration plays an critical role in outcome of various diseases including Lung adenocarcinoma (LUAD).

OBJECTIVES:

To understand the genetic and epigenetic factors affecting leukocyte infiltration and identification and validation of immune based biomarkers.

METHOD:

Correlation analysis was done to get the associations of the factors. CIBERSORT analysis was done for immune cell infiltration. Genetic and epigenetic analysis were performed. Cox regression was carried out for survival.

RESULTS:

We categorized the TCGA-LUAD patients based on Leukocyte fraction (LF) and performed extensive immunogenomic analysis. Interestingly, we showed that LF has a negative correlation with copy number variation (CNV) but not with mutational load. However, several individual genetic mutations, including KRAS and KEAP1, were significantly linked with LF. Also, as expected, patients with high LF showed significantly increased expression of genes involved in leukocyte migration and activation. DNA methylation changes also showed a strong association with LF and regulated a significant proportion of genes associated with LF. We also developed and validated an independent prognostic immune signature using the top six prognostic genes associated with LF.

CONCLUSION:

Together, we have identified clinical, genetic, and epigenetic variations associated with LUAD LF and developed an immune gene-based signature for disease prognostication.

Keywords

Lung adenocarcinoma Immune landscape leukocyte fraction prognosis

Abbreviations

LUAD	Lung adenocarcinoma
LF	Leukocyte fraction
CNV	Copy number variation
HR	Hazard ratio
IRS	Immune risk score

1. Introduction

Lung adenocarcinoma (LUAD) is a significant cause of cancer related fatalities globally. Even with improvements in therapy protocols, the 5-year survival rate remains abysmal [1, 2]. Late diagnosis of patients ( $\sim$ 66 percent) and high relapse rate are majorly responsible for the high mortality of LUAD [2]. Recent developments in immunotherapy have shown encouraging results, and it has been integrated into the standard care regimen of LUAD treatment [3]. Nevertheless, the response to immunotherapy varies from patient to patient, and only a small population of patients are benefited [3, 4, 5]. While lung cancer has traditionally been described as low immunogenic, emerging evidence suggests that inadequate immune response is triggered by active immune evasion mechanisms [6]. The response to immunotherapy is multidimensional and relies on a more complex immune microenvironment than just the expression of checkpoint genes, mutation burden, and antigen presentation [7, 8]. Recently, several genetic correlations, such as NF1 loss in GBM and increased infiltration of macrophages, TP53 mutation, and decreased infiltration of T cells, have been associated with the immune composition of tumors [9, 10]. The leukocyte fraction (LF) is the total immune cell content of the tumors. The presence of different types of immune cells has also been correlated with the prognosis and reaction of patients to therapy [11]. It is therefore, crucial to understand the genetic and epigenetic makeup of the tumor cells in order to understand the microenvironment of the tumor. Here, in this manuscript, we have performed mutation, copy number aberration, and methylation analysis to understand the genetic and epigenetic changes associated with LF of the LUAD patients. We have also used the statistical methods to develop an immune risk score (IRS) for prognostic prediction of LUAD patients.

2. Material and methods

2.1 Patients sample LF and CIBERSORT analysis

The clinical, expression, DNA methylation, CNV, and mutation data of LUAD patients were downloaded from TCGA (https://portal.gdc.cancer.gov/). Second cohort of LUAD RNA-seq data with clinical follow-up was downloaded from GSE66728. Another dataset with 172 LUAD patients set with count value and survival data was downloaded from East Asian LUAD study (EGAD00001004421) and used as second testing cohort. For survival analysis, TCGA data was used as the training set, and GSE66728 and EGAD00001004421 data was used as the testing set. The LF data was used as described earlier(Supplementary Table 1) [12]. In short, the DNA methylation probes, which were most differentially methylated in leukocyte and healthy lung tissue, were used to calculate the LF in each patient’s sample. To calculate the fraction of the immune cell type, CIBERSORT (https://cibersort.stanford.edu/) analysis was done as described earlier [13]. The LM22, 22 immune cell signatures were used as the signature file, and absolute and relative modes with no quantile normalization were used for analysis.

2.2 Expression analysis

To identify the genes differentially expressed between tumor and normal, a $t$ -test was carried out. Genes with abs (log2(FC)) $\geqslant$ 2 and FDR $\leqslant$ 0.05 were considered as differentially expressed. Also, Pearson correlation analysis was carried out to identify the genes whose expression was significantly correlated with the LF.

2.3 Mutation analysis

MAF file for LUAD samples were obtained from the GDC data portal (https://portal.gdc.cancer.gov/). To identify the genes whose mutation was correlating with the LF, Point-biserial correlation was computed between the mutation data and LF for all the samples. All non-silent mutations were used for analysis.

2.4 Methylation analysis

The standard method of preprocessing Illumina 450K array DNA methylation data was used [14, 15, 16]. Additionally, probes were filtered based on standard deviation (sd) and those with sd $>=$ 0.2 were used for further analysis. To identify the probes correlating with LF, Pearson’s correlation analysis was carried out. And to identify the differentially methylated (DM) probes, Limma was used on quantile normalized and log transformed beta values.

2.5 Copy number analysis

Masked segment mean for all the LUAD samples was obtained from the GDC data portal. To recognize the regions that are significantly amplified or deleted across the samples in low and high LF cohorts, GISTIC 2.0 algorithm (https://www.genepattern.org/modules/docs/GISTIC_2.0) was applied [17]. Besides the arm-level amplification and deletion $q$ -values obtained from GISTIC for each cohort, Fisher’s exact test was carried out to identify significant amplification and deletion for each arm at the cohort-level with FDR correction. All the arms with $q$ -value $\leqslant$ 0.1 were considered significant.

2.6 Survival analysis

Cox regression analysis was performed to identify the genes associated with survival prediction. The genes with at least $\pm$ 0.3 correlation coefficient and significantly differentially expressed in LUAD compared to normal were used as input in Cox regression analysis, and genes with a significant hazard ratio ( $p$ -value $\leqslant$ 0.01) were selected for further analysis. The immune risk score was calculated using the expression and Cox regression coefficient, as described earlier [18]. Cox regression analysis was performed in the survival package of R, and Kaplan-Meier graphs were plotted in GraphPad 8.3.

2.7 Other statistical analysis

For comparison of two groups, the two-sided t-test was performed in GraphPad8.3, and tests with $p$ -value $\leqslant$ 0.05 were considered significant. For multiple tests, FDR correction was done, and FDR $\leqslant$ 0.05 was considered significant. For three or more groups comparison Kruskal-Wallis test was performed, and $p$ -value $\leqslant$ 0.05 was considered significant.

3. Results

3.1 Characterization of immune nature of LUAD with high and low LF

Infiltration of leukocytes in the tumor environment plays a key role in the development and progression of cancers, including LUAD [12, 19]. While LUAD’s LF is the highest of all cancers, LUAD is also a rather heterogeneous tumor [12]. To understand the immune nature of LUAD, patients were divided into low (LF $<$ 0.2), moderate (LF 0.2–0.4), and high (LF $>$ 0.4) groups. (Fig. 1A and Supplementary Table 1). The LF includes all immune cells in the tumor milieu, suggesting that the higher LF should also have higher infiltration of all immune cells. Interestingly, however, lymphocyte cells (CD4, CD8, Treg, etc. cells) and macrophages of type M1 showed higher infiltration in high-LF samples, while APC cells such as Dendritic cells and M2 macrophages showed reduced infiltration in high-LF tumors (Fig. 1B, Supplementary Table 2). Furthermore, Kaplan-Meier analysis showed that the LF of LUAD tumors is not associated with the outcome of the disease (Fig. 1C and Supplementary Table 1).

Figure 1.

LF of LUAD patients: A. LF is plotted against TCGA-LUAD patient numbers. The patients with LF $<$ 0.02 were considered as low LF, Patients with 0.2–0.4 were considered as moderate LF, and patients with $>$ 0.4 LF were considered as high LF. B. The Spearman correlation between LF and relative infiltration of immune cells (from CIBERSORT analysis) was performed, and the coefficient of cell types with a significant correlation was plotted in a heatmap. ${}^{**}<$ 0.005, ${}^{***}<$ 0.0001. C. A log-rank test was performed to check the survival difference among the three groups of patients. KM plot is plotted to show that there was no significant association of LF with survival. D. The association of LF of LUAD patients with gender was tested with the $\chi^{2}$ test and found significant. A bar plot was developed to show the percent of male and female patients in low, moderate, and high LF groups. E. The normalized fraction of Th1, Th2, and Th17 was calculated and plotted as a boxplot. Kruskal-Wallis test was performed to obtain the $p$ -value. F. The infiltration of Treg cells from CIBERSORT analysis was plotted in three LF groups and plotted as a boxplot. Kruskal-Wallis test was performed to obtain the $p$ -value. G. Intratumoral heterogeneity was found to be significantly associated with LF in LUAD. A Box plot was plotted to show the association and $p$ -value obtained from the Kruskal-Wallis test. H. The expression of immune checkpoint genes was compared between low and high LF groups and plotted in the form of a volcano plot. I. A heatmap representing the expression of immune checkpoint genes associated with high LF. The expression of these genes was also compared with normal lung and shown. DE Differentially expressed (FDR $<$ 0.05). J. Three clinically most important immune checkpoint genes expression was plotted and in normal, low, moderate, and high LF patients. Kruskal-Wallis test was performed to obtain the $p$ -value.

Figure 2.

Identification of genetic aberration with LF in LUAD: A. The correlation matrix showing correlation coefficient of mentioned genomic elements. B. The total mutation burden and CNV fraction of LUAD patients was plotted in low, moderate, and high LF patients and shown in the box plot. Kruskal-Wallis test was performed to obtain the $p$ -value. C. The genes with at least five percent mutation rate were correlated with LF, and gene mutations with significant correlation were plotted as a heatmap. Red bars represent the mutation in a patient in given genes. D. Bar plot showing the mutation frequency of most commonly mutated genes in low and high LF samples. E. Bar plot showing significantly enriched GO terms for genes whose mutation was significantly correlated with LF. Darker colour represents a lower $p$ -value.

Next, we looked at the effect of different clinicopathological features on the LF of LUAD patients. Interestingly, we found that the gender of the patient has a significant association with the LF. A higher proportion of female patients were associated with a higher leukocytes fraction (Fig. 1D, Supplementary Table 3). However, other clinical parameters, such as age and stage of patients, did not show significant association with LF(Supplementary Table 3). CD4 $+$ T cells (Th cells) play a crucial role in the innate and adaptive immune response of the tumor microenvironment. The primary role of the Th cells depends on the differentiation of the Th cells in the Th1, Th2, Th17, and Treg cells. We checked the enrichment of the different Th cell lineages in patients with low, moderate, and high LFs. Interestingly, the infiltration of Th1 and Treg cells was significantly higher in high LF patients, and the infiltration of Th17 cells was higher in low LF samples (Fig. 1E and F). In addition, intratumoral heterogeneity was significantly positively associated with the LF of LUAD patients (Fig. 1G).

More recently, immune checkpoint genes have become a crucial determinant of the therapy outcome due to advances in immunotherapy protocols [20]. The association of stimulatory and inhibitory checkpoint molecules with LF of patients was calculated, and the expression of these genes were compared in high vs. low LF using the $t$ -test (Fig. 1H). The analysis showed that 30 out of 41 examined immune checkpoint regulators, which included both co-stimulatory ( $n=$ 16) and co-inhibitory ( $n=$ 14) genes, were positively associated with the LF (Fig. 1H). Interestingly, we also found 7 down-regulated and 11 up-regulated immune checkpoint genes in LUAD compared to normal (Fig. 1I). More importantly, the checkpoint genes currently being targeted for therapy, PDL1 (CD274), PDCD1, and CTLA4 were significantly up regulated in LUAD and showed a strong association with LF (Fig. 1J).

Taken together, the LUAD samples are heterogeneous in terms of the LF, and the enrichment of cells varies by cell type. Gender and intertumoral heterogeneity also play an essential role in the LF of LUAD patients.

3.2 LF of LUAD is associated with specific genetic aberrations and genomic instability

The genetic makeup of cancer is not only a determining factor in the treatment option and the response to therapy, but it also affects the immune contexture [21]. We performed a correlation analysis to understand the relationship between genetic aberrations and the LF. There was a significant positive correlation between the LF and the number of segments and the modified fraction, suggesting that copy number variation is associated with LUAD immune cell infiltration (Fig. 2A). However, there was no correlation with LF in other genetic aberrations. To validate the finding, we compared the mean mutation and CNV fraction in LUAD patients with low, moderate, and high LFs. The analysis showed that there was no significant difference in the overall mutation burden, although the CNV fraction was significantly lower in the high LF (Fig. 2B).

Genetic aberrations were associated with the presence or absence of immune cells. For example, VHL and STK11 mutations are associated with lower macrophages infiltration [22]. In the same way, PI3KCA, MET mutations are associated with CTL infiltration [22]. We have therefore identified recurrent gene mutations (mutation frequency $>$ 5 percent) that are significantly associated with the LF. Interestingly, most mutations were associated with a lower LF. Only 7 genes mutations (PLCL1, SLITRK4, CDH6, LYST, FRY, KMT2C, FLG2) had a significant association with a high LF (Fig. 2C). Importantly, the KRAS, KEAP1, key genes in LUAD, were associated with a lower LF (Fig. 2C). The mutation frequency of the most commonly mutated genes in low and high LF group is compared in Fig. 2D. Gene ontology analysis revealed that gene mutations associated with the LF regulate different developmental processes (Fig. 2E).

Significant number of CNV in low and high LFs have been identified. At chromosome arm level, the low LF group showed significantly higher copy number variations in 1p, 1q, 6p, 7p, 8q and 17q and lower copy number of 3p, 3q, 9p, 19p and 19q (Fig. 3A). Many focal level differences between high LF and low LF of CNVs were also identified (Fig. 3B and C). The common genes associated with cancer with significant focal copy number aberration are given in Fig. 3B and C.

Figure 3.

Copy number changes associated with LF: A. Arm level copy number analysis was performed and plotted. ${}^{*}$ indicates the significantly amplified/deleted changes from the standard copy number. Green dots: significant amplification in low LF vs high LF group. Red: Significant deletions in low vs high LF group. B. Significant focal amplification of the various chromosomes. C. Significant focal deletion of various chromosomes.

Figure 4.

Methylation and expression landscape of LF in LUAD: A. The CpGs methylation was correlated with LF and CpGs with $\pm$ 0.2 correlation coefficient and FDR $<$ 0.05 were selected and plotted as a heatmap. The methylation of CpGs was also compared with Normal and plotted. DM: Differentially methylated (FDR $<$ 0.05). B. The number of hyper and hypomethylated CpGs, CpGs associated with PcGs, LncRNA, and all other genomic locations were plotted in the form of a bar plot. C. The PcGs whose CpGs methylation was associated with LF were used for GO analysis, and significantly enriched GO terms were plotted. D. Expression of genes was compared between high and low LF LUAD patients and plotted as a volcano plot. FDR cut off $<$ 0.05, Fold change- at least 2-fold. E. The overlap of genes with differentially expressed between high and low LF and normal and LUAD was plotted as Venn diagram. F. A heatmap was plotted to show the expression pattern of the genes whose expression was associated with LF in LAUD. The expression of these genes was also compared with normal and shown. DE: Differentially expressed, FDR $<$ 0.05. G. The genes associated with LF in LUAD were used for Metascape analysis. The obtained GO network is shown, and important GO terms are mentioned. H. The relationship of methylation and expression of genes associated with LF was calculated and plotted as a starburst plot.

Figure 5.

Development of immune risk signature for LUAD: A. Cox regression analysis was performed using the genes associated with LF and plotted as a volcano plot. The Red genes are negatively associated with survival, and blue genes are positively associated with survival. B. The six selected survival genes were used for GO analysis, and enriched terms are shown in the form of a bar plot. C. A heatmap of six prognostic immune genes in normal and LUAD. The LUAD patients are arranged in increasing order of LF. D. The IRS was calculated and plotted. The patients were divided into three groups based on the IRS. E. The LUAD patients were divided into low and high risk at the median, and the Kaplan-Meier plot was made. The survival was compared using log-rank analysis. F. The patients were divided into low, intermediate, and high immune risk groups, as shown in Fig. 4D, and the Kaplan-Meier plot was made. The survival of the groups were compared using the log-rank test. G. The LUAD patients with stage I were taken and divided at the median of the IRS, and the Kaplan-Meier plot was obtained. The survival was compared using the log-rank test. H. The univariate and multivariate Cox regression analysis was performed with the given variable, and a forest plot was prepared to show the results. $P$ -value is given. In the plot dot represents the HR and error bars represent 95% CI. I. The CIBERSORT analysis of patients with high and low IRS was performed, and a $t$ -test was performed to identify the significantly enriched cell type between two groups. The significantly enriched cell types proportion is plotted as a bar diagram.

Figure 6.

Validation of IRS in both the independent cohorts: A to E for GSE66728 and F to H for EGAD00001004421 cohort. A. Heatmap of six prognostic immune genes in the testing cohort and normal lung. B. The IRS was calculated using the expression of six genes in the testing cohort. The patients were divided into low, intermediate, and high immune risk groups. C. The testing cohort patients were divided into a high and low immune risk group, and the Kaplan-Meier plot was made. A Log-rank test was performed to compare survival. D. The training group patients were divided into three groups, and the Kaplan-Meier plot was made. A log-rank test was done to compare survival. E. Univariate and multivariate Cox regression analysis was performed in the testing cohort, and a forest plot was made to show the HR (CI). $P$ -values are given. In the plot dot represents the HR and error bars represent 95% CI. F. The heatmap of the expression of the six IRS component genes. G. IRS was calculated in second training cohort and plotted. H. Second testing cohort patients were divided into low ( $n=$ 133) and high ( $n=$ 39) immune risk group and KM plot was plotted.

3.3 Epigenetic regulation of LF in LUAD

DNA methylation changes have a significant impact on the behaviour of the tumor, including the immune microenvironment [23]. In order to understand the role of DNA Methylation in the LF, the beta-values of CpG from TCGA-LUAD 450K-methylation data were correlated with the LF. A selection criterion (described in the method section) was used to identify 787 CpGs, the methylation level of which was strongly correlated with the LF (303 negative correlation and 484 positive correlation) (Fig. 4A and B). A major fraction (95.42%) of these CpGs have also been significantly differentially methylated in LUAD compared to normal lung samples (Fig. 4A) suggesting that changes in methylation pattern during carcinogenesis is an important factor in determining the immune cell infiltration of tumors. GO analysis showed that genes whose methylation was associated with a higher LF, regulated metabolism pathway associated with transcription, and genes whose methylation was correlated with a lower LF, were involved in the regulation of developmental pathways (Fig. 4C).

Further, we compared the expression pattern of the genes in low and high leukocyte using selection criteria of the correlation coefficient, and differential expression as explained in the method section. A total of 631 genes were identified as associated with LF in LUAD (Fig. 4D). Interestingly, the majority of genes were positively associated with LF, and more than 80% of these genes ( $n=$ 518) were also differentially expressed in LUAD compared to normal lungs (Fig. 4E and F). These observations suggest that gene expression variation is also a vital determinant of the LF of LUAD. In addition, network analysis using genes associated with LF revealed enrichment of pathways associated with leukocyte activation, cytokine production, lymphocyte activation, and leukocyte migration (Fig. 4G). DNA methylation is a common mechanism for gene regulation. Therefore, we tested the relationship between the LF associated with DNA methylation and gene expression. Combining the expression and methylation gene list identified 57 hypermethylated and downregulated and 41 hypomethylated and up-regulated genes in high LF compared to low LF LUAD patients (Fig. 4H).

3.4 Development of LUAD prognostic model using LF associated genes

In our initial analysis, we have shown that the LF is not associated with survival of the LUAD patients (Fig. 1C). However, further investigation of genes whose expression was correlated with LF, led to the identification of prognostic genes in LUAD. To identify the genes which were strongly correlated with survival, unbiased Cox-regression analysis was performed as detailed in method section. The study identified 6 genes (ARNTL2, ITRIP-risk genes), (RP11-750H9.5, TSPAN32, ARHGEF6, and C11orf21-protective genes) with a strong association with survival (Fig. 5A, Supplementary Table 1). Further analysis showed that these prognostic genes regulate various immune cell regulatory pathways (Fig. 5B). Interestingly, ARNTL2 was downregulated, and all other genes were up regulated in LUAD compared to normal (Fig. 5C). Next, the expression value and Cox-regression correlation coefficient were used to calculate an immune risk score (IRS) (Material and method, and Fig. 5D). The patients were divided at the median, and the Log-rank test was performed. The analysis showed that patients with high IRS are at 1.75 times higher hazard than the patients with low IRS (Fig. 5E). Furthermore, the division of patients into three groups (high IRS $>$ 0, intermediate immune score – 1.75–0 and low immune score $<-$ 1.75) showed extremely different patients’ survival between high and low IRS patients (Fig. 5F). Prognostication of early-stage LUAD patients may play an important role in therapy decisions as opposed to later-stage patients who show poor outcomes even with multimodal therapy. Thus, we evaluated the performance of the IRS in stage I patients. We divided the stage I patients at the median of the IRS and performed Kaplan-Meier (KM) analysis and showed that stage I LUAD patients with high IRS had significantly worse survival compared to patients with low IRS (Fig. 5G). In addition, we also conducted a univariate Cox-regression analysis and showed that the IRS was a significant predictor of prognosis (Fig. 5H). In order to understand the independence of the IRS, we performed a multivariate Cox-regression analysis. The analysis established the IRS as an independent survival predictor (Fig. 5H). Further, to check the immune cells’ infiltration difference in high and low IRS patients, CIBERSORT analysis was performed, and a $t$ -test was done to identify the differentially infiltrated immune cells in both the groups. The analysis showed high enrichment of CD4 T-cells, Plasma cells, M0 Macrophages, and NK cells and low enrichment of CD8 T cells, B cells, Monocytes, M1 macrophage in High risk group compared to low risk group (Fig. 5I).

3.5 Validation of prognostic markers

It is crucial to validate the performance of prognostic markers in an independent dataset. Hence, we used another set of data set comprising 67 LUAD patients to test the prognostic strength of the IRS. First, we checked the expression of all the genes in normal ( $n=$ 6) and LUAD samples of the validation cohort. As in the training cohort, the expression of ARNTL2 was higher, and other genes were lower in LUAD compared to normal (Fig. 6A). Further, we calculated the risk score of each LUAD patient from the validation cohort similar to the training cohort (Material method section, Fig. 6B, Supplementary Table 4). Like the training cohort, the patients were divided at the median. The Kaplan-Meier analysis showed that patients with high immune scores live significantly worse than patients with low immune scores (Fig. 6C). We also divided the patients into three groups. As expected, in the validation cohort, also patients with the highest immune scores survived significant poor (Fig. 6D). Next, we performed univariate and multivariate cox regression analysis to show that IRS is independent prognosticator in validation cohorts as well (Fig. 6E, Supplementary Table 6).

Next, we used the second testing cohort of samples with 172 LUAD East Asian patients for the validation of immune risk score. First, we checked the expression of six genes in normal and LUAD samples and found that like other two cohort, expression of ARNTL2 was higher and other genes was lower in LUAD compared to normal samples (Fig. 6F). The risk score was calculated as described in method section and plotted (Fig. 6G). The Kaplan-Meier analysis was done at multiple cut off and the cut off at which low and high immune risk groups were found most significantly different was used for division (Fig. 6G and H). The analysis showed the significant poor survival of the high IRS patients compared to low IRS patients (Fig. 6H).

4. Discussion

The immune microenvironment of cancer plays a key role in patients’ behaviour towards therapy and overall outcomes [19, 24]. In our current research, multiple analyses were conducted to identify the genetic and epigenetic markers associated with the LF in LUAD tumors. Although leukocytes consist of all the different adaptive and innate immune response cells, our results have shown that the LF is positively associated with lymphocytes and M1 macrophages and negatively associated with dendritic cells, M2 macrophages, and mast cells. This result suggests that as more cells infiltrate the tumor, lymphocyte and tumor-promoting macrophages infiltrate, preferably. Other cells, such as dendritic cells, infiltrate more in tumors with a lower LF. The LF was not found to be associated with survival in LUAD tumors. Interestingly, we observed that significantly more female patients showed a high LF. This was a unique observation and needs further investigation. CD4 ${}^{+}$ T cells are essential to the regulation of the immune response. CD4 ${}^{+}$ T cells attack cancer cells either by destroying them directly or by modulating the microenvironment of the tumor [25, 26]. CD4 ${}^{+}$ T cells divide into at least seven subtypes, including Th1, Th2, Th17, and Treg cells, which may have a pro-or anti-tumor response [27]. Here, Th1 and Treg cells were found to be highly enriched in patients with a high fraction of leukocyte, while Th17 cells were substantially more concentrated in tumors with a low fraction of leukocyte. The tumor with a higher LF will require an increased Th17/Treg cell ratio in order to improve its immune response [27]. As expected, tumors with high LF also had higher expression of immune checkpoint genes, including PDL1 (CD274), PDCD1, and CTLA4.

The genetic makeup of the tumor is a major influencer of the immune microenvironment [12, 28]. Correlation analysis has shown that intratumoral heterogeneity and copy number aberrations are significantly associated with the LF. Patients with high copy number aberration are less infiltrated by leukocytes. Overall mutations did not show any correlation with the LF in contrast. However, several individual mutations were associated with the LF. Importantly, clinically relevant mutations such as KRAS and KEAP1 have been more frequently mutated in patients with lower LFs. KRAS mutation has been associated with lower T cell infiltration [29]. These results suggest that individual mutation may play a role in defining the immune microenvironment. The LF was also significantly associated with overall changes in the copy number variation. Many focal CNV were also associated with LF.

In the recent past, DNA methylation has emerged as a major mechanism for regulating genes [30]. We hypothesized that changes in DNA methylation must also regulate the immune microenvironment. We identified 787 CpGs with a strong correlation with the LF. Majority of these CpGs were also differentially methylated in LUAD compared to normal. Such observations suggest that tumor DNA methylation shifts during carcinogenesis, and that the immune microenvironment is regulated according to the specific changes. It is noteworthy, in comparison to genes with a negative correlation, which regulated developmental pathways, that genes that correlated positively with leukocyte infiltration were enriched with transcription-related metabolism. We also observed vast differences in gene expression in patients with high and low LFs. Most of these genes have been differentially regulated in tumors in comparison with normal samples and regulated immune pathways, as have methylation changes. It is possible that changes in gene expression during tumor development regulate the infiltration of immune cells. Next, we used the LF related genes to develop and validate an independent and strong immune risk score for the prognostication of LUAD patients. The signature member genes showed an association with immune pathways. ARNTL2 is a core gene associated with circadian clock [31]. ARNTL2 has been associated with LUAD patients outcome and shown to regulate the pro-metastatic secretome [32]. ITRIP is a transmembrane protein which may be involved in immune modulation [33]. Similarly, TSPAN32 is a T cell associated gene which may function as tumor suppressor [34]. ARHGEF6 is a guanine nucleotide exchange factor for Rho protein which is required for medulloblastoma genesis [35]. In this analysis we also identified two novel prognostic genes RP11-7560H9.5 and C11orf21. The function of these genes is still unknown.

5. Conclusion

Taken together, we have identified genetic, epigenetic aberrations associated with immune cell infiltration in LUAD tumors. We showed that these changes in LUAD determine the genetic and epigenetic makeup during carcinogenesis and the immune microenvironment. Our data suggest that the genetic and epigenetic nature of the tumors is the significant determinant of immune surveillance in LUAD. We also developed and validated a strong and robust prognostic signature for survival prediction of LUAD patients.

Authors’ contribution

Conception: Sudhanshu Shukla

Interpretation or analysis of data: Sudhanshu Shukla and Seema Khadirnaikar

Preparation of the manuscript: Sudhanshu Shukla, Seema Khadirnaikar, Annesha Chatterjee

Revision for important intellectual content: Sudhanshu Shukla, Seema Khadirnaikar and Annesha Chatterjee

Supervision: Sudhanshu Shukla

All the authors have read and approved the manuscript.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-203071.

sj-pdf-1-cbm-10.3233_CBM-203071.pdf - Supplemental material

Supplemental material, sj-pdf-1-cbm-10.3233_CBM-203071.pdf

Footnotes

Acknowledgments

The Cancer genome atlas (TCGA) is acknowledged for the clinical and genomic datasets. Sudhanshu Shukla is supported by the grants from DBT, Govt of India (BT/PR27478/MED/30/1954/2018) and SERB, DST, Govt. Of India (ECR/2018/000528) and CSIR, Govt of India (27(0366)/20/EMR-II.

References

Bade

B.C.

and Dela Cruz

C.S.

, Lung cancer 2020: epidemiology, etiology, and prevention, Clin Chest Med 41 (2020), 1–24.

Siegel

R.L.

Miller

K.D.

and Jemal

, Cancer statistics, 2020, CA Cancer J Clin 70 (2020), 7–30.

Russo

McCusker

M.G.

Scilla

K.A.

Arensmeyer

K.E.

Mehra

et al., Immunotherapy in lung cancer: from a minor god to the olympus, Adv Exp Med Biol 1244 (2020), 69–92.

Mao

and Fang

, Immunotherapy for non-small cell lung cancers: biomarkers for predicting responses and strategies to overcome resistance, BMC Cancer 18 (2018), 1082.

Sharma

Hu-Lieskovan

Wargo

J.A.

and Ribas

, Primary, adaptive, and acquired resistance to cancer immunotherapy, Cell 168 (2017), 707–723.

Carbone

D.P.

Gandara

D.R.

Antonia

S.J.

Zielinski

and Paz-Ares

, Non-small-cell lung cancer: role of the immune system and potential for immunotherapy, J Thorac Oncol 10 (2015), 974–984.

Braun

D.A.

Hou

Bakouny

Ficial

Sant’ Angelo

et al., Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma, Nat Med 26 (2020), 909–918.

Borcoman

Kanjanapan

Champiat

Kato

Servois

et al., Novel patterns of response under immunotherapy, Ann Oncol 30 (2019), 385–396.

Pinto

Petriella

Lacalamita

Montrone

Catino

et al., KRAS-driven lung adenocarcinoma and B cell infiltration: novel insights for immunotherapy, Cancers 11 (2019).

10.

Chen

and Hambardzumyan

, Immune microenvironment in glioblastoma subtypes, Front Immunol 9 (2018), 1004.

11.

Barnes

T.A.

and Amir

, HYPE or HOPE: the prognostic value of infiltrating immune cells in cancer, Br J Cancer 117 (2017), 451–460.

12.

Thorsson

Gibbs

D.L.

Brown

S.D.

Wolf

Bortone

D.S.

et al., The immune landscape of cancer, Immunity 48 (2018), 812–830.e14.

13.

Newman

A.M.

Steen

C.B.

Liu

C.L.

Gentles

A.J.

Chaudhuri

A.A.

et al., Determining cell type abundance and expression from bulk tissues with digital cytometry, Nat Biotechnol 37 (2019), 773–782.

14.

Maksimovic

Phipson

and Oshlack

, A cross-package Bioconductor workflow for analysing methylation array data, F1000Res 5 (2016), 1281.

15.

Maros

M.E.

Capper

Jones

D.T.W.

Hovestadt

von Deimling

et al., Machine learning workflows to estimate class probabilities for precision cancer diagnostics on DNA methylation microarray data, Nat Protoc 15 (2020), 479–512.

16.

Capper

Jones

D.T.W.

Sill

Hovestadt

Schrimpf

et al., DNA methylation-based classification of central nervous system tumours, Nature 555 (2018), 469–474.

17.

Mermel

C.H.

Schumacher

S.E.

Hill

Meyerson

M.L.

Beroukhim

et al., GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers, Genome Biol 12 (2011), R41.

18.

Shukla

and Khadirnaikar

, RNA-sequencing analysis pipeline for prognostic marker identification in cancer, Methods Mol Biol 2174 (2021), 119–131.

19.

Whiteside

T.L.

, The tumor microenvironment and its role in promoting tumor growth, Oncogene 27 (2008), 5904–5912.

20.

Qin

et al., Novel immune checkpoint targets: moving beyond PD-1 and CTLA-4, Mol Cancer 18 (2019), 155.

21.

Bruni

Angell

H.K.

and Galon

, The immune contexture and Immunoscore in cancer prognosis and therapeutic efficacy, Nat Rev Cancer (2020).

22.

Wellenstein

M.D.

and de Visser

K.E.

, Cancer-cell-intrinsic mechanisms shaping the tumor immune landscape, Immunity 48 (2018), 399–416.

23.

Mitra

Lauss

Cabrita

Choi

Zhang

et al., Analysis of DNA methylation patterns in the tumor immune microenvironment of metastatic melanoma, Mol Oncol 14 (2020), 933–950.

24.

Pitt

J.M.

Marabelle

Eggermont

Soria

J.C.

Kroemer

et al., Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy, Ann Oncol 27 (2016), 1482–1492.

25.

Borst

Ahrends

Ba̧bała

Melief

C.J.M.

and Kastenmüller

, CD4+ T cell help in cancer immunology and immunotherapy, Nature Reviews Immunology 18 (2018), 635–647.

26.

Jiang

Wei

and Wang

, Regulatory T cells in tumor microenvironment: New mechanisms, potential therapeutic strategies and future prospects, in: Molecular Cancer, 2020.

27.

Knochelmann

H.M.

Dwyer

C.J.

Bailey

S.R.

Amaya

S.M.

Elston

D.M.

et al., When worlds collide: Th17 and Treg cells in cancer and autoimmunity, Cell Mol Immunol 15 (2018), 458–469.

28.

Wang

Zhang

Ou-Yang

et al., Multi-omics analysis of microenvironment characteristics and immune escape mechanisms of hepatocellular carcinoma, Front Oncol 9 (2019), 1019.

29.

Hanggi

and Ruffell

, Oncogenic KRAS drives immune suppression in colorectal cancer, Cancer Cell 35 (2019), 535–537.

30.

Greenberg

M.V.C.

and Bourc’his

, The diverse roles of DNA methylation in mammalian development and disease, Nat Rev Mol Cell Biol 20 (2019), 590–607.

31.

Tao

Zhang

Fan

and Mao

, Pan-cancer analysis reveals disrupted circadian clock associates with T cell exhaustion, Front Immunol 10 (2019), 2451.

32.

Brady

J.J.

Chuang

C.H.

Greenside

P.G.

Rogers

Z.N.

Murray

C.W.

et al., An arntl2-driven secretome enables lung adenocarcinoma metastatic self-sufficiency, Cancer Cell 29 (2016), 697–710.

33.

Cheng

L.C.

Baboo

Lindsay

Brusman

Martinez-Bartolome

et al., Identification of new transmembrane proteins concentrated at the nuclear envelope using organellar proteomics of mesenchymal cells, Nucleus 10 (2019), 126–143.

34.

Lombardo

S.D.

Mazzon

Basile

M.S.

Campo

Corsico

et al., Modulation of tetraspanin 32 (TSPAN32) expression in T cell-mediated immune responses and in multiple sclerosis, Int J Mol Sci 20 (2019).

35.

Hemmesi

Squadrito

M.L.

Mestdagh

Conti

Cominelli

et al., miR-135a inhibits cancer stem cell-driven medulloblastoma development by directly repressing arhgef6 expression, Stem Cells 33 (2015), 1377–1389.

Genetic and epigenetic landscape of leukocyte infiltration identifies an immune prognosticator in lung adenocarcinoma

Abstract

BACKGROUND:

OBJECTIVES:

METHOD:

RESULTS:

CONCLUSION:

Keywords

Abbreviations

1. Introduction

2. Material and methods

2.1 Patients sample LF and CIBERSORT analysis

2.2 Expression analysis

2.3 Mutation analysis

2.4 Methylation analysis

2.5 Copy number analysis

2.6 Survival analysis

2.7 Other statistical analysis

3. Results

3.1 Characterization of immune nature of LUAD with high and low LF

3.4 Development of LUAD prognostic model using LF associated genes

3.5 Validation of prognostic markers

4. Discussion

5. Conclusion

Authors’ contribution

Supplementary data

sj-pdf-1-cbm-10.3233_CBM-203071.pdf - Supplemental material

Footnotes

Acknowledgments

References