Simulation results

To evaluate the performance of our proposed method, simulations of mRNA and miRNA expression data in 100 samples were performed in different scenarios (S1 Table). Notably, we focused on the sensitivity and the positive predictive value (PPV) because the number of ceRNA triplets was small among all possible combinations of triplets.

In the beginning, we presumed the sample distribution of each gene follows normal distribution because about 98% of genes’ sample distributions passed the normality test based on 9,835 samples in The Cancer Genome Atlas (TCGA) pan-cancer atlas (S1 Fig). Therefore, we firstly presumed synthetic expression data were generated from a multivariate normal distribution with a mean value of 0 and a covariance matrix whose entries are 0.9. This ensured that the ceRNAR algorithm supports the hypothesis that pairs of ceRNA binding partners of a specific miRNA are highly correlated and that it works well to sensitively identify them among a pool of target pairs. However, it is unclear whether such an event between two target genes would occur at the lower or higher expression of a specific miRNA. Therefore, five scenarios were designed to capture how the molecular elements within a triplet interact with each other, and simulations of different parameters were performed. Notably, the number of identified ceRNAs dropped as the window size increased (S2 and S3 Tables, Fig 1). This is because more uncorrelated samples were included in the analysis when longer window sizes were used, especially larger than 30%. In other words, higher noise was included in the analysis, resulting in difficulty in identifying true ceRNA triplets. Next, we simulated different scenarios in which the ceRNA peak was located at different miRNA expression levels (scenarios 1 to 4). As shown in S2 Table and Fig 1A, the performance of our proposed algorithm was highly similar across the four scenarios with true ceRNA signals. Such performance also was shown when removing the permutation test (Fig 1B and S3 Table). Lastly, to reduce the calculation complexity and computation time, we evaluated the performance of the algorithm without the random walk step, which we called the “fast” version. Only minor differences were observed in these two versions (0.1 to 0.25 in terms of PPV value, S2 and S3 Tables), suggesting that both versions were able to identify true ceRNA triplets without reporting high false positive results (Fig 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Performance of our proposed method in four scenarios (1 to 4).

(A) The complete version. (B) The fast version. The numbers in blue represent the average number of identified ceRNA-miRNA triplets after 100 simulations.

https://doi.org/10.1371/journal.pcbi.1010497.g001

In addition to the analytical parameters for the algorithm, a simulation of different proportions of the correlated samples was performed. Because no major differences were observed in the simulated scenarios and running versions, only scenario 3 was evaluated in this round of testing. As shown in S4 Table, four settings for the proportion of correlated samples were analyzed (10%, 20%, 30%, and 40%). Notably, the proposed algorithm showed similar performance in the three highest settings (Fig 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Performance of our proposed method in three extended cases.

(A) Different proportions of correlated samples in scenario 3 using the fast version. (B) Comparison between complete and fast versions in the null scenario (i.e., scenario 5). (C) Different distances for merging peaks when window size is equal to 40 using the complete version. The numbers in blue represent the number of identified ceRNA-miRNA triplets.

https://doi.org/10.1371/journal.pcbi.1010497.g002

To further examine whether more false positive ceRNA triplets are reported in the fast version, a null scenario without a true ceRNA signal was simulated (scenario 5). The results showed that the negative predictive values of these versions were all higher than 0.9 (S1 and S2 Tables), suggesting there is a low chance of identifying false positive signals in both versions of the algorithm (Fig 2B). Subsequently, we examined how much time can be saved by using the fast version of the proposed algorithm, which omitted the random walk step. On average, the fast version accelerated the algorithm by approximately 76 times (283.967 seconds versus 21,600 seconds), and only slightly higher false positive rates were reported (Fig 2B). Lastly, to evaluate whether different window sizes for merging peaks are critical to the performance of the proposed algorithm, four settings of the proportion of correlated samples and merging window sizes were analyzed. As shown in Fig 2C and S5 Table, the performance of the proposed algorithm was not sensitive to the window sizes for merging peaks.

In addition, to ensure the algorithm only sensitively identified ceRNA events under higher correlation values between them, we also generated simulated data from a multivariate normal distribution with the same mean value but a covariance matrix whose entries are 0.6 (S2A Fig and S6 Table) or 0.3 (S2B and S7 Table). Although the ceRNAR algorithm was still able to detect ceRNA pairs with moderate sensitivity values (0.4 to 0.5) when the correlation values between target genes within a pair were 0.6, it did not work well (sensitivity values were all below 0.25) when the correlation values between target genes within a pair were 0.3. Together, these results support our correlation-based hypothesis that ceRNA events tend to occur when they are highly correlated. However, the above-mentioned results were based on simple and naïve simulations. We also mimicked the distribution of expression from pan-cancer data (9,835 samples), allowing the synthetic data to be generated from a multivariate normal distribution with a randomly selected mean value and a covariance matrix whose entries are also randomly selected. The simulation studies were also performed under different correlation levels. As shown in S3 Fig and S8–S11 Tables, the performance of the ceRNAR algorithm was determined by the level of the correlation values. Nevertheless, it still performed well (sensitivity values > 0.75) when the correlation between genes was high and when the window size was 10 throughout all scenarios (S3B Fig.), suggesting the ceRNAR algorithm is efficient to identify most potential ceRNA pairs when their correlation pattern is relatively high (0.8–0.9) to compete with a specific miRNA, and such a pattern exists in at least 20% of the sample. The optimal parameter settings for real data were also observed. For 100 samples, the best window size was 10, and the best cutoff correlation value for selecting the most significant event was 0.7.

Comparison with other tools

We have compared our tool with other state-of-the-art tools, including SPONGE, JAMI, GDCRNATools, and CERNIA, in terms of their performances using synthetic data and their running time using real data. Fig 3A and S12 Table illustrate that ceRNAR workflow generally outperformed the other tools in terms of sensitivity and PPV in all scenarios when the window size is set to 10 and the correlation cutoff is set to 0.7. It can be mainly observed when the correlation values among correlated genes are relatively high. However, all of the tools could identify valid ceRNA triplets without reporting high false-positive results except JAMI (Fig 3B and S12 Table). GDCRNATools generally had a high sensitivity compared to the other tools, and a lower specificity than ceRNAR, suggesting that it is good for catching ceRNAs but comes with a relatively high rate of false positives. Regarding running time, we used different sample numbers (250, 500, and 100) on a small subset of the pan-cancer dataset with 15 genes which form 105 triplets; we also used different triplet numbers (105, 1,225, and 4,950) on a small subset of the pan-cancer dataset with 250 samples (S4 Fig and S13 Table). Although the ceRNA algorithm was not fast with a large sample size and a large number of triplets compared with SPONGE, GDCRNATools and CERNIA, it was slightly faster than JAMI.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Performance comparisons between ceRNIA, JAMI, SPONGE, GDCRNATools and ceRNAR packages using synthetic data that mimic the real-world case.

(A) Sensitivity and PPV of various tools at different correlation levels for 100 samples of 105 pairs of target genes under four scenarios (1 to 4). (B) Specificity and NPV of various tools at different correlation levels for 100 sample of 105 pairs of target genes under null scenario (scenario 5). “cor_levels” represents four different correlation levels (all: 0.3–0.9; high: 0.8–0.9; low: 0.3–0.4; mid: 0.5–0.7).

https://doi.org/10.1371/journal.pcbi.1010497.g003

Application to TCGA cancer cohort datasets

To further validate the applicability and robustness of the ceRNAR algorithm, we also applied the algorithm to two TCGA-derived lung cancer cohorts–TCGA-LUAD and TCGA-LUSC (S14 Table). The top bridging miRNAs and the hub genes among ceRNA triplets are shown in S15 and S16 Tables. Intriguingly, the two cancer cohorts (LUAD and LUSC) shared some common triplets involving 53 miRNAs and 905 ceRNAs, which allowed construction of a miRNA-modulated ceRNA regulatory network (S5 Fig). Among them, PLEKHG6 had the largest number of co-expressed ceRNAs, and of the 53 common miRNAs, the top three miRNAs that bridged over 20 ceRNA pairs were hsa-miR-183-5p, hsa-miR-133a-3p, and hsa-miR-142-5p. MAP4K3 was another common hub gene in both datasets, and its bridging miRNA was hsa-let-7c-5p, around which a regulatory network of corresponding ceRNAs was built (S6A Fig), and the expression level related to the regulatory occurrence of its bridging ceRNAs is displayed in S6B Fig. Since lung cancers share some common molecular characteristics, such results demonstrate the applicability and robustness of the ceRNAR algorithm in multiple cancers or diseases.

In addition, we compared our findings against experimentally validated miRNA-gene pairs in the miRSponge database to endorse the potential ceRNAs identified by ceRNAR. As shown in S7 Fig, around 1% (ceRNA pairs that can be validated among total ceRNA founded based on ceRNAR algorithm) of experimentally validated ceRNA triplets were identified in LUAD and LUSC. The low proportion may be attributed to the fact that only 158 ceRNA triplets were analyzed and that those ceRNA triplets may not be expressed in lung tissues. Furthermore, approximately 13–14% of ceRNA triplets with at least one experimentally validated miRNA-target interaction were identified by using the ceRNAR algorithm. A Chi-square test was performed to examine whether the findings from the ceRNAR package were significantly enriched in the TCGA data, and the results showed that the P-values obtained from the LUAD and LUSC datasets were both less than 2.2e-16, suggesting our ceRNAR package can successfully identify previously reported ceRNA triplets.

Pipeline of ceRNAR

The ceRNAR package is written in R (version 4.0.5) and is available in the Github repository. The main pipeline of ceRNAR is illustrated in Fig 4 and contains three major components for the identification and analysis of ceRNA-miRNA triplets:

• Data preprocessing
• Identification of ceRNA-miRNA triplets
• Downstream analyses

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. An overview of our ceRNAR package.

This package includes three main parts. First, “Input data” has three functions: ceRNAputativePairs() for extracting curated miRNA-target lists, ceRNATCGA() for retrieving TCGA data, and ceRNAcustomize() for importing customized data. Second, “Identification of miRNA-ceRNA triplets” involves ceRNApairFiltering() and SegmentClusteringPlusPeakMerging() and is wrapped into a ceRNAmethod() function. Lastly, “Downstream analyses” contains six functions for different tasks, that is, ceRNAFunction(), ceRNAModule(), ceRNASurvival(), ceRNALocation(), ceRNAValidate() and ceRNAIntegrate().

https://doi.org/10.1371/journal.pcbi.1010497.g004

To reduce the computational complexity and time cost, the interactions between miRNAs and target genes were based on two experimentally validated miRNA-target databases (miRTarBase [35] and miRecords [36]) and seven computationally predicted miRNA-target databases (DIANA-micro T-CDS [37], EIMMO [38], miRDB [39], miRanda [40], PITA [41], RNA22 [42] and TargetScan [43]). In the default settings, only those interactions that were validated by experiments and/or predicted by more than half of the databases are retained as target miRNAs and target genes (S8A Fig). Conceptually, the ceRNAR algorithm iteratively goes over each miRNA-target list and runs through each mRNA pair in a list to evaluate the chance of the potential ceRNA event involved. For a specific triplet (i.e., a miRNA and its two targets), their expression vectors are extracted from the original expression matrix. Therefore, a miRNA expression vector, miRNA_m = [m_m1, m_m2, …], and two mRNA expression vectors, mRNA_i = [gene_i1, gene_i2, …] and mRNA_j = [gene_j1, gene_j2, …], are used as inputs into the ceRNAR algorithm to iteratively evaluate whether each mRNA pair is a potential ceRNA event (S8B Fig).

Data preprocessing

To prepare expression data for further analyses, ceRNAR can automatically retrieve TCGA data, including mRNA expression, miRNA expression, and survival data, by entering the cancer acronym [44], but it also supports the use of customized miRNA and mRNA expression matrices that are pre-normalized and formatted according to the instructions. In ceRNAR, we implement two functions to fulfill these two approaches: ceRNATCGA and ceRNACustomize.

Identification of ceRNA-miRNA triplets

To identify miRNA-ceRNA triplets (defined here as a miRNA and two target genes) from expression profiles at a specific miRNA expression level, the ceRNAMethod function can be used, and it contains three modules sequentially: ceRNApairFiltering, SegmentClustering, and PeakMerging (Fig 5). We have two assumptions in this study: (1) the expression levels of two target genes tend to be highly correlated when a possible ceRNA event occurs; (2) such events between target genes of a certain miRNA occur at a specific expression interval of that miRNA. However, for each ceRNA triplet from the real data, it is difficult to know which levels of miRNA expression (low, middle, or high expression intervals) will lead to a high correlation between a pair of target genes. Notably, for one miRNA, the high correlation values between two target genes can only be observed in a specific range of the miRNA expression. Definitely, the expression level of one miRNA can be regulated by many factors, such as compensation and/or other interactors. However, with our current understanding of all miRNAs, it is not feasible to consider all potential regulators of one specific miRNA at the same time. Therefore, for one single miRNA, we adapted the approach of utilizing its expression level as the final output instead of considering all possible confounding factors, which can be regarded as the hidden layers of the miRNA expression value.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Three main modules involved in the identification of ceRNA-miRNA triplets.

First, ‘ceRNA pair filtering’ contains sample sorting based on a miRNA expression vector, data augmentation in terms of correlation calculation through a sliding window, and permutation test by random walk. Next, ‘Segment clustering’ is based on circular binary segmentation to identify certain miRNA expression intervals (i.e., segments) with the highest correlation between potential ceRNA pairs. Short segments defined by our criteria are further merged and whether a segment is a peak or a trough is defined. Lastly, ‘Adjacent peak merging’ checks whether a triplet is a candidate ceRNA through two filtering conditions.

https://doi.org/10.1371/journal.pcbi.1010497.g005

The ceRNApairFiltering method

We adopted the sliding window approach to identify the correlation patterns of the two target genes within a specific range of expression values for the miRNA. That is the reason why we ranked the samples based on their miRNA expression value from each triplet to identify the correlation patterns and provided the number of samples that meet the criteria. Therefore, the purpose of this function is to identify ceRNA-miRNA triplets based on the Pearson correlation coefficient through the sliding window approach (S9 Fig) and a running sum statistic for such values by the random walk approach (S10 Fig). First, the samples are sorted based on their miRNA expression levels. For a putative triplet (gene_i, miRNA_m, gene_j) among N samples, the correlation coefficients of gene expression values between mRNA and miRNA are calculated within each window (i.e., a length of sample size that is always less than N) with predefined window size (w) as follows: (1) Here k is the number of windows and a predefined integer. Technically, the correlation value between two genes is calculated using the gene expression values of all samples. By using a sliding window approach [45], we can artificially create different varieties of a real dataset to increase its size, which is a sort of data augmentation technique [46].

Because the accumulated changes in terms of correlation values between two genes among samples may tend to increase the chance that a gene pair is highly correlated and, further, will affect a specific phenotype, we borrowed the concept of gene set enrichment analysis (GSEA) [47], which captures the accumulated changes in the expression of all genes within a pathway through a random-walk method, to identify a significant triplet according to the number of the samples enriched in such an event. The main idea of the GSEA algorithm is to understand whether the differentially expressed genes are significantly enriched in the samples belonging to the same phenotype (case or control). First, the differentially expressed genes were ranked by using the differences in expression level between the two phenotypes. Next, the GSEA algorithm gave a positive score to the genes located in the pathway, whereas a negative score was assigned to the genes outside of the pathway. One possible scenario to obtain a high score is that the differentially expressed genes were clustered and enriched in one phenotype instead of being randomly distributed. This scenario is the same as what we want to identify for the ceRNA triplet. That is, the highest correlation values were observed in a specific range of the miRNA expression and thus we used the two statistics, P_obey(k) and P_violate(k) to identify the range (S11 Fig). Conceptually, to calculate a score (S) to represent the enriched correlation levels among samples against a specific phenotype, we first rank ordered the k windows to form L = {miRNA₁, miRNA₂, …, miRNA_k} based on the average miRNA expression within each window: (2) Next, the S is computed by walking down the window list to evaluated P_obey(k), which is the proportion of samples whose gene correlations are over 0.3 (i.e., ), weighted by their corresponding correlation and P_violate(k), which is the proportion of the rest of the samples at a given ranking window position k across all samples. The formulas are as follows: (3) (4) S is then defined as the maximum distance from zero of the substations of P_obey(k) from P_violate(k). If the samples with similar correlation values are not enriched at a particular miRNA expression interval, it means those gene pairs are not biologically relevant to compete with miRNA in a miRNA-gene triplet; that is, there is no ceRNA event observed in this triplet (S11 Fig). Finally, we accessed the significance of an observed S by comparing it with the theoretical S computed by randomly permutating the expression of the candidate miRNA_m 1,000 times to provide assessment of significance. When an entire sample of potential triplets is evaluated, these p-values will be adjusted by the false discovery rate, which provides the estimated probability of whether a triplet within an entire set of potential triplets is a false positive finding or not. After that, we report the ceRNA pairs with statistical significance of the observed S (e.g., adjusted p-value <0.05) as significant ceRNA interactions. For one miRNA, a low score S will be observed if the two target genes have no correlation within a specific range of miRNA expression. On the other hand, the score S will be high if a large proportion of the samples showed high correlation among the two target genes (P_obey(k) is high). Consequently, we can use the score S to identify a ceRNA triplet when S is high, and a statistical approach was performed to ensure that S cannot be identified randomly.

The SegmentClustering method

The motivation of the segment clustering is to group the samples showing high correlation between the expression levels of the two target genes into one single cluster. In our analysis, we divided the samples into small groups (i.e., a window) to calculate the correlation values of the two target genes. Accordingly, we may identify several different groups showing high correlations that actually can be clustered into one single cluster because of their similar correlation values. Therefore, in this method, the concept of a circular binary segment algorithm, which was originally designed for change-point problems such as the identification of copy number variation, is used [48] to evaluate whether those small windows showing high correlation values should be clustered into a larger group. This method has been widely used in the analysis of copy number variations (CNVs) by many algorithms [49, 50]. Since several previous studies [18, 27, 28] indicated that ceRNA triplets may be observed in a specific range, the purpose of such an approach is to explore the clustering patterns of samples in terms of their gene-gene correlation values per triplet and then group samples with similar correlation values so that a certain miRNA expression interval with the highest correlation within a ceRNA pair can be observed. This algorithm starts with all segments (i.e., several intervals of rank-ordered miRNA expression) identified from the whole dataset. Similar to the previous method, the samples are sorted by the miRNA expression values. A recursive test for the change-points is calculated based on the correlation values of each gene-miRNA pair between each set of two neighboring regions of rank-ordered miRNA expression, and it stops when no significant changes can be found in any two segments. The maximal t-statistic with a permutation reference distribution is chosen to obtain the corresponding p-value. Let X₁, …,X_n be the correlation coefficients of two genes, which are indexed by the corresponding miRNA expression of the n samples being studied. The test statistic is given by , where Z_ij is the two-sample t-statistic to compare the mean of the correlation of two genes with the index, which refers to the position within the rank-ordered miRNA expression from i+1 to j, to the mean of the correlation of the rest of the genes. That is, (5) where are the partial sums. If the p-values are smaller than a threshold level α (typically 0.01), a change is declared to be statistically significant. The locations of the change-points as the i and j that maximize the test statistic are also estimated by using either Monte Carlo simulations [51] or the approximation approach [52]. After performing the segment clustering approach, only a few groups of samples showing high correlation remain for further analyses. We have a basic assumption here for the ceRNA triplet: the correlation values should be stable across the sliding window approach and thus the correlation values should not be bouncing up and down within a small sample size. Following this assumption, we performed the peak merging step to ensure two nearby peaks were merged into one under certain criteria.

The PeakMerging method

This method is designed to prevent smashed segments resulting from noise. First, two segments are merged when the difference of their correlation values is smaller than a predefined value, and whether the finalized segments are peaks or troughs is defined compared to the baseline. Then, the Fisher transformation is performed to test the difference between two adjacent peaks by comparing their differences against the null hypothesis. Two adjacent peaks are further merged if there is no statistically significant difference between them. Because we presumed it is less likely that two genes are highly correlated to compete for a miRNA at more than two miRNA expression intervals, the triplet was abandoned when more than two peaks occurred in such a relationship. Lastly, candidate ceRNAs are finally selected when their correlation pattern contains a peak with a correlation value over a predefined threshold value (0.3, 0.5, or 0.7).

The final output of the ceRNAMethod function provides information on each miRNA, its candidate ceRNA pairs, its peak position within the rank-ordered miRNA expression interval, and the number of samples involved in such ceRNA interactions (S12 Fig).

Downstream functional analyses

To characterize the biological functions and pathways of identified ceRNA pairs, we implemented the ceRNAModule function to generate ceRNA modules from their interaction networks. For the ceRNA modules, the ceRNAFunction function is used to perform functional enrichment analysis based on two ontology databases, the Gene Ontology database (GO, http://geneontology.org/) [53] and the Kyoto Encyclopedia of Genes and Genomes Pathway Database (KEGG, https://www.genome.jp/kegg/) [54]. Survival analysis has been widely used in biomedical fields to indicate whether the ceRNAs in the discovered interaction module are associated with the survival of cancer patients. Hence, in ceRNAR, we implemented the ceRNASurvival function to perform survival analysis. First, the risk score of each sample is calculated using a multivariate Cox model. All the samples (i.e., different patients) are divided into high or low risk groups based on their median risk scores. The Kaplan-Meier method with a log-rank test is used to test and visualize the difference between the high and low risk groups. Additionally, the ceRNALocation function allows users to further visualize the expression level of a specific miRNA when modulated by a specific ceRNA. We also implemented an integrated function (ceRNAIntegrate) to combine our results with other state-of-the-art tools, such as SPONGE [20] and JAMI [23], and a validation function (ceRNAValidate) based on the miRSponge database. The graphical outputs of the above-mentioned functions are presented in S13–S15 Figs.

Simulation study

In the simulation study, two types of expression data (ceRNA and miRNA) of 100 samples were generated. In each triplet, we presumed 10–40% of samples have correlated genes whereas the rest of the samples do not. The simulated expression profile of 100 samples of miRNA was distributed uniformly between 0 and 1 and was used to sort our samples. We simulated the gene expression profiles of these target gene pairs under two situations. The first is a naïve profile of the correlated genes simply generated from a multivariate normal distribution with a mean value of 0 and a covariance matrix whose entries are 0.9, while the expression distribution of the uncorrelated genes was specified by setting its mean and variance to zero. Another simulation scenario is mimicking a real-data profile (9,835 pan-cancer samples) generated from a multivariate normal distribution with a randomly selected mean value (±2, ±1 and 0), and a covariance matrix whose entries are also randomly selected from 0.3 to 0.9 (correlation level: all), while the expression distribution of the uncorrelated genes was specified by setting its mean (also randomly selected) at ±2, ±1, and 0, and its variance (randomly selected) at 0 to 0.2 (S16 Fig). We also specified three correlation levels (high = 0.8–0.9, mid = 0.5–0.7, and low = 0.3–0.4) and randomly selected values within them to construct the covariance matrix. The normality test for sample distribution of each gene was conducted by the Anderson-Darling method [55].

Five scenarios were considered and are summarized in S1 Table. The first three scenarios were designed for a single peak (i.e., the highest correlation value of the target genes compared to the baseline) occurring at different locations (i.e., different miRNA expression intervals): lower miRNA expression (left, scenario 3), higher miRNA expression (right, scenario 2), and average miRNA expression (center, scenario 1), which represents specific correlation coefficient values existing at different expression levels of miRNA. The fourth scenario was designed for considering the occurrence of two peaks. A null scenario was also compared to test the validity of the other conditions. Two parameters were set under each scenario: window size (10%, 20%, 30%, and 40%) and the threshold at which to define a peak (0.3, 0.5, and 0.7). Each scenario was simulated 1,000 times. We conducted these simulations using either the complete or fast versions of the algorithm, the latter of which reduced computational complexity by omitting the random walk step. The proportion of correlated samples (10%, 20%, 30%, and 40%) was also analyzed.

Real data application for tools comparison and validation

In order to compare many aspects, including the computational cost and the quality of the results, of ceRNAR with the other tools (SPONGE [20], CERNIA [21] and JAMI [23]), which similarly used miRNA and paired target gene expression to identify gene-miRNA-gene triplets, we used subsets of the TCGA pan-cancer atlas with different sample sizes and gene/triplet numbers. We ran these tools with default parameter settings in parallel processing with 8 cores with varying sample sizes and numbers of genes.

Regarding real data application, we retrieved two sequencing datasets (TCGA-LUAD and TCGA-LUSC) from TCGA [56] to demonstrate the potential application of our proposed algorithm. These datasets are all retrieved from public domains (GDG data portal: https://portal.gdc.cancer.gov/). The TCGA-LUAD dataset contained 510 samples from lung adenocarcinoma patients, and the TCGA-LUSC dataset contained 476 samples from squamous cell carcinoma patients. To validate the results, the potential ceRNAs identified by our tool were validated experimentally using the miRsponge database [57], and the ceRNAs shared by these two datasets were further analyzed.