Development of a novel embryonic germline gene-related prognostic model of lung adenocarcinoma

Data source and preprocessing

Data were collected and preprocessed as previously described (Yue, Ma & Zhou, 2019). Briefly, as the training dataset, the expression profile and clinical information of LUAD cohort were downloaded from TCGA database using GDC tool and summarized into an expression and a clinical matrix, respectively. The ensemble ids in the expression matrix were converted into gene symbols according to the annotation file and the multiple expression data of a single gene were replaced by an average value. Genes that are not expressed in all samples are also removed from the expression matrix. For testing sets, we used lung adenocarcinoma as a keyword to search in the GEO database. After removing the datasets which are devoid of clinical data or expression data of genes in the model, four independent GEO datasets GSE30219, GSE50081, GSE37745 and GSE72094 were selected and downloaded via GEOquery package (Davis & Meltzer, 2007). For pan-cancer analyses, the expression profiles and clinical data of 33 TCGA cancer types were downloaded from TCGA database, and preprocessed according to previously mentioned procedures. Tumor cohorts with normal samples <5 were not included in the pan-cancer expression analysis. The gene expression values were all converted to transcripts per million transcripts (TPM) to facilitate subsequent analysis. The protein expression data of LUAD were obtained from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) database (https://cptac-data-portal.georgetown.edu/).

Identification of the differentially expressed embryonic germline genes (DEGGs)

The embryonic germline genes were downloaded from the website https://venn.lodder.dev (Bruggeman et al., 2020), and their expression data were extracted into a matrix. Then, the differentially expressed embryonic germline genes (DEGGs) between tumor and normal samples were identified using methods mentioned previously (Yue, Ma & Zhou, 2019). The multiple comparisons were controlled using the false discovery rate (FDR) (Benjamini & Hochberg, 1995). The screening criteria were —log2Fold Chang— > 1 and FDR < 0.05.

Functional enrichment and pathway analysis

The Gene Ontology(GO) enrichment and Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway analysis were performed according to previously described methods (Yu et al., 2012; Yue, Ma & Zhou, 2019). P. adjust (FDR) <0.05 was considered statistically significant.

Construction of EGRS in training dataset

The EGRS model was established using methods described previously (Yue, Ma & Zhou, 2019). Briefly, the 454 LUAD samples were included in the univariate Cox regression analysis to determine the prognostic value of DEGGs. P < 0.05 was considered statistically significant. The prognostic DEGGs were further screened out using LASSO regression via glmnet package in R software for avoiding model overfitting (Friedman, Hastie & Tibshirani, 2010). Then, the resulting key DEGGs were included in multivariate Cox analysis, and risk score formula was constructed as follow: ${\displaystyle \text{Risk score}=\sum \beta \ast \text{gene}.}$ The β and gene represent the coefficient of multivariate Cox regression and the expression value of corresponding DEGGs in EGRS, respectively. The * represents the multiplication operator.

Validation of the EGRS in the training and testing sets

To further validate the expression pattern, the protein expression profiles of DEGGs in EGRS between normal and tumor samples were analyzed using proteomic data of LUAD cohort obtained from CPTAC database. Then, the predicted effect of the EGRS was evaluated in the training and testing sets, respectively. Samples in the training set were divided into high- and low-risk group according to the median of risk score values. The risk score distribution, survival status and the expression level of 11 DEGGs in EGRS between the high- and low-risk samples were evaluated. In addition, we divided the patients in training set into two groups according to commonly used clinicopathological factors and analyze the differentiation of the expression levels of 11 DEGGs between the two groups by student’s t-test. P < 0.05 was considered statistically significant. Kaplan–Meier (KM) survival analysis and receiver operating characteristic (ROC) curve were used to evaluate the predicting power of the EGRS, and the prognostic performance of other clinicopathological factors was also compared with that of EGRS model. The KM plot and ROC curve were used to assess the generalization power of EGRS model in four independent GEO datasets. Finally, a pan-cancer expression and prognostic analyses of 11 DEGGs were performed to extend our findings to other types of cancers.

Construction and evaluation of a nomogram

For predicting the OS of individual LUAD patient, Nomogram was generated based on the results of the multivariate Cox analysis to predict 1-, 3-, and 5-year OS. The performance of the nomogram was evaluated using calibration plots.

Identification the DEGGs in LUAD

The transcriptome data of LUAD cohort which contains 497 tumor and 54 normal samples were downloaded from TCGA database, and the corresponding expression profiles of embryonic germline genes were extracted according to the data published previously (Bruggeman et al., 2020). The DEGGs between normal and tumor samples were identified by limma package in R software (Ritchie et al., 2015). Totally, 269 DEGGs were screened out, among which, the upregulated DEGGs were 249, and the down-regulated DEGGs were 20 (Fig. 1A, Table S1). The representative DEGGs are shown by heat map (Fig. 1B) and box plot (Fig. 1C), respectively.

Functional enrichment and KEGG analysis of DEGGs

The GO enrichment analysis of the DEGGs revealed that the DEGGs were mainly enriched in biological process of cell proliferation and DNA damage repair, such as nuclear division (GO:0000280), DNA replication (GO:0006260), cell cycle checkpoint (GO:0000075), double-strand break repair (GO:0006302) and interstrand cross-link repair (GO:0036297) (Figs. 2A, 2B). The KEGG analysis showed that the enriched pathways were cell cycle (hsa04110), homologous recombination (hsa03440), and Fanconi anemia pathway (hsa03460) (Figs. 2C, 2D). The detailed GO terms of the three categories, biological process (BP), cellular component (CC), molecular function (MF) and KEGG results were presented in Tables S2 and S3, respectively.

The construction of embryonic germline gene related prognostic model

For EGRS construction, the univariate Cox regression analysis was conducted and the resulting 103 DEGGs with prognostic value (p < 0.05, Fig. 3A) were further screened by LASSO regression. Then, the sixteen key DEGGs were selected to performed the multivariate Cox regression analysis (Figs. 3B, 3C). Finally, a prognostic signature containing eleven DEGGs (RAD54L, ZNF322, CENPI, IGF2BP1, IGF2BP3, RAD51AP1, E2F7, HMMR, DNAJC5B, ADAMTS12, NANOS1) was constructed and the risk score formula was presented as follow: ${\displaystyle \text{Risk score}=\left(-0.4294\right)\ast \text{RAD54L}-0.2411\ast \text{ZNF322}-0.4645\ast \text{CENPI}+0.2243\ast \text{IGF2BP1}+0.2095\ast \text{IGF2BP3}-0.5093\ast \text{RAD51AP1}+0.7979\ast \text{E2F7}+0.7968\ast \text{HMMR}-0.3272\ast \text{DNAJC5B}+0.1997\ast ADAMTS12-0.5248\ast \text{NANOS1}.}$

The validation of the EGRS for survival prediction

Due to the potential inconsistencies in gene expression between the transcriptome and the proteome, we first verified the protein expression of the 11 DEGGs in EGRS in normal and tumor samples through the CPTAC database. Given the limitations of available data, only 4 genes of the 11 DEGGs were analyzed. Consistent with the transcriptome data, the protein expression of 4 genes (IGF2BP1, IGF2BP3, HMMR and ADAMTS12) were up-regulated in LUAD tumor samples. The log2FoldChanges of IGF2BP1, IGF2BP3, HMMR and ADAMTS12 were 0.49, 1.14, 0.50 and 0.87 and the p-values of the 4 genes were 2.7E−02, 9.22E−08, 1.32E−07 and 7.27E−18, respectively. In addition, the pan-cancer expression analysis also revealed that 11 DEGGs in EGRS showed differential expression in most types of tumors, but diverse expression trends between tumor and paired adjacent normal tissues of some types of cancers were observed, suggesting that the DEGGs may play different roles in the development of different types of tumors (Fig. S1). Then, Kaplan–Meier (KM) survival analysis indicated that the overall survival between the high- and low-risk groups in training set was significantly different, and the low-risk patients had significantly better outcomes than that of high-risk patients (P = 4.083e−09) (Fig. 4A). Besides, our model had the best predictive power, compared with other clinicopathological factors such as age, gender, TNM and stage. Its area under curve (AUC) of ROC was 0.790 (Fig. 4B). The DEGGs in EGRS indicated differential expression between high- and low-risk samples (Fig. 4C), and with the increase of risk score, the proportion of death cases was also increasing (Figs. 4D, 4E). Between different pathological groups, the expression levels of genes in EGRS also showed significant difference (Fig. 5). Contrary to CENPI, the expressions of DNAJC5B in elderly patients (>65 years old) were higher than those in younger patients (≤65 years old) (Figs. 5A, 5B). The HMMR, RAD51AP1 were highly expressed in men than in women (Figs. 5D, 5E). ZNF322 was related to T stage, and its expressions in T3-4 were significantly lower than those in T1-2 (Fig. 5G). The CENPI, DNAJC5B, HMMR, RAD51AP1 and ZNF322 were correlated to N stage. Except for ZNF322, other genes were highly expressed in N1-2 stage (Figs. 5I–5L). The expressions of ADAMTS12 in M1 stage were lower than those in M0 stage (Fig. 5N). The expressions of DNAJC5B, HMMR, IGF2BP1, IGF2BP3, RAD51AP1, ZNF322 were correlated to the pathological stage. Except for ZNF322 and DNAJC5B, the expressions of other genes in stage III-IV were significantly higher than those in stage I-II (Figs. 5P–5U). More importantly, the risk scores in the advanced LUAD were significantly higher than those in the early stage (Fig. 5H, 5M and 5V), although its distribution had no statistical difference in some pathological statuses such as age, gender and tumor metastasis (Figs. 5C, 5F and 5O).

The univariate and multivariate Cox analyses suggested the risk score was an independent prognostic factor (Figs. 6A, 6B). The pan-cancer prognostic analysis also indicated the 11 DEGGs had prognostic value in some cancer types. Likewise, different prognostic effects based on cancer type were observed in some DEGGs (the detail results were packaged into a Supplementary file). In addition, the predictive power of the EGRS was further validated in four GEO testing datasets. Consistent with the training set, our model also stratified the samples in the testing sets into low-risk and high-risk groups, and the low-risk patients had better outcomes (Figs. 6C–6F). The AUC results also showed that our model had stable predictive ability. The AUC in four independent GEO datasets were 0.726 (GSE30219), 0.764 (GSE50081), 0.657 (GSE37745) and 0.668 (GSE72094) (Figs. 6G–6J).

Construction and evaluation of a nomogram

In order to more straightforwardly predict the prognosis of LUAD patients, a nomogram was constructed (Fig. 7A). According to the expression level of 11 genes in an individual patient, a numeric value was calculated to directly predict the 1-, 3-year and 5-year prognostic survival rate. The calibration plots showed that the predictive values of OS probabilities in 1-, 3-, and 5-year fit well with the observation values, suggesting that the nomogram is suitable for predicting the prognostic survival rate of LUAD patients (Figs. 7B, 7C, 7D).