# Graph-Based Analysis of RNA Secondary Structure Similarity Comparison

Authors: Lina Yang; Yang Liu; Xiaochun Hu; Patrick Wang; Xichun Li; Jun Wu
Journal: Complexity
 (2021)

---

## Abstract

In organisms, ribonucleic acid (RNA) plays an essential role. Its function is being discovered more and more. Due to the conserved nature of RNA sequences, its function mainly depends on the RNA secondary structure. The discovery of an approximate relationship between two RNA secondary structures helps to understand their functional relationship better. It is an important and urgent task to explore structural similarities from the graphical representation of RNA secondary structures. In this paper, a novel graphical analysis method based on the triple vector curve representation of RNA secondary structures is proposed. A combinational method involving a discrete wavelet transform (DWT) and fractal dimension with sliding window is introduced to analyze and compare the graphs derived from feature extraction; after that, the distance matrix is generated. Then, the distance matrix is analyzed by clustering and visualized as a clustering tree. RNA virus and noncoding RNA datasets are applied to perform experiments and analyze the clustering tree. The results show that the proposed method yields more accurate results in the comparison of RNA secondary structures.


---

## Body

## 1. Introduction

The secondary structure of RNA is a double-stranded structure constructed according to the principle of complementary base pairing. In RNA, U (uracil) is complementary to A (adenine) andG (guanine) is complementary to C (cytosine). RNA is usually transcribed from DNA and acts as a bridge between DNA and proteins [1]. The versatility of an RNA molecule depends on its secondary structure. RNAs with similar structures tend to have similar functions or properties, but the opposite does not necessarily hold true. Many RNA molecules are conserved at the structural level, but they have little sequence similarity. Therefore, the comparison of RNA secondary structures is key to elucidating their functional and evolutionary relationships.Most recent studies have focused on RNA secondary structure prediction [2–4], and comparisons of RNA secondary structures have not yet been sufficiently studied. At the present stage, the comparison methods for RNA secondary structures are mainly divided into two types: alignment-based methods and alignment-free methods.Alignment-based methods mainly rely on an RNA secondary structure represented by a string or tree [5–9].The Sankoff algorithm used a free-energy minimization method to synchronize the folding and alignment of two or more RNA sequences [9]. Combining the information from folding and alignment optimizes the objective function when processing the problem separately, while generalizing the alignment method to the original sequence reconstruction problem and selecting appropriate, existing constraints for possible solutions. However, this suffers from a problem. The time cost is high and about O (n3kN) (N is the number of sequences, n is the maximum sequence length, and the value of k is small enough). In practice, it is impossible to solve and satisfy both the time complexity problem and the larger memory requirements, which, together with the impractical constraints, are all difficult to achieve.Therefore, to be more effective in RNA secondary structure alignment, a number of improved algorithms were subsequently developed, such as Consan [10] and, Dynalign [11], PMcomp [12], Stemloc [13], Foldalign [14], locARNA [15], SPARSE [16], MARNA [17], FoldAlignM [18], Murlet [19], CARNA [20], and RAF [21].Consan is based on Sankoff algorithm and was proposed by Dowell to address two problems not solved by Sankoff: the fact that it cannot be generalized to RNA structures of different lengths and the need to deal with the complexity of the algorithm. First, the method of modelling secondary structures and alignments using random context-independent grammars (SCFG) [22] is extended to construct pair-SCFG descriptions and build unconstrained algorithms. Second, the complete pair-SCFG algorithm is almost infeasible due to the O(N4) memory requirement and O(N6) time consumption. The structure algorithm is constrained by adding pins. Long sequence comparison is achieved using anchors fixed partial comparison. Also, the time and space complexity of the algorithm is effectively reduced. However, it is clear that Consan syntax lacks the handling of more complex statistical functions [10].PMcomp used McCaskill’s algorithm to obtain the base-pairing probability matrix and, next, the alignment probability matrix. This algorithm computes the input pairing matrix independently, unlike other combinatorial alignment and folding procedures (e.g., Dynalign and FoldAlign). Promising candidate sequences are preselected using a pattern search procedure, and an algebraic dynamic programming approach is used to scan the candidate sequences. It is folded by McCakill’s algorithm, and finally, each individual is compared with each other separately using familiar techniques. The method reached a lightweight energy computation and can be described as simplified Sankoff’s algorithm [12].LocARNA is a novel RNA local comparison tool that has achieved high sensitivity in multiple experiments based on different datasets, but its complexity is still high [15]. For this feature of LocARNA, a sparsified prediction and alignment of RNAs based on their structure ensembles (SRARSE) is introduced. It supports the operation of loop structures for the first time and retains the full flexibility of Sankoff’s algorithm. In the benchmark test, high accuracy and speed were both maintained [16].Non-Sankoff algorithms divide the process into two parts: folding and alignment. Among them, some tree-based methods, such as RNAforester [23], RNAdistance [8], RNAStrAt [24], and RNApdist [25], are widely used.The RNAforester tree comparison algorithm computes the local similarity in RNA secondary structures. The method provides a new representation of RNA secondary structure in forest representation and provides elaborate dynamic programming implementations using dense two-dimensional tables, which can greatly reduce space requirements in practice. RNAdistance calculates the similarity of the RNA secondary structures by measuring the editing distance of the tree. RNApdist compares the RNA secondary structures based on base-pairing probabilities. The Vienna RNA package can now be used to implement both the RNAdistance and RNApdist methods [26, 27]. Considering the structure with pseudoknots, Michela established ASPRAlign and obtained a method for comparing ASPRA distances [28]. A comparison of algebraic RNA trees and structural RNA trees [29] is made with ignoring the primary sequence and considering only the structure of the molecule. The execution time of the ASPRAlign workbench tool to process molecular pairs was compared experimentally. The experiments show that the time complexity of the algorithm is O(n2). Notably, these comparison algorithms have high time complexity. General alignment-free methods are usually based on the numerical representations of RNA secondary structures, followed by the development of the many graphical representations of RNA secondary structures. The visualization of an RNA secondary structure using a graphical representation is more intuitive, thus providing a new way of thinking about the comparison of RNA secondary structures. Previously, researchers used eight symbols to represent RNA, but this representation was accompanied by a loss of information [30]. A feature sequence may correspond to diverse RNA secondary structures. That is, the feature sequences obtained by such a method are not unique. To reduce the loss of information, based on a 4 × 4 matrix of RNA sequences proposed by Randic [31], Yao used four horizontal lines and eight symbols to represent the RNA secondary structure [32]. The method avoids the loss of information due to the crossover and overlap of curves, but the information loss caused by the limitations of feature-invariant extraction still exists. Randic proposed a method to visualize the secondary structure of RNA without loss of information [33]. The method is based on [30] labeling bases before and after hydrogen bonds, using 12 symbols to represent bases at different positions and, in turn, using numbers and graphs to represent RNA sequences with structural information. To illustrate the validity of the method, two RNAs that are similar were selected, and visualization of RNA sequences using the method enables visual observation of the sequence differences. In [34], Li proposed a novel graphical representation of RNA secondary structures (TV-Curve) and introduced a wavelet decomposition to compare RNA secondary structures. Based on the feature representation of RNA secondary structure with eight symbols proposed in [30], each character is represented by three vectors, and the vectors are connected sequentially to generate a unique graphical representation of RNA secondary structure. After feature extraction, RNA similarity is estimated using weighted correlations of wavelet domains at different scales. The method integrates the global and local structures. The effect of the program was verified using 100 RNA sequences from the RFAM database and 9 viruses to obtain a similarity metric closer to the actual data. Based on the approach of [34], Li further created a web server for multiscale similarity RNA secondary structure comparison based on the curve representation of triple vectors (RNA-TV-curve), including RNA visualization, mutation analysis, and multiple RNA structure comparison [35].Neutral components (NCs) in the genotype-phenotype (GP) maps were found in [36] to be able to influence the biological sequence function. Thus, based on a method that reveals the hierarchical community structure, which reveals the hierarchical community structure solely from the sequence constraints and composition of the genotypes that form a given NC, community detection has been used to analyze RNA secondary structure, and recent works on community detection have been studied in [37]. Some novel classifiers and learning models have also been proposed for multigraph analysis. Multigraph feature-based representation and learning methods for multigraph classification were used to achieve an effective learning model for multigraph classification [38]. Subsequently, to address the complexity of graph data generalization in multigraph classification, Wu proposed a boosting-based multigraph classification framework (bMGC). The framework uses a two-level weighting strategy to balance label ambiguity and introduces dynamic subgraph selection criteria to solve the structured data representation problem [39].Besides, the classification and comparison of noncoding RNAs is also gaining attention. The deep-learning-based model ncRDeep is a recently proposed tool to efficiently classify ncRNA families [40]. It uses a simple convolutional neural network and RNA sequence information only to predict the class of ncRNAs. The experiments in the article evaluate the performance of ncRDeep using a benchmark dataset. Ultimately, the method was effective in improving the accuracy of prediction.Wavelet analysis is a synthesis of ideas in pure mathematics, physics, and engineering. A commonly used windowed Fourier transform, such as the short-time Fourier transform (STFT), analyzes the signal with a fixed sliding window. Obviously, fixed-length sliding window processing is not suitable for all signals. The linear time-frequency analysis for nonstationary signals should have different resolutions at different positions in the time-frequency plain; in other words, it should be a multiresolution analysis method. Wavelet transform is a multiresolution analysis method. Furthermore, wavelet is a fairly simple and effective mathematical tool that has been widely used. A combination of wavelet transform and neural network has been proposed to achieve highly accurate machine fault diagnosis [41]. In [42], wavelet analysis was used to identify membrane protein types. It was also used to realize RNA secondary structure similarity analysis. In practical applications, especially when implementing the wavelet transform on a computer, the signal must be discretized and analyzed by the discrete wavelet transform. By detecting and analyzing protein secondary structures by discrete wavelet transforms, the correlation of the amino acid sequence and secondary structure was tested by the hydrophobic value of the amino acids [43].The fractal dimension (FD) reflects the validity of the space occupied by a complex form, which is a measure of the irregularity of a complex form [44]. FD acknowledges that various parts of the world may be similar in some way to the whole region under certain conditions or processes, and it recognizes that changes in spatial dimensions can be discrete and continuous [45]. Reference [46] combines empirical modal decomposition and multifractal detrended fluctuation analysis to study the fractal characteristics of harmonic signals. Additionally, the fractal dimension has been introduced in the studies of biological molecules. Yu exploited the hydrophobicity of amino acids and fractal analysis to classify the protein structure [47]. In [48], Yang applied the fractal dimension to protein sequence similarity analysis and obtained a high degree of similarity. Yang also integrated the fractal dimension with empirical mode decomposition to compare protein sequences and optimise the feature extraction process. Experiments show that the combination of the two methods extracts higher quality features [49]. The fractal dimension has also been recently used to determine the distribution of purines and pyrimidines in the miRNAs of humans, gorillas, chimpanzees, mice, and rats [50].In this paper, we propose a novel combined algorithm. Based on the characteristics of discrete wavelet transform (DWT) multiresolution analysis and the universality of FD, DWT is used to analyze the graphical representation of RNA secondary structures and FD is used to quantitatively characterize their geometric structure features. The purpose of this paper is to explore the application of DWT and fractal dimension in the comparison of RNA secondary structures.The paper is organized as follows. We outline the theory of the fractal dimensions and the graphical characterization of RNA secondary structure in Section2. Numerical experiments and the results are then given and discussed in Section 3. Section 4 is the conclusion.


## 2. Materials and Methods

In this paper, a combination of the discrete wavelet transform, detrended fluctuation analysis, and sliding window is proposed to evaluate the similarity of RNA secondary structures based on the TV-curve of RNA. The TV-curve graphical representation of RNA is used to construct the feature vector, which is analyzed by the discrete wavelet transform, and then, the sliding window is introduced. For a signal of fixed length within the window, the fractal dimension is calculated using the detrended fluctuation analysis (DFA) method to obtain the distance matrix. Finally, the clustering analysis is performed, and the corresponding cluster analysis graph is obtained. The algorithmic process is shown in Figure1.Figure 1
The flowchart of our method.


### 2.1. Basic Concepts of Fractal Dimension

#### 2.1.1. Theoretical Fractal Dimension

The German mathematician F. Hausdorff studied the properties and quantities of singular aggregates and proposed the concept of fractal dimension called the Hausdorff dimension in1918 [51]. It is one of the most important fractal dimensions. The Hausdorff dimension can be used for any set. The theory of the Hausdorff dimension is as follows:Suppose a nonempty subsetA of an n-dimensional Euclidean space ℝn, where the diameter of A is the maximum distance of any two points in A; then,(1)A=supx−y:x,y∈A,where sup ⋅ indicates the supremum of ⋅. If C⊂∪i=1∞Ai and 0≤Ai≤δ, for any i, Aiis called a δ-cover of C.AssumingC⊂ℝn and 0≤s≤∞, for any δ>0,(2)HδsC=inf∑i=1∞Uis:Aiis aδ−cover ofC.The abovementioned equation refers to the coverage ofC whose diameter does not exceed δ, and the sum of these diameters is minimized to the power of s. The symbol inf⋅ takes the minimum value of ⋅.Asδ decreases, the clusters covering C in equation (2) also decreases, and the infimum HδsC increases accordingly. Moreover, when δ⟶0, it tends to a limit, denoted as(3)HsC=limδ⟶0HδsC.For any subsetC in ℝn, the limit exists, and the limiting value is commonly 0 or ∞, and HsC is called the s-dimensional Hausdorff measure of C.In equation (2), for a given set C⊂ℝn and δ<1, HδsC does not increase with respect to s, so it can be shown from equation (3) that HsC does not increase either. Furthermore, if t > s and Ai is the δ-cover of C, then(4)∑iAit≤∑iAit−sAis≤δt−s∑iAis.Taking the infimum,(5)HδtC≤δt−sHδsC.In the case oft > s, let δ⟶0; if HsC<∞, thenHtC=0; thus, the existence of a critical point of s makes HsC “jump” from ∞ to 0. This critical value is known as the Hausdorff dimension of F, recorded as dimHC, and occasionally called the Hausdorff−Basicovitch dimension.The definition is expressed as(6)dimHC=infs≥0:HsC=0,=sups:HsC=∞.Therefore,(7)HsC=∞,if0≤s<dimHC,0,ifs>dimHC.In addition, ifs=dimHC, HsC satisfies(8)0≤HsC≤∞.


### 2.2. Algorithms of Fractal Dimension Calculation

Although the fractal dimension is widely defined, it is subject to some limitations in practical applications. Due to its high computational complexity and the discrete and finite nature of the scale factor, we usually use a number of algorithms to approximate the fractal dimension. To date, many methods have been developed to compute fractal dimensions, such as Katz’s algorithm [52], Petrosian’s algorithm [53], Higuchi’s algorithm [54], multifractal detrended fluctuation analysis (MFDFA) [55], detrended fluctuation analysis (DFA) [56], and fluctuation analysis (FA) [57]. In terms of processing time, Katz’s algorithm and the DFA method are relatively slow. However, in terms of processing accuracy, Higuchi’s algorithm, detrended fluctuation analysis (DFA), and fluctuation analysis (FA) are relatively accurate. In the following, we will introduce Higuchi’s algorithm, fluctuation analysis (FA), and detrended fluctuation analysis (DFA).


#### 2.2.1. Higuchi’s Algorithm

Treating an RNA sequence as a one-dimensional signals1,s2,…,sN, we construct a subsequence as(9)smk=sm,sm+k,sm+2k,…,sm+N−mkk,m=1,…,k,where m indicates the starting position of the signal and k is the measurement scale. N−m/k means the integer part of N−m/k, which is the number of terms in smk.We calculate the average lengthLmk of the curve with starting position m according to the measurement scale k.Form=1,…,k,(10)Lmk=∑i=1N−m/ksm+ik−sm+i−1kN−1N−m/k,where N is the length of the signal and N−1/N−m/kk is the normalization term. With the given measurement scale, the approximate length of the signal with starting position m is computed, m=1,…,k.Thus, the approximated signal length can be obtained by(11)Lk=1k∑m=1kLmk,k=1,…,kmax.Finally, the fractal dimensionh∗ of the signal is calculated by the least-squares method [58]:(12)h∗=argminh∑k=1Kh∗log1k−logLk+c2,where c is the bias.


#### 2.2.2. Fluctuation Analysis (FA) Method

The fluctuation analysis (FA) method is commonly used to calculate the Hurst parameter for time series. It works as follows.SupposeX=Xi,i=1,2,…,N is a random process with a mean value of μ and a variance of σ2. First, we remove the mean from the signal and represent the new signal as x=xi,i=1,2,…,N, where xi=Xi−μ.Then, we construct a new one-dimensional signaly=yn,n=1,2,…,N with the sum of the first n terms of x,(13)yn=∑i=1nxi.We test whether the following formula satisfies the power-law formula:(14)F2m=yn+m−yn21/2∼mH.That is, we test whetherlog2F2m and log2m are satisfied:(15)log2F2m=H∗log2m+c,where c is a constant.From the abovementioned equation, we can obtain the Hurst parameterH.In the end, we can obtain FD by the following formula:(16)FD=2−H.


#### 2.2.3. Detrended Fluctuation Analysis (DFA) Method

Owing to the complicated characterization of long-range correlations and power-law statistics that RNA TV-curves have, it follows that comprehensively describing and studying the internal features of nonstationary signals such as RNA TV-curves by traditional methods is difficult. For complex, highly nonstationary signals, time series analysis methods derived from statistical physics are currently used. The DFA method can effectively eliminate various unknown trends in the signal, thus avoiding the interference caused by noise and signal and instability.The DFA method has also been widely used in image processing and other fields [59]. The DFA method has been proven to be an effective method to analyze nonstationary signals. In summary, as the FA, the signal xr,r=1,…,N is preliminarily processed, and the mean value x is removed. The summation of i terms before the series xr as the term i in the new sequence is(17)Yi=∑r=1ixr−<x>,i=1,2,…,N.Next, the signalY(i) is equally split into data boxes of length k. k, which is the measurement scale, is the number of points per data box.In each box of lengthk, a least-squares line is fit to the data. As a result, fitting data box ykj is obtained. Then, we detrend the signal Y.(18)Ykj˜=Yj−ykj.In the end, the different data obtained for each box are integrated and averaged.(19)F2k=1N∑j=1NFDFA2k1/2,where FDFA2k=1/k∑j=1kYk˜2j.The scaling exponentα is the slope of the line of fit of F2k and k.The fractal dimension is calculated by(20)FD=2−α.


### 2.3. Graphical Characterization of RNA Secondary Structure

#### 2.3.1. Dot-Bracket Representation

The sequence of RNA consists of the basesC, G, U, and A. The primary sequence does not contain structural information. During the exploration of the RNA function, two bases combine with each other to form base pairs. Typically, the base pairs are A-U, G-C. Thus, the secondary structure of RNA is formed. Currently, “the dot-bracket” representation is commonly used to signify the secondary structure of RNA, and the dissociative bases and base pairs are indicated by dots and parentheses, respectively. Among them, open brackets express the base pairs near the 50-terminal of the RNA chain, while closed brackets indicate the bases close to the 30-terminal. This representation can be obtained via the Vienna RNA package. However, this linear representation is degenerate in some cases in which different RNA secondary structures have the same characteristic sequence. Nevertheless, this circumstance usually occurs in short RNA sequences, which can be ignored in the analysis of long and complex RNA secondary structures.


#### 2.3.2. The TV-Curve Representation of RNA Secondary Structure

Based on the “dot-bracket” representation introduced previously, in this part, we will show a method for the feature extraction of RNA secondary structures, called “RNA triple vector curve” representation. This method was proposed by Li in 2012, as a 2-D graphical representation of RNA secondary structures, in which RNA sequence information and structural information are considered [34]. The bases in the RNA secondary structure are divided into two types: nucleotide bases paired by hydrogen bonds and unpaired nucleotide bases. The four unpaired nucleotide bases are denoted by C (cytosine), G (guanine), U (uracil), and A (adenine); in addition, C′, G′, U′, and A′ denote paired nucleotide bases. As shown in Figure 2, each of the eight symbols is described by three vectors.(21)1,−1,1,1,1,−1⇒C,1,1,1,−1,1,−1⇒C′,1,−1,1,−1,1,−1⇒G,1,1,1,1,1,−1⇒G′,1,1,1,−1,1,1⇒U,1,−1,1,1,1,1⇒U′,1,1,1,1,1,1⇒A,1,−1,1,−1,1,1⇒A′.Figure 2
Graphical representation of the eight nucleotides.
(a)(b)(c)(d)(e)(f)(g)(h)The feature sequence is read from 50-terminal to 30-terminal, and the vectors are sequentially connected to obtain the TV-curves. For example, Figure3 demonstrates the secondary structure of tRNA_Y00055.1 : 4327–4494 and the corresponding TV-curve. An RNA sequence of length N generates a TV-curve with an X-axis length of 3N. An RNA sequence and its secondary structure may produce only one TV-Curve, and similarly, a TV-curve may represent only one RNA secondary structure. Moreover, the TV-curve has no information missing from the process of feature extraction of the RNA.Figure 3
The TV-curve of tRNA_Y00055.1 : 4327–4494.
(a)(b)


### 2.4. Comparison of RNA Secondary Structures Based on the RNA TV-Curve by Detrended Fluctuation Analysis

#### 2.4.1. Discrete Wavelet Transform

The wavelet transform decomposes data, functions, or operators into components of different frequencies and then studies each component on the corresponding scale [60, 61]. We now study the wavelet transform within the range of signal analysis. Similar to the discrete Fourier transform, the discrete wavelet transform is also a time-frequency description method. The signal is decomposed into two characters: the approximation coefficients (AC) and the detailed coefficients (DC).For the signalx(t), we use two functions at the same time, wavelet function Ψx and scaling function Φx. Φx for the s(t) overviews approximation, and Ψx details the approximate function of x(t). Φx can be formulated as(22)Φx=2∑nhnΦ2x−n,where hn is the low-pass filter. The wavelet functionΨx can be calculated by(23)Ψx=2∑ngnΦ2x−n,where gn is the high-pass filterTwo sets of orthogonal functions can be extracted from the wavelet transform: shifted wavelet functionsΨx−k and scaling functions Φx−k.For any signalx=akj, the approximation coefficients (AC) and detailed coefficients (DC) in the next level can be quickly calculated as follows:(24)akj+1=∑i∈Zaijh¯i−2k,j=0,1,2,…,dkj+1=∑i∈Zaijg¯i−2k,j=0,1,2,….The decomposition is taken level by level. The lengths of AC and DC obtained after the progressive decomposition are always half the length of the level sequence antecedently. The processing of RNA secondary structures using DFA allows for the description of the overall trends and general characteristics of the signal, as well as the local details.


#### 2.4.2. Windowed Detrended Fluctuation Analysis (DFA) Method

If the fractal dimension is calculated by the DFA of the whole signal, only one scalar can be obtained. The sliding window only calculates the signal values in a fixed-length window [62]. We introduce a sliding window that moves along the signal and calculates the fractal dimension within the window. Calculating the fractal dimensions along a set sliding window yields a feature vector, not just a number. The length of this feature vector is N−k−1, where N is the signal length and k is the window width. M indicates the starting position of the signal. A least-squares line is fit to the data in each window of length k. As a result, a fitting data box ykmj is obtained. Then, we detrend signal Y:(25)Y˜kj=Yj−ykmj.For each window, we calculate the following formula:(26)FDFAm2k=1k∑j=1kYk˜2j.Next, the average is calculated:(27)F2k=1N∑j=1NFDFAm2k1/2.Transform the value ofk, plot the lb-lb graph of F2k−k, fit the straight line, calculate the slope, and obtain the scaling exponent α. The fractal dimension is calculated by(28)FD=2−α.


## 2.1. Basic Concepts of Fractal Dimension

### 2.1.1. Theoretical Fractal Dimension

The German mathematician F. Hausdorff studied the properties and quantities of singular aggregates and proposed the concept of fractal dimension called the Hausdorff dimension in1918 [51]. It is one of the most important fractal dimensions. The Hausdorff dimension can be used for any set. The theory of the Hausdorff dimension is as follows:Suppose a nonempty subsetA of an n-dimensional Euclidean space ℝn, where the diameter of A is the maximum distance of any two points in A; then,(1)A=supx−y:x,y∈A,where sup ⋅ indicates the supremum of ⋅. If C⊂∪i=1∞Ai and 0≤Ai≤δ, for any i, Aiis called a δ-cover of C.AssumingC⊂ℝn and 0≤s≤∞, for any δ>0,(2)HδsC=inf∑i=1∞Uis:Aiis aδ−cover ofC.The abovementioned equation refers to the coverage ofC whose diameter does not exceed δ, and the sum of these diameters is minimized to the power of s. The symbol inf⋅ takes the minimum value of ⋅.Asδ decreases, the clusters covering C in equation (2) also decreases, and the infimum HδsC increases accordingly. Moreover, when δ⟶0, it tends to a limit, denoted as(3)HsC=limδ⟶0HδsC.For any subsetC in ℝn, the limit exists, and the limiting value is commonly 0 or ∞, and HsC is called the s-dimensional Hausdorff measure of C.In equation (2), for a given set C⊂ℝn and δ<1, HδsC does not increase with respect to s, so it can be shown from equation (3) that HsC does not increase either. Furthermore, if t > s and Ai is the δ-cover of C, then(4)∑iAit≤∑iAit−sAis≤δt−s∑iAis.Taking the infimum,(5)HδtC≤δt−sHδsC.In the case oft > s, let δ⟶0; if HsC<∞, thenHtC=0; thus, the existence of a critical point of s makes HsC “jump” from ∞ to 0. This critical value is known as the Hausdorff dimension of F, recorded as dimHC, and occasionally called the Hausdorff−Basicovitch dimension.The definition is expressed as(6)dimHC=infs≥0:HsC=0,=sups:HsC=∞.Therefore,(7)HsC=∞,if0≤s<dimHC,0,ifs>dimHC.In addition, ifs=dimHC, HsC satisfies(8)0≤HsC≤∞.


## 2.1.1. Theoretical Fractal Dimension

The German mathematician F. Hausdorff studied the properties and quantities of singular aggregates and proposed the concept of fractal dimension called the Hausdorff dimension in1918 [51]. It is one of the most important fractal dimensions. The Hausdorff dimension can be used for any set. The theory of the Hausdorff dimension is as follows:Suppose a nonempty subsetA of an n-dimensional Euclidean space ℝn, where the diameter of A is the maximum distance of any two points in A; then,(1)A=supx−y:x,y∈A,where sup ⋅ indicates the supremum of ⋅. If C⊂∪i=1∞Ai and 0≤Ai≤δ, for any i, Aiis called a δ-cover of C.AssumingC⊂ℝn and 0≤s≤∞, for any δ>0,(2)HδsC=inf∑i=1∞Uis:Aiis aδ−cover ofC.The abovementioned equation refers to the coverage ofC whose diameter does not exceed δ, and the sum of these diameters is minimized to the power of s. The symbol inf⋅ takes the minimum value of ⋅.Asδ decreases, the clusters covering C in equation (2) also decreases, and the infimum HδsC increases accordingly. Moreover, when δ⟶0, it tends to a limit, denoted as(3)HsC=limδ⟶0HδsC.For any subsetC in ℝn, the limit exists, and the limiting value is commonly 0 or ∞, and HsC is called the s-dimensional Hausdorff measure of C.In equation (2), for a given set C⊂ℝn and δ<1, HδsC does not increase with respect to s, so it can be shown from equation (3) that HsC does not increase either. Furthermore, if t > s and Ai is the δ-cover of C, then(4)∑iAit≤∑iAit−sAis≤δt−s∑iAis.Taking the infimum,(5)HδtC≤δt−sHδsC.In the case oft > s, let δ⟶0; if HsC<∞, thenHtC=0; thus, the existence of a critical point of s makes HsC “jump” from ∞ to 0. This critical value is known as the Hausdorff dimension of F, recorded as dimHC, and occasionally called the Hausdorff−Basicovitch dimension.The definition is expressed as(6)dimHC=infs≥0:HsC=0,=sups:HsC=∞.Therefore,(7)HsC=∞,if0≤s<dimHC,0,ifs>dimHC.In addition, ifs=dimHC, HsC satisfies(8)0≤HsC≤∞.


## 2.2. Algorithms of Fractal Dimension Calculation

Although the fractal dimension is widely defined, it is subject to some limitations in practical applications. Due to its high computational complexity and the discrete and finite nature of the scale factor, we usually use a number of algorithms to approximate the fractal dimension. To date, many methods have been developed to compute fractal dimensions, such as Katz’s algorithm [52], Petrosian’s algorithm [53], Higuchi’s algorithm [54], multifractal detrended fluctuation analysis (MFDFA) [55], detrended fluctuation analysis (DFA) [56], and fluctuation analysis (FA) [57]. In terms of processing time, Katz’s algorithm and the DFA method are relatively slow. However, in terms of processing accuracy, Higuchi’s algorithm, detrended fluctuation analysis (DFA), and fluctuation analysis (FA) are relatively accurate. In the following, we will introduce Higuchi’s algorithm, fluctuation analysis (FA), and detrended fluctuation analysis (DFA).


### 2.2.1. Higuchi’s Algorithm

Treating an RNA sequence as a one-dimensional signals1,s2,…,sN, we construct a subsequence as(9)smk=sm,sm+k,sm+2k,…,sm+N−mkk,m=1,…,k,where m indicates the starting position of the signal and k is the measurement scale. N−m/k means the integer part of N−m/k, which is the number of terms in smk.We calculate the average lengthLmk of the curve with starting position m according to the measurement scale k.Form=1,…,k,(10)Lmk=∑i=1N−m/ksm+ik−sm+i−1kN−1N−m/k,where N is the length of the signal and N−1/N−m/kk is the normalization term. With the given measurement scale, the approximate length of the signal with starting position m is computed, m=1,…,k.Thus, the approximated signal length can be obtained by(11)Lk=1k∑m=1kLmk,k=1,…,kmax.Finally, the fractal dimensionh∗ of the signal is calculated by the least-squares method [58]:(12)h∗=argminh∑k=1Kh∗log1k−logLk+c2,where c is the bias.


### 2.2.2. Fluctuation Analysis (FA) Method

The fluctuation analysis (FA) method is commonly used to calculate the Hurst parameter for time series. It works as follows.SupposeX=Xi,i=1,2,…,N is a random process with a mean value of μ and a variance of σ2. First, we remove the mean from the signal and represent the new signal as x=xi,i=1,2,…,N, where xi=Xi−μ.Then, we construct a new one-dimensional signaly=yn,n=1,2,…,N with the sum of the first n terms of x,(13)yn=∑i=1nxi.We test whether the following formula satisfies the power-law formula:(14)F2m=yn+m−yn21/2∼mH.That is, we test whetherlog2F2m and log2m are satisfied:(15)log2F2m=H∗log2m+c,where c is a constant.From the abovementioned equation, we can obtain the Hurst parameterH.In the end, we can obtain FD by the following formula:(16)FD=2−H.


### 2.2.3. Detrended Fluctuation Analysis (DFA) Method

Owing to the complicated characterization of long-range correlations and power-law statistics that RNA TV-curves have, it follows that comprehensively describing and studying the internal features of nonstationary signals such as RNA TV-curves by traditional methods is difficult. For complex, highly nonstationary signals, time series analysis methods derived from statistical physics are currently used. The DFA method can effectively eliminate various unknown trends in the signal, thus avoiding the interference caused by noise and signal and instability.The DFA method has also been widely used in image processing and other fields [59]. The DFA method has been proven to be an effective method to analyze nonstationary signals. In summary, as the FA, the signal xr,r=1,…,N is preliminarily processed, and the mean value x is removed. The summation of i terms before the series xr as the term i in the new sequence is(17)Yi=∑r=1ixr−<x>,i=1,2,…,N.Next, the signalY(i) is equally split into data boxes of length k. k, which is the measurement scale, is the number of points per data box.In each box of lengthk, a least-squares line is fit to the data. As a result, fitting data box ykj is obtained. Then, we detrend the signal Y.(18)Ykj˜=Yj−ykj.In the end, the different data obtained for each box are integrated and averaged.(19)F2k=1N∑j=1NFDFA2k1/2,where FDFA2k=1/k∑j=1kYk˜2j.The scaling exponentα is the slope of the line of fit of F2k and k.The fractal dimension is calculated by(20)FD=2−α.


## 2.2.1. Higuchi’s Algorithm

Treating an RNA sequence as a one-dimensional signals1,s2,…,sN, we construct a subsequence as(9)smk=sm,sm+k,sm+2k,…,sm+N−mkk,m=1,…,k,where m indicates the starting position of the signal and k is the measurement scale. N−m/k means the integer part of N−m/k, which is the number of terms in smk.We calculate the average lengthLmk of the curve with starting position m according to the measurement scale k.Form=1,…,k,(10)Lmk=∑i=1N−m/ksm+ik−sm+i−1kN−1N−m/k,where N is the length of the signal and N−1/N−m/kk is the normalization term. With the given measurement scale, the approximate length of the signal with starting position m is computed, m=1,…,k.Thus, the approximated signal length can be obtained by(11)Lk=1k∑m=1kLmk,k=1,…,kmax.Finally, the fractal dimensionh∗ of the signal is calculated by the least-squares method [58]:(12)h∗=argminh∑k=1Kh∗log1k−logLk+c2,where c is the bias.


## 2.2.2. Fluctuation Analysis (FA) Method

The fluctuation analysis (FA) method is commonly used to calculate the Hurst parameter for time series. It works as follows.SupposeX=Xi,i=1,2,…,N is a random process with a mean value of μ and a variance of σ2. First, we remove the mean from the signal and represent the new signal as x=xi,i=1,2,…,N, where xi=Xi−μ.Then, we construct a new one-dimensional signaly=yn,n=1,2,…,N with the sum of the first n terms of x,(13)yn=∑i=1nxi.We test whether the following formula satisfies the power-law formula:(14)F2m=yn+m−yn21/2∼mH.That is, we test whetherlog2F2m and log2m are satisfied:(15)log2F2m=H∗log2m+c,where c is a constant.From the abovementioned equation, we can obtain the Hurst parameterH.In the end, we can obtain FD by the following formula:(16)FD=2−H.


## 2.2.3. Detrended Fluctuation Analysis (DFA) Method

Owing to the complicated characterization of long-range correlations and power-law statistics that RNA TV-curves have, it follows that comprehensively describing and studying the internal features of nonstationary signals such as RNA TV-curves by traditional methods is difficult. For complex, highly nonstationary signals, time series analysis methods derived from statistical physics are currently used. The DFA method can effectively eliminate various unknown trends in the signal, thus avoiding the interference caused by noise and signal and instability.The DFA method has also been widely used in image processing and other fields [59]. The DFA method has been proven to be an effective method to analyze nonstationary signals. In summary, as the FA, the signal xr,r=1,…,N is preliminarily processed, and the mean value x is removed. The summation of i terms before the series xr as the term i in the new sequence is(17)Yi=∑r=1ixr−<x>,i=1,2,…,N.Next, the signalY(i) is equally split into data boxes of length k. k, which is the measurement scale, is the number of points per data box.In each box of lengthk, a least-squares line is fit to the data. As a result, fitting data box ykj is obtained. Then, we detrend the signal Y.(18)Ykj˜=Yj−ykj.In the end, the different data obtained for each box are integrated and averaged.(19)F2k=1N∑j=1NFDFA2k1/2,where FDFA2k=1/k∑j=1kYk˜2j.The scaling exponentα is the slope of the line of fit of F2k and k.The fractal dimension is calculated by(20)FD=2−α.


## 2.3. Graphical Characterization of RNA Secondary Structure

### 2.3.1. Dot-Bracket Representation

The sequence of RNA consists of the basesC, G, U, and A. The primary sequence does not contain structural information. During the exploration of the RNA function, two bases combine with each other to form base pairs. Typically, the base pairs are A-U, G-C. Thus, the secondary structure of RNA is formed. Currently, “the dot-bracket” representation is commonly used to signify the secondary structure of RNA, and the dissociative bases and base pairs are indicated by dots and parentheses, respectively. Among them, open brackets express the base pairs near the 50-terminal of the RNA chain, while closed brackets indicate the bases close to the 30-terminal. This representation can be obtained via the Vienna RNA package. However, this linear representation is degenerate in some cases in which different RNA secondary structures have the same characteristic sequence. Nevertheless, this circumstance usually occurs in short RNA sequences, which can be ignored in the analysis of long and complex RNA secondary structures.


### 2.3.2. The TV-Curve Representation of RNA Secondary Structure

Based on the “dot-bracket” representation introduced previously, in this part, we will show a method for the feature extraction of RNA secondary structures, called “RNA triple vector curve” representation. This method was proposed by Li in 2012, as a 2-D graphical representation of RNA secondary structures, in which RNA sequence information and structural information are considered [34]. The bases in the RNA secondary structure are divided into two types: nucleotide bases paired by hydrogen bonds and unpaired nucleotide bases. The four unpaired nucleotide bases are denoted by C (cytosine), G (guanine), U (uracil), and A (adenine); in addition, C′, G′, U′, and A′ denote paired nucleotide bases. As shown in Figure 2, each of the eight symbols is described by three vectors.(21)1,−1,1,1,1,−1⇒C,1,1,1,−1,1,−1⇒C′,1,−1,1,−1,1,−1⇒G,1,1,1,1,1,−1⇒G′,1,1,1,−1,1,1⇒U,1,−1,1,1,1,1⇒U′,1,1,1,1,1,1⇒A,1,−1,1,−1,1,1⇒A′.Figure 2
Graphical representation of the eight nucleotides.
(a)(b)(c)(d)(e)(f)(g)(h)The feature sequence is read from 50-terminal to 30-terminal, and the vectors are sequentially connected to obtain the TV-curves. For example, Figure3 demonstrates the secondary structure of tRNA_Y00055.1 : 4327–4494 and the corresponding TV-curve. An RNA sequence of length N generates a TV-curve with an X-axis length of 3N. An RNA sequence and its secondary structure may produce only one TV-Curve, and similarly, a TV-curve may represent only one RNA secondary structure. Moreover, the TV-curve has no information missing from the process of feature extraction of the RNA.Figure 3
The TV-curve of tRNA_Y00055.1 : 4327–4494.
(a)(b)


## 2.3.1. Dot-Bracket Representation

The sequence of RNA consists of the basesC, G, U, and A. The primary sequence does not contain structural information. During the exploration of the RNA function, two bases combine with each other to form base pairs. Typically, the base pairs are A-U, G-C. Thus, the secondary structure of RNA is formed. Currently, “the dot-bracket” representation is commonly used to signify the secondary structure of RNA, and the dissociative bases and base pairs are indicated by dots and parentheses, respectively. Among them, open brackets express the base pairs near the 50-terminal of the RNA chain, while closed brackets indicate the bases close to the 30-terminal. This representation can be obtained via the Vienna RNA package. However, this linear representation is degenerate in some cases in which different RNA secondary structures have the same characteristic sequence. Nevertheless, this circumstance usually occurs in short RNA sequences, which can be ignored in the analysis of long and complex RNA secondary structures.


## 2.3.2. The TV-Curve Representation of RNA Secondary Structure

Based on the “dot-bracket” representation introduced previously, in this part, we will show a method for the feature extraction of RNA secondary structures, called “RNA triple vector curve” representation. This method was proposed by Li in 2012, as a 2-D graphical representation of RNA secondary structures, in which RNA sequence information and structural information are considered [34]. The bases in the RNA secondary structure are divided into two types: nucleotide bases paired by hydrogen bonds and unpaired nucleotide bases. The four unpaired nucleotide bases are denoted by C (cytosine), G (guanine), U (uracil), and A (adenine); in addition, C′, G′, U′, and A′ denote paired nucleotide bases. As shown in Figure 2, each of the eight symbols is described by three vectors.(21)1,−1,1,1,1,−1⇒C,1,1,1,−1,1,−1⇒C′,1,−1,1,−1,1,−1⇒G,1,1,1,1,1,−1⇒G′,1,1,1,−1,1,1⇒U,1,−1,1,1,1,1⇒U′,1,1,1,1,1,1⇒A,1,−1,1,−1,1,1⇒A′.Figure 2
Graphical representation of the eight nucleotides.
(a)(b)(c)(d)(e)(f)(g)(h)The feature sequence is read from 50-terminal to 30-terminal, and the vectors are sequentially connected to obtain the TV-curves. For example, Figure3 demonstrates the secondary structure of tRNA_Y00055.1 : 4327–4494 and the corresponding TV-curve. An RNA sequence of length N generates a TV-curve with an X-axis length of 3N. An RNA sequence and its secondary structure may produce only one TV-Curve, and similarly, a TV-curve may represent only one RNA secondary structure. Moreover, the TV-curve has no information missing from the process of feature extraction of the RNA.Figure 3
The TV-curve of tRNA_Y00055.1 : 4327–4494.
(a)(b)


## 2.4. Comparison of RNA Secondary Structures Based on the RNA TV-Curve by Detrended Fluctuation Analysis

### 2.4.1. Discrete Wavelet Transform

The wavelet transform decomposes data, functions, or operators into components of different frequencies and then studies each component on the corresponding scale [60, 61]. We now study the wavelet transform within the range of signal analysis. Similar to the discrete Fourier transform, the discrete wavelet transform is also a time-frequency description method. The signal is decomposed into two characters: the approximation coefficients (AC) and the detailed coefficients (DC).For the signalx(t), we use two functions at the same time, wavelet function Ψx and scaling function Φx. Φx for the s(t) overviews approximation, and Ψx details the approximate function of x(t). Φx can be formulated as(22)Φx=2∑nhnΦ2x−n,where hn is the low-pass filter. The wavelet functionΨx can be calculated by(23)Ψx=2∑ngnΦ2x−n,where gn is the high-pass filterTwo sets of orthogonal functions can be extracted from the wavelet transform: shifted wavelet functionsΨx−k and scaling functions Φx−k.For any signalx=akj, the approximation coefficients (AC) and detailed coefficients (DC) in the next level can be quickly calculated as follows:(24)akj+1=∑i∈Zaijh¯i−2k,j=0,1,2,…,dkj+1=∑i∈Zaijg¯i−2k,j=0,1,2,….The decomposition is taken level by level. The lengths of AC and DC obtained after the progressive decomposition are always half the length of the level sequence antecedently. The processing of RNA secondary structures using DFA allows for the description of the overall trends and general characteristics of the signal, as well as the local details.


### 2.4.2. Windowed Detrended Fluctuation Analysis (DFA) Method

If the fractal dimension is calculated by the DFA of the whole signal, only one scalar can be obtained. The sliding window only calculates the signal values in a fixed-length window [62]. We introduce a sliding window that moves along the signal and calculates the fractal dimension within the window. Calculating the fractal dimensions along a set sliding window yields a feature vector, not just a number. The length of this feature vector is N−k−1, where N is the signal length and k is the window width. M indicates the starting position of the signal. A least-squares line is fit to the data in each window of length k. As a result, a fitting data box ykmj is obtained. Then, we detrend signal Y:(25)Y˜kj=Yj−ykmj.For each window, we calculate the following formula:(26)FDFAm2k=1k∑j=1kYk˜2j.Next, the average is calculated:(27)F2k=1N∑j=1NFDFAm2k1/2.Transform the value ofk, plot the lb-lb graph of F2k−k, fit the straight line, calculate the slope, and obtain the scaling exponent α. The fractal dimension is calculated by(28)FD=2−α.


## 2.4.1. Discrete Wavelet Transform

The wavelet transform decomposes data, functions, or operators into components of different frequencies and then studies each component on the corresponding scale [60, 61]. We now study the wavelet transform within the range of signal analysis. Similar to the discrete Fourier transform, the discrete wavelet transform is also a time-frequency description method. The signal is decomposed into two characters: the approximation coefficients (AC) and the detailed coefficients (DC).For the signalx(t), we use two functions at the same time, wavelet function Ψx and scaling function Φx. Φx for the s(t) overviews approximation, and Ψx details the approximate function of x(t). Φx can be formulated as(22)Φx=2∑nhnΦ2x−n,where hn is the low-pass filter. The wavelet functionΨx can be calculated by(23)Ψx=2∑ngnΦ2x−n,where gn is the high-pass filterTwo sets of orthogonal functions can be extracted from the wavelet transform: shifted wavelet functionsΨx−k and scaling functions Φx−k.For any signalx=akj, the approximation coefficients (AC) and detailed coefficients (DC) in the next level can be quickly calculated as follows:(24)akj+1=∑i∈Zaijh¯i−2k,j=0,1,2,…,dkj+1=∑i∈Zaijg¯i−2k,j=0,1,2,….The decomposition is taken level by level. The lengths of AC and DC obtained after the progressive decomposition are always half the length of the level sequence antecedently. The processing of RNA secondary structures using DFA allows for the description of the overall trends and general characteristics of the signal, as well as the local details.


## 2.4.2. Windowed Detrended Fluctuation Analysis (DFA) Method

If the fractal dimension is calculated by the DFA of the whole signal, only one scalar can be obtained. The sliding window only calculates the signal values in a fixed-length window [62]. We introduce a sliding window that moves along the signal and calculates the fractal dimension within the window. Calculating the fractal dimensions along a set sliding window yields a feature vector, not just a number. The length of this feature vector is N−k−1, where N is the signal length and k is the window width. M indicates the starting position of the signal. A least-squares line is fit to the data in each window of length k. As a result, a fitting data box ykmj is obtained. Then, we detrend signal Y:(25)Y˜kj=Yj−ykmj.For each window, we calculate the following formula:(26)FDFAm2k=1k∑j=1kYk˜2j.Next, the average is calculated:(27)F2k=1N∑j=1NFDFAm2k1/2.Transform the value ofk, plot the lb-lb graph of F2k−k, fit the straight line, calculate the slope, and obtain the scaling exponent α. The fractal dimension is calculated by(28)FD=2−α.


## 3. Results and Discussion

### 3.1. Results

In this section, the optimal experimental setting was obtained by the trial-and-error method: the level of wavelet decomposition was set to 3 and the window width was set to 17. Figure4 shows the DWT and DFA on the TV-curve graphical representation of a RNA secondary structure using our method. Nine commonly used viruses were selected in this experiment. These viruses have conserved secondary structures, and the differences between them are subtle enough to verify the effectiveness of the experimental method in the classification of secondary structures with slight differences. Information on these nine viruses is presented in Figure 5. The nine viruses include the alfalfa mosaic virus (ALMV), apple mosaic virus (APMV), citrus leaf rugose virus (CiLRV), citrus variegation virus (CVV), elm mottle virus (EMV), lilac ring mottle virus (LRMV), prune dwarf ilarvirus (PDV), asparagus virus II (AVII), and tobacco streak virus (TSV) [63].Figure 4
(a) Graphical representation of RNA secondary structure. (b) The third-level discrete wavelet transform. (c) The DFA fractal dimension of (b) obtained with a sliding window width of 17.
(a)(b)(c)Figure 5
Secondary structure information of nine viruses.
(a)(b)(c)(d)(e)(f)(g)(h)(i)The distance matrix of the nine viruses obtained by the method of this paper is shown in Table1. Next, as shown in Figure 6, the proposed method is compared with the results obtained by MEGA software [64], while the classification effect was compared with classical methods such as RNApdist and the innovative methods ASPRA distance and Li’s method. The cluster analysis graph generated by the detrended fluctuation analysis method is based on the distance matrix calculated in MATLAB. In MEGA software, the maximum likelihood method is used to generate the cluster analysis graph.Table 1
The distance matrix of nine viruses generated by our method.
ALMVAPMVCiLRVTSVCVVPDVLRMVEMVAVIIALMV00.7180.6410.7690.6670.5380.6920.8460.692APMV00.7180.7180.6670.8460.7180.8210.641CiLRV00.6920.4620.7180.4870.7180.590TSV00.7440.7950.7440.7690.769CVV00.7690.2820.6150.718PDV00.7180.7440.718LRMV00.5900.692EMV00.846AVII0Figure 6
The cluster analysis graph of nine virus sequences generated by (a) our method, (b) Li’s method, (c) RNApdist, (d) ASPRAlign, and (e) MEGA software.
(a)(b)(c)(d)(e)In addition, an experiment containing 11 noncoding RNAs is performed to test the applicability of our method in comparing the similarities between noncoding RNAs. Eight sequences are randomly selected from the ncRNAs in the RFAM database; additionally, three human ncRNA sequences from the NONCODE database are randomly sampled for the experiment. The information of these ncRNAs is provided in Table2, and the distance matrix generated by the experiment is shown in Table 3. Based on this dataset, the results of the classification using ncRDeep are shown in Table 4.Table 2
11 noncoding RNA sequences from the RFAM and NONCODEv5 databases.
NameLengthFamilyNONHSAT000002.21653NONCODEv5_humanNONHSAT000003.21483NONHSAT000004.2632MF489813.1849tRNAMF489812.1535MF489811.1549MF489810.1566MF489809.1549MF489808.1552MF489806.1554NT_033777.3299RNase P RNATable 3
The distance matrix of eleven ncRNA sequences produced by our method.
NONHSAT000002.2NONHSAT000003.2NONHSAT000004.2MF489813.1MF489812.1MF489811.1MF489810.1MF489809.1MF489808.1MF489806.1NT_033777.3NONHSAT000002.200.38820.48430.51230.51780.51530.52160.52120.52790.53940.5510NONHSAT000003.200.45710.53370.55050.55040.51480.51730.51600.52520.5636NONHSAT000004.200.51840.41990.42140.46980.52210.51990.49990.5348MF489813.100.49890.48020.49050.48060.48240.53340.5929MF489812.100.02040.36030.36340.42060.45570.4776MF489811.100.35970.35190.41330.46160.4736MF489810.100.30490.41570.46320.5165MF489809.100.37310.44450.4902MF489808.100.41390.4933MF489806.100.4591NT_033777.30Table 4
The classification results of ncRDeep.
Sequence IDResultNONHSAT000002.2RibozymesNONHSAT000003.2RibozymesNONHSAT000004.2RibozymesMF489813.1Intron_gplMF489812.1Intron_gplMF489811.1Intron_gplMF489810.1Intron_gplMF489809.1Intron_gplMF489808.1Intron_gplMF489806.1Intron_gplNT_033777.3RibozymesUnfortunately, ncRDeep does not provide the corresponding distance matrix as well as the cluster analysis graph like the other methods, which would not provide a complete comparison in this paper. The remaining experimental results, including the proposed method, RNApdist, ASPRAlign, MEGA software, and Li’s method, are shown in Figure7.Figure 7
The cluster analysis graphs of 11 noncoding RNA sequences generated by (a) our method, (b) Li’s method, (c) RNApdist, (d) ASPRAlign, and (e) MEGA software.
(a)(b)(c)(d)(e)To get rid of the specificity, we randomly extracted 100 ncRNA sequences from four families in the RFAM database, including tRNA, RNase P arch, miRNA, and 5s rRNA. We constructed three cluster analysis graphs using our method, Li’s method, RNApdist, and ASPRAlign (see Figure8: Figure S1, Figure S2, and Figure S3). As shown in Figure 8, the new dataset was applied to build a cluster analysis graph using the proposed method. Four branches can be clearly shown in Figure 8, where the 5 s rRNA family is thoroughly distinguished, while four ncRNAs in tRNA failed to successfully assign with most of the same family. As well, in the miRNA and RNase P arch families, there are small numbers of ncRNAs that are distant from their major branches. However, among the method proposed by us, RNApdist, Li’s method, and ASPRAlign, the classification effect of the proposed method is clearly superior to that of the other three methods.Figure 8
The cluster analysis graph of 100 noncoding RNA sequences generated by our method.


### 3.2. Discussion

A detailed explanation of the cluster analysis graphs is given below; it is apparent from Figure6 that the cluster analysis graph shown in Figure 6(a) is more similar to the standard cluster analysis graph (Figure 6(e)) generated by MEGA software. First, AVII and APMV and ALMV and PDV have a high similarity, which is also clearly reflected in Figure 6(a). However, it is not reflected in the other cluster analysis graphs, such as Figures 6(b)–6(d), and they deviate from the standard cluster analysis graph generated by MEGA software. Obviously, there are three branches in Figure 6(e). CiLRV is closer to LRMV, EMV, and TSV than the remaining four, and Figures 6(a) and 6(b) have similar representations which is not reflected in Figures 6(c) and 6(d). In addition, EMV is closer to LRMV and CVV in all viruses, which is visible at a glance in Figure 6(e), while Figure 6(a) has the same branch. They are, however, separated in Figures 6(b) and 6(c), as well as in Figure 6(d). In conclusion, compared with other methods, the proposed method can obtain the similarities of RNA secondary structures more accurately.The experimental results of the same dataset, the cluster analysis graphs, and the distance between them are measured using the classical Robinson–Foulds (RF) metric [65]. In [65], the smallest number of steps to transform the topology of one of the trees into the other tree by the necessary operation (α or α−1) is noted as the distance between the two trees.α is the contraction of disjoint, possibly empty set, or zero distance nodes with other normal nodes into a single node while the label changes.(29)T2=αT1,prps,where trees T1 and T2 have the same set of species, prps is the edge, and T2 is obtained from T1 by shrinking the edge prps.In turn,α−1 is the inverse operation of α, which means dispersing a node into two nodes.In this paper, the metric was used to calculate the distance between evolutionary trees. Thus, the classification efficiency of the method is verified.The distance between evolutionary trees generated by the different methods and those generated by MEGA software is shown in Table5. Clearly, the cluster analysis graph generated by our approach is closer to the cluster analysis graph generated by MEGA software than the other approaches.Table 5
The Pearson correlation efficiency between the cluster analysis graphs generated by the different methods and those generated by MEGA software.
MethodRNApdistOur methodLi’s methodASPRAlignDistance0.15490.12880.29730.4709In the second experiment, as indicated in Figure7, the experimental sequences are accurately divided into three families in Figure 7(e), tRNA, RNase P RNA, and a group of human ncRNA in the NONCODE database. Figure 7(e) shows that the NT_033777.3 sequence has a high similarity to the tRNA family, which is also reflected in Figure 7(a) and Table 3. In Figure 7(b), NT_033777.3 and MF489813.1 are classified into a group, with NT_033777.3 being the most distant from the other sequences in the tRNA family, which is unreasonable. NT_033777.3 is mistakenly placed in a group close to the NONCODE human noncoding RNA in Figure 7(c). More unreliably, the grouping in Figure 7(d) does not completely separate the three groups of RNAs. On analysis of Table 4, it is apparent that ncRDeep did not succeed in accurately distinguishing ncRNAs. Based on the abovementioned analysis, the method based on DFA and DWT is an effective algorithm for RNA secondary structure similarity analysis. The distance between evolutionary trees generated by the method proposed by us, RNApdist, and Li’s method, ASPRAlign, and the cluster analysis graph generated by MEGA software are shown in Table 6. Evidently, the results obtained by our method in this experiment are more accurate.Table 6
The Pearson correlation efficiency between the cluster analysis graphs generated by the different methods and those generated by MEGA software.
MethodRNApdistOur methodLi’s methodASPRAlignDistance0.24940.17480.20480.2722


## 3.1. Results

In this section, the optimal experimental setting was obtained by the trial-and-error method: the level of wavelet decomposition was set to 3 and the window width was set to 17. Figure4 shows the DWT and DFA on the TV-curve graphical representation of a RNA secondary structure using our method. Nine commonly used viruses were selected in this experiment. These viruses have conserved secondary structures, and the differences between them are subtle enough to verify the effectiveness of the experimental method in the classification of secondary structures with slight differences. Information on these nine viruses is presented in Figure 5. The nine viruses include the alfalfa mosaic virus (ALMV), apple mosaic virus (APMV), citrus leaf rugose virus (CiLRV), citrus variegation virus (CVV), elm mottle virus (EMV), lilac ring mottle virus (LRMV), prune dwarf ilarvirus (PDV), asparagus virus II (AVII), and tobacco streak virus (TSV) [63].Figure 4
(a) Graphical representation of RNA secondary structure. (b) The third-level discrete wavelet transform. (c) The DFA fractal dimension of (b) obtained with a sliding window width of 17.
(a)(b)(c)Figure 5
Secondary structure information of nine viruses.
(a)(b)(c)(d)(e)(f)(g)(h)(i)The distance matrix of the nine viruses obtained by the method of this paper is shown in Table1. Next, as shown in Figure 6, the proposed method is compared with the results obtained by MEGA software [64], while the classification effect was compared with classical methods such as RNApdist and the innovative methods ASPRA distance and Li’s method. The cluster analysis graph generated by the detrended fluctuation analysis method is based on the distance matrix calculated in MATLAB. In MEGA software, the maximum likelihood method is used to generate the cluster analysis graph.Table 1
The distance matrix of nine viruses generated by our method.
ALMVAPMVCiLRVTSVCVVPDVLRMVEMVAVIIALMV00.7180.6410.7690.6670.5380.6920.8460.692APMV00.7180.7180.6670.8460.7180.8210.641CiLRV00.6920.4620.7180.4870.7180.590TSV00.7440.7950.7440.7690.769CVV00.7690.2820.6150.718PDV00.7180.7440.718LRMV00.5900.692EMV00.846AVII0Figure 6
The cluster analysis graph of nine virus sequences generated by (a) our method, (b) Li’s method, (c) RNApdist, (d) ASPRAlign, and (e) MEGA software.
(a)(b)(c)(d)(e)In addition, an experiment containing 11 noncoding RNAs is performed to test the applicability of our method in comparing the similarities between noncoding RNAs. Eight sequences are randomly selected from the ncRNAs in the RFAM database; additionally, three human ncRNA sequences from the NONCODE database are randomly sampled for the experiment. The information of these ncRNAs is provided in Table2, and the distance matrix generated by the experiment is shown in Table 3. Based on this dataset, the results of the classification using ncRDeep are shown in Table 4.Table 2
11 noncoding RNA sequences from the RFAM and NONCODEv5 databases.
NameLengthFamilyNONHSAT000002.21653NONCODEv5_humanNONHSAT000003.21483NONHSAT000004.2632MF489813.1849tRNAMF489812.1535MF489811.1549MF489810.1566MF489809.1549MF489808.1552MF489806.1554NT_033777.3299RNase P RNATable 3
The distance matrix of eleven ncRNA sequences produced by our method.
NONHSAT000002.2NONHSAT000003.2NONHSAT000004.2MF489813.1MF489812.1MF489811.1MF489810.1MF489809.1MF489808.1MF489806.1NT_033777.3NONHSAT000002.200.38820.48430.51230.51780.51530.52160.52120.52790.53940.5510NONHSAT000003.200.45710.53370.55050.55040.51480.51730.51600.52520.5636NONHSAT000004.200.51840.41990.42140.46980.52210.51990.49990.5348MF489813.100.49890.48020.49050.48060.48240.53340.5929MF489812.100.02040.36030.36340.42060.45570.4776MF489811.100.35970.35190.41330.46160.4736MF489810.100.30490.41570.46320.5165MF489809.100.37310.44450.4902MF489808.100.41390.4933MF489806.100.4591NT_033777.30Table 4
The classification results of ncRDeep.
Sequence IDResultNONHSAT000002.2RibozymesNONHSAT000003.2RibozymesNONHSAT000004.2RibozymesMF489813.1Intron_gplMF489812.1Intron_gplMF489811.1Intron_gplMF489810.1Intron_gplMF489809.1Intron_gplMF489808.1Intron_gplMF489806.1Intron_gplNT_033777.3RibozymesUnfortunately, ncRDeep does not provide the corresponding distance matrix as well as the cluster analysis graph like the other methods, which would not provide a complete comparison in this paper. The remaining experimental results, including the proposed method, RNApdist, ASPRAlign, MEGA software, and Li’s method, are shown in Figure7.Figure 7
The cluster analysis graphs of 11 noncoding RNA sequences generated by (a) our method, (b) Li’s method, (c) RNApdist, (d) ASPRAlign, and (e) MEGA software.
(a)(b)(c)(d)(e)To get rid of the specificity, we randomly extracted 100 ncRNA sequences from four families in the RFAM database, including tRNA, RNase P arch, miRNA, and 5s rRNA. We constructed three cluster analysis graphs using our method, Li’s method, RNApdist, and ASPRAlign (see Figure8: Figure S1, Figure S2, and Figure S3). As shown in Figure 8, the new dataset was applied to build a cluster analysis graph using the proposed method. Four branches can be clearly shown in Figure 8, where the 5 s rRNA family is thoroughly distinguished, while four ncRNAs in tRNA failed to successfully assign with most of the same family. As well, in the miRNA and RNase P arch families, there are small numbers of ncRNAs that are distant from their major branches. However, among the method proposed by us, RNApdist, Li’s method, and ASPRAlign, the classification effect of the proposed method is clearly superior to that of the other three methods.Figure 8
The cluster analysis graph of 100 noncoding RNA sequences generated by our method.


## 3.2. Discussion

A detailed explanation of the cluster analysis graphs is given below; it is apparent from Figure6 that the cluster analysis graph shown in Figure 6(a) is more similar to the standard cluster analysis graph (Figure 6(e)) generated by MEGA software. First, AVII and APMV and ALMV and PDV have a high similarity, which is also clearly reflected in Figure 6(a). However, it is not reflected in the other cluster analysis graphs, such as Figures 6(b)–6(d), and they deviate from the standard cluster analysis graph generated by MEGA software. Obviously, there are three branches in Figure 6(e). CiLRV is closer to LRMV, EMV, and TSV than the remaining four, and Figures 6(a) and 6(b) have similar representations which is not reflected in Figures 6(c) and 6(d). In addition, EMV is closer to LRMV and CVV in all viruses, which is visible at a glance in Figure 6(e), while Figure 6(a) has the same branch. They are, however, separated in Figures 6(b) and 6(c), as well as in Figure 6(d). In conclusion, compared with other methods, the proposed method can obtain the similarities of RNA secondary structures more accurately.The experimental results of the same dataset, the cluster analysis graphs, and the distance between them are measured using the classical Robinson–Foulds (RF) metric [65]. In [65], the smallest number of steps to transform the topology of one of the trees into the other tree by the necessary operation (α or α−1) is noted as the distance between the two trees.α is the contraction of disjoint, possibly empty set, or zero distance nodes with other normal nodes into a single node while the label changes.(29)T2=αT1,prps,where trees T1 and T2 have the same set of species, prps is the edge, and T2 is obtained from T1 by shrinking the edge prps.In turn,α−1 is the inverse operation of α, which means dispersing a node into two nodes.In this paper, the metric was used to calculate the distance between evolutionary trees. Thus, the classification efficiency of the method is verified.The distance between evolutionary trees generated by the different methods and those generated by MEGA software is shown in Table5. Clearly, the cluster analysis graph generated by our approach is closer to the cluster analysis graph generated by MEGA software than the other approaches.Table 5
The Pearson correlation efficiency between the cluster analysis graphs generated by the different methods and those generated by MEGA software.
MethodRNApdistOur methodLi’s methodASPRAlignDistance0.15490.12880.29730.4709In the second experiment, as indicated in Figure7, the experimental sequences are accurately divided into three families in Figure 7(e), tRNA, RNase P RNA, and a group of human ncRNA in the NONCODE database. Figure 7(e) shows that the NT_033777.3 sequence has a high similarity to the tRNA family, which is also reflected in Figure 7(a) and Table 3. In Figure 7(b), NT_033777.3 and MF489813.1 are classified into a group, with NT_033777.3 being the most distant from the other sequences in the tRNA family, which is unreasonable. NT_033777.3 is mistakenly placed in a group close to the NONCODE human noncoding RNA in Figure 7(c). More unreliably, the grouping in Figure 7(d) does not completely separate the three groups of RNAs. On analysis of Table 4, it is apparent that ncRDeep did not succeed in accurately distinguishing ncRNAs. Based on the abovementioned analysis, the method based on DFA and DWT is an effective algorithm for RNA secondary structure similarity analysis. The distance between evolutionary trees generated by the method proposed by us, RNApdist, and Li’s method, ASPRAlign, and the cluster analysis graph generated by MEGA software are shown in Table 6. Evidently, the results obtained by our method in this experiment are more accurate.Table 6
The Pearson correlation efficiency between the cluster analysis graphs generated by the different methods and those generated by MEGA software.
MethodRNApdistOur methodLi’s methodASPRAlignDistance0.24940.17480.20480.2722


## 4. Conclusions

This paper proposes a hybrid method for the similarity comparison of RNA secondary structures. The algorithm is based on the existing RNA triple vector graphical representation, uses DWT to process the feature sequences, and captures the fractal characteristics using the DFA method. Compared with several commonly used RNA comparison methods, the approximate relationships between the RNA secondary structures obtained by the DFA and wavelet transform method are close to the actual relationships. However, the secondary structures predicted by the minimum free energy in the Vienna RNA package are not optimal, and finding the optimal secondary structure for RNA rapidly and efficiently remains a challenging problem. In addition, the method is not yet excellent for analyzing shorter RNAs, and systematic studies will be hoped to be carried out on RNAs of different lengths and characteristics.



---

*Source: 8841822-2021-03-23.xml*