Integration and deconvolution methodology deciphering prognosis-related signatures in lung adenocarcinoma

doi:10.21203/rs.3.rs-3233906/v1

Download PDF

Research Article

Integration and deconvolution methodology deciphering prognosis-related signatures in lung adenocarcinoma

https://doi.org/10.21203/rs.3.rs-3233906/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 14 Sep, 2023

Read the published version in Journal of Cancer Research and Clinical Oncology →

You are reading this latest preprint version

Purpose

To develop a machine learning classifier using 16 prognosis-related genes (PRGs) to stratify lung adenocarcinoma (LUAD) patients according to their risk score. Validate the model's performance, interpretability and generalizability. Investigate the role of PRGs in LUAD prognosis through single-cell and bulk transcriptome analysis.

Methods

We constructed a machine learning classifier based on 16 PRGs to categorize LUAD patients by risk score. We combined the risk score with clinicopathological factors to predict overall survival. We assembled a single-cell atlas of LUAD comprising over 50,000 cells to identify cell types linked to prognosis and studied PRG expression in different cell types. We analyzed PRG involvement in tumor pathways and interactions using gene set variation analysis (GSVA).

Results

The machine learning classifier combining risk score and clinicopathological factors showed strong performance in predicting overall survival. The single-cell atlas revealed that epithelial cells were primarily associated with prognosis. PRGs were predominantly expressed in malignant epithelial cells and influenced epithelial cell growth and progression, especially in tumor states. PRGs were involved in tumor pathways like epithelial-mesenchymal transition, hypoxia and KRAS_UP. High PRG GSVA scores correlated with worse outcomes in LUAD patients.

Conclusions

The model provides a valuable tool for clinicians to personalize LUAD treatment based on risk stratification. The study elucidated the biological basis of PRG signatures in LUAD through integrated single-cell and bulk transcriptome analysis, contributing to a better understanding of LUAD prognosis and guiding targeted therapy development.

Lung adenocarcinoma

Prognosis-related signatures

Single cell transcriptome

Deconvolution

Lung adenocarcinoma (LUAD), a particularly virulent subtype of non-small cell lung cancer, remains a formidable global health adversary largely due to its dismal patient prognosis (Denisenko et al. 2018; Siegel et al. 2022). The advent of single-cell RNA sequencing has revolutionized our understanding of this disease, enabling a high-resolution exploration of tumor heterogeneity and prognostic cell subpopulations (Suvà and Tirosh 2019; Gohil et al. 2021). In this complex landscape, machine learning emerges as a powerful tool to integrate multidimensional omics datasets for prognostic prediction (Cruz and Wishart 2007; Komorowski et al. 2022). A diverse array of machine learning algorithms, including artificial neural networks, Bayesian networks, support vector machines, and decision trees, have been harnessed to tackle various cancer prediction tasks (Akay 2009; Lg and At 2013). The unveiling of a plethora of gene expression signatures tethered to LUAD prognosis has further amplified the potential of machine learning models (Langerød et al. 2007; Corrêa and Augsburger 2016; Stelloo et al. 2016). The integration of these signatures into machine learning models holds tantalizing promise for patient risk-stratification and the tailoring of treatments.

However, rigorous validation, model interpretability, and generalizability remain significant hurdles. The “black-box” nature of machine learning further complicates matters, obscuring the biological interpretation of model features. Confronting these limitations is paramount for the clinical translation of robust, personalized prediction models that harness the bounty of omics data to enhance lung adenocarcinoma patient outcomes. In this study, a machine learning pipeline was employed to thoroughly investigate prognosis-related signatures across several molecular profiles and develop a reliable prognostic model for LUAD. In addition, we elucidated the biological activities of these predictive risk gene signatures through the implementation of an integrated and deconvolution analysis of both single-cell and bulk transcriptome data of LUAD.

Publicly available data collection and processing

The single-cell RNA-Seq data used in this study consist of nine samples from patients with lung adenocarcinoma (LUAD) obtained from the GSE171145 dataset. Additionally, eight samples (GSM4058906, GSM4058911, GSM4058912, GSM4058913, GSM4058914 GSM4058917, GSM4058921, and GSM4058926) from individuals with healthy lungs were included, sourced from the GSE136831 dataset. The bulk RNA-seq data utilized in this study were obtained from the following Gene Expression Omnibus (GEO) datasets: GSE13120, GSE68465, GSE31210, and GSE72094. The GEO datasets can be accessed at http://www.ncbi.nlm.nih.gov/geo. The additional bulk RNA-Seq datasets for OncoSG (Chen et al. 2020) and the Cancer Genome Atlas (TCGA) were obtained from the cBioPortal database (http://www.cbioportal.org) and the UCSC Xena (http://xena.ucsc.edu), respectively. All datasets underwent batch effect correction. The mRNA expression data from the bulk lung adenocarcinoma cohorts were converted to z-scores to improve comparability across different databases. The “ComBat” function from the R package sva (version 3.46.0) was used to eliminate batch effects.

Univariate Cox regression and correlation analysis

The univariate Cox regression analysis was performed on both the training and testing sets using the R package glmnet (version 4.1-7), focusing on the intersection genes. The selection of consensus prognosis-related genes for further investigations was based on the criteria of statistical significance (p < 0.05) and hazard ratios (HRs) greater than 1 or less than 1. The correlation analysis was performed using the “cor” function in R software (version 4.2.3).

Developing a robust machine learning-based classifier

A combination of eight machine learning algorithms, including random survival forest (RSF), elastic network (Enet), stepwise Cox regression (StepCox), CoxBoost, partial least squares regression for Cox (plsRcox), and supervised principal components (SuperPC), was utilized to develop a classifier to generalize boosted regression modeling (GBM) and survival support vector machine (survival-SVM). The assessment of the classifiers’ performance was carried out utilizing the Concordance index (C-index). The assignment of GBM risk score for each sample was performed using the “gbm.perf” function from the R package gbm (version 2.1.8.1). The Kaplan-Meier (KM) curve analysis was performed using the R package survival (version 3.5-5) and survminer (version 0.4.9). The receiver operating characteristic (ROC) analysis was conducted using the R package timeROC (version 0.4) and survivalROC (version 1.0.3.1).

Multivariate Cox regression analysis

To assess the statistical independence of the risk score, a multivariate Cox regression analysis was performed using the R package glmnet (version 4.1-7). This analysis offers estimates of the contribution of the risk score while controlling for other relevant factors.

Nomogram model development and evaluation

The PRGs risk score and clinical pathological factors were enrolled into multivariate Cox regression analysis to confirm prognostic factors associated with overall survival (OS) of LUAD patients. To predict the survival probability at 1, 3, and 5 years for LUAD patients, a nomogram model was created using the R package regplot (version 1.1). Calibration curves were plotted to evaluate the agreement between the predicted survival probabilities and the actual observed survival rates, using the R package rms (version 6.7-0). In addition, decision curve analysis (DCA) was conducted using the R package caret (version 6.0–94) and ggDCA (version 1.2) in order to assess the clinical utility of a predictive model by quantifying the net benefit of using the model for decision-making compared to other strategies.

Single-cell RNA-seq data generation, quality control and across-sample integration

The samples were converted into a Seurat object using the "CreateSeuratObject" function from the Seurat package (version 4.3.0.1) in R software (version 4.2.3) (Satija et al. 2015). We removed cells with either lower than 500 or higher than 6000 expressed genes (Supplementary Figs. 1A). In addition, we discarded cells with mitochondria content higher than 20% (Supplementary Figs. 1B). Finally, a total of 55,426 cells were acquired for the subsequent analysis. The R package harmony (version 0.1.1) was used as the batch effect removal method (Korsunsky et al. 2019). We used the R package Seurat to normalize expression matrices by the functions of “NormalizeData” and “ScaleData”. Then “FindVariable” function of Seurat package was applied to select the top 2000 variable genes (Supplementary Figs. 1C).

Dimension reduction and unsupervised clustering

Principal components analysis (PCA) was performed on the integration-transformed expression matrix using the “RunPCA” function from R package Seurat, and the first 15 principal components (PCs) were used in the “FindNeighbors” function from R package seurat (Supplementary Figs. 1B). The resolution parameters of the “FindClusters” function were different for different cell types, with 0.5 for all cells. Uniform manifold approximation and projection (UMAP) was performed for visualization in two dimensions using the “RunUMAP” function with the same PCs and other default parameters (Supplementary Figs. 1D)

Cell type identification and differential gene expression analysis

The R package SingleR (version 2.0.0) was used to evaluate the cell type scores (Aran et al. 2019), and “FindAllmakers” function from R packsge Seurat was used to identify the top10 genes for each cluster. Besides, the R package COSG (version 0.9.0) was applied to evaluate the top 100 genes in each cluster (Dai et al. 2022). The differentially expressed genes (DEGs) of each cluster were generated by Seurat “FindAllMarkers” function.

Single cell deconvolution using Scissor and Gene sets enrichment analysis

The R package Scissor (version 2.0.0) was used to associate phenotypic data from bulk RNA-seq with the single-cell data (Sun et al. 2022). To determine the primary effector subgroups of PRGs, the enrichment scores of PRGs in different cell types was assessed using five enrichment score algorithms conducted by the R package irGSEA (version 2.1.4) and Seurat (version 4.3.0.1), namely “UCell”, “AUCell”, “Singscore”, “Jasmine” and “AddModuleScore”.

CopyKAT analysis

The R package copykat (version 1.1.0) was used to infer copy number variations (CNVs) based on the gene expression of all the epithelial cells in our single-cell RNA data (Gao et al. 2021). The parameters containing ngene.chr = 5, win.size = 25, KS.cut = 0.1, distance= ‘euclidean’ and id.type= ‘S’ were used in function “copycat”. Three samples with more than 100 epithelial cells in healthy lungs were used as reference, and CNV inference analysis on epithelial cells from LUAD was performed.

Cell-cell interplay and xCell infer cell composition

To investigate the potential communication between epithelial cells and other types of cells, the ligand-receptor interactions between different cells was analyzed using the R package CellChat (version 1.6.1) (Jin et al. 2021). The function “xCellAnalysis” in package xCell (version 1.1.0) was used to infer cell composition from bulk RNA-seq.

Pseudo-time trajectories analysis of epithelial cells

To further investigate the potential role of PRGs in epithelial cell development, we extracted all the epithelial cells and performed cell trajectory analysis using the R package monocle (version 2.26.0) (Qiu et al. 2017). The function “reductionDimension” in package monocle was used to reduce dimension, with the parameters of ‘max_components = 2’ and ‘method = ‘DDRTree’’. The function “orderCells” and “plot_cell_trajectory” was used to order cells and plot cell development trajectory.

Gene set enrichment analysis (GSEA) and Gene set variation analysis (GSVA)

Gene set variation analysis was conducted by function “gsva” in the R package GSVA (version 1.46.0) (Hänzelmann et al. 2013). The samples were divided into high PRGs GSVA score (HPGS) groups and low PRGs GSVA score (LPGS) groups based on the GSVA scores of the PRGs gene set. Differential expression genes analysis between the two groups was performed by using the R package limma (version 3.54.2). All differential expression genes were ranked in descending order according to log2FoldChange (log2FC). The gene set of 50 hallmark pathways investigated in this study was downloaded from GSEA website (https://www.gsea-msigdb.org/gsea/index.jsp). Then, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enriched pathways was identified by the R package clusterProfiler (version 4.6.2) (Yu et al. 2012).

Statistical analysis and graphical visualization

All data processing, statistical analysis, and plotting were conducted using the R software (version 4.2.3). The optimal cutoff value was determined using the R package cutoff (version 1.3). The plots displayed in this study was generated using R package ggplot2 (version 3.4.2), corrplot (vesion 0.92), survival (version 3.5-5), circlize (version 0.4.15 ), ComplexHeatmap (version 4.2.3), ggsci (3.0.0), ggpubr (0.6.0), timeROC (version 0.4), survivalROC (version 1.0.3.1), ggstats (version 0.3.0), regplot (version 1.1), rms (version 6.7-0), ggDCA (version 1.2), enrichplot (version 1.18.4), glmnet (version 4.1-7), clusterProfiler (version 4.6.2), cowplot (version 1.1.1), forestplot (version 3.1.1), ggplotify (version 0.1.1), ggrepel (version 0.9.3), ggunchull (version 0.1.0), gplots (version 3.1.3), jjAnno (version 0.0.3), monocle (version 2.26.0), Cellchat (version 1.6.1), patchwork (version 1.1.2), pROC (version 1.18.4), Scissor (version 2.0.0), and scRNAtoolVis (version 0.0.5).

Building a robust and generalizable LUAD classifier

To explore prognosis-related signatures in LUAD, we performed an integrating analysis of different RNA-Seq data sets from multiple cohorts (Fig. 1). Firstly, 17 consensus prognosis-related genes (PRGs) were identified using a univariate Cox regression analysis performed on intersection genes of six independent cohorts (GSE13210, GSE68465, GSE31210, GSE72094, TCGA, and OncoSG). Then, a subset of 16 PRGs with correlation coefficient < 0.75 (Supplementary Figs. 2A) was selected for further investigation after correlation analysis.

The study utilized the 16-PRG set as a basis for developing a total of 26 classifiers. These classifiers were constructed by employing a mixture of eight distinct machine learning techniques (Fig. 2A). Upon evaluating the performance of 26 classifiers, we utilized the Concordance index (C-index) derived from both training and test sets. Our findings revealed that the classifier generated by employing the Gradient Boosting Machine (GBM) algorithm, known as the generalized boosted-regression-modeling classifier, attained a C-index of 0.69 in the training set. In addition, the consistency indices in the test sets (GSE13120, GSE31210, GSE68465, GSE72094, and OncoSG) ranged from 0.61 to 0.69 (Fig. 2A). All of these results collectively demonstrate that the GBM classifier exhibits excellent stability and accuracy. Subsequently, the scoring system of the GBM classifier was employed to allocate risk scores, and the samples were segregated into high-risk and low-risk categories based on a cutoff of GBM risk scores.

Additionally, the results of ROC analysis revealed the following area under the curve (AUC) values for 1-, 3-, and 5-year predictions, respectively: 0.725, 0.714, and 0.723 in TCGA; 0.924, 0.719, and 0.693 in GSE13120; 0.70, 0.686, and 0.706 in GSE31210; 0.685, 0.653, and 0.589 in GSE68465; 0.613, 0.656, and 0.601 in GSE72094; and 0.625, 0.638, and 0.756 in OncoSG (Fig. 2B-G). The analysis of Kaplan-Meier (KM) curves revealed that LUAD patients in the greater GBM-risk group exhibited a considerably worse prognosis compared to those in the lower GBM-risk group (Fig. 2H-M). Collectively, our findings indicated that the GBM classifier, which consisted of 16 PRGs, had a high level of accuracy in both patient classification and prognosis prediction for LUAD.

Integrating GBM risk score and clinicopathological factors in a prognostic model

Next, we examined the autonomy of the GBM-risk score as a prognostic model. The independence of GBM risk score was validated through multivariate Cox regression analysis, with a statistically significant p-value of less than 0.001, log (HR) = 5.75 (95%CI, 4.39–7.11) (Fig. 3A). Considering the inconsistent clinical information categories across distinct cohorts, we have opted to exclude non-shared clinical information. Consequently, we selected the TCGA dataset as the training set, while the GSE13120, GSE31210, GSE72094, and OncoSG datasets as the test sets. Subsequently, a prognostic model was constructed utilizing the Cox proportional hazards (Coxph) approach, and shared clinical data was integrated from several cohorts, including variables such as gender, age, tumor grade, and risk score. The results indicated that the Coxph model demonstrated promising predictive accuracy in both the training set of TCGA and the independent external test sets of GSE13120, GSE31210, GSE72094, and OncoSG (Fig. 3B-F). Specifically, in the training set of TCGA, the 1-, 3-, and 5-year AUC values were 0.783, 0.754, and 0.753, separately (Fig. 3B). In the independent external test set of GSE13120, the 1-, 3-, and 5-year AUC values were 0.904, 0.792, and 0.752, separately. Similarly, in the independent external test sets of GSE31210, GSE72094, and OncoSG, the 1-, 3-, and 5-year AUC values were as follows: GSE31210 (0.843, 0.780, and 0.732), GSE72094 (0.659, 0.727, and 0.667), and OncoSG (0.666, 0.591, and 0.606) (Fig. 3C-F). Consequently, a visual representation of the Cox proportional hazards model, incorporating the GBM risk score, age, sex, and grade, was created utilizing the R package 'regplot' (Fig. 4A). The constructed nomogram showed a visual depiction of the prognostic model, enabling healthcare practitioners to approximate the likelihood of survival for particular patients by considering the prognostic markers that have been established.

In comparison to other clinical features, the GBM-risk score had a greater area under the curve (AUC) value as determined by ROC analysis (Fig. 4B-D). The calibration curves offered empirical evidences that support the reliable prediction ability of the nomogram in determining the survival rates at 1-year, 3-year, and 5-year intervals for LUAD patients (Fig. 4E-G). The DCA showed that the nomogram established in this study demonstrated enhanced precision in prognosticating survival outcomes in patients with LUAD as compared to conventional methods of either administering all available treatments or refraining from treatment altogether (Fig. 4H). These findings indicated that the nomogram has a potential to be a useful tool in the clinical decision-making process for LUAD patients. Collectively, the utilization of the GBM-based risk score, in conjunction with clinical characteristics, demonstrated encouraging efficacy in prognosticating the prognosis of LUAD patients.

Single-cell deconvolution revealing epithelial cells as a predominant cell subgroup associated with worse survival

Additionally, we conducted a reanalysis of the single-cell transcriptomes of a total of 17 samples, consisting of nine samples from patients with LUAD and eight samples from individuals with healthy lung sample. Following a comprehensive set of stringent quality control and filtering procedures, we conducted an integrative analysis of the transcriptomes from a collective of 55,426 cells. Eleven major cell types were identified based on characteristic canonical cell markers (Supplementary Figs. 2B, Supplementary Figs. 3A, C), including epithelial cells and various immune cell types such as T cells, B lymphocytes, plasma cells, monocytes, macrophages, mast cells, and dendritic cells (Supplementary Figs. 3A, C). Furthermore, stromal cell types such as fibroblasts and endothelial cells, as well as cycling cells, were also identified (Supplementary Figs. 3A, C). Subsequently, we excluded cycling cells and undefined cells from our examination and proceeded to develop a comprehensive single-cell atlas for LUAD (Fig. 5A-B). This atlas revealed notable heterogeneity in the proportions of stromal cells and immune cells across the samples (Fig. 5C), which aligns with findings from prior studies (Kim et al. 2020; Wu et al. 2021).

Leveraging the deconvolution algorithm, scissor, we projected the prognostic information extracted from LUAD BULK-seq data onto the single-cell atlas of LUAD. The findings revealed epithelial cells as a predominant cell subgroup associated with worse survival across all of six cohorts (Fig. 6A-F). This suggests that the cellular composition of epithelial cells plays crucial role in influencing the prognosis of LUAD.

Epithelial cells with a significant enrichment of PRGs, and malignant epithelial cells with higher levels of PRGs than non-cancerous counterparts

To gain a more comprehensive understanding of the operational mechanisms of consensus PRGs at a subcellular level, we conducted a thorough examination focusing on individual cellular samples obtained from clinical sources. Multiple methods were employed to evaluate the enrichment of PRGs within certain subpopulations of cells. It was revealed that PRGs had a prominent enrichment in epithelial cells (Fig. 7A), and there was a statistically significant increase in PRGs in LUAD when compared to healthy lung tissues (Fig. 7B). These findings provide additional evidence that epithelial cells are the primary subgroup of cells associated with prognosis in LUAD. In order to determine if malignant epithelial cells maintain a significant enrichment of PRGs, we conducted copy number karyotyping (CopyKAT) analysis. The epithelial cells were divided into two groups, high-grade malignancy group and normal group, followed by the scoring of PRGs in both cell categories utilizing five distinct algorithms. A notably elevated PRGs score was seen in malignant cells compared to their normal counterparts (Fig. 7C-D). In summary, these findings indicate that PRGs may exert a substantial influence on tumor formation at a subcellular level.

Correlation between PRGs and tumor microenvironment in LUAD

In order to examine the potential interactions among epithelial cells and other biological components, we employed CellChat as a tool for analyzing cell-cell communication. Our study revealed an intricate interplay among epithelial cells, macrophages, and T cells (Fig. 8A-B). Based on the findings obtained from the CopyKAT analysis, we proceeded to conduct a more detailed investigation into the communication between epithelial cells in both malignant and normal states. The results revealed a noteworthy contact between malignant cells and macrophages, as well as between malignant cells and T cells (Fig. 8C-D).

Moreover, in order to further investigate the correlation between PRGs and tumor microenvironment, we utilized the deconvolution method known as xCell to calculate the cellular makeup of each bulk RNA-seq data sample. This was done in combination with the Gene Set Variation Analysis (GSVA) score for each sample. By utilizing these metrics, we performed the correlation analysis between PRGs and the invasion of non-neoplastic cells. A positive correlation was identified between the levels of PRGs and M1 macrophages, as well as between PRGs and Th2 cells (Fig. 9A-C, Supplementary Figs. 4A-D). Conversely, a negative correlation was found between PRGs and CD4 Tcm cells, as well as between PRGs and M2 macrophages (Fig. 9A-C, Supplementary Figs. 4A-D). Furthermore, a notable disparity was observed in samples with varying tumor purities based on PRGs GSVA score. Specifically, samples with elevated tumor purity demonstrated markedly higher PRGs GSVA scores (Fig. 9D). In summary, these findings established a connection between PRGs and tumor microenvironment in LUAD, specifically in relation to the infiltration of non-neoplastic cells, including macrophages and T lymphocytes.

PRGs exhibiting higher enrichment in tumor-associated cell states

To investigate the potential function of PRGs in epithelial cell development, we extracted all of the epithelial cells and analyzed their cell trajectory using monocle2. The results demonstrated that the epithelial cells were separated into three states (Fig. 10A), with the majority of tumor cells being enriched in states 1 and 3 (Fig. 10B). The PRGs enrichment score in various states were evaluated by five aforementioned enrichment score algorithms, ‘AddModuleScore’, ‘UCell’, ‘AUCell’, ‘Singscore’, and ‘Jasmine’. The results showed the cells of states 1 and 3 exhibited significantly higher PRGs score enrichment than state 2 cells (Fig. 10C). These results indicate that PRGs may have a substantial effect on the growth and progression of epithelial cells, particularly in cell states associated with tumors.

High GSVA score of PRGs in LUAD experiences worse outcomes

To delineate the disparate biological functions of PRGs in LUAD, we further performed GSVA and GSEA analysis on previously aforementioned six cohorts. We found that the high PRGs GSVA score (HPGS) groups showed significantly enrichment of signal pathways including epithelial-mesenchymal transition (EMT), hypoxia, IL6/JAK/STAT3 and KRAS_UP (Fig. 11A-F, Supplementary Figs. 5A-F). Meanwhile, we also observed that the HPGS groups had a markedly reduced survival span (Fig. 11G-L). EMT plays a pivotal role in the advancement of cancer, as it enables the migration, invasion, and metastasis of cells (Craene and Berx 2013; Yuan et al. 2014; Pastushenko and Blanpain 2019). The increased prevalence of EMT markers in the HPGS groups may provide a plausible explanation for their poorer clinical outcomes of LUAD. Hypoxia is a widely recognized catalyst for the advancement of tumors and the development of resistance against therapeutic interventions (Eltzschig and Carmeliet 2011; Gilkes et al. 2014). The activation of the IL6/JAK/STAT3 pathway, involved in inflammation and immunity, has the potential to contribute to the remodeling of the tumor microenvironment, hence facilitating tumor growth and progression (Johnson et al. 2018; Choy et al. 2020). The KRAS_UP pathway, which is commonly shown to be mutated in LUAD, is closely linked to the manifestation of aggressive tumor characteristics and the development of resistance towards therapeutic interventions (Buscail et al. 2020; Drosten and Barbacid 2020). Individuals exhibiting elevated PRGs GSVA scores may potentially derive advantages from the implementation of more assertive treatment approaches or the utilization of specialized therapeutic interventions specifically designed to counteract the activation of these implicated pathways. Taken together, our findings indicated that LUAD patients with high PRGs GSVA score experienced worse outcomes, which was associated with the activation of cancer-related pathways.

This study demonstrates the robust prognostic utility of the 16-gene classifier and integrated prognostic model in LUAD. Through validation in multiple independent cohorts, these signatures exhibit generalizability without overfitting to a single dataset. The 16-gene classifier, comprised of prognosis-related genes (PRGs), provides additive and independent prognostic information beyond conventional clinical and pathological factors.

Multivariate and nomogram analyses revealed that the classifier's risk score significantly improved prognostic prediction even after accounting for variables like age, sex, and grade. This indicates that the biological processes represented by these PRGs, including glutathione metabolism (GCLM), tumor suppression (HTATIP2), RNA splicing (SNRPD1), apoptosis inhibition (BIRC5), and potential roles in metastasis (NUP37), confer substantial prognostic signals that complement clinical variables. The results of cell trajectory analysis showed that PRGs are essential to epithelial cell growth and development, especially in tumor-associated cells. Monocle2 (Qiu et al. 2017) has helped to clarify epithelial cell states, with most tumor cells in states 1 and 3. Since malignancy is characterized by increased cell proliferation, tumor cells may be more abundant in certain settings. State 2 cells with lower PRGs content may exhibit reduced proliferative potential, enhanced differentiation potential, and a decreased ability to produce tumors. Notably, Some PRGs have established mechanistic roles in LUAD biology governing proliferation, apoptosis, and metastasis. For instance, GCLM overexpression promotes tumor progression and decreases survival rate (Harris et al. 2015). The HTATIP2 protein acts as a tumor suppressor inhibiting proliferation and metastasis (Tong et al. 2009; Li et al. 2013). Elevated SNRPD1 correlates with an aggressive cancer phenotype (Schümperli and Pillai 2004; Quidville et al. 2013). NUP37 upregulation may enhance metastasis (Bilokapic and Schwartz 2012; Liu et al. 2012). BIRC5 inhibits apoptosis during carcinogenesis (Li et al. 2019; Frazzi 2021). Overall, these signatures provide clinically relevant prognostic information reflecting the biological processes driving LUAD pathogenesis. The generalizability and additive prognostic value of the 16-gene classifier underscores its potential clinical utility for risk stratification and precision medicine approaches in LUAD. Further research into the functional roles of these PRGs may yield biological insights and novel therapeutic targets.

Thorough examination of single-cell transcriptome data from both LUAD and healthy lung tissues, our findings highlight the notable involvement of epithelial cells and PRGs in the prognosis of LUAD. Firstly, the investigation of single-cell deconvolution of this study revealed that epithelial cells were the predominant cell subgroup that had a favorable correlation with worse survival in patients with LUAD. This discovery aligns with the established understanding of the involvement of epithelial cells in the development of lung cancer, given that these cells are the main site of initiation in LUAD (The Cancer Genome Atlas Research Network 2014; Kim et al. 2020; Wu et al. 2021). Our study also underscored the significance of the cellular composition of epithelial cells in exerting an impact on the prognosis of LUAD. This implies that implementing techniques focused on regulating the makeup of epithelial cells may have the potential to enhance the prognosis of individuals with LUAD. Additionally, our investigation unveiled a substantial augmentation of PRGs in epithelial cells, along with a statistically significant elevation of PRGs in LUAD in comparison to lung tissues in a healthy state. This observation implies that PRGs have a significant impact on the advancement and evolution of LUAD. The significant abundance of PRGs in malignant epithelial cells provided additional evidences in favor of this concept, indicating that the pursuit of targeting PRGs may be a feasible treatment approach for LUAD. However, additional research is required to clarify the precise functions of these PRGs in LUAD and to verify their potential as targets for therapeutic intervention.

The examination of intercellular communication and the relationship between PRGs and immune cell infiltration yields significant insights into the tumor microenvironment of LUAD. The observed contact between cancer cells and macrophages/T cells demonstrates a notable intercommunication that has the potential to facilitate the advancement of tumors. The noteworthy observation lies in the positive correlation between PRGs and M1 macrophages/Th2 cells, given that both cell types have the potential to display pro-tumoral activities. M1 macrophages, which have historically been regarded as having antitumoral properties, can also contribute to the advancement of cancer by promoting angiogenesis, modifying the extracellular matrix, and suppressing the adaptive immune response (Kessenbrock et al. 2010; Quail and Joyce 2013). Similarly, the pro-cancer role of Th2 cells is largely due to their ability to suppress Th1/CTL responses (Kidd 2003; Walker and McKenzie 2018). The elevated levels of PRGs observed in these cells indicate a potential involvement of PRGs in modulating an immunosuppressive milieu. On the other hand, the intriguing aspect is in the negative connection seen between PRGs and M2 macrophages/CD4 Tcm cells. M2 macrophages and CD4 Tcm cells are commonly associated with the provision of antitumor immunity (Mueller et al. 2013; Yunna et al. 2020). The observed reduction in cell infiltration, along with elevated levels of PRGs, suggests that PRGs might play a role in modifying the microenvironment to avoid immune surveillance. Additionally, our findings that samples exhibiting high tumor purity displayed elevated levels of PRGs suggest that the upregulation of PRGs is associated with the increased dominance of cancer cells inside the tumor. All of these provide evidences for PRGs involved in the formation of an immunosuppressive microenvironment, which facilitates uncontrolled tumor growth.

One potential constraint of this study pertains to the comparatively limited sample size, particularly in relation to the analysis conducted using single-cell RNA sequencing data. Although the bulk RNA sequencing data yielded valuable insights into the patterns of gene expression across a substantial cohort of patients, the single-cell analysis was limited to a small subset of patients. In order to validate the findings and achieve a more comprehensive understanding of the variability of cell types within lung tumors and their microenvironment, it would be necessary to conduct single-cell RNA sequencing on a larger cohort. Furthermore, the examination solely focused on the analysis of gene expression data, neglecting to explore other potential determinants of patient prognosis, such as gene mutations, copy number variations, and protein expression levels. The integration of multi-omics data from bigger cohorts is expected to enhance the robustness of prognostic models and uncover novel prognostic signals that extend beyond gene expression profiles. Notwithstanding these constraints, the utilization of both single-cell and bulk RNA sequencing in this investigation yields significant findings and underscores the potential of an integrated multi-omics strategy in analyzing larger groups of individuals. This approach has the capacity to enhance the development of more accurate prognostic tools and the identification of actionable therapeutic targets for personalized treatment of lung cancer.

Our investigation employed a combination of single-cell and bulk RNA sequencing methodologies to acquire significant knowledge regarding the microenvironment of lung tumors and to ascertain prospective targets for therapeutic intervention. A prognostic model was developed utilizing gene signatures to assign risk scores, thereby equipping clinicians with a valuable tool to facilitate informed treatment decisions for patients. The utilization of both single-cell and bulk analysis techniques has provided insights into the functional significance of specific genes that have elevated expression in tumor epithelial cells. These genes have been found to be crucial in facilitating intercellular communication, driving malignant transformation, and influencing the overall prognosis of the tumor. Our study found that individuals who had high-risk scores experienced worse outcomes, which was associated with an increased presence of cancer-related pathways. The results of our study emphasize potential targets for immunotherapy and alternative treatment strategies that have the potential to disrupt critical processes that promote tumor growth. Our findings can aid healthcare professionals in developing more precise and individualized treatment plans. In brief, the integration of single-cell and bulk RNA-seq data holds significant potential for elucidating the intricate dynamics of the tumor microenvironment and identifying actionable targets that can inform more precise treatment strategies tailored to individual patients' prognoses.

AUC Area under the curve

GO Gene ontology

GSEA Gene set enrichment analysis

GSVA Gene set variation analysis

KEGG Kyoto encyclopedia of genes and genomes

OS Overall survival

scRNA-seq Single-cell RNA sequencing

TME Tumor microenvironment

LUAD Lung adenocarcinoma

TCGA The Cancer Genome Atlas

GEO Gene Expression Omnibus

UMAP Uniform manifold approximation and projection

PRGs Prognosis-related genes

GBM Generalized boosted regression modeling

EMT Epithelial–mesenchymal transition

DCA Decision curve analysis

Acknowledgements

We would like to express our gratitude to the researchers who generously shared the original data of open single-cell and bulk RNA sequencing in public databases. Additionally, we extend our appreciation to the providers of open-source R packages, whose contributions were instrumental in our analysis.

Author contributions

The study was designed by YM, SJY, and LYQ. YM and SJY conducted the analysis. YM, SJY, TXL, ZXY, and TDC performed the validation in the independent cohort. YM, SJY, TXL, and ZXY contributed to the preparation of figures and tables. YM and SJY drafted the manuscript. LYQ, SJY, and YY provided revisions to the manuscript. All authors reviewed and approved the final version of the manuscript.

Data accessibility

The datasets utilized in this study are available in online repositories. The specific repository/repositories and corresponding accession number(s) can be found in the article or Supplementary Material.

Declaration of competing interest

Conflict of interest The authors have no conflict of interest

Ethic approval and consent to participate This study does not involve ethical approval.

Funding

This work was supported by Grants from the National Natural Science Foundation of China (81773159 and 81871203).

Consent for publication

Not applicable.

Appendix A. Supplementary data

Supplementary data to this article can be found online at

Akay MF (2009) Support vector machines combined with feature selection for breast cancer diagnosis. Expert Systems with Applications 36:3240–3247. https://doi.org/10.1016/j.eswa.2008.01.009
Aran D, Looney AP, Liu L, et al (2019) Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20:163–172. https://doi.org/10.1038/s41590-018-0276-y
Bilokapic S, Schwartz TU (2012) Molecular basis for Nup37 and ELY5/ELYS recruitment to the nuclear pore complex. Proceedings of the National Academy of Sciences 109:15241–15246. https://doi.org/10.1073/pnas.1205151109
Buscail L, Bournet B, Cordelier P (2020) Role of oncogenic KRAS in the diagnosis, prognosis and treatment of pancreatic cancer. Nat Rev Gastroenterol Hepatol 17:153–168. https://doi.org/10.1038/s41575-019-0245-4
Chen J, Yang H, Teo ASM, et al (2020) Genomic landscape of lung adenocarcinoma in East Asians. Nat Genet 52:177–186. https://doi.org/10.1038/s41588-019-0569-6
Choy EH, De Benedetti F, Takeuchi T, et al (2020) Translating IL-6 biology into effective treatments. Nat Rev Rheumatol 16:335–345. https://doi.org/10.1038/s41584-020-0419-z
Corrêa ZM, Augsburger JJ (2016) Independent Prognostic Significance of Gene Expression Profile Class and Largest Basal Diameter of Posterior Uveal Melanomas. Am J Ophthalmol 162:20-27.e1. https://doi.org/10.1016/j.ajo.2015.11.019
Craene BD, Berx G (2013) Regulatory networks defining EMT during cancer initiation and progression. Nat Rev Cancer 13:97–110. https://doi.org/10.1038/nrc3447
Cruz JA, Wishart DS (2007) Applications of machine learning in cancer prediction and prognosis. Cancer Inform 2:59–77
Dai M, Pei X, Wang X-J (2022) Accurate and fast cell marker gene identification with COSG. Brief Bioinform 23:bbab579. https://doi.org/10.1093/bib/bbab579
Denisenko TV, Budkevich IN, Zhivotovsky B (2018) Cell death-based treatment of lung adenocarcinoma. Cell Death Dis 9:117. https://doi.org/10.1038/s41419-017-0063-y
Drosten M, Barbacid M (2020) Targeting the MAPK Pathway in KRAS-Driven Tumors. Cancer Cell 37:543–550. https://doi.org/10.1016/j.ccell.2020.03.013
Eltzschig HK, Carmeliet P (2011) Hypoxia and inflammation. N Engl J Med 364:656–665. https://doi.org/10.1056/NEJMra0910283
Frazzi R (2021) BIRC3 and BIRC5: multi‐faceted inhibitors in cancer. Cell Biosci 11:8. https://doi.org/10.1186/s13578-020-00521-0
Gao R, Bai S, Henderson YC, et al (2021) Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol 39:599–608. https://doi.org/10.1038/s41587-020-00795-2
Gilkes DM, Semenza GL, Wirtz D (2014) Hypoxia and the extracellular matrix: drivers of tumour metastasis. Nat Rev Cancer 14:430–439. https://doi.org/10.1038/nrc3726
Gohil SH, Iorgulescu JB, Braun DA, et al (2021) Applying high-dimensional single-cell technologies to the analysis of cancer immunotherapy. Nat Rev Clin Oncol 18:244–256. https://doi.org/10.1038/s41571-020-00449-x
Hänzelmann S, Castelo R, Guinney J (2013) GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. https://doi.org/10.1186/1471-2105-14-7
Harris IS, Treloar AE, Inoue S, et al (2015) Glutathione and Thioredoxin Antioxidant Pathways Synergize to Drive Cancer Initiation and Progression. Cancer Cell 27:211–222. https://doi.org/10.1016/j.ccell.2014.11.019
Jin S, Guerrero-Juarez CF, Zhang L, et al (2021) Inference and analysis of cell-cell communication using CellChat. Nat Commun 12:1088. https://doi.org/10.1038/s41467-021-21246-9
Johnson DE, O’Keefe RA, Grandis JR (2018) Targeting the IL-6/JAK/STAT3 signalling axis in cancer. Nat Rev Clin Oncol 15:234–248. https://doi.org/10.1038/nrclinonc.2018.8
Kessenbrock K, Plaks V, Werb Z (2010) Matrix metalloproteinases: regulators of the tumor microenvironment. Cell 141:52–67. https://doi.org/10.1016/j.cell.2010.03.015
Kidd P (2003) Th1/Th2 balance: the hypothesis, its limitations, and implications for health and disease. Altern Med Rev 8:223–246
Kim N, Kim HK, Lee K, et al (2020) Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun 11:2285. https://doi.org/10.1038/s41467-020-16164-1
Komorowski M, Green A, Tatham KC, et al (2022) Sepsis biomarkers and diagnostic tools with a focus on machine learning. EBioMedicine 86:104394. https://doi.org/10.1016/j.ebiom.2022.104394
Korsunsky I, Millard N, Fan J, et al (2019) Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 16:1289–1296. https://doi.org/10.1038/s41592-019-0619-0
Langerød A, Zhao H, Borgan Ø, et al (2007) TP53 mutation status and gene expression profiles are powerful prognostic markers of breast cancer. Breast Cancer Res 9:R30. https://doi.org/10.1186/bcr1675
Lg A, At E (2013) Using Three Machine Learning Techniques for Predicting Breast Cancer Recurrence. J Health Med Inform 04: https://doi.org/10.4172/2157-7420.1000124
Li A, Zhang C, Gao S, et al (2013) TIP30 loss enhances cytoplasmic and nuclear EGFR signaling and promotes lung adenocarcinogenesis in mice. Oncogene 32:2273–2281, 2281e.1–12. https://doi.org/10.1038/onc.2012.253
Li F, Aljahdali I, Ling X (2019) Cancer therapeutics using survivin BIRC5 as a target: what can we do after over two decades of study? J Exp Clin Cancer Res 38:368. https://doi.org/10.1186/s13046-019-1362-1
Liu X, Mitchell JM, Wozniak RW, et al (2012) Structural evolution of the membrane-coating module of the nuclear pore complex. Proc Natl Acad Sci USA 109:16498–16503. https://doi.org/10.1073/pnas.1214557109
Mueller SN, Gebhardt T, Carbone FR, Heath WR (2013) Memory T cell subsets, migration patterns, and tissue residence. Annu Rev Immunol 31:137–161. https://doi.org/10.1146/annurev-immunol-032712-095954
Pastushenko I, Blanpain C (2019) EMT Transition States during Tumor Progression and Metastasis. Trends Cell Biol 29:212–226. https://doi.org/10.1016/j.tcb.2018.12.001
Qiu X, Mao Q, Tang Y, et al (2017) Reversed graph embedding resolves complex single-cell trajectories. Nat Methods 14:979–982. https://doi.org/10.1038/nmeth.4402
Quail DF, Joyce JA (2013) Microenvironmental regulation of tumor progression and metastasis. Nat Med 19:1423–1437. https://doi.org/10.1038/nm.3394
Quidville V, Alsafadi S, Goubar A, et al (2013) Targeting the Deregulated Spliceosome Core Machinery in Cancer Cells Triggers mTOR Blockade and Autophagy. Cancer Research 73:2247–2258. https://doi.org/10.1158/0008-5472.CAN-12-2501
Satija R, Farrell JA, Gennert DG, et al (2015) Spatial reconstruction of single-cell gene expression data
Schümperli D, Pillai RS (2004) The special Sm core structure of the U7 snRNP: far-reaching significance of a small nuclear ribonucleoprotein. Cell Mol Life Sci 61:2560–2570. https://doi.org/10.1007/s00018-004-4190-0
Siegel RL, Miller KD, Fuchs HE, Jemal A (2022) Cancer statistics, 2022. CA Cancer J Clin 72:7–33. https://doi.org/10.3322/caac.21708
Stelloo E, Nout RA, Osse EM, et al (2016) Improved Risk Assessment by Integrating Molecular and Clinicopathological Factors in Early-stage Endometrial Cancer-Combined Analysis of the PORTEC Cohorts. Clin Cancer Res 22:4215–4224. https://doi.org/10.1158/1078-0432.CCR-15-2878
Sun D, Guan X, Moran AE, et al (2022) Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol 40:527–538. https://doi.org/10.1038/s41587-021-01091-3
Suvà ML, Tirosh I (2019) Single-Cell RNA Sequencing in Cancer: Lessons Learned and Emerging Challenges. Mol Cell 75:7–12. https://doi.org/10.1016/j.molcel.2019.05.003
The Cancer Genome Atlas Research Network (2014) Comprehensive molecular profiling of lung adenocarcinoma. Nature 511:543–550. https://doi.org/10.1038/nature13385
Tong X, Li K, Luo Z, et al (2009) Decreased TIP30 expression promotes tumor metastasis in lung cancer. Am J Pathol 174:1931–1939. https://doi.org/10.2353/ajpath.2009.080846
Walker JA, McKenzie ANJ (2018) TH2 cell development and function. Nat Rev Immunol 18:121–133. https://doi.org/10.1038/nri.2017.118
Wu F, Fan J, He Y, et al (2021) Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun 12:2540. https://doi.org/10.1038/s41467-021-22801-0
Yu G, Wang L-G, Han Y, He Q-Y (2012) clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16:284–287. https://doi.org/10.1089/omi.2011.0118
Yuan X, Wu H, Han N, et al (2014) Notch signaling and EMT in non-small cell lung cancer: biological significance and therapeutic application. J Hematol Oncol 7:87. https://doi.org/10.1186/s13045-014-0087-z
Yunna C, Mengru H, Lei W, Weidong C (2020) Macrophage M1/M2 polarization. Eur J Pharmacol 877:173090. https://doi.org/10.1016/j.ejphar.2020.173090

No competing interests reported.

Supplementalfiles.docx

Download PDF

Journal Publication

published 14 Sep, 2023

Read the published version in Journal of Cancer Research and Clinical Oncology →

Editorial decision: Major revision
28 Aug, 2023
Reviews received at journal
28 Aug, 2023
Reviewers agreed at journal
07 Aug, 2023
Reviewers invited by journal
07 Aug, 2023
Editor assigned by journal
07 Aug, 2023
Submission checks completed at journal
04 Aug, 2023
First submitted to journal
04 Aug, 2023

You are reading this latest preprint version

Integration and deconvolution methodology deciphering prognosis-related signatures in lung adenocarcinoma

Status:

Journal Publication

Version 1

Abstract

Purpose

Methods

Results

Conclusions

Figures

Introduction

Materials and method

Publicly available data collection and processing

Univariate Cox regression and correlation analysis

Developing a robust machine learning-based classifier

Multivariate Cox regression analysis

Nomogram model development and evaluation

Single-cell RNA-seq data generation, quality control and across-sample integration

Dimension reduction and unsupervised clustering

Cell type identification and differential gene expression analysis

Single cell deconvolution using Scissor and Gene sets enrichment analysis

CopyKAT analysis

Cell-cell interplay and xCell infer cell composition

Pseudo-time trajectories analysis of epithelial cells

Gene set enrichment analysis (GSEA) and Gene set variation analysis (GSVA)

Statistical analysis and graphical visualization

Results

Building a robust and generalizable LUAD classifier

Integrating GBM risk score and clinicopathological factors in a prognostic model

Single-cell deconvolution revealing epithelial cells as a predominant cell subgroup associated with worse survival

Correlation between PRGs and tumor microenvironment in LUAD

PRGs exhibiting higher enrichment in tumor-associated cell states

High GSVA score of PRGs in LUAD experiences worse outcomes

Discussion

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1