Recognition of driver genes with potential prognostic implications in lung adenocarcinoma based on H3K79me2

Abstract

Lung adenocarcinoma is a malignancy with a low overall survival and a poor prognosis. Studies have shown that lung adenocarcinoma progression relates to locus-specific/global changes in histone modifications. To explore the relationship between histone modification and gene expression changes, we focused on 11 histone modifications and quantitatively analyzed their influences on gene expression. We found that, among the studied histone modifications, H3K79me2 displayed the greatest impact on gene expression regulation. Based on the Shannon entropy, 867 genes with differential H3K79me2 levels during tumorigenesis were identified. Enrichment analyses showed that these genes were involved in 16 common cancer pathways and 11 tumors and were target-regulated by trans-regulatory elements, such as Tp53 and WT1. Then, an open-source computational framework was presented (https://github.com/zlq-imu/Identification-of-potential-LUND-driver-genes). Twelve potential driver genes were extracted from the genes with differential H3K79me2 levels during tumorigenesis. The expression levels of these potential driver genes were significantly increased/decreased in tumor cells, as assayed by RT-qPCR. A risk score model comprising these driver genes was further constructed, and this model was strongly negatively associated with the overall survival of patients in different datasets. The proportional hazards assumption and outlier test indicated that this model could robustly distinguish patients with different survival rates. Immune analyses and responses to immunotherapeutic and chemotherapeutic agents showed that patients in the high and low-risk groups may have distinct tendencies for clinical selection. Finally, the regions with clear H3K79me2 signal changes on these driver genes were accurately identified. Our research may offer potential molecular biomarkers for lung adenocarcinoma treatment.

Keywords: Driver genes; Gene expression; H3K79me2; Lung adenocarcinoma.

Graphical abstract

Fig. 1

Correlation analysis of HM signals and gene expression. (A, C) Up-DEGs, (B, D) Down-DEGs. (A, B) HM signal changes (purple bars) and significant differences (pink bars) during LUAD tumorigenesis. The red lines in the 1st and 2nd circles correspond to the change ratio of HM signal = 1 and –log10(P value) = 2. (C, D) Spearman correlation coefficients between gene expression changes and HM signal changes in each of the 100 bins. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

H3K79me2 signals are strongly correlated with gene expression levels in LUAD tested via random forest algorithm. (A) ROC curves for the top five HMs. (B) The prediction results tested by the RF, GLM and SVM algorithms and the top five vital bins for each HM obtained by the RF algorithm. (C) AUC distributions for HM combinations with the same number of HMs. (D) and (E) HM occurrence times in the studied 5-HM-based models and 6-HM-based models, respectively. (F) AUC distributions for 11 HMs in the same bin across the 100 bins around the TSS. (G) The ranks of 11 HMs in the same bin. Top-ranked HMs indicate higher ‘IncMSE’ values and greater effects on gene expression.

Fig. 3

Enrichment analyses of the genes with differential H3K79me2 levels during oncogenesis. (A) The top 20 significant pathways in which the genes with differential H3K79me2 levels participated. (B) The interconnections of the 20 pathways described in A. (C) The top 20 diseases related to the genes with differential H3K79me2 levels and the corresponding number of genes in each disease. (D) The top 20 trans-regulatory elements and the number of genes regulated by these trans-regulatory elements.

Fig. 4

Construction, evaluation and validation of the risk score model. (A) RT–qPCR analyses for the 12 PLDGs. The relative expression levels were calculated via$2^{-\mathrm{\Delta }\mathrm{\Delta }{C}_t}$. (B and D) Kaplan–Meier survival curves for all-, high- and low-risk groups in the TCGA and GEO cohorts. (C and E) ROC curves for the risk score model in the TCGA and GEO databases. The distributions of the risk scores in LUAD patients with different (F) T-categories, (G) N-categories, (H) M-categories and (I) tumor stages. (J) Univariate and (K) multivariate Cox analyses for the risk score model and clinical characteristics.

Fig. 5

The risk score model predicts the sensitivities of drug therapies. (A) Landscape of immune cell infiltrations in the high- and low-risk groups. Red and blue represent cells with higher and lower infiltration levels, respectively. *P < 0.05; **P < 0.01; ***P < 0.001; ns not significant. (B) Association between risk scores and immune cell infiltration (all P < 0.001). Immunophenoscores comparison between the high- and low-risk groups for LUAD patients treated with (C) CTLA4_neg_PD1_neg, (D) CTLA4_neg_PD‐1_pos, (E) CTLA4_ pos_PD‐1_neg, and (F) CTLA4_ pos_PD‐1_pos. PD‐1_pos or CTLA4_pos indicates anti‐PD‐1 or anti‐CTLA4 therapy, respectively. (G) The ratios of normalized IC50 values of the 112 drugs between the high- and low-risk groups. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

The distributions of H3K79me2 signals on the 12 potential LUAD driver genes. The green and yellow histograms represent the average distributions of H3K79me2 in LUAD tumor cells and normal cells, respectively. The grey lines show the ratios of H3K79me2 signal changes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)