A innovative prognostic symbol based on neutrophil extracellular traps (NETs)-related lncRNA signature in non-small-cell lung cancer

Identification of the candidate NETs-related genes

DEGs between NSCLC patients and healthy people from TCGA were analyzed (Figure 2A, 2B). Simultaneously, we obtained the results of DEGs from the transcriptome RNA microarray in the NETs treated and untreated group (Figure 2C, 2D). Then, with | log2 (fold change) | ≥ 1 and FDR < 0.05, we performed a comprehensive analysis on 1115 prognostic NSCLC-related DEGs and 834 NETs-related DEGs, ultimately obtaining 119 genes that were mainly were correlated with NSCLC and NETs (Figure 3A).

Figure 2. Identification of the candidate prognostic NSCLC-related differentially expressed genes (DEGs) and NETs-related DEGs. (A, B) Heatmap and volcano plot of DEGs between NSCLC patients and healthy people that from TCGA datasets in the training cohort. (C, D) Heatmap and volcano plot showing DEGs in the transcriptome RNA microarray of A549 cells treated with or without NETs for12 h.

Figure 3. Identification and analyses of lncRNAs correlated with NETs. (A) Venn diagram identifying 119 overlapping NETs-related genes between TCGA and transcriptome RNA microarrays. (B) Forest plots demonstrating the results of the univariate Cox regression analysis between NETs-related lncRNA expression and overall survival (OS). (C) Heatmap revealing the NETs-related lncRNAs that were associated with OS by univariate Cox regression analysis in NSCLC patients and healthy people. (D) The coexpression network between 12 NETs-related lncRNAs and prognostic DEGs. Red indicates NETs-related lncRNA, and the purple indicates prognostic DEGs. (E) Sankey diagram showing the relationships among 12 NETs-related lncRNAs, prognostic DEGs and risk types.

Identification of the prognostic NETs-related lncRNA signature

A total of 119 NSCLC and NETs-related genes were identified, as summarized above, and 13413 lncRNAs were referred to in TCGA. A NETs-related gene lncRNA co-expression network was constructed to identify the NETs-related lncRNAs, and 1039 NETs-related lncRNAs were selected (|R²| ≥ 0.3 and P ≤ 0.001). Subsequently, according to the univariate Cox regression analysis, 19 NETs-related lncRNAs had prognostic value for patients with lung cancer (P < 0.05, Figure 3B, 3C). As a result of Lasso Cox regression, 12 NETs-related lncRNAs were identified (Supplementary Figure 1), among which three lncRNAs (AP000695.2, MYOSLID, and NKILA) were considered as harmful prognostic factors, while the others (AC020765.2, AC090152.1, AF131215.5, AL035587.1, AP004608.1, MMP2.AS1, SCAMP1.AS1, SNHG10, and TRG.AS) were considered as favorable prognostic factors (Figure 4A–4L). Ultimately, a NETs-related co-expression network consisting of 12 lncRNAs revealed the correlations between the genes and lncRNAs referenced above (Figure 3D, 3E).

Figure 4. Survival data analyses of 12 NETs-related lncRNAs in the training cohort. (A–L) Kaplan-Meier survival analysis of 12 nominated lncRNAs. Three NETs-related lncRNAs (AP000695.2, MYOSLID, NKILA) were as independent adverse factors for NSCLC. The other lncRNAs, that represented independent favorable factors, were as follows: AC020765.2, AC090152.1, AF131215.5, AL035587.1, AP004608.1, MMP2.AS1, SCAMP1.AS1, SNHG10, and TRG.AS.

Establishment of a prognostic model in the training cohort

A prognostic model based on the 19 lncRNAs referenced above was built following univariate Cox regression analysis, and 12 lncRNAs were determined through Lasso Cox regression based on 19 lncRNAs. The formula, risk score = 0.044 * expression level of AC020765.2 + 0.079 * expression level of AC090152.1 + 0.133 * expression level of AF131215.5 + 0.067 * expression level of AL035587.1 + 0.044 * expression level of AP000695.2 + 0.059 * expression level of AP004608.1 + 0.216 * expression level of MMP2.AS1 + 0.006 * expression level of MYOSLID + 0.076 * expression level of NKILA + 0.018 * expression level of SCAMP1.AS1+ 0.027 * expression level of SNHG10+ 0.063 * expression level of TRG.AS, was utilized to calculate the risk score. By evaluating these 12 NETs-related lncRNAs, we were able to acquire the risk score of each patient. The patients were divided into a high-risk group (n=241) and a low-risk group (n=241) in accordance with the risk score (Figure 5C). Next, the Kaplan-Meier curve suggested that patients in the high-risk group had a significantly worse OS than their low-risk counterparts (Figure 5A, P < 0.05). Additionally, PCA indicated that the patients in different risk groups were distributed in two directions (Figure 5B). Patients with high risk had a higher probability of earlier death than those with low risk (Figure 5D, 5E).

Figure 5. Prognostic analyses of the 12 NETs-related lncRNA signature model in the training cohort. (A) Kaplan–Meier survival analysis based on 12 NETs-related lncRNAs in the high-risk group and low-risk groups. (B) Principal component analysis (PCA) plot of the training cohort. (C) The distribution and median values of the risk scores in different groups. (D) The distributions of survival time in the training cohort. (E) Heatmap of the expression levels of the 12 lncRNAs related to OS.

Verification of the 12 NETs-related lncRNA signature in the validation cohort

To test the strength of the model built from the training cohort, the patients in the validation cohort were also categorized into high- and low-risk groups by the median value calculated with the same formula as that of the training cohort. Similar to the results obtained from the training cohort, patients in the high-risk group, had a reduced survival time, compared with those in the low-risk group (Figure 6A, P < 0.05) and were more likely to die earlier (Figure 6D, 6E). Additionally, PCA confirmed that patients in the two subgroups were distributed in discrete directions (Figure 6B).

Figure 6. Prognostic analyses of the validation cohort. (A) Kaplan–Meier survival analysis. (B) PCA plot. (C) The risk scores in the different groups. (D) The survival time in the training cohort. (E) Heatmap of the expression levels of the 12 lncRNAs related to OS.

Independent prognostic value of the 12 NETs-related lncRNA signature in both two cohorts

Univariate and multivariate Cox regression analyses were carried out among the available variables to determine whether the risk score was an independent prognostic predictor for OS. In univariate Cox regression analyses, the risk score was significantly associated with OS in both the training and validation cohorts (HR = 3.715, 95% CI = 2.622-5.263, P < 0.05; HR = 3.539, 95% CI = 1.841-6.802, P < 0.05, respectively) (Figure 7A, 7C). After correction for other confounding factors, the risk score still proved to be an independent predictor for OS in the multivariate Cox regression analysis (training cohort: HR = 3.316, 95% CI = 2.307-4.767, P < 0.05; test cohort: HR = 2.908, 95% CI = 1.496-5.652, P < 0.05; Figure 7B, 7D).

Figure 7. Univariate and multivariate Cox regression analyses for OS. (A, B) Analyses of the training cohort. (C, D) Analyses of validation cohort.

Validation of the expression of the 3 adverse prognostic lncRNAs at the NSCLC cell level

As mentioned before, survival data analyses clearly revealed that AP000695.2, MYOSLID, and NKILA were unfavorable prognostic factors. Here, we conducted further verification to understand the characteristics of these 3 lncRNAs at the cell level. As shown in Figure 8A–8C, the levels of AP000695.2 and NKILA in A549 cell were much higher than those in human normal epithelial lung cells (BEAS-2B), whereas the expression of MYOSLID in the A549 cells was lower. The expression levels of 3 lncRNAs, namely, AP000695.2, MYOSLID, and NKILA, in NSCLC cells (H1299, SK-MES-1, H1703) was distinctly higher than those in BEAS-2B, as detected by RT-qPCR (Figure 8D–8F). Therefore, we tested the 3 lncRNAs that were expressed at higher levers in A549 and SK-MES-1 cells following NETs treatment or not for 12 h (Figure 9A–9F), and the results showed that the expression levels of MYOSLID and NKILA were elevated in both NSCLC cells lines after NETs treatment for 12 h (P < 0.05).

Figure 8. Validation of the expression levels of 3 adverse prognostic lncRNAs in different NSCLC cells by RT-qPCR. (A–F) Expression levels of lncRNA AP000695.2, MYOSLID, and NKILA in human normal lung epithelial cells (BEAS-2B) and NSCLC cells (A549, H1299, SK-MES-1, and H1703), and presented as bar charts. ^*P < 0.05.

Figure 9. Validation of the expression levels of 3 adverse prognostic lncRNAs treated with NETs by RT-qPCR. (A–F) Expression levels of lncRNA AP000695.2, MYOSLID, and NKILA in NSCLC cells (A549, and SK-MES-1) treated with or without NETs for 12h, and presented as bar charts. ^*P < 0.05.

Cell culture

The human bronchial epithelial cells line (BEAS-2B), and human lung adenocarcinoma cells (A549 and H1299), human squamous cell carcinoma cell lines (SK-MES-1 and H1703), purchased from the Type Culture Collection of the Chinese Academy of Sciences, Shanghai, China, were all grown in high glucose DMEM medium (BI, Israel) supplemented with 10% fetal bovine serum (Gibco, Grand Island, USA), penicillin (100 U/ml), and streptomycin (100 μg/ml). The cells were incubated in a humidified atmosphere containing 5% CO² at 37° C.

NETs isolation and treatment

NETs were obtained from rat neutrophils following phorbol 12-myristate-13-acetate (PMA) 100 uM stimulated 4h as previous studies [34, 35]. A549 cells were randomly divided into two groups (each group n=3): without NETs-stimulated group, NETs-stimulated 12h group. TRIzol reagent (Invitrogen, Carlsbad, USA) was used to isolate and purify the total RNA of A549 cells. Then, RNA samples, A549 cells with or without NETs treated 12h, were sent to the laboratory of the OE Biotech Company (Shanghai, China) and the microarray profiling was carried out.

Data and sample collection

The RNA sequencing (RNA-seq) data of 1037 NSCLC patients and 108 healthy people samples were received from the TCGA database (https://portal.gdc.cancer.gov/repository). Additionally, the detailed clinical materials of partial patients were obtained from TCGA. Patients without clinical data recorded were eliminated when the related clinical prognostic analysis was performed. Both NSCLC patients and healthy people samples fell into two groups at random equally and served as training and validation data sets, respectively. The gene expression profiles data from TCGA were normalized, and normalization was not needed. The data from TCGA is openly getable, so that the present study was not required the permission of local ethics committees. The current study abides by the access policies and publication guidelines of TCGA.

The NETs-related genes were obtained from the transcriptome RNA microarray data of A549 cells with or without NETs treated 12h.

Identification of NETs-related genes and lncRNAs

"limma" R package was exploited to identify the differentially expressed genes (DEGs) between NSCLC and healthy people with | log2(Fold Change) | ≥ 1, and FDR < 0.05, in the training cohort, and the differentially expressed NETs-related genes from the transcriptome RNA microarray were screened in the same way. Then, prognostic NETs-related genes were selected by using Venny 2.1.0 based on further analysis of genes aforesaid. Subsequently, Pearson correlation was utilized to measure the relationship between lncRNAs and NETs-related genes. NETs-related lncRNAs were ensured according to the square of correlation coefficient |R²| > 0.3 and P < 0.05.

Establishment and test of a prognostic NETs-related lncRNA signature

NETs-related lncRNAs with prognostic values were picked out via univariate Cox analysis of overall survival (OS). P values were adjusted by Benjamini and Hochberg (BH) correction. The Lasso Cox regression analysis [33, 36] constructed an NSCLC prognostic model to lessen the chance of overfitting as much as possible. The Lasso algorithm and "glmnet" R package were utilized for selecting variable and shrinkage together. The normalized expression array of prognostic NETs-related lncRNAs was the independent variable in the regression, and the dependent variables were the overall survival and status of patients in the training cohort. Then, the Co-expression networks were depicted through Cytoscape software 3.7.2.

Judging by the normalized expression level of each lncRNA and its homologous regression coefficients, the risk scores of the NSCLC patients were gained. The formula was operated: score = sum (each gene’s expression × corresponding coefficient). The patients were divided into high- and low-risk groups severally according to the median value of the risk score. PCA was performed with the "ggplot2" function in the "Rtsne " R package in light of the expression of lncRNAs in the signature. The ideal cut-off expression value was calculated via the "surv_cutpoint" function of the "survminer" R package to achieve survival analysis for each lncRNA.

Real-time quantitative polymerase chain reaction

Total RNA from the BEAS-2B, A549, H1299, SK-MES-1, and H1703 cell, according to the manufacturer’s specification, was separated using Trizol (Invitrogen, Carlsbad, USA). Then, the first-strand cDNA at a volume of 20μL was compound by using the reverse transcription kit and finally diluted to 100uL. Quantitative real-time polymerase chain reaction (RT-qPCR) analyses were processed utilizing an SYBR Green mix in the Real-Time PCR System (TAKARA, Kyoto, Japan). 20μL of PCR reaction was prepared as follows: 10μL of 2×SYBR Green PCR master mix, 1 μL of 10μM of suitable forward and reverse primers, 7μL of RNase-free water, and 2μL of cDNA template. RT-qPCR was performed for 10 minutes at 95° C, followed by 40 cycles of 5 seconds at 95° C and 30 seconds at 60° C. A dissociation melting curve was generated using thermal conditions from 60° C to 95° C. Human GAPDH was employed as an internal housekeeping reference. The sequences for all primers were as follows: NKILA (forward primer): GGCTAGTCTGGCTGGGAGAAGTC; (reverse primer): AGCGTTGTGGGTAGGTTTGGTTTC. MYOSLID (forward primer): TCTGCCTAGTCCTGCTGCCTTC; (reverse primer): ATGGGAAGCTGTGTTCACTTTGGG; AP000695.2 (forward primer): CGGAAGCCACCACATGACCTTG; (reverse primer): TTCCAACCGCATGGGTGAAAGTC. GAPDH (forward primer): GTCAGTGGTGGACTGACCT; (reverse primer): TGCTGTAGCCAAATTCGTTG.

Statistical analysis

The OS between different groups was compared by Kaplan-Meier analysis with the log-rank test. Univariate and multivariate Cox regression analyses were implemented to identify independent predictors of OS. All statistical analyses were finished with R software (Version 3.5.3) or GraphPad Prism 8 (Version 8.0). If not particular remark, P < 0.05 was considered statistically significant, and all P values were two-tailed.