SaPt-CNN-LSTM-AR-EA: a hybrid ensemble learning framework for time series-based multivariate DNA sequence prediction

Transformation of biological sequences

A set of DNA sequences with the length n, are denoted as $Seq=S_1{S}_2\cdots S_i\cdots S_n\left(S_i\in \left(A,C,T,G\right)\right)$. The following three methods are used to transform biological sequences into time series (Chou, Chen & Huang, 2022; Li, Dai & He, 2022)_ENREF_49.

Model framework and construction

BTS $x_t$ is defined as:

(4) $$x_t=Merge\left(funi\left(x_{t-T:t-1}\right),fmulti\left(x_{t-T:t-1}\right)\right),$$where $x_{t-T}$ indicates the sequence with a length of $\mathit{T}$ before time $\mathit{t}$; $Merge(\cdot )$ fuses single-sequence pattern $funi(\cdot )$ with multivariate-sequence pattern $fmulti(\cdot )$. In this article, by modelling $funi(\cdot )$ and $fmulti(\cdot )$, the ensemble learning module $Merge(\cdot )$ is constructed and then the hybrid ensemble learning model (SaPt-CNN-LSTM-AR-EA) is obtained.

The model framework is shown in Fig. 1.

SaPt-CNN-LSTM-AR-EA has two modules: single-sequence module and multi-sequence module. The input data of SaPt-CNN-LSTM-AR-EA are $N$ biological sequences $x^{\left(1\right)},\cdots x^{\left(N\right)}$. In the single-sequence module, ARFIMA method is used to obtain N single-sequence output vectors. The weighted fusion is carried out to obtain the output of single-sequence module $O_t^{AR}$. In the multi-sequence module, N biological sequences are combined together to form the input of multi-sequence vectors, which are input into LSTM for modelling. Finally, the weighted fusion of multivariate output sequences is carried out to obtain the output of multi-sequence module.

In the single-sequence module, the classical statistical method ARFIMA (Bhardwaj, Gadre & Chandrasekhar, 2020) is used to establish the model. The three parameters of ARFIMA model are optimized with the classical particle swarm evolution algorithm. The particle structure is illustrated as:

(5) $$y_{pso1}=f_{pso1}\left(p_i,d_i,q_i\right)$$

In the multi-sequence module, the temporal convolution of one-dimensional CNN is used to represent the dependence between multi-ple biological sequences. The one-dimensional convolution operation is illustrated as:

(6) $$C^{\left(k\right)}=W^{\left(k\right)}\times X_{t-T:t-1},$$where $C^{\left(k\right)}$ is the convolution result and $W^{\left(k\right)}$ is the k-th convolution kernel.

After the execution of the one-dimensional convolution feature extraction, the results are used as the input of the model LSTM (Karim et al., 2019; Singaravel, Suykens & Geyer, 2018) to obtain the output of the multi-sequence pattern $O_t^{\mathrm{L}\mathrm{S}\mathrm{T}\mathrm{M}}$. The three parameters of LSTM model are also optimized by particle swarm optimization algorithm (Bi et al., 2023; Dias et al., 2023; Li, Liu & Wu, 2022). The particle structure is illustrated as:

(7) $$y_{pso2}=f_{pso2}\left(N_{tr},N_{hn},D_{tr}\right).$$

Finally, the output of the frame A is obtained by fusing the output of single sequence and multivariate sequence. The procedure is illustrated as:

(8) $$C_t^{\left(k\right)}=W_t^{AR}\otimes O_t^{AR}+W_t^{LSTM}\otimes O_t^{LSTM}.$$

Evaluation and testing

Experimental data

Source of experimental data

The experimental data were DNA sequences of six viruses (Table 1) downloaded from NCBI (https://www.ncbi.nlm.nih.gov).

Table 1:

Source of experimental data.

Label	Source	Accession	Length (bp)
A	Human adenovirus C	NC_001405	35,937
B	Dubowvirus MR25	NC_010808	44,342
C	Infectious bronchitis virus	NC_048213	27,464
D	Phietavirus MR11	NC_010147	43,011
E	Abalone shriveling syndrome-associated virus	NC_011646	34,952
F	Clostridium phage phiCD505	NC_028764	49,316

DOI: 10.7717/peerj.16192/table-1

Biological sequences are transformed into time series so as to obtain the corresponding BTS. According to the length of each BTS, the sequence is divided into several data to form a BTS database, which is used for model training, verification, and testing. The procedure is illustrated in Fig. 2.

After sequence division, each BTS consists of several groups of subsequences. To improve the generalization ability of the model, each subsequence is divided into three subsets: training set (70%), validation set (15%), and test set (15%), as shown in Table 2.

Table 2:

Experimental data distribution.

Label	Training set	Validation set	Testing set	Data points
A	420	89	89	598
B	519	110	110	739
C	321	68	68	457
D	502	107	107	716
E	408	87	87	582
F	575	123	123	821

DOI: 10.7717/peerj.16192/table-2

Stationary analysis

Taking the BTS generated with sequences A, C, and E as an example, the BTS corresponding to the first 3,000 bp is shown in Fig. 3.

In Fig. 3, the time series curves obtained with sequences A, C, and E are characterized by large fluctuations, unequal amplitudes, and unequal position intervals, indicating that BTS has typical non-stationary characteristics. Other series have similar non-stationary characteristics.

Results

The operating environment is Windows 10 64-bit OS (16 GB of memory and Intel (R) Core™ i7-12700F processor). The deep learning framework was constructed with Matlab2020a. The transformed and preprocessed 6 BTSs were used for model testing and prediction (see the architecture of the model and corresponding parameters is in the Supporting Information). Finally, the output results of the model were reversely normalized to make the model more explanatory.

Firstly, SaPt-CNN-LSTM-AR-EA was trained with the training set. Secondly, various parameters of the model were adjusted to minimize the training error in the training process (see Supporting Information). The relationship between the predicted value and the actual value of each data sample in the training set is shown in Fig. 4.

The closer the predicted data points are to the experimental line, the smaller the prediction error is. As shown in Fig. 4, the predicted data points of the six biological sequences are basically distributed on the experimental line, indicating that the training effect of SaPt-CNN-LSTM-AR-EA was good on each sequence data set and that the model had been fully trained. The validation set was used to verify the reliability of the model. The parameters of the model were fine-tuned in the validation process so as to reduce the output error of the model. The distribution relationship between the predicted value and the actual value of the model in the validation set is shown in Fig. 5.

In the validation set, the data points predicted by the model are basically distributed near the straight line, indicating that the predicted value was in good agreement with the experimental value. The predicted value of SaPt-CNN-LSTM-AR-EA was more consistent with the experimental value, indicating that the trained model had the reliable and accurate prediction ability and could be used to predict BTS.

After training and validation, the prediction experiment was performed with the trained model. With the test set, the model was tested against each biological sequence. The test data of each biological sequence were used for model testing. The prediction results and relevant data statistics of the model are shown in Fig. 6.

The distributions of predicted values and experimental values and the mean distribution of errors are shown in Fig. 6A. The error statistics and error distribution of the model in each biological sequence test are respectively shown in Figs. 6B and 6C.

The predicted values of each sequence were basically consistent with the actual values except the results of some loci. In the error analysis results, most of the error points were distributed around the value of 0 and the number of sequences with large errors accounted for a small proportion. The predicted values were highly consistent with the experimental values in six biological sequences, indicating that the model had good prediction performance. The performance indexes of the model in the training set, validation set, and test set are shown in Table 3.

The values of performance indexes fully reflect the comprehensive performance of the model in the three data sets. The error in the validation set was small and the correlation was high. The prediction ability of training and validation sets was better than that of the test set. From the perspective of prediction mechanism, model training and validation aims to reduce the output error, so the correlation coefficient of the training and validation sets should be larger than that of the new samples in the test set.

Comparison with other benchmark models

In order to verify the performance of the model, several models with better performances in time series prediction were selected as benchmark comparison models in this article. The theories and parameters of all the models are shown in Table 4.

One hundred samples were randomly selected from six BTS data sets to form the test data set of the comparison model. Each benchmark model was used to predict each sample in the test data set. The prediction results of each model on Sequence A are shown in Fig. 7. The results on other sequences are similar to those of Sequence A.

The correlation between the predicted and experimental values is shown in Figs. 7A and 7B. The prediction results of the SaPt-CNN-LSTM-AR-EA model were closer to the experimental line, indicating that the predicted results were more consistent with experimental results. The statistical distribution of predicted values shown in Figs. 7C–7E indicated that the SaPt-CNN-LSTM-AR-EA model performed better than other models. The errors of each model are shown in Fig. 8.

In the error curve (Fig. 8A), the error bar of the SaPt-CNN-LSTM-AR-EA model was closest to the origin and the prediction error was smaller, indicating that the model had the highest prediction accuracy. According to the error statistics (Fig. 8B), most of the error points of the SaPt-CNN-LSTM-AR-EA model were distributed between 1.4 and 2.0 and the average error was also the smallest, indicating that the accuracy of the model was relatively high. The correlation coefficient and calculation time of each benchmark model are shown in Fig. 9.

According to Fig. 9A, the correlation curve of SaPt-CNN-LSTM-AR-EA model is at the top of the coordinate and the coordinate value is closest to 1. In addition, the statistics of correlation data points also showed that the predicted values of the model in this article had the most significant correlation with experimental values. The BI-ARFIMA model had the shortest computation time and CONV-LSTM model had the longest computation time (Fig. 9B). The SAPT-CNN-LSTM-AR-EA model also had the acceptable computation time. Table 5 shows the performance statistics of each model. The SaPt-CNN-LSTM-AR-EA model performed better in terms of both prediction accuracy and correlation and its computation time was also acceptable.

The SaPt-CNN-LSTM-AR-EA model had obvious advantages over other models in terms of prediction accuracy and correlation due to the following factors. Firstly, the characteristics of multivariate sequences were fully utilized in the model. Secondly, adaptive pre-training mechanism improved the training performance. Thirdly, the advantages of CNN-LSTM in feature extraction had been fully utilized. The SaPt-CNN-LSTM-AR-EA model belonged to a multi-layer cyclic deep learning framework and fused ARFIMA, so its computation time was long. Compared with CTS-LSTM and CONV-LSTM models, SaPt-CNN-LSTM-AR-EA had a slight advantage in computation time because it adopted one-dimensional convolution operation.

Discussion of different evaluation metrics

Accuracy growth rate

In this article, the accuracy growth rate of each benchmark model was calculated in order to verify the prediction accuracy of the SaPt-CNN-LSTM-AR-EA model on the test dataset. The accuracy growth rate statistics of the SaPt-CNN-LSTM-AR-EA compared with other five benchmark models are shown in Fig. 10.

Compared with the five benchmark models (BI-ARFIMA, ARFIMA-LSTM, EA-LSTM, CTS-LSTM, and Conv-LSTM), the SaPt-CNN-LSTM-AR-EA model had different MAPE growth rates in six different biological sequences. The highest accuracy was increased by nearly 50% and the average accuracy was increased by about 30%. The accuracy improvement was obvious.

The output of the SaPt-CNN-LSTM-AR-EA model had obvious advantages in accuracy due to the comprehensive results of various algorithms. In addition, the performance of the model was significantly improved.

VAR

The proposed model is to improve the efficiency, accuracy, and stability of prediction. We used variance of residuals (VAR) to test the stability of the model. The VAR change curves of each model are shown in Fig. 11. The VAR curve of the SaPt-CNN-LSTM-AR-EA model was close to the abscissa, indicating that the model was more stable. In addition, BI-ARFIMA and ARFEIMA-LSTM had similar stability and CTS-LSTM was more stable than EA-LSTM.

The results of VAR stability test indicated that the stability of the SaPt-CNN-LSTM-AR-EA model was basically the same among all the BTS because the fusion of multivariate sequences in the model avoided the defect of single sequence output error and improved the generalization ability. The results also proved the good scalability of the model.

Discussion of ablation study

To validate the contribution of these modules to the framework, we performed ablation study by respectively removing one module from the four modules in four times. Then, four models were obtained: CNN-LSTM-AR-EA, SaPt-LSTM-AR-EA, SaPt-CNN-LSTM-EA, and SaPt-CNN-LSTM-AR. The performance indexes of these models are shown in Table 7.

After any module was removed from SaPt, CNN, AR, and EA, the error was significantly larger and the correlation decreased, indicating that each module contributed to the SaPt CNN LSTM AR-EA model. The SaPt-LSTM-AR-EA model without the CNN module had the worst experimental results, indicating that the CNN module had the greatest contribution. The contributions of the other three modules ranked in the following order: AR > SaPt > EA. The contribution rate is defined as the ratio of the contribution of one module to the total contribution (assuming that the total contribution of the four modules is 100%). Figure 12 shows the accuracy contribution and correlation contribution of each module.

The contributions of the four modules ranked in the following order: CNN > AR > SaPt > EA and all modules have a contribution rate of over 17% (Figs. 12A and 12B). The ablation study indicated that the four modules were critical and indispensable and jointly contributed to the SaPt-CNN-LSTM-AR-EA model.