Integrated bulk and single-cell RNA sequencing identifies an aneuploidy-based gene signature to predict sensitivity of lung adenocarcinoma to traditional chemotherapy drugs and patients’ prognosis

Study source and data preprocessing

A pan-cancer genomic aneuploidy score was downloaded from a previous study (Taylor et al., 2018). After excluding the aneuploidy score of other cancer samples, The Cancer Genome Atlas Lung Adenocarcinoma (TCGA-LUAD) samples with aneuploidy score were retained for further study (Table S1). Subsequently, the RNA-Seq data of TCGA-LUAD were acquired from the The Cancer Genoma Atlas GDC application programming interface (API) website. For validation, GSE31210 (Okayama et al., 2012) and GSE30219 (Rousseaux et al., 2013) data were downloaded from the GEO database. For TCGA-LUAD samples, samples with clinical information and survival time >0 days were retained but those without survival status were removed. Finally, a total of 483 LUAD samples were compiled for further analysis (Table S2). For GEO data, samples without survival status were deleted, and mean expression value of the genes was adopted when multiple probes matched to one gene. All the data used in this article were derived from publicly available databases without access restriction and therefore did not require additional ethical approval. Data collection and analysis were performed strictly following the relevant rules and regulations.

Selection of differentially expressed genes (DEGs) related to aneuploidy and functional annotation

Differential gene analysis between tumor and adjacent tissue samples was conducted using the TCGA-LUAD database. False discovery rate (FDR) <0.05 and |log2FC| > 1 were the screening criteria. At the same time, the correlation between DEGs and aneuploidy was analyzed using Spearman method (R > 0.3 and p < 0.05). Finally, the clusterprofiler package (Yu et al., 2012) was employed to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on these aneuploidy-related DEGs.

Developing a riskscore model based on aneuploidy-related genes

The TCGA dataset and the GEO dataset served as a training and a validation set, respectively. Prognostic genes for LUAD were screened from all the aneuploidy-related DEGs using univariate Cox analysis, followed by performing Lasso regression to shrink the number of genes in the risk model. The risk score for each patient was calculated using the following Eq. (1):

(1) $$ARS=\sum \beta i\times Expi$$where Expi represented the expression of prognosis-related gene, β represented the Cox regression coefficient of the corresponding gene.

Samples with ARS score > 0 after zscore transformation were assigned into high-ARS group, while those with ARS score < 0 were assigned into low-ARS group. Survival curves were plotted using Kaplan-Meier for prognosis analysis. Additionally, the log-rank test was used to assess the significance of the differences.

Immune landscape analysis

To comprehensively assess the immune cell composition and activity in LUAD, we used two complementary computational approaches (ssGSEA and CIBERSORT). Based on the reference dataset from Charoentong et al.’s (2017) study ssGSEA was used to calculate the score of 28 immune cells in each sample. In addition, CIBERSORT is a more precise quantification method that estimates the proportion of 22 immune cells using a predefined set of immune cell-specific marker genes (Chen et al., 2018).

Immunotherapy and drug sensitivity evaluation

An immunotherapy cohort (IMvigor210) was included to evaluate the influence of ARS on predicting immunotherapy sensitivity. The IMvigor210 cohort of patients demonstrated different degrees of responses to PD-L1 receptor blockers, including progressive disease (PD), stable disease (SD), partial response (PR), and complete response (CR). In addition, pRRophic package (Geeleher, Cox & Huang, 2014) was adopted to estimate potential sensitive chemotherapy drugs between two ARS groups in the TCGA queue. A lower half maximal inhibitory concentration (IC50) value indicated a higher sensitivity to drugs.

Single-cell RNA sequencing analysis of LUAD samples

The GSE203360 dataset (He et al., 2022) including six LUAD samples was used for scRNA sequencing analysis. The R package Seurat (Stuart et al., 2019) (v4.0.4) was adopted for following analysis. Cells were primarily filtered based on the criteria that each gene expressed at least in three cells and each cell contained more than 200 expressed genes. Cells with mitochondrial genes >25% were excluded to delete aging cells with limited biological significance. Subsequently, dimensionality reduction clustering analysis was performed employing principal component analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE), FindNeighbor and FindCluster function. Based on cell markers from the CellMarker2.0 website, initial subpopulations were reclassified according to high-expressed marker genes. The DEGs in annotated subgroups were screened by FindAllMarkers function, followed by KEGG enrichment analysis.

Screening of differential regulons and cell communications in tumor cells divided by ARS

Cancer cells were extracted from scRNA-seq data. The ARS score of each cell was calculated using the Eq. (2). Thereafter, these cancer cells were stratified into high and low ARS cell groups under the same grouping criteria used in TCGA-LUAD queue. The human RcisTarget database was downloaded from website (https://resources.aertslab.org/cistarget/). R package “SCENIC” (Van de Sande et al., 2020) was employed for the development of TF regulatory networks. A potential TF-target regulatory relationship was established first and then motif analysis was used to eliminate genes indirectly regulated by TF. A regulon is a set of TFs and their directly regulated target genes. The connections among the regulons were analyzed using connection specific indicators (CSI). CellChat software (Jin et al., 2021), a free online R package (https://GitHub.com/sqjin/CellChat), was applied to infer intercellular communications and select predominant receptor–ligand pairs.

Cell culture and transient transfection

Human LUAD cell lines including H1299, A549 and human normal bronchial epithelial cell BEAS-2B were commercially obtained from BNCC (Beijing, China). DMEM F12 with 10% FBS (Gibco, Thermo Fisher, Massachusetts, USA) was used for culturing A549 and BEAS-2B cells in a humidified environment at 37 °C with 5% CO₂. The H1299, A549 cell lines were transfected with CKS1B siRNA (Dharmacon, Lafayette, CO, USA) utilizing Lipofectamine RNAi Max (Invitrogen, Waltham, MA, USA). The target sequence for CKS1B siRNA was CGCACAAACAAATTTACTATTCG (5′–3′).

QRT-PCR

Total RNA was separated using TRIzol (Thermo Fisher, Waltham, MA, USA) reagent. QRT-PCR was conducted on the extracted RNA from each sample (2 μg) using FastStart Universal SYBR Green Master on a LightCycler 480 PCR System (Roche, Switzerland, USA). The cDNA consisted of a total reaction volume of 20 μl (an appropriate water volume, 0.5 μl of forward and reverse primers, 2 μl of cDNA template, and 10 μl of PCR mixture) and served as a template. PCR reaction was reacted following the cycling conditions that began with initial DNA denaturation phase at 95 °C for 30 s, followed by 45 cycles at 94 °C for 15 s, at 56 °C for 30 s, and at 72 °C for 20 s. Each sample was run in triplicate. Threshold cycle (CT) data were obtained and standardized to the level of GAPDH applying the 2^−ΔΔCT method. The mRNA expression in both tumor and normal control tissues was compared. The sequences of primer pairs for the genes were listed in Table 1.

Transwell assay

Migration and invasion of H1299 and A549 cell lines were detected by Transwell assays. Briefly, chambers coated (for invasion) or uncoated with Matrigel (BD Biosciences, Franklin Lakes, NJ, USA) (for migration) were first inoculated with cells (5 × 10⁴). Serum-free medium and complete DMEM medium were added to the upper and lower layer, respectively. After incubation for 24 h, 4% paraformaldehyde was used for fixing migrating or invading cells, followed by cell dyeing using crystalline violet (0.1%).

Cell viability

Cell viability was detected by performing cell counting kit-8 assay (Beyotime, Beijing, China). Specifically, the cells were different treatments were cultured in 96-well plates at a density of 1 × 10³ cells/well and added with CCK-8 solution at indicated time points. After 2-h incubation at 37 °C, a microplate reader (Thermo Fisher, Waltham, MA, USA) was utilized to measure the O.D 450 values of each well.

Statistical analysis

All analysis were conducted utilizing R package software (version 4). Survival data were summarized using Kaplan–Meier approach, and statistical significance was detected by the log-rank test. The time receiver operating characteristic (ROC) curves were adopted to evaluate the prediction ability of the model. Student t-tests were used to determine statistical significance with a one-way ANOVA used for comparative testing. A p < 0.05 was thought to be statistically significant. SangerBox (http://sangerbox.com/) provided analytical assistance (Shen et al., 2022).

DEGs related to aneuploidy and functional annotation

Firstly, RNA-seq data of the TCGA-LUAD dataset were used to screen differential genes (FDR < 0.05 and |log2FC| > 1) between tumor and para-cancerous tissue samples and a total of 3,316 DEGs were obtained for subsequent analysis (Fig. 1A). Correlation analysis detected a significant correlation between 182 DEGs and AS, and these 182 DEGs were high-expressed in tumor tissue (Fig. 1B). Further gene functional enrichment analysis showed that these 182 DEGs were enriched mainly in DNA replication and cell cycle-related pathways (Figs. 1C–1F).

Establishment and assessment of the riskscore model based on the aneuploidy-related genes

Univariate Cox analysis identified 125 genes with significant prognostic impact (p < 0.01) (Table S3). Lasso cox regression in R software package glmnet (Simon et al., 2011) was used to reduce the gene number. Ten-fold cross validation was applied to develop the model. The model was the optimal when penalty parameter lambda was 0.038 (Figs. 2A and 2B), therefore, five genes at lambda = 0.038 were chosen as key prognostic genes. The coefficient of the corresponding gene was shown in Fig. 2C. Each patient’s riskscore was calculated following Eq. (2).

(2) $$\begin{array}{rl}ARS=& -0.004\times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{P}\mathrm{L}\mathrm{K}1+0.208\times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{A}\mathrm{N}\mathrm{L}\mathrm{N}+0.059\\ & \times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{H}\mathrm{M}\mathrm{M}\mathrm{R}+0.033\times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{C}\mathrm{K}\mathrm{S}1\mathrm{B}+0.073\times \mathrm{E}\mathrm{x}\mathrm{p}\mathrm{U}\mathrm{B}\mathrm{E}2\mathrm{S}\end{array}$$

Based on Eq. (2), samples in TCGA-LUAD queue were separated into high and low ARS groups. It could be observed that the OS time in low ARS group was longer than that in high-ARS group (p < 0.0001, Fig. 2D). We observed that in the TCGA cohort, the area under the ROC curve (AUC) was 0.68, 0.65, and 0.67 for 1-, 3-, and 5-year OS, respectively (Fig. 2D). The five key genes were all overexpressed in tumor tissues (Fig. 2D). The model was confirmed using GSE31210 and GSE30219 datasets (Figs. 2E and 2F). Specifically, patients in the GSE31210 cohort had AUC values of 0.63, 0.68, and 0.74 for 1-, 3-, and 5-year OS, respectively, while those in the GSE30219 cohort had AUC values of 0.77, 0.74, and 0.75 for 1-, 3-, and 5-year OS, respectively.

The ARS model was an independent indicator for LUAD prognosis

To verify the applicability of the five-gene ARS signature in clinical setting, COX regression was performed. In the TCGA cohort, staging and risk scores and patients’ OS were significantly correlated (p < 0.001) (Fig. S1A). Multivariate analysis showed that staging and risk score were also independent predictors (Fig. S1B). Subsequently, similar results were found in the GSE30219 cohort, and age, staging and risk score were also the independent predictors of LUAD prognosis (p < 0.001) (Figs. S1D, S1E). Importantly, the nomogram all showed the same predictive function (Figs. S1C, S1F). These results suggested that the risk model developed based on these two datasets was a favorable tool for predicting the prognostic outcomes of LUAD patients, and that both staging and risk scores were independent prognostic factors.

CKS1B was a characteristic gene in the ARS to promote the malignant phenotypes of LUAD

To validate the reliability of the ARS model in predicting the prognosis of LUAD and the correlation between the model genes and LUAD progression, we measured the expression of PLK1, ANLN, HMMR, CKS1B and UBE2S in LUAD cell lines H1299 and A549 and bronchial normal epithelial cells BEAS-2B using PCR. Consistent with the bioinformatics results, the level of PLK1, ANLN, HMMR, CKS1B and UBE2S was elevated in the H1299 and A549 cell lines, in particular CKS1B was significantly upregulated in tumor cell lines (Figs. 3A–3E). Next, changes in migration and invasive of H1299 and A549 after inhibiting the expression of CKS1B in the two cell lines were examined. We found that CKS1B expression promoted the development of LUAD, and that the migration and invasion abilities of the two cell lines were significantly suppressed after inhibition of CKS1B (Figs. 3F and 3G). In addition, silencing CKS1B expression affected the viability of H1299 and A549 cell lines relative to controls (Figs. 3H and 3I). The above results demonstrated the reliability of using ARS model to predict the prognostic survival of LUAD, and in particular, the model genes significantly contributed to the malignant phenotype of LUAD.

Analysis of immune cell infiltration in the two cohorts

As shown in Fig. 4A, most of immune cells including T follicular helper cell, central memory CD4 T cell, immature B cell, activated B cell, immature, activated and plasmacytoid dendritic cell, eosinophil cell, macrophage, type 17 T helper cell, mast cell, monocyte and natural killer cell had higher score in low ARS group (p < 0.05), while few immunocytes such as activated, effector memory CD4 T cell, memory B cell, Type 2 T helper cell, gamma delta T cell, natural killer T cell and neutrophil were enriched in high ARS group. CIBERSORT was used to perform a more detailed assessment on different states of the immune cells. It could be found that low ARS group had obviously a higher proportion of T cells CD4 memory resting, Tregs, dendritic cells resting, and mast cells resting (p < 0.01), whereas that of T cells CD4 memory activated and macrophages at M0 and M1 stage were noticeably reduced (p < 0.0001, Fig. 4B).

The responsiveness of ARS model to immunotherapy and drug sensibility

Immunotherapy is a treatment that fights against cancer with synergistic survival benefits (Curran et al., 2010). Herein, we tested the prognostic significance of ARS model in immunotherapy dataset IMvigor210 cohort. Kaplan-Meier curves displayed that patient with a higher ARS had shorter OS (p < 0.05, Fig. 5A). Among them, SD/PD patients exhibited a higher ARS than CR/PR patients (p < 0.05, Fig. 5B). In high ARS group, the proportion of SD/PD was higher than that in the low ARS group (Fig. 5C). Applying pRRophic package, we further predicted potential sensitive drugs for different patients and found that patients in high ARS group were more sensitive to cosplatin, BI-2536, pyrimethamine and GW843682X, whereas those in low ARS group were more responsive to erlotinib and roscovitine (Figs. 5D and 5E).

Single-cell landscape of LUAD samples

He et al. (2022) generated scRNA-seq count matrix. PCA and t-SNE dimension reduction were performed based on six specimens derived from invasive adenocarcinoma patients. A total of 17 subpopulations were obtained applying FindNeighbor and FindCluster (Dim = 30, Resolution = 0.5) and reclassified according to the expression of cell markers in the subpopulations (Fig. S2). Finally, seven cell types were determined and denoted as alveolar epithelial cell, B cells, cancer cells, mast cells, monocytes/macrophages, myeloid dendritic cells and T cells (Figs. 6A and 6B). By using the FindAllMarkers function under the screening criteria of logfc = 0.5 (multiple of differences) and Minpct = 0.25 (smallest expression ratio of differential genes), the DEGs in the seven cell types were retained and the top five genes were shown in a heatmap (Fig. 6C). Subsequently, to explore the function of these cell subpopulations, we selected the top five genes for enrichment analysis. Figure 6D showed the enrichment of marker genes in different pathways in different cell types. Among them, monocytes/macrophages were mainly involved in pathways such as phagosome and lysosome. B and T cells were also predominantly enriched in the pathway of protein processing in endoplasmic reticulum. These results may provide insight into the complex biology within LUAD and offer potential targets for further improve personalized treatment for LUAD patients.

TF-target connection and function analysis

Applying SCENIC algorithm and CSI, key regulons were identified and clustered. As seen in Fig. 7A, regulons with similar CSI were clustered together. Regulons with analogous high CSI may jointly regulate downstream genes and were involved in cell functions. Next, we calculated the regulon specificity score (RSS) for the two ARS cell groups. The horizontal and vertical axes represented the ranking of the regulons and the RSS, respectively. The top five TFs in the ranking were represented by red dots. A higher RSS indicated a closer relationship between the regulons and cell type specificity. As shown in Figs. 7B and 7C, regulons such as CREB3L1, JUN_extended, CEBPD, FOS_extended and JUNB_extended were enriched in high ARS cell group. Nevertheless, the regulons including BCLAF1_extended, TAF7_extended, RAD21_extended, BCLAF1 and UQCRB ranked the top 5 in low ARS cell group.

CellChat identified communication patterns between immune cells and LUAD cells

Cell communication analysis helps to analyze the interaction between cells, explore the tumor immune microenvironment and dig potential therapeutic targets for diseases (Jin et al., 2021). As mentioned above, high and low ARS groups had different infiltration status of immune cells. Herein, we investigated the communication between immune cells and LUAD cells using scRNA-seq analysis. The results demonstrated that the cancer cells in high ARS group could receive signals from immune cells through the SPP1 signaling pathway via ligand receptor pair ITGA4-ITGB1 (Fig. 8A). As shown in Fig. 8B, tumor cells communicated with other immune cells through the MK signaling pathway via ligand receptor pair MDK-NCl. Moreover, high ARS group had thicker lines, which indicated stronger interaction signals. Communications through GDF (GDF15-TGFBR2) and ECF (AREG-EGRF) signaling pathway between high and low ARS cell groups were also found (Figs. 8C and 8D). These results demonstrated that ARS helped to detect an enhanced tumor-immune system interplay in LUAD.