Tumor microenvironment related novel signature predict lung adenocarcinoma survival

A total of 576 samples, gene expression RNAseq and clinical information such as pathological stage, survival of LUAD were downloaded from the TCGA database and prepared for the analysis of differential expression. The data set was submitted by University of North Carolina TCGA genome characterization center based on the Illumina HiSeq 2000 RNA Sequencing platform (Oct 13, 2017). We downloaded the estimate, immune and stromal score of 517 LUAD samples from TCGA in ESTIMATE website. Stromal score captures the presence of stroma in tumor tissue, immune score represents the infiltration of immune cells in tumor tissue, and estimate score that infers tumor purity. The estimate score is equal to the immune score plus the stromal score. Fifty-nine samples were deleted for lack of immune and stromal score in ESTIMATE website. Gene expression profiles and LUAD clinical covariates (stage, age, sex) were obtained from GEO database with the accession number, consisting of 442 samples (). After using the Robust Multichip Averaging (RMA) method to normalize the sequence matrix data, the data was used to validate the differential expression gene and test the prognostic model. The data set was submitted by Mervi Heiskanen based on Affymetrix HG-U133A (May 01, 2015) (). A total of 133 non-small cell lung cancer samples were downloaded from GEO database with the accession number, were to compare the prognostic model with other gene expression signature. GSE68465 GSE14814 [Shedden et al., 2008] [Marti et al., 2003]

The raw data of TCGA was pre-processed by limma algorithm (). The adjusted-values (adj. p) < 0.05 and |Log2 (FC)| > 1 were set as the cut-offs to screen for differentially expressed genes (DEGs). [Ritchie et al., 2015] p

ClustVis web tool was used to create heatmaps and clustering (). [Metsalu & Vilo, 2015]

DAVID conducted functional enrichment analysis, including biological process, molecular function and cell component to analyze DEGs (), an essential foundation for visualization, annotation, and integrated discovery. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways was also performed by the database of DAVID. Whole genome was the background using DAVID. False discovery rate (FDR) < 0.05 was used as the cut-off. [Huang, Sherman & Lempicki, 2009]

Kaplan–Meier plots and univariate Cox regression were used to visualize the association between the genes expression and overall survival of patients to explore its prognostic value. Statistical significance was examined using the Log-rank test.< 0.05 were considered statistically significant. p

Search Tool for the Retrieval of Interacting Genes (STRING) online database was used to retrieve predicted PPIs. Only interactions with a combined score >0.4 of all associations obtained in STRING were selected to construct the PPI network using Cytoscape software (). The Molecular Complex Detection (MCODE) plugin in Cytoscape was then utilized to find clusters of the PPI network. Top 8 nodes ranked by degree were identified as hub genes. [Bader & Hogue, 2003]

From Molecular Signatures Database, Hallmark and C7 gene sets v6.2 collections were downloaded as the target sets. GSEA performed using the software gsea-3.0. 517 samples from TCGA were used for GSEA, and only gene sets with nominal (NOM)-value < 0.05 and FDR < 0.25 were considered as significant. p

CIBERSORT computational method was applied for estimating the tumor-infiltrating immune subsets (TICs) abundance profile in all tumor samples, which followed by quality fifiltering that only 517 tumor samples with< 0.05 were selected for the following analysis. p

To reduce the dimensionality and eliminate collinearity, factor analysis was performed on 106 prognostic genes. Then we used the Step-wise multivariate Cox regression analysis to obtain factors correlated with overall survival significantly (SurFactors). By combining the score coefficient of factors weighted by their regression coefficients and standard deviation, the risk index of each gene was calculated as follows: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$${ {\rm \alpha} _i} = \displaystyle{{\sum\limits_{j = 1}^m {{{\rm \beta} _j} \times fa{c_{ij}}} } \over {st{d_i}}}$$\end{document} α i = ∑ j = 1 m β j c i j × f a s t d i

wherewas the risk index of gene,the number of SurFactors,the regression coefficient of factorin multivariate Cox regression analysis,the factor score coefficient of-th factor over geneandthe standard deviation of gene. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${ {\rm \alpha} _i}$\end{document} α i \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$i$\end{document} i \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$m$\end{document} m \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${{\rm \beta} _j}$\end{document} β j \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$j$\end{document} j \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$fa{c_{ij}}$\end{document} f a c i j \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$j$\end{document} j \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$i$\end{document} i \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$st{d_i}$\end{document} s t d i \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$i$\end{document} i

By combining the expression values of prognostic genes weighted by their risk index, the following risk scores can be established for each patient: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$Risk\_score = \sum\limits_{i = 1}^n {{{\exp }_i}} * {\alpha _i}$$\end{document} R i s k _ s c o r e = ∗ ∑ i = 1 n exp i α i

wherewas the number of prognostic genes, andthe expression value of gene. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$n$\end{document} n \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}${\exp _i}$\end{document} exp i \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$i$\end{document} i

We calculated risk score (RS) for every sample from testing data set based on the prognostic model and divided the specimens into two different groups based on the median RS. Kaplan–Meier survival curves were drawn and compared the two subgroups via log-rank tests. We divided the samples into two group according to their overall survival, receiver operating characteristic (ROC) curves were drawn by IBM SPSS statistics 22, and area under the curve (AUC) was calculated. Moreover, the univariate and multivariate analyses of survival were conducted to identify the prognostic factors for LUAD patients from TCGA data set. A nomogram and Calibration plots were established based on the TCGA LUAD cohort. The nomogram and calibration plot analysis were conducted by using the R package “rms” and “rmda.” All tests were two-tailed, and< 0.5 was considered statistically significant. All statistical analyses were conducted using R software version 3.4.2. p

Transcriptional expression profiles and clinical information of 517 LUAD patients were collected from TCGA database. Among them, 277 (53.6%) patients were female, 240 (46.4%) were male. Pathological stage included 277 patients (53.6%) of stage I, 122 (23.6%) stage II, 84 (16.3%) stage III, 26 (5%) stage IV and 8 (1.5%) unknown. Based on ESTIMATE algorithm, the scores of stromal and immune were calculated, ranging from −1355.85 to 3286.67 and −1959.31 to 2098.77, respectively. The Estimate score was significantly associated with Pathological stage (,= 0.0436). The lowest Estimate score was at the most advanced pathological stage, the stage IV. Immune score showed remarkable prognostic potential, correlated with the pathological stage (,= 0.026), while stromal score showed no correlation with the pathological stage (,= 0.145). Fig. 1A Fig. 1B Fig. 1C p p p

To explore whether the potential correlation existing between survival benefits and immune/stromal score, 517 LUAD patients were divided into high and low score groups based on their scores. Kaplan–Meier survival curves () showed that immune score significantly correlated with overall survival. Compared to the cases in the low score group, patients with high immune score have longer median overall survival (1725 d vs 1235 d,= 0.0152 in log-rank test). Although there was no statistically significant correlation between the stromal score and overall survival, patients with higher stromal score had longer median overall survival (, 1830 d vs 1293 d,= 0.0599 in log-rank test). Fig. 1D Fig. 1E p p

To identify DEGs profiles according to immune/stromal score, we obtained the gene expression array of 517 LUAD patients obtained TCGA database. Based on the comparison of immune score of high group and low group, after analysis with the limma software package algorithm, there were 903 genes up-regulated and 56 genes down-regulated. We draw the heatmap with the top highly variate gene to (). Similarly, 1,007 up-regulated genes and 30 down regulated genes were obtained for high stromal score compared with low score (). It is obvious that up-regulated DEGs, no matter for comparison based on immune score nor stromal score, take the major part in total DEGs. Therefore, we decide to focus on up-regulated DEGs in further analysis. Fig. 2A Fig. 2B

To reveal the biological function of the DEGs, GO and KEGG pathway enrichment analysis were performed for the up-regulated DEGs. GO analysis showed that immune up-regulated DEGs were remarkably enriched in immune response by biological processes, receptor activity by molecular functions, and plasma membrane by cellular component, respectively (–). Similarly, stromal up-regulated DEGs mainly enriched in cell adhesion by biological processes, extracellular matrix binding by molecular functions, and extracellular space by cellular component, respectively (–). KEGG pathway were both mainly enriched in cytokine-cytokine receptor interaction and cell adhesion molecules (CAMs) pathway (;). Figs. 3A 3C Figs. S1A S1C Fig. 3D Fig. S1D

To make further efforts to elucidate if up-regulated DEGs could give benefits to the LUAD patient survival, Kaplan–Meier survival and univariate Cox proportional hazards regression analyses were performed on the up-regulated DEGs. The cases whose overall survival <30 days were excluded. The results showed that 446 immune related DEGs correlated with patient survival (< 0.05), 387 stromal related DEGs correlated with patient survival (< 0.05). Among all those genes, there were 291 duplicates of immune and stromal genes, a total of 542 genes (prognostic DEGs) associated with overall survival in patients with LUAD (). p p Table S1

Based on the STRING database, PPI network was obtained by using Cytoscape software to clarify the interaction between prognostic DEGs. The network was constructed by 13 modules, and we choose the top three significant modules ranked by mcode score for further study (). A total of 527 nodes and 4,661 edges were screened from the PPI network. According to the criteria above, eight nodes (PTPRC, ITGAM, LCP2, CTLA4, CD80, ITGAX, CD19, CCR5) were identified as hub genes. Consistently, these hub genes serve crucial roles in maintaining the top three modules. Fig. 4

To test whether these prognostic genes have prognostic significance in other LUAD cases. Microarray expression profile datasetfrom the GEO database was downloaded and the data were subjected to prognostic gene selection. There are 176 genes among 542 prognostic genes cannot find indataset because of platform differences. Of the remaining 366 genes, 106 genes were confirmed involving in LUAD patient survival. By intersection analysis the 106 prognostic genes with hub genes of PPI network, only PTPRC and CD19 were associated with overall survival in patients with LUAD. GSE68465 GSE68465

Based on the study above, PTPRC and CD19 expression level were correlated with overall survival of LUAD patient, and PTPRC and CD19 high-expression group with longer survival. To explore the correlation of PTPRC and CD19 expression with clinical characteristic in LUAD patient, we analysis the PTPRC and CD19 expression level with TNM stage. The result indicated that the PTPRC and CD19 expression were negative correlated with the TNM stage of LUAD patient, and with the TNM stage rising, the expression of PTPRC and CD19 declined (). Fig. 5

To further elaborate the role of PTPRC and CD19 in LUAD, 517 samples from TCGA were divined into high/low expression group based on median expression. We implement GSEA in the high and low expression group of PTPRC and CD19. The result display that, in the high expression level whether PTPRC or CD19, the genes were basically enriched in immune relative activities including allograft rejection, complement and inflammatory response. In the low expression level of PTPRC, the genes were basically enriched in glycolysis, and typical tumor pathway including MYC-targets-V1 and MYC-targets-V2. As to the low expression group of CD19, genes were gathered in metabolic pathway. In C7 collection, many genes were concentrated in the high expression group of no matter PTPRC nor CD19, and both low expression group enrich few gene (). Fig. 6

Further, we performed correlation analysis between PTPRC and CD19 and 22 kind of immune infiltration cell. The result reveal that 16 kinds of immune infiltration cell were related with the expression of PTPRC (), 15 kinds associated with CD19 (). Fig. 7 Fig. S2

Factor analysis was used to determine common axes (or dimensions) of patterns and structures, which were measured by a reduced set of 106 predicted genes. The first twelve factors explained about 81.808% of the variation in the 106 prognostic genes. After multivariate Cox regression analysis, five factors correlated with survival were obtained (< 0.05). Through the formula mentioned before, the final prognostic model combined with 106 candidate gene were constructed. Further, to test the predictive ability of the prognostic model, the microarray expression profile datasetwas downloaded and we calculate the risk score for each patient. Based on the threshold risk score, we divided the patients into high- and low-risk group. We choose the log-rank test to identify the differences of OS between subgroups. The results showed that the overall survival of the high-risk patients is shorter than those in the low-risk ones in both TCGA dataset and(< 0.0001,= 0.0002 respectively,and). Then, to estimate whether the prognostic model is predictive of relapse free survival (RFS). We divided the patients into high- and low-risk group compared with the median risk score, log-rank test indicate that the prognostic model we build could predict the RFS of LUAD. ROC curves were also applied to evaluate the sensitivity of the prognostic model, and the ROC curves showed that the AUC value of the prognostic model reached 0.649, 0.617 in TCGA dataset andrespectively. The result indicates a substantially effective performance of the prognostic model for overall survival prediction (and). GSE68465 GSE68465 GSE68465 p p p Figs. 8A 8B Figs. 8C 8D

We compared the prognostic model constructed by 106 candidate gene with 9-gene biomarker, as it also use estimate score on TCGA datasets to build the prognostic model ().included all 106 genes of our prognostic model and 9 genes needed to constructed the 9-gene biomarker. Fordata set, the 106-gene prognostic model achived a higher C-index (C-index, 0.56; 95% CI [0.48–0.64];< 0.01) compared with the 9-gene biomarker (C-index, 0.50; 95% CI [0.44–0.56];= 0.7) (). GSE14814 GSE14814 [Chen et al., 2020] Fig. S3 p p

Further, univariate and multivariate Cox regression analysis have been executed between some of clinical pathological parameters and risk score of the prognostic model (). The result showed that the risk assessment model was an independent prognostic factor for prognostic. As the constructed risk assessment model with great prognostic value, we intended to improve the prognostic accuracy by intergrating with LUAD clinicopathological factors. We designed a nomogram to predict the survival of LUAD patient by combing T stage, lymph nodes metastasis, recurrence and risk score (). The AUC of the model was achieved 0.716 for 3 year and 0.699 for 5 year (and).andshow the nomogram calibration plots for predicting the overall survival of 3 years and 5 years of patient. Table 1 Fig. 9A Figs. 9D 9E Figures 9C 9D