Multiple Linear Regression Analysis of lncRNA–Disease Association Prediction Based on Clinical Prognosis Data

Long noncoding RNAs (lncRNAs) have an important role in various life processes of the body, especially cancer. The analysis of disease prognosis is ignored in current prediction on lncRNA–disease associations. In this study, a multiple linear regression model was constructed for lncRNA–disease association prediction based on clinical prognosis data (MlrLDAcp), which integrated the cancer data of clinical prognosis and the expression quantity of lncRNA transcript. MlrLDAcp could realize not only cancer survival prediction but also lncRNA–disease association prediction. Ultimately, 60 lncRNAs most closely related to prostate cancer survival were selected from 481 alternative lncRNAs. Then, the multiple linear regression relationship between the prognosis survival of 176 patients with prostate cancer and 60 lncRNAs was also given. Compared with previous studies, MlrLDAcp had a predominant survival predictive ability and could effectively predict lncRNA–disease associations. MlrLDAcp had an area under the curve (AUC) value of 0.875 for survival prediction and an AUC value of 0.872 for lncRNA–disease association prediction. It could be an effective biological method for biomedical research.


Introduction
Long noncoding RNAs (lncRNAs) are noncoding RNA molecules, including miRNAs [1], lncRNAs [2], tRNAs [3], piRNAs [4], and more than 200 nucleotides. They were initially thought to be nonfunctional RNA fragments and the only by-product of massive transcription [5][6][7][8]. A large number of recent studies have shown that lncRNAs have abundant biological functions, including the silencing of X chromosome [9] and activation and interference of transcription [10][11][12]. At the same time, the abnormal expression of lncRNAs leads to various diseases [13][14][15]. Therefore the investigation on lncRNA-disease associations is of great significance at the molecular level to radically cure the disease.
Many computational methods have been applied to human lncRNA-disease association prediction in recent years. These methods have two prominent features: machinelearning-based feature and network-based feature.
The machine-learning-based feature of lncRNA-disease association prediction is to establish a learning model in the training dataset and then to perform tests in the test dataset using this learning model. For instance, Zhao et al. [16] developed a learning model based on the Bayesian classifier for lncRNA-disease association prediction. The key issue was that the learning model regarded unknown lncRNA-disease associations as negative sets, restricting the performance of the learning model. In fact, the negative sets for lncRNA-disease association prediction were difficult to achieve. To avoid this problem, Chen et al. [17] put forward the method of LRLSLDA to predict lncRNA-disease associations. It was based only on positive sets, not on negative sets. It adopted the strategy of Laplacian regularized least squares, was a semisupervised learning model, and needed selected optimal parameters to obtain optimal prediction results. LRLSLDA had two limitations: (a) they were under the assumption that functionally similar lncRNAs were related 2 BioMed Research International to similar diseases and (b) they were restricted to selecting optimal parameters. Moreover, Chen et al. [18] developed another method, KATZLDA, which incorporated known lncRNA-disease associations, lncRNA expression profiles, lncRNA functional similarity, disease semantic similarity, and Gaussian interaction profile kernel similarity. KATZLDA could adapt to new diseases and lncRNAs without any known associations. However, KATZLDA still was built on lncRNAs with more known associated diseases or/and miRNA interaction partners.
The network-based feature of lncRNAs-disease association prediction is to establish a learning network of lncRNA-disease associations using known associations. For instance, Yang et al. [19] constructed the bipartite network about lncRNA-disease associations and predicted these associations by the method of transmission in the network, which was the first prediction method based on the network model. A coding-noncoding gene-disease bipartite network was constructed to improve the prediction results in which better prediction results were obtained. However, the approach did not take into account the interaction between lncRNAs and coding genes, and the forecast range was limited. The results were relatively general. Sun et al. [20] proposed a new method based on the network model: RWRlncD. It was used to construct the functional similarity network of lncR-NAs using the known lncRNA-disease association network and the similarity network. Subsequently, the reactivated random walk was performed in the functional similarity network of lncRNAs to predict the potential lncRNA-disease associations. However, the edge of the test set was used to calculate the functional similarity of lncRNAs before cross validation, which overestimated the verification results. Zhou et al. [21] presented a novel method (named RWRHLD) to distinguish candidate lncRNA-disease associations using the hybrid network and then performed the random walk algorithm on this hybrid network. RWRHLD was used only to predict lncRNAs in known lncRNA-miRNA associations, where an incomplete coverage of lncRNAs cross talk network and lncRNA-disease association network might lead to inaccurate prediction results. Chen et al. [22] improved the traditional random walk with restart and proposed the method of improved random walk with restart for lncRNA-disease association prediction (IRWRLDA). But the existing problem of IRWRLDA was how to obtain integrated lncRNA similarity based on lncRNA functional similarity and lncRNA Gaussian interaction profile kernel similarity. Chen et al. [23] developed two novel LNCSIMs and proposed a new method LRLSLDA-LNCSIM that could improve the predictive ability of LRLSLDA. LRLSLDA-LNCSIM still had the limit that a semantic contribution decay factor was not well solved. Yu et al. [24] employed multidimensional heterogeneous data to construct an lncRNA function similarity network, employed the disease ontology to construct a disease network, and then proposed the BRWLDA to predict lncRNA-disease associations. Although the prediction performance was improved by BRWLDA, the defect of random walk algorithm still existed. Chen et al. [25] developed a method of HGLDA by integrating miRNA-disease associations and lncRNA-miRNA interactions. However, HGLDA could not be used in the lncRNAs without any known miRNA interaction partners. Ganegoda et al. [26] developed the computational model of KRWRH network, which was a heterogeneous network formed by integrating a disease-disease similarity network, lincRNA-lincRNA similarity network, and known lincRNA-disease association network.
The reviews of Chen et al. [27] and the aforementioned discussions showed that few references were made to the combination of clinical prognosis data with lncRNA-disease association in the existing studies on lncRNA-disease association prediction. In the present study, the analysis of disease prognosis was ignored, and the existing prediction model was limited to a single lncRNA prediction. The prognosis information of a disease associated with lncRNAs was rarely involved (such as the survival time of patients, current state of disease, and family history of genetic diseases). In fact, an analysis related to the prognosis of lncRNA-disease association has more realistic meaning and value.
To overcome the aforementioned issues, a multiple linear regression model for lncRNA-disease association prediction based on clinical prognosis data (MlrLDAcp) was constructed to predict the potential associations between lncRNAs and diseases. At the same time, the survival time of patients with prostate cancer was also predicted in MlrL-DAcp. The concepts of predictive correlation factor Θ, decay coefficient , Γ operation, and Γ correction were proposed in this study to construct the multiple linear regression model. An algorithm for developing the multiple linear regression model was designed, in which 481 lncRNA transcripts with values less than 0.001 were cut back to 60 most closely related to the survival time of patients with prostate cancer. Finally, the potential multiple linear regression relationship between the prognosis survival time of 176 patients with prostate cancer and 60 lncRNAs was proposed. MlrLDAcp could realize not only cancer survival prediction but also lncRNA-disease association prediction. Compared with previous findings, MlrLDAcp had a predominant survival predictive ability and could effectively predict lncRNA-disease associations.

LncRNA Expression Data.
The lncRNA expression data of prostate cancer was obtained from the lncRNAtor database (http://lncrnator.ewha.ac.kr). A total of 44 normal samples and 176 prostate cancer samples were obtained, and the prostate cancer samples were denoted in an ascending order according to sample ID (denoted by 1 , 2 , ⋅ ⋅ ⋅ , 176 ). The expression level of each lncRNA transcript in normal and prostate cancer samples was calculated. The differential expressions between normal and prostate cancer samples were calculated using the aforementioned expression quantities. A total of 481 lncRNA transcripts with a significant difference were obtained (denoted by 1 , 2 , ⋅ ⋅ ⋅ , 481 using a P value in the ascending order) by selecting a P value less than 0.001, and the 481 transcript expression quantities on were denoted by ∼1 , ∼2 , ⋅ ⋅ ⋅ , ∼481 (1 ≤ ≤ 176). The lncRNA expression training data matrix (denoted by ) was constructed based on the aforementioned data.

Abstracting the Issue.
This study on lncRNA-disease associations was conducted from the following two aspects: (a) A part of lncRNAs in the 481 lncRNAs, which were most closely related to prostate cancer, was screened out through the analysis of prognosis survival. Hence, a subset of (denoted by , sub ⊂ ) was obtained, and set Δ = ⟨ ⟩ (⟨ ⟩ was the number of elements in , 1 ≤ Δ < 481). contained Δ lncRNAs (denoted by 1 , 2 , ⋅ ⋅ ⋅ , Δ ). The expression quantity of on was denoted by ∼1 , ∼2 , ⋅ ⋅ ⋅ , ∼Δ . Therefore, = ].
(b) The potential relationship between and was given using multiple statistical methods, and finally the prognosis prediction of lncRNA-disease associations was realized using predicting .
Definition 1 (predictive correlation factor Θ). Θ was defined , where Θ corresponded to in . The value of predictive correlation factor Θ was the coefficient of multiple linear regression . is shown in (1). For the prognosis prediction of lncRNA-disease associations, the formal definition was as follows. Two tasks needed to be completed while establishing : (1) to calculate and (2) to calculate Θ. Definition 3 (decay coefficient ). denoted the decay duration in Γ operation. ( ), which was the value of the ℎ iteration in Γ operation, is shown in (2).
In (2), is the maximal iterations of Γ operation and is the current iterations of Γ operation. When takes the minimal value 1, is close to 1. When takes the maximal value , is 0.1. When takes the value between 1 and , is the reduction value between 1 and 0.1. The greater is, the greater the decay efficiency is.
Definition 4 (Γ operation). Γ operation was divided into three stages. (a) The variation center ( ) of the ℎ iteration on was calculated according to (3). (b) The variation range ( ) of the ℎ iteration on was calculated according to (4). (c) Bitwise inversion was performed within the variation range ( ( ) ∩ [1, ( )]). Figure 1 shows the schematic of Γ operation.
In (3) and (4), is the maximal iterations of Γ operation and is the current iterations of Γ operation.
Definition 5 (Γ correction). The components (their value was 1) in remained valid only within the top period in Γ operation, and the rest were set to 0.

Stepwise.
( ) was performed on . Step 2. Whether a new component was added to { } was judged; if true, then both + and − were set to 0 and Step 1 was performed; otherwise, + was set to 1 and Step 3 was performed.
Step 4. A component of insignificant effect on was removed from { }.
Step Stepwise (X d ) AIC( −X d + ) was better than AIC(X optimal ) X optimal = -X d + Led sub = X optimal on the previous  components of X d set to 0 and Step 4 was performed; otherwise, − was set to 1 and Step 1 was performed.   current best bulletin board. 0 , the value of which was less than 10 −12 in , was set to 1, and the rest was set to 0.
Step 2. Γ operation was executed. If top components of were the same as top components of −1 and > 0, bitwise inversion was performed on the previous components of . After that, Γ correction was performed.
Step 4. Whether the maximal iterations were reached was judged; if true, then Step 5 was performed; otherwise, Step 2 was performed.
Step 5. was set to .
Step 6. The multivariate linear regression analysis on was performed using . Θ was set to the multiple linear regression coefficients.
was the output.

3.1.
( ). The calculation results of variation center ( ) in Γ operation are shown in Figure 3. The figure shows that the variation center value of Γ operation was evenly distributed in 10 iterations. Furthermore, the variation center value was scattered around each interval from 1 to 481. The variation center value had 5 points on each side of the center point (240). This ensured that the components of all 481 LncRNAs had an equal opportunity to perform variations. It was more beneficial to obtain the global optimal solution. Figure 4 shows the calculation results of decay coefficient . The variation operation proposed in this study aimed to enrich the diversity of sample space. Meanwhile, the variation operation should be dynamic rather than fixed. For the aforementioned issue, the decay coefficient was proposed to control the strength of the variation operation. Figure 4 shows that the decay coefficient decreased with the increase in the number of iterations. This was because the variation operation should be strong at the incipient iteration to obtain the global optimization ability. On the contrary, the variation operation should be weak at the late iteration to obtain the local development ability. It was not hard to see that the decay coefficient was the value of decreasing change between 1 and 0.1, and it controlled the lncRNA components that performed the bitwise inversion in 481 lncRNA components. Table 1 shows the detailed calculation process of , including the variation position, interval, and result in the ℎ iteration. As shown in Table 1, the variation position and the interval distribution were relatively discrete, reducing the blind area of Γ operation. In addition, the two factors proposed in the calculation to ensure optimal performance were significant differences in (denoted by ) and of ( ) (denoted by ( ) ), both of which were characterized by the better performance with a smaller value. Two factors should not be considered unilaterally but comprehensively. For example, if only was taken into account, of 0 + , which was the set of less than 10 −12 , was equal to 2249.24, as shown in Figure 5. Obviously, it was not an optimal solution but a poorer solution. Based on the aforementioned considerations, Γ operation, which was proposed in this study, was a combination of and ( ) . In each variation process, the smaller part of performed the variation, rather than the whole. Finally, 1 + ( =2208.47) was the optimal solution. Therefore, was equal to 1 + (i.e., the matrix built using the expression quantity of 1 + on , 1 ≤ ≤ 176 ). Then, a multivariate linear regression analysis was performed on , and the regression coefficients were obtained as required ( Table 2). The ensemble transcript ID of X14 was ENST00000559477 in the Table 2. The ensemble transcript ID of X16 was ENST00000560882 in the Table 2.  3, 9, 10, 11, 13, 14, 16, 18, 19, 20, 22, 23, 25, 26, 27, 29, 30, 33, 34, 35, 36, 108, 110, 112, 113, 114,  116, 118, 120, 122, 125, 126, 127, 130, 131, 134, 135,  137, 138, 139, 140, 145, 146, 149, 150, 151, 152, 153,  154, 157, 163, 164, 165, 169, 171, 172 190,192,193,194,198,199,200,201,202,203,204,205,208 190,193,194,197,198,200,203,206,208 {7, 9,11,12,14,15,16,18,20,23,24,108,110,111,112,114,116,120,121,122,123,124,125,184,185,187,190,193,194,196,197,199,201,203 The test was performed on the regression model (the results are shown in Table 3). As was less than 0.0001, the regression model had the statistical significance. It also indicated that MlrLDAcp was feasible and effective.

Computing .
The prediction model (MlrLDAcp) proposed in this study had two potential aspects: (a) The survival of cancer patients was predicted by combining with the multiple linear regression model of MlrLDAcp.
(b) The association between lncRNAs and diseases was predicted using MlrLDAcp.
The performance of evaluation was expanded from the two aforementioned aspects.

Survival Predictive Ability.
Receiver operating characteristic (ROC) analyses were performed to compare the predictive accuracies of prostate cancer samples between MlrLDAcp and Huang's method [28] (the state-of-the-art method), to evaluate the survival predictive ability. The 5-year biochemical recurrence survivals of the two methods were compared between TCGA and lncRNAtor databases. Figure 6 shows the experimental results. The value of the area under 8 BioMed Research International  the curve (AUC) was calculated from the corresponding area under the ROC curve. As shown in Figure 6, MlrLDAcp with an AUC value of 0.875 was better than Huang with an AUC value of 0.833. As a result, the prediction accuracy  associations in { } were defined as -(1 ≤ ≤ , and was the number of known lncRNA-disease associations). In each step of the LOOCV, eachwas implemented on { − -} and { + -}, and then the model learning was carried out on { − -}. The ROC curve plotted the sensitivity (that was true-positive rate = /( + )) versus the 1-specificity (that was false-positive rate = /( + )), where TP denoted true positives, FP denoted false positives, TN denoted true negatives, and FN denoted false negatives. The sensitivity was the ratio of positive samples which could be accurately distinguished, and the specificity represented the percentage of negative samples which could be correctly predicted. Figure 7 shows the experimental results. The value of AUC was calculated from the corresponding area under the ROC curve. As shown in Figure 7, MlrLDAcp with an AUC value of 0.872 was better than KRWRH with an AUC value of 0.838 and LRLSLDA with an AUC value of 0.822. As a result, the prediction accuracy of lncRNA-disease associations in MlrLDAcp increased by 3.4% (versus KRWRH) and 5.0% (versus LRLSLDA). These results suggested that MlrLDAcp might have a preferable ability to predict lncRNA-disease associations.

Conclusions
In this study, a model of MlrLDAcp was constructed. MlrL-DAcp took the expression quantity of lncRNAs transcript as an independent variable and the clinical prognosis data as a dependent variable. Using MlrLDAcp, 60 lncRNAs, which were most closely related to cancer prognosis information (survival time), were selected from 481 alternative lncR-NAs. MlrLDAcp could realize not only the cancer survival prediction but also the lncRNA-disease association prediction.
Further research directions about lncRNA-disease association prediction are as follows.
(a) The lncRNA-disease association prediction should take into account clinical prognostic data in future investigations. The lncRNAs associated with diseases may have a clinical value as therapeutic targets. Hence, the clinical prognostic data is quite valuable to lncRNA-disease association prediction. The clinical implications and the mechanism underlying the association of lncRNAs with diseases are definitely worth exploring further.
(b) How to build an effective computational model to construct an lncRNA similarity function, which can reasonably integrate the similarity scores of different biological information, is worthy of further research.
(c) With the increase in lncRNA-disease correlation, the prediction accuracy can be further improved. Furthermore, most computing models rely heavily on unobtainable negative samples, which is an urgent problem to be solved.
(d) The new network-based computing model should be implemented on heterogeneous networks instead of single networks. Hence, more heterogeneous networks, such as lncRNA-disease network, disease similarity network, lncRNA functional similarity network, and lncRNA interactive networks, should be integrated in the future.

Data Availability
The data used to support the findings of this study are included within the article.