LncRNACNVIntegrateR: a novel framework for correlating long non-coding RNAs with copy number variation abnormalities and disease progression

Design and implementation

The structure and design of the lncRNACNVIntegrateR package is modular and user-centric, ensuring ease of use and flexibility for integrating and analyzing multi-omics data. The package is organized into several core functions, each performing specific tasks such as data input, preprocessing, correlation analysis, risk score model development, and functional enrichment analysis. This modular approach allows users to utilize each function independently or as a comprehensive workflow. Figure 1 illustrates the end-to-end workflow implemented in lncRNACNVIntegrateR, highlighting the major data integration and analysis steps. The analysis begins with three essential input files consisting of the following datasets:

• (1)

  Expression profiles: Provided as a raw count matrix or a tab-delimited file. Rows represent genes (including lncRNAs), and columns represent samples. The first row should include sample identifiers (header), followed by gene identifiers and their corresponding expression values.

• (2)

  Absolute CNV profiles: A matrix or tab-delimited file where rows correspond to CNV features and columns to samples. The header includes sample IDs, followed by rows with CNV identifiers and their absolute copy number values (indicating amplifications or deletions).

• (3)

  Clinical data: A tab-delimited file with clinical annotations. Rows represent samples and columns represent clinical variables, including at minimum: Sample_ID, Survival_Time, and Status.

These datasets are processed through sequential stages of data harmonization, correlation analysis, survival modelling, and downstream functional exploration.

Core functionalities of the package

The lncRNACNVIntegrateR package provides a comprehensive suite of functions, each corresponding to a key step in the analysis pipeline (Fig. 1). These include:

Data preprocessing: The function automatically harmonizes sample IDs across the expression, CNV, and clinical datasets by identifying the common samples shared among all three. It then proceeds with downstream analysis using only these matched samples, ensuring consistency and compatibility across data types.

Correlation analysis: Computes statistical associations between lncRNA expression and CNV status to identify CNV-driven lncRNAs.

Survival association testing: Evaluates the prognostic significance of candidate lncRNAs using clinical survival data.

Signature identification: Filters for the most informative CNV-deregulated lncRNAs based on statistical significance and biological relevance.

Risk score modelling: Constructs and validates a prognostic model using selected lncRNA signatures.

Functional enrichment analysis: Investigates the downstream mRNA targets of identified lncRNAs to interpret their indirect roles in biological pathways and processes.

Each of these functionalities is implemented as a standalone function, allowing users to use each function in a standalone mode in their data analysis (as shown in Fig. 2). The structure and design of the lncRNACNVIntegrateR package is modular and user-centric, ensuring ease of use and flexibility for integrating and analyzing multi-omics data.

The lncRNACNVIntegrateR package processes the three primary input files: expression profiles, absolute CNV profiles, and clinical data (as discussed above) through various dedicated functions. It facilitates comprehensive analysis and integration of the data, encompassing several key tasks. These tasks include data preprocessing and preparation, evaluating the correlation between lncRNA expression and CNVs, and identifying survival-associated lncRNAs. The package further identifies significant CNV-deregulated lncRNA signatures, constructs a risk score model for these signature lncRNAs, and performs functional enrichment analysis on the target genes of the identified lncRNA signatures. Detailed descriptions of the functions provided by this package are outlined below.

Input functions

To streamline data handling, the package provides dedicated input functions for each data type:

• input_expr for loading raw gene expression count matrices (whole transcriptome data),

• input_cnv for importing absolute CNV profiles, and

• input_clin for clinical data, including survival time and event status.

In addition, the package supports annotation using the GTF file via the input_gtf function (sourced from GENCODE version 22) and retrieves a curated list of lncRNAs using the input_lncRNAs function. The ‘input_lncRNA_positions’ function helps incorporate genomic positions of lncRNAs for downstream analyses and is available in the example extdata directory of the package.

These datasets serve as inputs for the analysis and are passed into the package using the corresponding input functions. For the GBM and COAD case study, gene expression and clinical data were obtained from TCGA via the GDC portal (https://gdc.cancer.gov/), while absolute CNV data were downloaded from the UCSC Xena Browser (http://xena.ucsc.edu/). These datasets serve as inputs for the analysis and are passed into the package using the corresponding input functions.

Pre-processing

The ‘preprocessing_and_preparation’ function conducts preprocessing on individual data types. The first step involves identifying common samples across the three data types, extracting the relevant information from each data type for these common samples, and then proceeding with the preprocessing. For gene expression data, first, data normalization is performed using the DESeq2 Variance Stabilizing Transformation (VST) (Love, Huber & Anders, 2014). Further, the lncRNAs and protein-coding genes (PCGs) are mapped using GTF information (obtained from gencode version 22, https://www.encodeproject.org/files/gencode.v22.annotation). Genes categorized in the GTF file under labels such as “processed_transcripts,” “3prime_overlapping_ncRNA,” “sense_intronic,” “sense_overlapping,” “antisense,” and “lincRNA” are recognized as lncRNAs, whereas those labeled as “protein_coding” are identified as PCGs. Following the extraction of lncRNAs and PCGs expression profiles, the CNV profiles for the lncRNAs are then extracted from the total absolute CNV data. Important parameters, including survival-related information, are extracted from the clinical data.

Expression data undergoes lncRNA and PCG annotation based on the GTF file. Genes annotated as “processed_transcript,” “3prime_overlapping_ncRNA,” “sense_intronic,” “sense_overlapping,” “antisense,” and “lincRNA” are considered lncRNAs, while those labeled “protein_coding” are designated as PCGs. Expression data is normalized using the VST method from the DESeq2 package (Love, Huber & Anders, 2014). Simultaneously, lncRNA-specific CNV profiles are extracted from the absolute CNV dataset. Clinical parameters such as survival time and status are also extracted at this stage.

lncRNA-CNV correlation

After the data preprocessing and preparation step, the processed lncRNA expression and CNV profiles are utilized to compute the Pearson correlation coefficient (PCC) for each lncRNA using the ‘Correlation_lncRNAs_and_CNVs’ function. The correlation distribution is then plotted against a random distribution to evaluate any correlation. Subsequently, prognosis-related lncRNA signatures are extracted using a univariate Cox proportional risk regression model with a p-value cut-off <0.05 under function ‘extract_survival_related_significant_lncRNA’. Significant lncRNAs exhibiting amplification or deletion are further identified using a violin-box plot based on the Kruskal-Wallis test (cut-off p-value < 0.05).

Risk score model construction

The list of screened lncRNAs displaying significant CNV amplification/deletion, identified through the Kruskal-Wallis test in the ‘visualization_of_lncRNAs_and_CNVs_association’ function, is refined using a p-value cut-off <0.05. Key prognostic lncRNAs (pr-lncRNAs) are selected to construct a risk score model using the PredictABEL (a R package available at https://cran.r-project.org/package=PredictABEL, version 1.2-4). With the identified key lncRNAs, the model can stratify patients into two risk groups using their respective risk scores. This functionality is implemented through the ‘Risk_score_model_development’ function. A survival plot is created to determine if these two groups differ in survival probabilities (p-value < 0.05).

The list of screened lncRNAs showing significant CNV amplification/deletion, determined by the Kruskal-Wallis test, was further refined using a multivariate Cox regression model with a p-value < 0.05. Key prognostic lncRNAs (pr-lncRNAs) were chosen to build a risk score model using the PredictABEL R package (https://cran.r-project.org/package=PredictABEL). Based on these key lncRNAs, this model can classify patients into two risk categories using their respective risk scores. A survival plot was created to determine if these two groups differ in survival probabilities (p-value < 0.05).

Functional enrichment analysis

To explore the function of the identified pr-lncRNAs, the tool extracts highly correlated genes associated with these pr-lncRNAs. To analyze the functional enrichment of PCGs correlated with these lncRNAs, we utilize the enrichR package (Kuleshov et al., 2016) within the ‘Correlation_with_PCGs_and_functional_enrichment’ function.

Case study 1: identification of CNV-driven lncRNA signatures in glioblastoma using lncRNACNVIntegrateR

To demonstrate the utility of lncRNACNVIntegrateR, we analyzed multi-omics data from GBM, the most aggressive primary brain tumor with a 5-year survival rate of 5.8% (Zhu et al., 2021). The lncRNACNVIntegrateR pipeline was applied to integrate gene expression, CNV, and clinical data from the GBM dataset, accessed via the ‘input_expr’, ‘input_cnv’, and ‘input_clin’ functions. Data preprocessing was performed using the ‘preprocessing_and_preparation’ function, which extracted lncRNA and PCG expression profiles mapped to a GTF reference file.

The ‘Correlation_lncRNAs_and_CNVs’ function was used to evaluate associations between lncRNA expression and CNV profiles. This analysis revealed statistically significant correlations (p < 0.05), which were visualized as correlation distribution plots. The Pearson correlation distribution is shown in Fig. 3A, and the genome-wide distribution of lncRNA CNVs is presented in Fig. 3B. Next, the ‘extract_survival_related_significant_lncRNA’ function was applied to identify lncRNAs significantly associated with patient survival. The ‘visualization_of_lncRNA_and_CNVs_association’ function further confirmed CNV alterations, amplifications and deletions for these lncRNAs, as illustrated by representative violin plots.

By applying the Kruskal-Wallis test to compare CNV groups (diploid vs. amplified/deleted), and selecting those with statistically significant differences, we identified 21 lncRNAs significantly associated with both CNV alterations and patient survival. These include C1orf229, C3orf35, C9orf147, FENDRR, FTX, HCG18, LINC00115, LINC00265, LINC00271, LINC00324, LINC00354, LINC00641, LINC00645, LINC00667, LINC00676, LINC00858, LINC00907, LINC00963, LINC00970, MIAT, NEAT1, PCA3, RNU6ATAC35P, SNHG6, SNHG7, and SNX29P2.

A risk score model was then developed using the ‘Risk_score_model_development’ function. The model achieved a strong predictive performance with an AUC of 0.794 (95% CI [0.722–0.867]), as shown in Fig. 3C. Further, Kaplan-Meier survival analysis based on the expression of these 21 lncRNA signatures stratified GBM patients into high-risk and low-risk groups, demonstrating a statistically significant difference in overall survival (p < 0.05, Fig. 3D).

Using the full functionality of lncRNACNVIntegrateR, integration analysis of GBM datasets identified 21 lncRNA signatures for GBM as described above. Among these, LINC00115, LINC00324, LINC00641, LINC00963, MIAT, NEAT1, SNHG6, and SNHG7 play significant roles in glioma progression and therapy, suggesting their potential as biomarkers and therapeutic targets for glioblastoma (Chen, Lin & Li, 2023; Chen et al., 2020; Liang et al., 2022; Nie et al., 2021; Rao et al., 2024; Tang et al., 2019; Yang et al., 2020; Ye et al., 2020). Thereby validating the findings and demonstrating the tool’s reliability in this study. In addition, several other lncRNAs identified in this study, including LINC00271, LINC00265, LINC00858, and LINC00667, have also been implicated in the pathogenesis of various cancer types beyond GBM. The risk score model achieved an exceptional AUC of 0.794 with a 95% Confidence Interval [0.722–0.867]. The survival analysis also effectively distinguished high-risk and low-risk groups and verified the model’s accuracy.

Functional enrichment analysis (using function ‘Correlation_with_PCGs_and_functional_enrichment’) of the top correlated PCGs with the identified lncRNAs revealed significant enrichment in Gene Ontology (GO) categories, providing valuable insights into their potential biological functions. Among these, several key Biological Processes (BP) were identified to be enriched, which include ncRNA processing, p53-mediated signal transduction, nonsense-mediated mRNA decay, and regulation of gene expression, as illustrated in Fig. 3G. Corresponding Cellular Component (CC) and Molecular Function (MF) categories were also enriched, supporting the regulatory influence of these lncRNAs on diverse cellular mechanisms. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis demonstrated that these lncRNA-targeted genes are involved in critical signaling pathways, such as the Rap1 signaling pathway, regulation of the actin cytoskeleton, and ABC transporter pathways (Fig. 3H). Collectively, these results highlight the involvement of lncRNA-associated PCGs in essential oncogenic processes, reinforcing their potential as biomarkers and therapeutic targets in GBM.

Case study 2: identification of CNV-driven lncRNA signatures in colon adenocarcinoma using lncRNACNVIntegrateR

To further evaluate and demonstrate the utility of lncRNACNVIntegrateR, we analyzed multi-omics data for COAD dataset, the predominant subtype of colorectal cancer, which ranks as the fourth most diagnosed cancer and second leading cause of cancer-related deaths in the United States (Benson et al., 2018). Cancer prevention, screening, and treatment advances have significantly reduced CRC incidence and mortality rates (Edwards et al., 2010). However, the prognosis for patients with advanced colon cancer remains poor (Sadanandam et al., 2013), with 90% of these cases being COAD (Hajmanoochehri et al., 2014).

Effective prognostic stratification is essential for improving patient outcomes in COAD. Using lncRNACNVIntegrateR, we integrated gene expression, CNV, and clinical data from COAD dataset. Significant correlations between lncRNA expression and CNV profiles were observed (p < 0.05), as shown in the correlation distribution (Fig. 4A) and the genome-wide CNV distribution plot (Fig. 4B). Applying this framework, the analysis identified 14 distinct CNV-driven lncRNA signatures, including SNHG7, TMEM191A, PCAT6, MIR210HG, MIAT, LINC00847, LINC00486, LINC00115, GAS5, FENDRR, FAM85B, C7orf66, C10orf95-AS1, and BCYRN1. A risk score model achieved an AUC of 0.707 (95% CI [0.631–0.782]), demonstrating reliable prognostic performance (Fig. 4C). Further, based on these 14 lncRNA signatures, the stratification of COAD patients into high- and low-risk groups (Fig. 4D), with Kaplan-Meier survival analysis confirming significant survival differences.

Eight lncRNAs, including SNHG7, PCAT6, MIR210HG, MIAT, LINC00115, GAS5, FENDRR, and BCYRN1 have been previously linked to colorectal cancer progression, prognostic implications, and therapeutic response (Pan et al., 2023; Han et al., 2024; Ruan et al., 2019; Liu et al., 2018). Feng et al. (2020), Xie et al. (2022a), Yang et al. (2023), Gu et al. (2018) corroborating the tool’s ability to identify clinically relevant biomarkers. These results underscore lncRNACNVIntegrateR’s efficacy in uncovering CNV-driven lncRNA signatures with prognostic and therapeutic potential in COAD.

Figures 4E–4H showcases the results of functional enrichment analysis for the top correlated genes associated with these signature lncRNAs, emphasizing their roles in COAD tumorigenesis. Functional enrichment analysis revealed significant involvement of the correlated genes in key biological processes such as ncRNA processing, nonsense-mediated mRNA decay, translation, and gene expression. KEGG pathway analysis further highlighted enrichment in critical cancer-related pathways, including JAK-STAT signaling, necroptosis, viral infection, and central carbon metabolism in cancer. These findings demonstrate the utility of lncRNACNVIntegrateR in identifying prognostic lncRNA signatures and uncovering their functional roles in colorectal adenocarcinoma, supporting their relevance in disease progression and potential therapeutic targeting.