Sequence-Based Predicting Bacterial Essential ncRNAs Algorithm by Machine Learning

1  Introduction

Bacterial non-coding RNA (ncRNA) is a newly discovered gene expression regulator in bacterial genome in recent years [1]. The molecular size of ncRNA is 40–500 nucleotides. It plays an important role in RNA transcription regulation, chromosome replication, RNA processing and modification, mRNA translation and stability, protein degradation and transport, bacterial infection and other biological processes [2]. With the rapid increase in the number of bacterial ncRNAs discovered and their important role in organisms, bacterial ncRNAs have become one of the hot spots in microbial research. Since ncRNAs play important roles in organisms, the identification of new ncRNAs is of great scientific significance and commercial value [3].

Some ncRNAs are indispensable for the survival of organisms. Thus, they are called “essential non-coding RNA” (essential ncRNA) [4]. Although essential ncRNAs cannot encode proteins, they are as important as essential coding genes in biology. For example, bacterial essential ncRNAs play vital role ranging from structural regulation to catalysis, affecting a variety of processes, such as bacterial toxicity, development control, mRNA stability and protein degradation. Therefore, these groups of ncRNA could be used as a potential target for drug development [5]. Furthermore, significant number of theoretical researches revealed that essential ncRNAs are helpful to understand and determine the composition and lifestyle of minimum genome. For example, Gil and Christen et al. explored that a complete minimum genome should include regulatory and structural elements in addition to encoding proteins, such as 5′-UTRs and ncRNA [6,7]. Similarly, Maria et al., reported a minimal cell containing essential ncRNA [8]. In our previous work, we also constructed a minimum bacterial gene set, which also contain the minimum non-coding RNA set besides the minimum coding gene sets [9].

With the development of high-throughput sequencing and genome-scale inactivation technologies, it has been possible to identify the essential ncRNAs contained in the genome. At present, the determination of essential ncRNAs mainly relies on experimental methods. However, these traditional methods are time-consuming, laborious and costly. In addition, traditional experimental methods are limited by the organisms, as less than 1% of bacteria can be cultured in the laboratory. In this context, it is urgent to establish effective new algorithms for theoretical analysis, and to mine the information and knowledge of essential ncRNA [10]. Pattern recognition and informatics theory play an extremely important role in biological data analysis, especially in sequence recognition research. With the development of machine learning and artificial intelligence [11,12], more and more research groups are applying machine learning to the prediction of essential genes and essential ncRNAs [13]. Zhang integrated features derived from sequence data and protein-protein interaction network and developed a machine learn method called DeepHE to predict human essential genes [14]. In terms of recognition of essential genes of bacteria, we have developed general recognition software Geptop for genome-wide recognition of essential genes, with recognition accuracy is the best so far [15]. The intrinsic features of sequences are the easiest to access for constructing a predicting essential model [16]. Using the intrinsic features of sequences, the accuracy of the identification of essential genes could be significantly increased. Machine learning is an effective predictive method [17], which has been used in classification of medical field and image processing [11,17–23]. At present, there are few researches on the theoretical recognition of “essential ncRNA”. Therefore, it is important and necessary to develop theories and methods for the recognition of essential non-coding RNA.

At present, the number of species required for the determination of essential ncRNAs in the whole genome scale is still very few. The essential ncRNAs have important scientific significance and application in biomedical research. Facing the urgent need to obtain essential ncRNAs, we have conducted in-depth research on the theoretical identification of bacterial essential ncRNAs. The purpose of this study was to design a universal model, which is less dependent on training sets, and to develop a computer based automatic method for identification of essential ncRNA of bacteria with high accuracy and universality. The findings of the study can be utilized to extend the existing drug target library, greatly accelerate the process of new drug research and shorten the research cycle, reduce costs, promote the development of new rapid detection of pathogenic bacteria system. The findings can also assist development of new pathogenic bacteria prevention and treatment methods, and the minimal genome chassis in synthetic biology research.

2  Materials and Methods

2.1 Training and Dataset Constructed

The quality of training dataset has a great influence on the predicted performance in machine learning field. The cross-organism validation result of Geptop shows that not all species achieve a high area under the curve (AUC). The essential gene dataset of species with low AUC may be unreliable due to differences in wet-lab experimental conditions. To build a reliable training dataset for developing a predicted essential ncRNA model in bacteria, we analyze the bacterial essential genes in database of essential gene (DEG) by Geptop and CEG_Match [24,25]. There are nine bacteria with their essential ncRNAs determined through wet-lab experiment. Their essential ncRNAs are obtained from DEG, and the non-essential ncRNAs were downloaded from RNAcentral (https://rnacentral.org), as shown in Table 1 [26]. The ratio of the number of essential ncRNA to the number of non-essential ncRNA in each bacterium is calculated in order to construct a reliable training dataset for proposing a bacterial essential ncRNA recognition algorithm. Six datasets with ratio in the range from 0.1 to 0.7 are extracted for further analyses. Data of Sphingomonas wittichii is used as test dataset for assessing the accuracy of model of essential ncRNA recognition. Data of the remaining bacteria is used as the training dataset.

2.2 Training and Data Marking

In this paper, we use the support vector machine learning approach, a supervised machine-learning method, to predict the essential ncRNA. Thus, every piece of data in the training dataset must be marked. In the training dataset, the essential ncRNA is marked as 1 and non-essential ncRNA is marked as 0.

2.3 Features Extraction

2.3.1 Homology Feature

For a DNA sequence, the sequence homology is often used in prediction of protein function [10]. Because negative selection acts more strongly on essential ncRNA, essential ncRNAs are more conserved than non-essential ncRNAs. Thus, this is an important feature that distinguishes between essential and non-essential ncRNA. For a predicted gene or ncRNA, we compared it with both essential and non-essential ncRNA in the reference library by BLASTN with e < 10−10 [27]. After sequence alignment, we use the average value of 10 sequences with the highest BLASTN score as the indicator for homology feature. For those alignments less than 10 hits, we use the average value of all hit sequences as the as the indicator for homology feature.

2.3.2 Sequence Related Feature

The first sequence feature is the length of sequence, which is defined as L. The GC content is a main factor shaping the amino acid usage during bacterial evolution process [28]. Thus, the GC content is used as a feature for perdition of essential ncRNAs. The GC content (GC%) is calculated using the following formula:

GC%=count(G)+count(C)L (1)

where count (G) is the number of guanine in a DNA sequence, and count (C) is the number of cytosine in a DNA sequence. The L indicates the length of a DNA sequence, including guanine, cytosine, thymine and adenine.

2.3.3 Nucleotide Concatemer Feature

We obtain the contents of single-nucleotide, dimeric-nucleotide, triplet-nucleotide, quadruplexes-nucleotide, pentamers-nucleotide and hexon-nucleotide as the features in the prediction models [29]. The above-mentioned features were calculated considering phase-independent. We calculate the content of triplet-nucleotide as an example in following paragraph.

We set window length of three and step size of one when intercepting and calculating the counts of triplet-nucleotide in a sequence. For triplet-nucleotide, sixty-four (43) combinations were computed, such as ATG, ACG, GTC and so on. For each combination, we calculate its frequency in following formula:

fi=Ci∑i=164Ci,i=1,2,...64 (2)

where i is the ith triplet-nucleotide. Ci indicates the number of the ith triplet-nucleotide in a sequence. The frequencies of other nucleotide concatermers are calculated using similar approached used to calculate the triplet-nucleotide using Eq. (3).

fik=Ci∑i=14kCi,i=1,2,…4k (3)

2.3.4 Markov Remotely Related Base Feature

The Markov chain model in stochastic process theory assumes that the state of the next character is determined by successive characters adjacent to it [30]. Thus, we used this theory to obtain the remote correlation information of nucleotide in gene or ncRNA sequences. It means that a nucleotide in a gene or ncRNA is determined by successive nucleotides adjacent to it. DNA is not linear in organisms. Gene and ncRNA have a specific structure due to the curling of DNA and the principle of nucleotide complementarity among nucleotides at different locations.

By the same analogy, we extend Markov chain model that nucleotides spaced some distance apart (lambda nucleotides) are also correlated. This correlation is not the actual correlation of a Markov chain but the correlation between sequence and tertiary structure. Experimental evidence for this correlation is three-dimensional genomics studies. The formula used for calculating the correction is presented in Eq. (4).

{xkx=(pk(XA)+pk(XG))−(pk(XC)+pk(XT))ykx=(pk(XA)+Pk(XC))−(pk(XG)+pk(XT)),X=A,C,G,T, k=1:λ+1,2:λ+2,3:λ+3zkx=(pk(XA)+pk(XI))−(pk(XG)+pk(XC)) λ=1,2,3…L (4)

where A, C, G and T indicate the four nucleotides in DNA sequence, respectively. The x, y and z mean the effects among the four nucleotides (guanine, cytosine, thymine and adenine) depending on the bond properties among them, respectively [31]. k represents the position of two relevant nucleotides. The first nucleotide of k respectively falls into first, second and third phases. λ denotes the spacing of positions between two nucleotides in sequences level. And L represents the maximum spacing in a sequence. Thus, a DNA sequence could be transformed into 3 × 4 × λ variables.

2.3.5 Mutual Information Feature

The mutual information is used for measuring the correlation of two random variables [32,33]. Based on the mutual information, the information of a random variable could be provided by another random variable. In predicting essential gene and ncRNA model, the mutual information is used to measure the information between successive nucleotides X and Y of a DNA sequence, which is defined as:

I(X,Y)=∑x∈Ω∑y∈ΩP(x,y)log2P(x,y)P(x)P(y) (5)

where Ω is the set of,{A, T, C, G} and P(x,y) is the joint probability of X and Y. P(x) and P(y) correspond to the probability of X nucleotide and Y nucleotide, respectively. The mutual information of a DNA sequence is calculated with the Eq. (5). For each dimeric-nucleotide (x,y), its mutual information is defined by the below formula as a feature:

P(x,y)log2P(x,y)P(x)P(y) (6)

For a DNA sequence, there are 42 = 16 dimeric-nucleotides and hence generating 16 features. Thus, considering the mutual information of a DNA sequence, we obutain 17 mutual informations features in total.

2.3.6 Conditional Mutual Information Feature

Conditional mutual information refers to the mutual information between the random variables X and Y in the case of a third random variable Z [34]. It is defined as:

I(X;Y|Z)=∑z∈ΩP(z)∑x∈Ω∑y∈ΩP(x,y|z)log2P(x,y|z)P(x|z)P(y|z)=∑z∈Ω∑x∈Ω∑y∈ΩP(x,y,z)log2P(z)P(z,y,z)P(x,z)P(y,z) (7)

where P(x,y,z), P(x,z) and P(y,z) are the joint probability of the random variables in brackets. The three positions of triplet-nucleotide are considered to the random variables of X, Y and Z, in sequence. According to the formula (7), the mutual information between the first base and the third base in the condition of the second base is calculated, and is used as a feature. In a DNA sequence, we could obtain 43 = 64 triplet-nucleotide. The conditional mutual information of each triplet-nucleotide is calculated by the formula given in Eq. (8):

∑z∈Ω∑x∈Ω∑y∈ΩP(x,y,z)log2P(z)P(x,y,z)P(x,z)P(y,z) (8)

In total, there are 65 conditional mutual information features.

2.3.7 Kullback-Leibler Divergence Feature

Kullback-Leibler divergence is used for mensuration of the similarity between probability distribution P(x) and model distribution Q(x), which is calculated by the formula in [35] :

KLD=∑iP(x)log2P(x)Q(x) (9)

The frequencies of nucleotide, dinucleotides and trinucleotides in a given sequence to be predicted are compared with the corresponding frequencies of the genomes used to train the model. For a given dataset, one feature is extracted by Kullback-Leibler Divergence method. In total, we get three features in this step.

2.4 Feature Normalization

After obtaining a large number of features generated by the sequences and information theory of essential ncRNA and non-essential ncRNA, all features need to be normalized to the same scale. Accordingly,, we used the Min-Max Normalization approach [36]. The formula is used for normalizing each feature.

X∗=X−minmax−min (10)

2.5 Positive and Negative Samples Balanced

Table 1 shows that the number of essential ncRNA and non-essential ncRNA are differ greatly. It is possible to misevaluate the performance of a machine learning model. For example, Bacillus subtilis 168 has 271 essential genes and 4175 non-essential genes. If all genes are predicted as non-essential genes in a model, then the accuracy rate of the model is still as high as 94% (4175/(4175+271) = 0.939). However, this model is worthless.

To avoid this situation, balancing the number of positive and negative samples is crucial process. The conventional method to solve the imbalance of positive and negative samples is the subset sampling method, which randomly select the same amount of negative sample as the positive sample [37]. However, only part of the known information from the negative sample is used with subset sampling method [38]. Thus, a large amount of potentially helpful information can be lost, which in turn affects the accuracy of the model. In this paper, we attempt to introduce principal component regression to solve this problem. As a new multivariate statistical method, principal component regression is mainly used in multivariate correlation analysis. Principal component regression has hardly been applied to the classification of biomolecules. Principal component regression uses regression technique to make classification prediction, which is not affected by sample size. The principal component regression includes steps of principal component analysis and correlation analysis at the same time, which can be used to eliminate the redundancy of variables and screen significant variables at a time.

In general, the number of negative samples (non-essential coding gene and ncRNAs) is tens of times that of positive samples (essential coding gene and ncRNAs). For example, the number of non-essential gene is 14 times more than essential gene in Bacteroides thetaiotaomicron VPI 5482, (Table 1). And the number of non-essential ncRNA is 6 times more than essential ncRNA in Mycobacterium tuberculosis, (Table 2). In order to use information as much as possible of negative samples, we implemented a K-means clustering center substitution method. In K-means clustering method, if there are K sequences of positive samples, then the negative samples will be divided into K categories. In K-mean clusters, the center of each category of negative samples will be taken as the characteristic of negative samples. The information regarding each non-essential ncRNA is reservation, and the number of negative samples obtained should be consistent with that of positive samples. For instance, B. subtilis 168 has 271 essential genes and 4175 non-essential genes. The 4175 non-essential genes would be divided in to 271 categories with K-means clustering method and this gives balanced sample for building the predicting model.

2.6 Data Dimension Reduction

The features collected through interval correlation and phase heterogeneous adjacency correlation of sequence composition created a characteristic variable data matrix with huge dimensions. Thus, this requires a computational tool to reduce the dimensions of the matrix.

The least absolute shrinkage and selection operator (Lasso regression) can be used for feature extraction, which has been successfully applied in the classification of cancer [39,40]. Thus, we have applied the Lasso regression method to reduce the dimensions. In the end, we obtained 231 variables to construct the essential ncRAN prediction model by SVM.

2.7 Machine Learning Algorithm

In this work, we use sequence features to characterize essential and non-essential ncRNAs. The essential and non-essential ncRNAs were computationally identified and analyzed as the characteristic variables of machine learning.

The first step was analysis of sequence homology and corresponding features are extracted. And then, the support vector machine (SVM) is used to build essential ncRNA models.

2.7.1 SVM Model Training

The 70% of the data were selected as training data. The remaining 30% was used as a test data. The prediction model was created by executing the SVMLIB training program---Svm-train.exe with the default parameters.

2.7.2 Prediction Model

We marked two different labels (1 and 0) to essential ncRNAs and non-essential ncRNAs sequences respectively and combined the 231 parameters into a training data. After optimizing the SVM parameters, C and Gamma by running grid.py program while the training, an accuracy of 89% was obtained under 5-fold cross-validation, which further indicated that the model could significantly separate the positive and negative sets. Finally, we obtained the classification model with the optimized parameters by executing the command: svm-train -b l -C xx - g xx training set file.

3  Results and Discussions

To verify the performance of the algorithm, we used the across-species validation to test the accuracy of the algorithm.

Accordingly, we used sensitivity (SN), specificity (SP), precision (precision) and accuracy (Acc) to measure the predictive performance of algorithm, respectively. These indicators are calculated as follows:

SN=TPTP+FN (11)

SP=TNTN+FP (12)

precision=TPTP+FP (13)

Acc=TP+TNTP+FP+TN+FN (14)

where TP indicates the number of the essential genes predicted as essential genes, FP indicates the number of the non-essential genes predicted as essential genes, TN indicates the number of the non-essential genes predicted as non-essential genes, FN indicates the number of the essential genes predicted as non-essential genes. The results of predictive performance are shown in Table 2.

Only the six datasets with ratio in the range from 0.1 to 0.7 are used for construction of essential ncRNA prediction model. The accuracy was found to be between 0.69 and 0.81 through cross-species cross-validation. It shows that the features we selected have good performance for prediction of essential ncRNA. And there is a potential application for discovering essential ncRNAs in bacteria for further analysis, such as antimicrobial target discovery, minimal genome construction and evolution analysis. However, the prediction accuracy fluctuated greatly for different bacteria. In order to further analyze the reasons for the different accuracy, we analyzed the relationship between the above four indicators (SN, SP, Precision and Acc) and the ratio of the number of essential ncRNA to the number of essential ncRNA in each organism (ratio), which is shown in Fig. 1.

Figure 1: The relationship between SN, SP, Precision, Acc and ratio

Fig. 1 shows that the accuracy of the algorithm is related to the ratio. It had been mentioned that accuracy of an algorithm would be unreliable due to differences in wet-lab experimental conditions. It was further found that when the ratio was about 0.35, the prediction effect was the best.

Acknowledgement: The authors are grateful to anonymous reviewers for their valuable suggestions and comments, which have led to the improvement of this paper.

Funding Statement: This study was jointly funded by the National Natural Science Foundation of China (61803112, 32160151), the Science and Technology Foundation of Guizhou Province (2019-2811).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.