A methylation‐related lncRNA‐based prediction model in lung adenocarcinomas

We acquire raw LUAD and normal control RNA sequencing and clinical data fetched from The Cancer Genome Atlas (TCGA) (n = 566).2 Clinical outcomes and lncRNA expression values were gathered for 594 patients. According to the cBioPortal database (https://www.cbioportal.org/↗), 566 LUAD patients' disease‐free survival (DFS) data were obtained. The 165 methylated genes have been analysed in our study and are selected from a literature review. Six genes were methylated in m1A, and eight m5C methylated genes, 35 m6A methylated genes. There 118 m7G methylated genes were selected from gene cards (https://www.genecards.org↗).

The data on RNA‐methylation‐related gene expression levels were acquired from TCGA‐LUAD for the analysis. R software was used to carry out the analysis. Transcripts with log2FC| > 1 and FDR < 0.05 are selected to be differentially expressed genes (DEGs).

We performed Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, to explore the biological function of DEGs and the biological pathway enriched by the DEGs.3

The ‘limma’ software package was performed to analyse the correlation between methylation‐related genes and IncRNAs obtained from the TCGA database. Under screening standard that |R²| > 0.3 and P < 0.001, 1595 methylation‐related lncRNA was gained. Following univariate COX regression analysis was used to analyse the correlation between the survival data of patients in LUAD and 1595 methylation‐related lncRNAs to screen out the lncRNAs with prognostic value. Multivariate COX regression analysis identified 18 methylation‐related lncRNAs which were used to construct the predictive signature of LUAD. In addition to calculating methylation‐related risk scores for patients in LUAD using the following formula.

The validation efficiency of the predictive signature of LUAD was evaluated according to the area under the receiver operating characteristic (ROC) curve we drew. Univariate and multivariate Cox regression analyses were performed to evaluate the independent prognostic value of the Risk model. Cluster heatmaps were further applied to evaluate the distribution relationship between clinicopathological indicators and the predictive signal was analysed to test the accuracy of the prediction model based on methylated lncRNA.

Download ‘term set’ from the LncSEA website and use the ‘Transcription_Factor.gmt’ and ‘ceRNA.gmt’ files to analyse the TFs of these 18 lncRNAs and the ceRNA network.4 The file used by TF network is ‘Transcription_Factor.gmt’, the file used by ceRNA network is ‘ceRNA.gmt’ and the ‘RNA_Protein_Interaction.gmt’ file in ‘term set’ and the go package of R language are used for go enrichment analysis.

We used two websites, TRlnc and lncRNASNP, to analyse the SNP sites of lncRNAs, and a total of two of 18 lncRNAs were retrieved from TRlnc: AF186192.1↗ and LINC00996. Among the 18 lncrnas, only LncRNASNP was retrieved: AF186192.1↗.

For patients with LUAD, the risk score based on the methylation correlation model, age, sex, stage and N stage was incorporated into the analysis. After the analysis, we constructed a nomogram that predicts LUAD patients 1, 3 and 5 years OS after diagnosis. The calibration curve was used to describe the accuracy of predicted OS rates.3

To prove the prognostic ability of methylated‐related‐lnRNA signature independent of clinicopathological parameters, we divided LUAD patients into different subgroups according to age, gender, TNM stage and stage. Kaplan–Meier survival analysis of the subgroups. In addition, internal validation was performed to verify the prognostic value of lncRNA‐signature.

The patients with LUAD obtained from the TCGA database were divided into high‐risk and low‐risk groups according to the risk scores based on methylation‐related lncRNA signature. GSEA was used to analyse the functional enrichment differences of immune pathways in patients with different risk groups.ssGSEA programme was performed by the ‘GSVA’ R package to quantify the differences in 16 immune cell infiltration levels and 13 immune pathway activities between patients in different groups. Expression levels of Immune checkpoint genes in different risk groups were analysed using the Wilcoxon test.

The ‘prophetic’ software package was used to analyse the commonly used immunotherapy drugs in the clinical treatment of LUAD. Wilcoxon signed‐rank test was performed to compare the half maximal inhibitory concentration (IC50) values between the high‐risk group and the low‐risk group to evaluate the effectiveness of the predictive signature in predicting the prognosis of LUAD.

Throughout this article, R software is carried out for all analyses. Using the Wilcoxon test, the difference between tumour and normal tissues was analysed. Risk prognostic models were constructed using the Cox regression algorithm. Multivariate Cox regression algorithms were performed to show the independent prognostic value of the Risk model. The Chi‐square test was applied to draw a heatmap, which reveal the distribution and expression of 18 methylation‐related lncRNAs and clinicopathological factors in different risk groups.5 Using Kaplan–Meier analyses estimates the OS time of patients in different risk groups. R package ‘survivalROC’ was used to generate the receiver operating characteristic (ROC) curves for calculating the area under the curve (AUC values). For ssGSEA, the ‘GSVA’ package was utilized.

From TCGA‐LUAD, we downloaded 594 LUAD and paracancerous tissue data. We screening out 165RNA‐methylation‐related genes, including 6 m1C, 8 m5C, 35 m6A and 116 m7G. Subsequently, the different expression levels of the above four groups of RNA methylation‐related genes in cancer tissues and adjacent tissues were analysed, resulting in a total of 71 DEGs being screened out (Figure 1A–D). In addition, we examined the relationship between the methylation‐related RNA DEGs in each group (Figure 1E–H).

Among the 71 methylation‐related RNA DEGs (DEGs), seven were down‐regulated and 64 were up‐regulated (Figure 1I). KEGG pathway analysis indicated the significant enrichment pathways associated with the 71 RNA methylation‐associated DEGs are RNA polymerase, basal transcription factors, mRNA surveillance pathway, nucleocytoplasmic transport, nucleotide excision repair, and so on (Figure 1J). GO analysis revealed that the function of DEGs was primarily enriched in the regulation of translation, RNA capping, RNA methylation, 7‐methylguanosine RNA capping, and so on (Figure 1K).

The ‘LimmaR’ package was used to identify 1598 methylation‐related lncRNAs. In a univariate regression analysis, it was found that 83 of these variables were related to the prognosis of patients with LUAD. Multivariate Cox regression analysis was performed to screen out 18 methylation‐related lncRNAs (AP002761.1↗, AL118558.3↗, CH17‐340M24.3, AL353150.1↗, AC004687.1↗, LINC00996, AF186192.1↗, HSPC324, AC087752.3↗, FAM30A, AC106047.1↗, AC026355.1↗, ABALON, LINC01843, AL606489.1↗, NKILA, AP001453.2↗ and GSEC). We used these 18 lncRNAs independent predictors to build a prediction model. Figure 3A‐i illustrates the expression levels of 18 lncRNA signatures in the tissue of LUAD. Using the ‘gg alluvial’ R package and Cytoscape, we found 18 pairs of lncRNA‐mRNA co‐expressed in the network (Figure 3A‐ii, |R²| > 0.4 and P < 0.001). Table 1 illustrates the co‐expression relationship. Based on this, we further visualized the corresponding Sankey diagram showing the relationship between prognostic lncRNAs and risk type (Figure 3A‐iii). In total methylation‐related lncRNAs, ABALON, AC004687.1↗, AC087752.3↗, AL606489.1↗, AP001453.2↗, CH17‐340M24.3, GSEC and NKILA were risk factors. AP002761.1↗, AL118558.3↗, AL353150.1↗, LINC00996, AF186192.1↗, HSPC324, FAM30A, AC106047.1↗, AC026355.1↗ and LINC01843 were identified as protective lncRNAs. The risk scores were analysed as follows: riskscore = (−0.730145377 × AC004687.1↗ expression) + (1.411577652 × AC087752.3↗ expression) + (0.763757779 × AC034236.2↗ expression) + (0.516345501 × NKILA expression) + (−0.922753388 × AL118558.3↗ expression) + (0.78048419 × ABALON expression) + (0.493697701 × CH17‐7340 M24.3 expression) + (−0.933963185 × FAM30 expression) + (−0.297718694 × LINC01843 expression) + (−0.924673401 × LINC00996 expression) + (0.751538891 × AP001453.2↗ expression) + (−0.457662074 × AL353150.1↗ expression) + (0.60638866 × AL606489.1↗ expression) + (−0.548720538 × AC026355.1↗ expression) + (−0.666023466 × AF186192.1↗ expression) + (−0.854079483 × AP002761.1↗ expression) + (−0.539050891 × HSPC324 expression) + (0.536675573 × GSEC expression).

In order to analyse the role of these 18 lncRNAs in the process of tumourigenesis and development, we downloaded term set from LncSEA website and then used the clusterProfiler package of R language to analyse which term set in each file contains the 18 lncRNAs we need to identify and used the file ‘RNA_Protein_Interaction.gmt’ for enrichment analysis. It was found that the enriched pathways enriched by proteins interacting with 18 lncRNAs were mainly EMT, hypoxia, stemness and proliferation (Figure 2A), among which LINC00996 and AF186192.1↗ were regulated by multiple tumour‐related transcription factors, such as TP53 and TP63 (Figure 2B), and FAM30A and mRNA formed a ceRNA network (Figure 2C‐i,ii). Then TRlnc and LncRNASNP were used to analyse the SNP sites of lncRNAs, and 2319 SNP sites were identified in LILNC00996, 36 of which were risk SNP sites and 205 SNP sites in AF186192.1↗. Conservation indicates that AF186192.1↗ affects 95 conserved miRNAs and 123 non‐conserved miRNAs. Gain/loss indicates that af186192.1 promotes the binding of 149 pairs of miRNA:lncRNA and inhibits the binding of 95 pairs of miRNA:lncRNA (Figure 2D).

According to the cutoff value based on the median risk score we set, LUAD patients with different risk scores counted by the formula were assigned to high‐ and low‐risk groups. The Kaplan–Meier curve was used to estimate the overall survival (OS) time of LUAD patients in different risk score groups. As excepted, the OS time of the high‐risk group was shorter than the low‐risk group (Figure 3B‐i, P < 0.001). As shown in the curve of risk score and scatter plot, LUAD patients with significantly lower survival time had an increased risk score (Figure 3B‐ii,iii). Through the univariate and multivariate COX regression analysis. The lncRNA signature‐based risk score was verified the same independent prognostic value as some common clinicopathological parameters (Figure 3C‐i,ii). The AUC value for risk score, age, gender, stage, T stage, N stage and M stage were 0.790, 0.481, 0.449, 0.651, 0.604, 0.530 and 0.617, respectively (Figure 3C‐iii), which indicated the prognostic value of risk score in the prognosis of LUAD patients was superior to clinicopathological variables. This signature has well predictive value within 3 years. The AUCs value of risk score at 1, 3 and 5 years were 0.699, 0.715 and 0.807, respectively (Figure 3C‐iv). A heat map was constructed to analyse the expression level of 18 methylation‐related lncRNAs and the distribution of clinicopathological variables in different risk groups (Figure 4A). We also plotted a nomogram, incorporating clinicopathological variables and risk scores, which was used to predict the prognosis of patients with LUAD at 1, 3 and 5 years (Figure 4B‐i, P = 0.05). According to the calibration curve for the nomogram (Figure 4B‐ii, iii, iv), we further confirmed the predictive value of the nomogram for the prognosis of patients with LUAD.

To examine whether the predictive signature has an independent prognostic value for the patients in LUAD, first, we grouped LUAD patients according to different clinicopathological factors, and then performed Kaplan–Meier analyses to plot the survival curve for LUAD patients in different risk groups. Kaplan–Meier results show the OS rate of the patients in low‐risk is higher than the patients in high‐risk groups (Figure 5A–K). Predictive signatures can predict the prognosis of LUAD patients independently regardless of clinicopathological variables.

All patients (477)with LUAD acquired from the TCGA database were sorted into two cohorts, namely, the training group and the Test group (n₁ = 240, n₂ = 237), to test the applicability of the predictive signature. In Table 1, we present the demographic characteristics of both groups. As we excepted, The AUC of the 1‐, 3‐ and 5‐year ROC curves of the training group were respectively 0.788, 0.809 and 0.822(Figure 6A). Accordingly, test group survival rates of 1‐, 3‐ and 5 years were 0.708, 0.659 and 0.767 (Figure 6B). Both cohorts showed good predictive performance on ROC curves. The overall survival rate of the high‐risk group is higher than the low‐risk group both in the training group (Figure 6C, P = 5.574e−10) and the test group (Figure 6D, P = 67577e−5).

GSEA was performed to analyse the different expression levels of tumour‐ and immune‐related pathways in the high‐ and low‐risk groups. Results show that tumourigenesis and immune‐related pathways, such as pentose‐phosphate‐pathway, cell cycle pathway, N‐glycan‐biosynthesis, glycolysis‐gluconeogenesis, p53‐signalling‐pathway, notch signalling pathway, TGF‐beta‐signalling‐pathway, pancreatic cancer, thyroid cancer and prostate cancer were significantly enriched in the high‐risk group (Figure 7A).

We also performed ssGSEA for revealing the immune landscape of the LUAD patients with low or high‐risk scores. The infiltration difference of 16 immune cells was analysed, and the results showed in Figure 7B‐i that eight immune cell infiltration differences between the high‐risk group and the low‐risk group. Among these eight immune cells, in the low‐risk group, immature dendritic cells (iDCs), mast cells and T helper cells were significantly infiltrated more than in the high‐risk group. Next, the expression level of representative immune signalling pathways in LUAD patients within different risk groups was evaluated. The results show that the immune function score of human leukocyte antigen (HLA), type IIIFN response was significantly higher in the low‐risk group (Figure 7B‐ii). These results indicate that Immune function activation seems to be more prone for High‐risk individuals. Cancer is associated with immune checkpoint genes. Based on our analysis, some immune checkpoint genes were expressed differently in the two risk groups (Figure 7C).

We performed principal component analysis to visualize the spatial distribution of LUAD patients in high‐ and low‐risk. We used PCA profiles to depict patients based on total methylation‐related genes, total methylation‐related lncRNA, and lncRNA signature. According to the prediction signal constructed using lncRNA‐signature, the lowest and highest risk groups showed the most evident differences (Figure 8A‐i, ii, iii). This indicated that our prediction signal based on total methylation‐related lncRNA could distinguish well between high‐risk and low‐risk groups.

An analysis of the correlation between the predicted signals and conventional LUAD chemotherapy drugs was performed. The high‐risk group IC50s for A.443654, ErIotinib, and S.Trityl.L. Steine was higher than the low‐risk group (Figures 8B‐i, ii, iii). Nevertheless, low‐risk patients received less ATRA, Methotrexate, and Lenalidomide (Figures 8B‐iv, v, vi). Data from these studies can be used to reveal individualized treatment options for LUAD patients.

Disease‐free survival (DFS) is useful when screening LUAD patients. Also, methylation‐related lncRNAs for total lncRNAs were developed as a predictor of DFS. From the cBioPortal database, we downloaded DFS data for 567 patients with LUAD. Based on univariate analysis, we identified 22 RNA methylation‐related lncRNAs that were considerably associated with DFS. As a result of multivariate COX analysis, we obtained four lncRNAs associated with RNA methylation and used these four lncRNAs to construct predictive features. Calculate the risk score as follows:riskscore = (1.127506988 × AC127502.2↗ expression) + (0.595996496 × LINC00511 expression) + (−0.661724435 × AL139289↗.2expression) + (0.55340284 × UBE2Q1‐AS1 expression). LUAD patients with different risk scores, calculated with the formula, were divided into high‐ and low‐risk groups based on the median risk. Kaplan–Meier survival curve analysis shows the results that DFS in a high‐risk group is shorter (Figure 8C‐i). The AUC of 1‐, 3‐ and 5‐year survival were 0.686, 0.74 and 0.734 (Figure 8D‐i). Two internal cohorts(n₁ = 148; n₂ = 150) of 298 LUAD patients were randomly assigned to verify the predictability of DFS. The DFS and AUC trends results, calculated with the Kaplan–Meier, in the firster internal cohorts (Figure 8B‐ii, iv) or second internal cohort (Figure 8D‐iii, vi) were consistent with the whole LUAD patients dataset.