Construction and validation of a prognostic signature based on anoikis-related lncRNAs in lung adenocarcinoma

Acquisition and processing of data

Clinical data and RNA sequencing (RNA-seq) results from 539 patients with LUAD and 59 normal controls were retrieved from The Cancer Genome Atlas (TCGA) database (last assessed on September 2022) in the FPKM format (https://portal.gdc.cancer.gov/). Overall, 37 patients were excluded from the analysis due to the absence of survival time or survival status data. By setting a relevance score of >2, 65 anoikis-related genes (ARGs) were retrieved from GeneCards [13]. We employed limma in R (|R| > 0.4 and p-value < 0.001) to perform Pearson correlation analysis on the coexpression of lncRNA and ARGs to discover potential ARLRs. Since our research relied on publicly available data, there were no ethical concerns.

Generation and assessment of the ARLR signature

Two distinct groups were formed from the LUAD patients at random, with a ratio of 1:1. Table 1 displays the patients’ clinical characteristics in both the training and the testing groups. To begin, we screened the training cohort for ARLRs using differential expression and Pearson correlation analyses. Next, genes were selected using univariate Cox and Lasso-Cox regression analyses, based on the ARLRs, to reduce the probability of overfitting. Thereafter, a prognostic signature was developed utilizing multivariate Cox regression analysis, which incorporated ARLRs with p < 0.05. The following calculation determined the anoikis-related lncRNA signature (ARLSig) risk score: risk factor = (ARLRs 1 expression × coefficient) + (ARLRs 2 expression × coefficient) … + (ARLRs n expression × coefficient).

We classified the LUAD patients in the training, test, and entire cohorts as high- or low-risk as determined by the median risk score value. The pheatmap package in R was utilized to generate heatmaps, risk curves, and survival statuses for various patient sets depending on the risk score. We utilized the timeROC package to conduct a time-dependent receiver operating characteristic (ROC) curve analysis to ascertain the ARLSig’s accuracy and sensitivity. The survival package was employed to execute K-M analysis to determine overall survival (OS), progression-free survival (PFS), and disease-specific survival (DSS) in patients with varying risk scores. Univariate and multivariate analyses were conducted to establish the signature’s independent prognostic significance. Clinicopathological characteristics were also evaluated for their predictive significance in each group via K-M analysis. To determine how the patients with varied scores were distributed, principal component analysis (PCA) was carried out.

The independent parameters were subsequently utilized to develop a nomogram. To gauge the ARLSig score’s reliability, we compared the expected and actual values in the TCGA-LUAD cohort utilizing calibration and time-dependent ROC curves.

Biological function of the risk scores

The high- and low-risk groups were compared via differential gene expression (DEGs) analysis with criteria of |log2FC| > 1 and p < 0.05, and volcano plots were applied to display the results. The clusterProfiler package was adopted for gene ontology (GO) functional analysis. The GSEA software v4.2.1 was also utilized for gene set enrichment analysis (GSEA), and functional enrichment was considered significant when the FDR < 0.05.

Analysis of the risk score’s relationship with TMB

To better understand how tumor mutation burden (TMB) correlated with risk score, we developed a waterfall plot and subsequently assessed the varying TMB in the groups at high and low risk. Subsequently, the survival package was employed to generate survival curves.

Infiltration of immune cells and anticancer treatment

A heatmap was used to display the results of the immune-related activities identified by the limma program in LUAD patients. The statuses of immune infiltration were examined utilizing the following tools: XCELL, CIBERSORT, MCP-counter, QUANTISEQ, EPIC, TIMER, and CIBERSORT-ABS. The results are displayed as a bubble diagram and are based on the infiltration estimate profile in the TCGA database. Additionally, the estimate package was utilized to ascertain the link between the TME score and the risk score.

An assessment of the response to immunotherapy was conducted by retrieving each LUAD patient’s immunophenoscore (IPS) from the TCIA database (https://tcia.at/home), which was determined by the expression of genes associated with immunity and denoting four distinct categories of immune (selected immunomodulators, effector cells, MHC molecules, and immunosuppressive cells). We employed the pRRophetic package to determine the half-maximal inhibitory concentration (IC50) to evaluate anticancer drugs for LUAD.

Cell culture and quantitative real-time PCR (qPCR)

The human LUAD cell line and one human embryonic lung cell line, MRC-5, were both stored in our laboratory. The experimental parameters for cell culture included the following: complete RPMI-1640 medium that contained 10% fetal bovine serum (FBS) and 1% antibiotics (100 ng/mL streptomycin and 100 U/mL penicillin) in a humidified incubator at 37°C with 5% CO2. The cells were treated with TRI reagent to isolate total RNA. PrimeScript RT was employed to reverse transcribe the total RNA to acquire cDNA, followed by qPCR (Takara Bio Company). The primers used for qPCR included the following: AC026355.2, forward primer: 5′-CTGGATGCTTCCTGCCCTTGAAC-3′, reverse primer: 5′-CCAACAGCCCCTGCCAAACC-3′.

Cell transfection

Lipofectamine 2000 was used for cell transfection in compliance with established protocols. After 48 to 72 hours, the transfected cells were harvested to use for additional experiments. In this study, a negative control siRNA (si-NC) and si-AC026355.2 were chemically synthesized by GenePharma (Shanghai, China). The sequence of si-AC026355.2 was 5′-GUGACAGGCAACACCUAUATT-3′.

Transwell experiment

The transfected A549 cells (6 × 10⁵ cells/mL) were placed in the upper chamber in a serum-free culture medium. Thereafter, the bottom chamber of the transwell plate was gradually filled with 500 μL of a culture medium that contained 10% FBS. Following 24 hours, the infiltrating cells were fixed for a half-hour at room temperature utilizing 4% polyoxymethylene, after which 0.5% crystal violet was introduced to stain them for 10 minutes. The cell count was performed in four fields of view that were selected at random utilizing a microscope.

Wound healing assay

Approximately 4 × 10⁵ cells per well were inoculated in 6-well plates. The next day, a sterile pipette tip was employed to create a scratch at the bottom of each well. At both the 0 and 48-hour post-scratch intervals, the area was microscopically imaged. The following equation was employed to compute the scratch healing rate: (initial scratch width - observed scratch width at the designated time point)/initial scratch width × 100%.

Statistical analysis

R software was employed for all statistical analyses undertaken in this study. We evaluated the model’s performance via ROC analyses. The survival rates of the high- and low-risk groups were compared via K-M analysis. The tests were all two-tailed, and p < 0.05 was established as the significance criterion.

Data availability statement

The data used for our analysis in this study are openly available in a public database (https://portal.gdc.cancer.gov/).

ARLRs identification of and creation of the prognostic signature

Figure 1 depicts a comprehensive flow overview of our study. We enrolled 59 healthy controls and 502 LUAD patients who had accessible RNA seq data, and an aggregate of 3,202 differentially expressed lncRNAs was discovered. As illustrated in Figure 2, coexpression relationships between 65 ARGs from GeneCards and differentially expressed lncRNAs are shown in the Sankey diagram. The next step was to randomly assign database patients to either a training or testing group. Table 1 displays the essential clinical characteristics of the LUAD patients at baseline. As determined by univariate and LASSO regression analyses, the prognostic signature was composed of 9 ARLRs (Figure 3). The ARLSig formula was expressed as indicated: risk score = (−0.345342588082988 × AC090912.1) + (0.251214414934825 × LINC00707) + (−0.18011074824835 × AC026355.2) + (−0.493695647381226 × FOCAD-AS1) + (0.172008983013553 × LINC00460) + (0.366834641287433 × LINC01117) + (0.4332718020797 × AC068228.1) + (−1.04234276424881 × AP000346.1) + (0.428355116775672 × LINC01537).

ARLSig validation

Each patient was assigned a risk score depending on a previously indicated formula. Subsequently, they were categorized into high- and low-risk groups using the median score as the dividing line. (Figure 4A–4C) illustrates the distributions of patients categorized as low- or high-risk for LUAD in the training, testing, and entire cohorts, and (Figure 4D–4F) illustrates the survival states. The results indicated that the high-risk patients exhibited a higher mortality rate than low-risk patients. lncRNAs characteristics in various cohorts are depicted on the heatmap (Figure 4G–4I). We identified LINC01117, AC068228.1, LINC00460, LINC00707, and LINC01537 as lncRNAs associated with high risk, while the remaining four were identified as lncRNAs associated with low risk.

A negative relationship between the risk score and the overall survival (OS) rate was noted in all three groups, as illustrated in (Figure 4J–4L) (p < 0.001, p = 0.013, and p < 0.001, respectively). Additionally, the signature-based area under the ROC curves (AUC) (Figure 4M–4O) for 1, 3, and 5 years, was > 0.600, with respective values of 0.722, 0.704, and 0.709 for the entire cohort, signifying strong discrimination performance across the training, testing, and entire cohorts.

Next, we investigated how LUAD patients’ risk scores were correlated with their clinical characteristics. Figure 5A shows a heatmap visualization of the associations between high-risk scores and tumor stages (stages III–IV vs. I–II, p = 0.001), T stages (T3+4 vs. T1+2, p = 0.023), and N stages (N+ vs. N0, p = 0.001). In addition, the risk score, age, sex, and stage were employed to generate time-dependent ROC curves for OS prediction at 1, 3, and 5 years (Figure 5B–5D). Notably, the risk score’s concordance index (C-index) was close to 0.7 (Figure 5E), demonstrating strong predictive ability when contrasted with conventional clinicopathological markers. PCA was implemented for all genes, anoikis-related genes, anoikis-related lncRNAs and risk-related lncRNAs (Figure 5F–5I), and the findings proved that the signature was distributed clearly, which is a positive indicator that the ARLR signature is effective.

Comprehensive univariate and multivariate analyses of the entire cohort showed a robust link between the signature risk score and PFS, DSS, and OS (Figure 6A–6F). Notably, patients categorized as high-risk had a lower PFS compared to those categorized as low-risk (p = 0.005, Figure 6G). There was a substantial reduction in the DSS rate among individuals at high risk (p < 0.001, Figure 6H). This provided more evidence that the risk signature is an accurate indicator of LUAD patients’ prognoses.

Nomogram development and assessment using the risk score

The nomogram was developed by evaluating the risk score in conjunction with other clinical parameters, which include age, sex, and T, N, and M stages (Figure 7A). The model’s reliability was evaluated by creating time-dependent ROC curves at 1, 3, and 5 years of follow-up, with corresponding AUCs of 0.742, 0.725, and 0.738 (Figure 7B). Furthermore, the nomogram forecasts and actual values were found to be in excellent agreement according to the calibration curves (Figure 7C). This suggests that the nomogram model is more useful in clinical settings for patient prognostic prediction than the risk score alone.

Functional enrichment analysis

We performed coexpression analysis between protein-coding genes and lncRNAs from the LUAD cohort and visualized the results in a volcano plot to learn more about the biology of the signature’s anoikis-related lncRNAs (Figure 8A). The GO study also showed that genes associated with anoikis could be involved in several of biological processes (BPs), cellular components (CCs) and molecular functions (MFs), which include positive regulation of secretion, receptor-ligand activity, collagen-containing extracellular matrix (ECM) organization, and signaling receptor activator activity (Figure 8B). To examine the pathways enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) across the various categories, GSEA was implemented. The results illustrated that the PROTEASOME GLYCOSAMINGLYCAN BIOSYNTHESIS CHONDROTIN SULFATE, and ECM RECEPTOR INTERACTION, were predominantly enriched in the high-risk group (FDR < 0.05), whereas the ABC TRANSPORTERS pathway showed significant enrichment in the low-risk patients (FDR < 0.05) (Figure 8C).

TMB correlation with the risk score

In addition, 15 tumor mutations were identified in the high- and low-risk groups, and we demonstrated that the frequency of TP53 mutations was higher in the low-risk patients (high-risk: 40%, low-risk: 51%), whereas the high-risk patients had a greater frequency of KAS mutations (high-risk: 33%, low-risk: 22%, Figure 9A, 9B). According to survival analysis, patients who have a high TMB may experience a reduced duration of survival in comparison to those with a low TMB (p = 0.026, Figure 9C). We subsequently examined survival in terms of TMB and risk scores, and the groups with high TMB and low risk exhibited the best OS (Figure 9D).

Immune function and screening for anticancer drugs

An analysis of immune-related functions revealed, via heatmap, that the low-risk patients had significantly more active type II interferon (IFN) and HLA as opposed to those at high risk (Figure 10A). Furthermore, an extensive array of immune cells exhibited a negative association with the risk score, as determined by the evaluation conducted by various platforms (for example, B cells and CD8+ T cells in EPIC) (Figure 10B). Furthermore, it was observed that the immune score of the low-risk patients exceeded that of the high-risk patients. However, the stromal and ESTIMATE scores showed no significant variation between the two groups (Figure 10C). The findings suggest that immune-associated functions may be more prevalent in low-risk patients than those in high-risk individuals.

An investigation was conducted in LUAD patients to assess the link between risk scores and immunotherapy efficacy to examine the impact of the risk score on antitumor immune response. As demonstrated by violin plots derived from the IPS, low-risk score patients exhibited improved responsiveness to PD1 inhibitor therapy alone (p = 0.027, Figure 11A), CTLA4 inhibitor monotherapy (p = 0.00024, Figure 11B), and combined PD1 and CTLA4 inhibitors (p = 0.021, Figure 11C). Then, we determined the IC50 of several anticancer medications utilizing the pRRophetic package. These drugs included doxorubicin, talazoparib, palbociclib, phenformin, rapamycin and naviroclax. The findings demonstrated that doxorubicin, talazoparib, and palbociclib exhibited greater IC50 values in the low-risk patients contrasted with the other three medicines (Figure 11D–11I). Based on these results, the risk signature may be useful for guiding clinical anticancer treatment.

Further exploration of AC026355.2

To ascertain the role of these lncRNAs in LUAD development, we selected AC026355.2, which has rarely been reported, for further study. We first evaluated the anoikis-related lncRNA AC026355.2 expression in cell lines. Consistent with our hypothesis, AC026355.2 was more highly expressed in adenocarcinoma cell lines (HCC827, PC9, NCI-H1650, A549, and NCI-H1975) (Figure 12A) than in the human embryonic lung cell line MRC-5. Afterward, we transfected AC026355.2 siRNA into the human lung cancer cell line A549. The qRT-PCR results demonstrated that the AC026355.2 expression was substantially decreased in the treated group than in the negative control group (Figure 12B). Then, we used transwell and scratch assays to assess the role of AC026355.2 in LUAD cell migration and invasion. According to the results of the transwell assay, AC026355.2 silencing markedly increased the number of A549 cells that infiltrated the bottom chamber (Figure 12C). Similarly, the scratches persisted in their healing process, and silencing AC026355.2 enhanced the cells’ capacity to migrate into the scratch region (Figure 12D).