We downloaded human and mouse lncRNA-protein interactions from the NPInter database v4.0 (Hao et al., 2016). Self-interactions and duplicates were removed. The mouse lncRNA-protein interaction network involves 33255 lncRNAs, 182 proteins, and 102051 interactions. The human lncRNA-protein interaction network contains 41589 lncRNAs, 3237 proteins, and 394895 interactions. We downloaded human and mouse protein-protein interaction data from BioGrid database version 4.4 (Oughtred et al., 2021). The mouse protein-protein interaction network includes 9744 proteins and 52342 interactions. The human protein-protein interaction network includes 19106 proteins and 644235 interactions.

We combine the lncRNA-protein interactions and protein-protein interactions by matching the name of the proteins in both datasets, producing a heterogeneous network with two types of interactions. The mouse network was composed by 9845 proteins and 33255 lncRNAs with 102051 lncRNA-protein interactions and 52342 protein-protein interactions. The human network was composed by 19553 proteins and 41589 lncRNAs with 394895 lncRNA-protein interactions and 644235 protein-protein interactions. The sequences of all lncRNAs in both human and mouse interaction networks were obtained from the NONCODE database version 5 (Zhao et al., 2016).

According to literatures (Penny et al., 1996; Marahrens et al., 1997; Lee, 2000; Sado et al., 2001; Grote et al., 2013; Klattenhoff et al., 2013; Sauvageau et al., 2013; Yildirim et al., 2013; Zeng et al., 2018), eight mouse lncRNAs, including Xist, Gas5, Meg3, Tsix, Gt (ROSA) 26Sor, Dnm3os, Fendrr, and Braveheart, were identified as essential lncRNAs. The remaining 33247 lncRNAs in the mouse network were marked with unknown status. For human lncRNAs, we curated a set of lncRNAs that are reported to be essential in various conditions from literatures (Supplementary Table S1). This set contains 63 lncRNAs. The names of these lncRNAs and the conditions that they are reported to be essential, are listed in Supplementary Table S1, along with literatures of the original reports. In addition, 11 mouse lncRNAs, which are homologous of human essential lncRNAs, were also collected for validation purpose, as homologous usually have similar essentiality (Georgi et al., 2013).

Gene Importance Calculator (GIC) (Zeng et al., 2018) is a useful essentiality indicator for both protein-coding genes and noncoding genes. It is based solely on sequence information. The GIC score (g) is defined as follows: g = 1 1 + exp [ - θ ( p ) ] , (1) where θ(p) is derived from a logistic regression model. θ(p) can be defined as θ ( p ) = ln p 1 − p = β 0 + β 1 L + β 2 1 L e + ∑ i = 1 5 α i f i , (2) where α ₁, α ₂, …, α ₅, β ₀, β ₁ and β ₂ are regression coefficients, L the length of RNA sequence, e the minimum free energy of RNA secondary structure, p the conditional probability that a gene is essential, and f _i the occurrence frequency of a triplet in the sequence. The five types of triplets, which are considered in the GIC, are CGA, GCG, TCG, ACG and TCA (Zeng et al., 2018).

When calculating the GIC score, we need to use the external program RNAfold (Lorenz et al., 2011), which requires a sequence length less than 20000 nt. Therefore, only 24450 mouse lncRNAs and 29481 human lncRNAs can be calculated for GIC. All other lncRNAs have lengths too long for the RNAfold to work.

We formulate the heterogeneous graph as G = ( V , E ), where V is the set of all nodes, including lncRNAs and proteins, and E the set of all interactions, including lncRNA-protein and protein-protein interactions. Without losing generality, we note the number of all nodes as n. The network can be represented as an adjacency matrix A∈{0.1} ^n×n . The element on the ith row and the jth column of A can be denoted as a _i,j . If a _i,j = 1, the ith node and the jth node have interactions between them. If a _i,j = 0, there is no interaction between the ith node and the jth node. Given a _i,j , we can define four different centrality measures, including degree centrality (DC), betweenness centrality (BC), closeness centrality (CC), and eigenvector centrality (EC) for each node in the network.

The degree centrality of the ith node can be defined as follows: d i = 1 n − 1 ∑ j = 1 n a i , j (3)

The betweenness centrality of the ith node can be defined as follows: b i = 1 ( n − 1 ) ( n − 2 ) ∑ u ≠ i ≠ v ∈ V σ u , v ( i ) σ u , v (4) where σ _u,v is the number of shortest paths between the uth node and the vth node, and σ _u,v (i) the number of shortest paths between the uth node and the vth node that pass the ith node.

The closeness centrality of the ith node is defined as follows: c i = [ | R ( i ) | − 1 ] 2 ( n - 1 ) ∑ j ∈ R ( i ) d i , j (5) where R (i) is the set of nodes that can reach the ith node, d _i,j the length of the shortest path between the ith node and the jth node, and |.| cardinal operator of a set.

The eigenvector centrality of the ith node is defined as follows: e i = x max ( i ) (6) where x _max(i) is the ith dimension of the normalized eigen vector x that corresponds to the largest eigen value of adjacency matrix A. Let λ _max be the largest eigen value of A, the following relationships are satisfied in finding x: Ax = λ max x ,   and (7) ‖ x ‖ = 1 (8) where ||.|| is the vector norm operator.

Our network model contains two types of nodes, lncRNAs, and proteins. It also involves two types of interactions, the lncRNA-protein interactions and protein-protein interactions. Essentially, it is a lncRNA-protein-protein interaction (LPPI) heterogeneous network. Figure 1 illustrates a part of the LPPI network for human and mouse respectively.

We propose the Systematic Gene Importance Index (SGII) as a comprehensive measure of gene essentiality, particularly for non-coding genes. SGII is a combination of the sequence-based GIC score and centrality measures, which have been elaborated as above.

For the ith node in the LPPI network, we compute its BC, CC, DC and EC, which can be noted as b _i , c _i , d _i and e _i , respectively. Its GIC score is noted as g _i . We sort all nodes according to their BC, CC, DC, EC and GIC in a descending order, respectively. The rank of the ith node after sorting according to BC, CC, DC, EC and GIC can be noted as r _b (i), r _c (i), r _d (i), r _e (i) and r _g (i), respectively.

Let s _i be the degree of the ith node, which can be computed as follows: s i = ∑ j = 1 n a i , j (9)

Given a threshold z, if s _i ≥ z, the centrality measures will determine the essentiality of a gene directly. For convenience, we define the centrality-based essentiality indicator function for the ith node according to BC, CC, DC, and EC respectively as follows: I b ( i ) = { 1 r b ( i ) n < k % , 0 o t h e r w i s e (10) I c ( i ) = { 1 r c ( i ) n < k % , 0 o t h e r w i s e (11) I d ( i ) = { 1 r d ( i ) n < k % , 0 o t h e r w i s e , and (12) I e ( i ) = { 1 r e ( i ) n < k % , 0 o t h e r w i s e (13) where k is a rank threshold parameter. The ith node is identified as essential when I b ( i ) I c ( i ) I d ( i ) I e ( i ) = 1 (14) is satisfied.

If s _i < z, we rely on the GIC score to determine the essentiality of a gene. Similarly, we can define the indicator function for GIC ranking, as follows: I g ( i ) = { 1 r g ( i ) n < t % , 0 o t h e r w i s e (15) where t is another rank threshold parameter. The ith node is essential if I g ( i ) = 1 (16) is satisfied.

The whole flowchart of SGII is illustrated in Figure 2.

In evaluating SGII, we use three statistics to describe its predictive performance. These statistics include sensitivity (s), false positive rate (r), and Fisher’s exact test score (f), which are defined as follows: s = n t n + (17) r = n f n − , and (18) f = − log 10 p (19) where n _t is the number of known essential lncRNAs that are identified as essential lncRNAs, n ₊ the total number known essential lncRNAs, n _f the number of lncRNAs with unknown essentiality that are identified as essential, n _- the total number of lncRNAs with unknown essentiality and p the p-value of Fisher’s exact test. Since SGII is a direct scoring method with manually configurable cutoff values, no training procedure is involved in the whole process. This is different to machine learning based methods. We cannot treat the above sensitivity and false positive rate as comparable to those in evaluating machine learning methods, as the knowledge of essential lncRNAs is too limited to perform any kind of cross-validations. This is also why we introduced the Fisher’s exact test to further quantifying the quality of our results. It will measure how likely a result in whole is random or not. The bigger f value is, the results are less likely to be random.

There are eight parameters in the GIC, which represent all the coefficients in the model built by GIC method. We took all the parameter values from literature (Zeng et al., 2018). The values for the mouse model are β ₀ = 0.1625, β ₁ = 2.638 × 10^–4, β ₂ = 2.194, α ₁ = 19.88 (for CGA), α ₂ = 37.59 (for GCG), α ₃ = 50.37 (for TCG), α ₄ = 35.44 (for ACG), and α ₅ = -64.66 (for TCA). The values for human model are β ₀ = 0.7417, β ₁ = 2.612 × 10^–4, β ₂ = 4.295, α ₁ = 48.66, α ₂ = 15.64, α ₃ = 76.23, α ₄ = -1.113, and α ₅ = -60.29.

Three parameters are introduced in combining centralities and GIC, which are noted as z, k and t. We first perform a grid search of k and t with a given value of z. The pairs of k and t, which maximize the score f, are recorded for every different z. These values are further sorted to find the best z, k and t combination. When performing the grid search on the mouse dataset, k = 1, 3, 5, 7, 9, and t = 1, 3, 5, 7, 9. When performing the grid search on the human dataset, k = 5, 10, 15, 20, 25 and t = 5, 10, 15, 20, 25. For both datasets, z = 5, 10, 15, 20. Finally, we set z = 15, k = 5, t = 9 for mouse dataset, and z = 5, k = 20, t = 5 for human dataset. All results for different parameters are provided in supplementary materials, as Supplementary Table S2.

We first explore the basic statistical characters of the LPPI network. We plot the degree distribution of the mouse and human network respectively in Figure 3. It is intuitively that the distribution of the degree follows the common power law distribution, which is similar to the PPI networks (Jeong et al., 2001). Since in the PPI network, essential proteins are usually rare and with high degrees, we assume that in our LPPI network, the essential lncRNAs have similar properties.

As we have mentioned in the method section, several lncRNAs with a length too long to calculate its secondary structure were not counted in our analysis. It becomes a question whether these lncRNAs have preferences to large or small amounts of interactions. We plot the degree distribution with and without those over-length lncRNAs for mouse and human datasets, respectively, in Figure 4. It is hard to find differences on the degree distributions. We therefore believe that, for a lncRNA, its length alone is not a major contributing factor to its interactions in the LPPI network. This also implied that the essentiality, which we believe to be associated with local network structure, has no direct relationship with the length of the lncRNA. These over-length lncRNAs were kept in the network as dummy nodes, which means we did not compute their essentiality at all, regardless of whether they have a degree over the threshold or not.

Figure 5 gives scatter plots of GIC pairing with each of the four types of centralities on human and mouse datasets, respectively. For the mouse dataset, the red dots, which represent essential lncRNAs, tend to appear in the top-right part of the plots, while the blue dots, which denote all other lncRNAs, spread much wider. Although the red dots are relatively rare, but their top-right preference is still observable. For human dataset, this preference is not intuitively obvious.

This allows us to carry out further quantitative analysis on combining the centrality measures and the GIC scores. A primary challenge is that the number of known essential lncRNAs is too small for a machine learning algorithm to train on. In addition, some essential lncRNAs are only involved in a very limited number of interactions. For example, the Braveheart (Bvht) lncRNA, which is essential, has only one interaction record in the database. We think this may be due to the incomprehensive knowledge of the lncRNA-protein interaction network. As the estimation of centrality measures highly rely on the interaction enrichment of a node in the network, when dealing with a lncRNA with limited number of interactions, we turn to rely on the GIC score.

With the settings in the method section, we combined four types of centrality measures and the GIC scores. On the mouse dataset, we identified 2284 essential lncRNAs from altogether 24450 lncRNAs. Among the 2284 lncRNAs, eight lncRNAs are known to be essential, accounting for 100% of all known essential lncRNAs, resulting a p-value = 5.73 × 10^–9 (Fisher’s exact test). On the human dataset, we identified 5063 essential lncRNAs, from altogether 29481 lncRNAs, Among the 5063 essential lncRNAs, 41 lncRNAs are reported to be essential in various conditions in literatures, accounting for 65% of all curated essential lncRNAs (p-value = 3.59 × 10^–17, Fisher’s exact test). This result clearly indicates that our method is effective to identify essential lncRNAs.

As SGII is the first attempt to combine the network information and sequence information to identify essential lncRNAs, we explore which kind of centrality measure is more capable to identify essential lncRNAs along with the GIC scores. We first plot the distribution density of different centralities and the GIC scores on mouse and human datasets respectively. As in Figure 6, BC and DC centrality measures along with the GIC scores appear to have much better separation than the CC and EC measures on the mouse dataset, while on the human dataset, only BC and DC present an intuitive separation.

However, considering the large differences on axis scale for essential lncRNAs and all lncRNAs, these intuitive observations may be misleading. Therefore, we performed a quantitative comparison using eight different conditions, GIC alone, GIC combined with each one of four types of centralities, GIC combined with BC and DC, GIC combined with CC and EC, and GIC combined with all four types of centralities. The parameters of all comparison are optimized as in method section (Table 1).

The first observation on Table 1 is that the best combination of centrality measure and the GIC is not the combination of all four types of centralities. For the mouse dataset, the BC + DC + GIC method has the best significance level and lowest FPR value. For the human dataset, the BC + GIC method reaches the highest significance level. A second to the best significance level is obtained again by BC + DC + GIC method, with the highest sensitivity value. Therefore, we think that the BC + DC + GIC may be a better way to identify essential lncRNAs than the current configuration of SGII. This consists with the impression from Figure 6. However, due to the limited number of available data and current results, it is possible that this observation does not reflect a comprehensive scene of identifying essential lncRNAs. Therefore, we keep the configuration of SGII to combine all four kinds of centralities and the GIC score, for an unbiased way of identifying essential lncRNAs.

At a closer look to Table 1, it appears that the Fisher’s exact test reports much more significant results on both datasets when GIC is combined with centralities, which proves that integration of centrality measures and GIC is effective. Another observation is that SGII gives under-expected sensitivity values on the human dataset. However, the significance levels on the human dataset are generally way higher than that of the mouse dataset. This may be the results of two differences between the mouse and the human datasets. First, the human dataset is collected from literatures of lncRNAs in various conditions, including tumor cell line experiments. Essential lncRNAs, which are identified by one type of cell line experiments, may be different to those from the original essential gene definitions. As direct essential gene experiments on human are not feasible, the quality of the dataset is not comparable to the mouse dataset. This also applies to the coding gene data (Austin et al., 2004). Secondly, the number of essential lncRNAs in the human dataset is roughly eight times of that of mouse dataset. Since the computation process of the significance level is affected by the raw counts, it is anticipated that systematic differences on significance levels exist.

To further confirm the above explanations, we performed the following analysis. We find homologous genes of human essential lncRNAs in mouse. According to the studies in coding genes, these genes are likely to also produce essential lncRNAs(Georgi et al., 2013). Altogether 11 homologous genes in mouse were identified as lncRNA genes in the mouse LPPI network. We used SGII to test if we can identify these homolog essential lncRNAs (Table 2).

Obviously, sensitivity is dropping in comparison to the mouse essential lncRNAs. However, it should be noted that the FPR is also dropping, which indicates much less false positives. The significance levels remain almost the same as the mouse essential lncRNAs. Again, the BC + GIC method obtained the best significance level, while the BC + DC + GIC method obtained a second to the best significance level with the highest sensitivity. This result confirmed that the significance level difference between human and mouse dataset is largely caused by the raw counts of the dataset. It also suggests that the BC + GIC or BC + DC + GIC method may be a better choice than combining all types of centralities and the GIC score.

The importance of BC can be understood intuitively. If we think the cellular system as a system composed of molecules. The interactions between molecules transfer information. A high BC value indicated that the node is critical as an information hub in many shortest paths between other nodes. Therefore, dropping such nodes will easily break many information channels simultaneously, which will eventually destroy the whole system. That makes it an essential node in the network.

For the DC measure, the intrinsic mechanism is similar. The DC measure is directly associated to the degree of a node. If a node with many edges is dropped, it is more likely that the whole network collapses. This consists with the observations in coding genes. In addition, although some other kinds of centralities, like the NC (new centrality) (Wang et al., 2012), can identify essential coding genes better, it does not work well in non-coding genes. This is an expected result. For NC to work in the LPPI network, it requires that dense interactions exist among the proteins that interacting the same lncRNAs. However, we did not observe this phenomenon in our dataset. The NC is difficult to be estimated for many lncRNAs, due to lacking such kind of interactions.

We took the essential lncRNA gene in mouse genome for functional analysis. For every lncRNA that was predicted as essential in mouse genome, we first map this lncRNA to the Ensembl database (Howe et al., 2021) using either gene name or sequence information. The mapped genes are then uploaded to the Gene Ontology online system for functional enrichment analysis. The top three enrichment of functions are “nucleic acid binding” (GO:0003676), “heterocyclic compound binding” (GO:1901363) and “organic cyclic compound binding” (GO:0097159). As we have mentioned, this is expected for lncRNAs. They realize their functions through bindings with other molecules.