6mAPred-MSFF: A Deep Learning Model for Predicting DNA N6-Methyladenine Sites across Species Based on a Multi-Scale Feature Fusion Mechanism

Abstract

DNA methylation is one of the most extensive epigenetic modifications. DNA N6-methyladenine (6mA) plays a key role in many biology regulation processes. An accurate and reliable genome-wide identification of 6mA sites is crucial for systematically understanding its biological functions. Some machine learning tools can identify 6mA sites, but their limited prediction accuracy and lack of robustness limit their usability in epigenetic studies, which implies the great need of developing new computational methods for this problem. In this paper, we developed a novel computational predictor, namely the 6mAPred-MSFF, which is a deep learning framework based on a multi-scale feature fusion mechanism to identify 6mA sites across different species. In the predictor, we integrate the inverted residual block and multi-scale attention mechanism to build lightweight and deep neural networks. As compared to existing predictors using traditional machine learning, our deep learning framework needs no prior knowledge of 6mA or manually crafted sequence features and sufficiently capture better characteristics of 6mA sites. By benchmarking comparison, our deep learning method outperforms the state-of-the-art methods on the 5-fold cross-validation test on the seven datasets of six species, demonstrating that the proposed 6mAPred-MSFF is more effective and generic. Specifically, our proposed 6mAPred-MSFF gives the sensitivity and specificity of the 5-fold cross-validation on the 6mA-rice-Lv dataset as 97.88% and 94.64%, respectively. Our model trained with the rice data predicts well the 6mA sites of other five species: Arabidopsis thaliana, Fragaria vesca, Rosa chinensis, Homo sapiens, and Drosophila melanogaster with a prediction accuracy 98.51%, 93.02%, and 91.53%, respectively. Moreover, via experimental comparison, we explored performance impact by training and testing our proposed model under different encoding schemes and feature descriptors.

1. Introduction

Epigenetics refers to the reversible and heritable changes in gene function when there is no change in the nuclear DNA sequence [1]. DNA methylation modifications play important roles in epigenetic regulation of gene expression without altering the sequence, and it is widely distributed in the genome of different species [2]. It can be divided into three categories according to the position of methylation modification: N6-methyladenine (6mA), 5-Methylcytosine (5mC), and N4-methylcytosine (4mC) [3,4]. DNA methylations at the 5th position of the pyrimidine ring of cytosine 5-Methylcytosine(5mC) and at the 6th position of the purine ring of adenine N6-methyladenine (6mA) are the most common DNA modifications in eukaryotes and prokaryotes, respectively [5]. Previous studies have shown that DNA N6-methyladenine (6mA) is associated with germ cell differentiation, stress response, embryonic development, nervous system, and other processes [6,7,8,9]. N6-Methyladenine (6mA) DNA methylation has recently been implicated as a potential new epigenetic marker in eukaryotes, including the Arabidopsis thaliana, Rice, Drosophila melanogaster, and so on [5,8,10]. Zhang et al. reveal that 6mA is a conserved DNA modification that is positively associated with gene expression and contributes to key agronomic traits in plants [11]. Some studies have found that N6-methyladenine DNA modification is also widely distributed in the Human Genome and plays important biological functions. Xiao et al. demonstrate that 6mA DNA modification is extensively present in human cells, and the decrease of genomic DNA 6mA promotes human tumorigenesis [12]. Zhou et al. proposed a 6mA DNA modification area as a new mechanism for the epigenetic regulation of stem cell differentiation [13]. Xie et al. report that N6-methyladenine DNA modifications are enriched in human glioblastoma, and targeting regulators of this modification can inhibit cancer growth by altering heterochromatin landscapes and downregulating oncogenic programs [14]. Therefore, how to quickly and accurately detect the DNA 6mA modification sites is also an important research topic in these epigenetic studies.

Due to the rapid development of high-throughput sequence technology, various experimental techniques have been proposed to detect DNA 6mA modifications and study protein function. Pormraning et al. developed a protocol using bisulfite sequencing and a methyl-DNA immunoprecipitation technique to analyze genome-wide DNA methylation in eukaryotes [15]. Krais et al. reported a fast and sensitive method for the quantification of global adenine methylation in DNA, using laser-induced fluorescence and capillary electrophoresis [16]. Flusberg et al. proposed the single-molecule real-time sequencing (SMRT) technology to detect the 4mC and 6mA sites from the whole genome [17]. Greer et al. used ultra-high performance liquid chromatography coupled with mass spectrometry technique to access DNA 6mA levels in Caenorhabditis elegans [18]. By performing mass spectrometry analysis and 6mA immunoprecipitation followed by sequencing (IP-seq), Zhou and his colleagues obtained the 6mA profile of the rice genome [10].

Although experimental methods indeed yielded encouraging results, the technology cannot detect m6A sites from the whole genome and the cost of the technique is high. Therefore, it is necessary to develop a computational model that can efficiently and accurately predict and identify 6mA sites. Recent studies focus more on the recognition of 6mA sites using machine learning [19], which is capable of predicting 6mA sites based on genome sequences, without any prior experimental knowledge. Chen et al. developed the first ML-based method named i6mA-Pred for identifying DNA 6mA sites and provided a benchmark 6mA dataset containing 880 6mA sites and 880 non-6mA sites in the rice genome. Their method used a Support Vector Machine (SVM) classifier based on chemical features of nucleotides and position-specific nucleotide frequencies [20]. The i6mA-Pred shows a good classification performance in rice 6mA data. However, the association information among nucleotides near 6mA sites is ignored. Pian et al. proposed a new classification method called MM-6mAPred based on a Markov model which makes use of the transition probability between adjacent nucleotides to identify 6mA sites [21]. Pian et al. built and evaluated their MM-6mAPred based on the 6mA-rice-Chen benchmark dataset. Their results show that MM-6mAPred outperformed i6mA-Pred in prediction of 6mA sites. Basith et al. developed a novel computational predictor, called SDM6A, which explores various features and five encoding methods to identify the DNA 6mA sites [22]. Basith et al. also trained and evaluated their SDM6A based on the 6mA-rice-Chen benchmark dataset, and they found that SDM6A outperformed i6mA-Pred on the 6mA-rice-Chen benchmark dataset. The above three prediction models are trained on the 6mA-rice-Chen benchmark dataset including only the 880 rice 6mA sites and 880 non-6mA sites. Even though the above methods have improved the performance for identifying 6mA sites, too few data sets have been adopted to fully reflect the whole genome and to build robust models. Lv et al. developed a machine learning method of predicting 6mA sites named iDNA6mA-rice which was trained and evaluated on the dataset 6mA-rice-Lv containing 154,000 6mA sites and 154,000 non-6mA sites in the rice genome [23]. iDNA6mA-rice utilized Random Forest to perform the classification for identifying 6mA sites after using mono-nucleotide binary encoding to formulate positive and negative samples.

In recent years, deep learning has not only been developed as a new research direction in machine learning, but has also made a lot of achievements in data mining, machine translation, natural language processing, and other related fields. In the field of computational biology [24,25,26,27,28], deep learning has been widely applied [29,30,31], especially in solving the problems of genome sequence-based by convolutional neural networks (CNN). Tahir et al. proposed an intelligent computational model called iDNA 6mA (5-step rule) which extracts the key features from DNA input sequences via the convolution neural network (CNN) to identify 6mA sites in the rice genome [32]. Yu et al. developed a simple and lightweight deep learning model named SNNRice6mA to identify DNA 6mA sites in the rice genome and showed its advantages over other methods [33]. Li et al. developed a deep learning framework named Deep6mA to identify DNA 6mA sites. Deep6mA, composed of a CNN and a bidirectional LSTM (BLSTM) module, is shown to have a better performance than other methods above on 6mA prediction [34]. Although the above methods have made great progress, their performance is still not satisfactory. Furthermore, most of the methods are designed for one specific species such as the rice genome. Although some methods provide the across species validation test, the results are not as good as the original species-specific model.

In this paper, we proposed a novel deep learning method called 6mAPred-MSFF to identify DNA 6mA sites. In this predictor, we establish the model by integrating the inverted residual block, the multi-scale channel attention module (MS-CAM), Bi-directional Long Short Term Memory Network (Bi-LSTM), and attentional feature fusion (AFF). To capture more feature information, we use the inverted residual block to expand the input features to a high dimension and filter them with a lightweight depthwise convolution. After that, we project them back to a low-dimensional representation as the input of the multi-scale channel attention module (MS-CAM). Bi-LSTM is designed to capture the long dependent information in sequences. The features generated from the Bi-LSTM are fed into the AFF to be fused. In the experimental part, we evaluate our predictor and other 6mA predictors with 5-fold cross validation on seven benchmark datasets. We also evaluate these 6mA predictors trained by the 6mA-rice-Lv dataset to predict the 6mA sites of other species. Finally, we explore the performance impact under the different encoding methods and feature descriptors.

2. Materials and Methods

2.1. Dataset Collection

Previous studies have demonstrated that a stringent dataset is essential for building a robust predictive model. We collected seven 6mA site datasets of six species including Rice, Arabidopsis thaliana (A. thaliana), Rosa chinensis (R. chinensis), Fragaria vesca (F. vesca), Homo sapiens (H. sapiens) and Drosophila melanogaster (D. melanogaster). There are two benchmark datasets 6mA-rice-Lv and 6mA-rice-Chen for Rice. The 6mA-rice-Lv dataset include 15,400 positive samples and 15,400 negative samples which are obtained from NCBI Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/genome/10 (accessed on 20 August 2021)) following the steps proposed by Lv et al. [23]. The 6mA-rice-Chen benchmark dataset includes the 880 rice 6mA sites and 880 non-6mA sites which are obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/ (accessed on 20 August 2021)) under the accession number GSE103145 following the steps proposed by Chen et al. [20].

To demonstrate that our model has the ability to predict the 6mA sites of other species and 6mA shares similar patterns across different species, we collected DNA 6mA sequences of Arabidopsis thaliana and Drosophila melanogaster from the MethSMRT database (http://sysbio.gzzoc.com/methsmrt/ (accessed on 20 August 2021)) [35] and also obtained the DNA 6mA sequences of Fragaria vesca and Rosa chinensis from the MDR database [36]. Preliminary trials indicated that, when the length of the segments is 41 bp with the 6mA in the center, the highest predictive results could be obtained. Thus, the sequences of all positive samples are 41 bp. In order to construct a high-quality benchmark dataset, the following two steps were performed. Firstly, according to the the Methylome Analysis Technical Note [35], the sequences with modQV no less than 30 are left for the subsequent analysis. It should be noted that, in order to obtain statistically significant results, if the raw data are too small, this step was ignored to get more samples. Secondly, to avoid redundancy and reduce homology bias, sequences with more than 80% sequence similarity were removed using the CD-HIT program [37]. After the above two steps, the objective and strict positive datasets for the above species were obtained. Negative samples for Arabidopsis thaliana, Drosophila melanogaster, Fragaria vesca, and Rosa chinensis were also collected. These samples are 41 nt long sequences with Adenine in the center and were proved to be unmethylated by experiments. For convenience, we used the Homo sapiens dataset (http://lin-group.cn/server/iDNA-MS/download.html (accessed on 20 August 2021)) [38] directly, including 18,335 positive samples and 18,335 negative samples. The details of the datasets are presented in Table 1.

2.2. Architecture of 6mAPred-MSFF

In this work, we proposed a new predictive model based on a multi-scale feature fusion mechanism. An overview of our model architecture is illustrated in Figure 1. There are four modules in our model including the sequence embedding module, the feature extraction module, feature fusion module, and prediction module. First, in the sequence embedding module, we use four different encoding schemes (1-gram, NAC, 2-grams, DNC) for the input of the embedding layer. Each nucleotide is represented as a dense vector of float point values. As a result, the genomic sequences can be represented by four feature matrices, respectively. We concatenate the feature matrices of the 1-gram and NAC and the feature matrices of the 2-grams and DNC matrices, respectively. Second, the resulting two feature matrices are fed into the feature extraction module which consists of the inverted residual block and the multi-scale channel attention mechanism module (MS-CAM). To extract more information from the global features and filter the features as the source of the MS-CAM, we use the inverted residual block including an expansion point-wise layer, a convolution layer, and a projection point-wise layer. The MS-CAM extract and combine the global and local features of the two feature matrices, respectively. Before fusing the two feature matrices, we use the bidirectional LSTM layer to learn the information about the long-distance dependence. Afterward, in the feature fusion module, we combine the features extracted from the bidirectional LSTM layer by adding operation and feed them into the MS-CAM module to calculate the fusion weight of the two feature matrices. We multiply the features matrices with the corresponding fusion weights and combine the results by adding operation. In this way, we obtain the features by aggregating the local and global contexts of the four different encoding schemes. Third, in the prediction module, the outputs of the features are passed through dropout layers for the prevention of overfitting during the training. After the dropout layer, to unstack the output, a flatten function was used to squash the feature vectors from the previous layers. Right after a flattened layer, a fully connected (FC) dense layer used as Dense(n), with the number of n neurons set as 32. At last, an FC layer was applied and used a sigmoid function for the binary classification. Sigmoid is used to squeeze the values between the range of 0 and 1 to represent the probability of having 6mA and non-6mA sites.

The 6mAPred-MSFF is trained by the Adam optimizer and the binary cross-entropy loss function. Via giving a prediction score, the 6mAPred-MSFF can determine whether the 6mA sites are detectable in DNA sequences. It is detectable if the prediction score is >0.5; otherwise, it is not. The model is implemented using Keras 2.4.3. Below, each of ourmodules will be described in detail.

2.2.1. Sequence Embedding Module

DNA sequences are consisting of four nucleotides: “A” (adenine), “G” (guanine), “C” (cytosine), and “T” (thymine). Undetermined bases are annotated as “N”. The DNA sequences can be expressed as $S=s_1,s_2,\dots s_L$ where $s_i$ is the i-th word and the L is the length of the sequence. We use four encoding methods to define ‘words’ in DNA sequences including 1-gram, NAC, 2-grams, and DNC.

1-gram encoding method. The n-grams are the set of all possible subsequences of nucleobases [39]. One-gram encoding sets the n-grams number n to 1. As a result, the 1-gram nucleobase sequences include ‘A’, ‘T’, ‘C’, ‘G’. We define a dictionary that maps the nucleotides to numbers, such as A (adenine) to 1, T (thymine) is 2, C (cytosine) to 3, and G (guanine) to 4. In this way, we map a DNA sequence into a real number vector of length 41, which is denoted as $V_{i,1-gram}$:

(1)$\begin{array}{cc}{V}_{i,1-gram}={[v_1{v}_2{v}_3\dots v_{39}{v}_{40}{v}_{41}]}_i,& v\in \left[1,2,3,4\right]\end{array}$

NAC encoding method. The Nucleic Acid Composition (NAC) encoding calculates the frequency of each nucleic acid type in a nucleotide sequence. The frequencies of all four natural nucleic acids (i.e., “ATCG”) can be calculated as:

(2)$\begin{array}{cc}{f}_{NAC}\left(t\right)=\frac{N\left(t\right)}{N},& t\in \left\{A,T,C,G\right\}\end{array}$

for a DNA sequence S of length L, and it can be converted into a vector of length 4, which is denoted as $V_{i,NAC}$:

(3)$V_{i,NAC}={[f_{NAC}\left(A\right),f_{NAC}\left(T\right),f_{NAC}\left(C\right),f_{NAC}\left(G\right)]}_i$

2-gram encoding method. A 2-gram encoding sets the n-grams number n to 2. We split the DNA sequences into overlapping n-gram nucleobases. The total possible number is 16, since there are four types of nucleobases. For example, we split a DNA sequence into overlapping 2-gram nucleobase sequences as follows: GTTGT…CTT $\to$ ‘GT’, ‘TT’, ‘TG’, ‘GT’, …, ‘CT’, ‘TT’. We define a dictionary that maps the 2-gram nucleobase sequences to the 1–16 numbers set. In this way, we map a DNA sequence into a real number vector of length 40, which donates $V_{i,2-gram}$:

(4)$\begin{array}{cc}{V}_{i,2-gram}={[v_1{v}_2{v}_3\dots v_{39}{v}_{40}]}_i,& v\in \left[1,2,3,\dots 14,15,16\right]\end{array}$

DNC encoding method. The Di-Nucleotide Composition (DNC) encoding calculates the frequency of every two nucleic acid types in a nucleotide sequence which gives 16 descriptors. It is defined as:

(5)$\begin{array}{cc}{f}_{DNC}\left(r,s\right)=\frac{N_{rs}}{N-1},& r,s\in \left\{A,T,C,G\right\}\end{array}$

where $N_{rs}$ is the number of di-nucleotide represented by nucleic acid types r and s. For a DNA sequence S of length L, it can be converted into a vector of length 16, which is donated as $V_{i,DNC}$:

(6)$V_{i,DNC}={[f_{DNC}\left(AA\right),f_{DNC}\left(AC\right),\dots ,f_{DNC}\left(TC\right),f_{DNC}\left(TG\right)]}_i$

Next, the vectors $V_{i,1-gram}$, $V_{i,NAC}$, $V_{i,2-gram},$ and $V_{i,DNC}$ are fed to an embedding layer, transforming into learnable embedding vectors. To establish the relationship between $V_{i,1-gram}$ and $V_{i,NAC}$, we concatenate the two embedding vectors of $V_{i,1-gram}$ and $V_{i,NAC }$ to a new feature vector that is denoted as $X_i$:

(7)$X_i=Concatenate\left(Embedding\left(V_{i,1-gram}\right),Embedding\left(V_{i,NAC}\right)\right)$

We do the same operation on the embedding vectors of $V_{i,2-gram}$ and $V_{i,DNC}$, generating a new feature vector which is denoted as $Y_i$:

(8)$Y_i=Concatenate\left(Embedding\left(V_{i,2-gram}\right),Embedding\left(V_{i,DNC}\right)\right)$

To extract more feature information, the vectors $X_i$ and $Y_i$ are fed into the feature extraction module.

2.2.2. Feature Extraction Module

The inverted residual block. To reduce the computation of the network, we always design the network with a low-dimensional space. However, the network can not capture enough information in low-dimensional space. To solve these problems, the inverted residual with linear bottleneck was proposed by Sandler et al. [40] in the MobileNetV2. This module takes an input as a low-dimensional compressed representation, which is first expanded to high dimensions and filtered with a lightweight depthwise convolution. Features are subsequently projected back to a low-dimensional representation with a linear convolution which can prevent the destruction of information in low-dimensional space. Our proposed method follows the idea to build an inverted residual block that consists of an expansion layer, a depthwise layer, and a projection layer. The expansion layer and the projection layer are the 1 × 1 convolution layer, called the pointwise convolution with a different expansion ratio. The depthwise layer is a depthwise separable convolution layer that is a key building block for many efficient neural network architectures [40,41,42,43,44]. The basic idea is to replace a full convolutional operator with a factorized version that splits convolution into two separate layers. The first layer is called a depthwise convolution, and it performs lightweight filtering by applying a single convolutional filter per input channel. The second layer is a point-wise convolution, which is responsible for building new features through computing linear combinations of the input channels. The inverted residual block operator can be expressed as a composition of three operators which is denoted as $F_{IRB}\left(X\right)$:

(9)$F_{IRB}\left(X\right)=\left[A\circ N\circ B\right]X$

where A is the expansion layer that is a linear transformation: $A:R^{s\times s\times k}\to R^{s\times s\times n}$, N is the depthwise layer which contains several non-linear per-channel transformation: $N:R^{s\times s\times n}\to \dots \to R^{s^{\prime }\times s^{\prime }\times n}$, and B is the projection layer which is again a linear transformation to the output domain: $B:R^{s\times s\times n}\to R^{s^{\prime }\times s^{\prime }\times k^{\prime }}$. For our inverted residual block $F_{IRB}\left(X\right)$, A and B are the linear transformation without any nonlinear activations. We use the ELU activation function in the operator $N=ELU\circ dwise\circ ELU$. which accelerates the learning rate of the network. The ELU function is shown in the following equation:

(10)$f\left(x\right)=\{\begin{array}{cc}a\left(e^x-1\right)& if\left(x<0\right)\\ x& if\left(0\le x\right)\end{array}$

Thus, we can obtain $k^{\prime }$ feature maps containing the 6mA site features captured by the inverted residual block.

Multi-scale attention mechanism block. The attention mechanism in deep learning, which mimics the human visual attention mechanism [45,46], is originally developed on a global scale. For example, the matrix multiplication in self-attention draws global dependencies of each word in a sentence [47] or each pixel in an image [48,49,50,51]. The Squeeze-and-Excitation Networks (SENet) squeeze global spatial information into a channel descriptor to capture channel-wise dependencies [52]. However, merely aggregating contextual information on a global scale is too biased and weakens the features of small objects. To solve the problem, the multi-scale channel attention module (MS-CAM) proposed by Yimian Dai et al. [53] aggregates global and local channel contexts given an intermediate feature $X\in R^{H\times W\times C}$ with C channels and feature maps of size $H\times W$. The attention channel weights $M\left(X\right)\in R^C$ can be computed as:

(11)$M\left(X\right)=L\left(X\right)\oplus G\left(X\right)$

where $L\left(X\right)$ denotes the local channel context, and $G\left(X\right)$ denotes the global channel context. The stacked architecture of MS-CAM is shown in Figure 1. $L\left(X\right)\in R^{C\times H\times W}$ is computed as follows:

(12)$L\left(X\right)=\beta \left(PWCon{v}_2\left(\beta \left(PWCon{v}_1\left(X\right)\right)\right)\right)$

where $\beta$ denotes the Batch Normalization (BN) [54]. The pointwise convolution (PWConv) is used as the local channel context aggregator, which only exploits pointwise channel interactions for each spatial position. The expansion ratios of $PWCon{v}_1$ and $PWCon{v}_2$ are $r$ and $1/r$, respectively. Different from the original MS-CAM, we remove the nonlinear activation function between the convolution layer to prevent the destruction of feature information. $G\left(X\right)\in R^C$ is computed as follows:

(13)$G\left(X\right)=\beta \left(PWCon{v}_2\left(\beta \left(PWCon{v}_1\left(g\left(X\right)\right)\right)\right)\right)$

where $g\left(X\right)$ is global average pooling (GAP) and computed as follows:

(14)$g\left(X\right)=\frac{1}{H\times W}{\displaystyle \sum }_{i=1}^H{\displaystyle \sum }_{j=1}^W{X}_{\left[:,i,j\right]}$

The expansion ratios of $PWCon{v}_1$ and $PWCon{v}_2$ are $r$ and $1/r$, respectively. $\beta$ denotes the Batch Normalization (BN). The refined $X^{\prime }\in R^{H\times W\times C}$ can be obtained as follows:

(15)$X^{\prime }=X\otimes \sigma \left(M\left(X\right)\right)=X\otimes \sigma \left(L\left(X\right)\oplus G\left(X\right)\right)$

where $\sigma$ is the sigmoid function. $\oplus$ denotes the broadcasting addition and $\otimes$ denotes the element-wise multiplication.

Bi-directional Long Short Term Memory Network. Recurrent Neural Network (RNN) is a powerful neural network for processing sequential data [55]. As a special type of RNN, Long Short Term Memory Network (LSTM) is not only designed to capture the long dependent information in sequence but also overcome the training difficulty of RNN due to gradient explosion and disappearance [56]; thus, it is the most widely used RNN in real applications. In LSTM, a storage mechanism is used to replace the hidden function in traditional RNN, with a purpose to enhance the learning ability of LSTM for long-distance dependency. Bi-directional LSTM (BLSTM), compared with unidirectional LSTM, better captures the information of sequence context. The LSTM components are updated by the following formulations:

(16)$\begin{array}{c}\begin{array}{c}{f}_t=\sigma \left(W_f·\left[h_{t-1},x_t\right]+b_f\right)\\ i_t=\sigma \left(W_i·\left[h_{t-1},x_t\right]+b_i\right)\\ {C^{\prime }}_t=\tanh (W_C·\left[h_{t-1},x_t\right]+b_C)\end{array}\\ \begin{array}{c}{C}_t=i_t{C}^{\prime }+f_t{C}_{t-1}\\ o_t=\sigma \left(W_x{h}_{t-1}+W_h{x}_t+W_c{C}_t+b_0\right)\\ h_t=o_t\tanh (C_t)\end{array}\end{array}$

where $i_t$,$f_t$, and $o_t$ represent the input, forget, and output gate, respectively; ${C^{\prime }}_t$ is temp value for calculating $C_t$; $t$ denotes the recurrent time step; $W_f$, $W_i$, $W_C$ and $b$ are the weights of each equation. $h_t$ is the output of LSTM at the $t$ time step. The bi-directional LSTM consists of two direction LSTM networks. Therefore, the $h_t$ is computed as follows:

(17)$h_t=\left[{\overrightarrow{h}}_t\oplus {\overleftarrow{h}}_t\right]$

2.2.3. Feature Fusion Module

Feature fusion, the combination of features from different layers or branches, is an omnipresent part of modern network architectures. It is often implemented via simple operations, such as summation or concatenation, but this might not be the best choice. We use a uniform and general scheme, namely attentional feature fusion (AFF) which is proposed by Yimian Dai et al. [53]. The architecture of AFF is shown in Figure 1. Based on the multi-scale channel attention module $M\left(X\right)$, the AFF can be expressed as follows:

(18)$Z=M\left(X\oplus Y\right)\otimes X+\left(1-M\left(X\oplus Y\right)\right)\otimes Y$

where $Z\in R^{H\times W\times C}$ is the fused feature, and the element-wise summation of $X\oplus Y$ is used as the initial integration for MS-CAM. $M\left(X\oplus Y\right)$ and $1-M\left(X\oplus Y\right)$ are the fusion weights for $X$ and $Y$, respectively. Please note that the fusion weights consist of real numbers between 0 and 1, which enable the network to conduct a soft selection or weighted averaging between $X$ and $Y$.

2.2.4. Prediction Module

The fused feature $Z$ generated from the AFF module is fed into the dropout layer which prevents overfitting and can improve the reliability and robustness of the model by discarding some intermediate features. We use a flatten function to integrate the feature vectors that are generated from the dropout layer. Afterwards, the flattened feature vectors are fed to the fully connected (FC) layer that has 32 hidden units. The activation function used in the FC layer is the ELU activation function. Finally, a sigmoid activation function is used to combine the outputs from the FC layer to make the final decision.

6mAPred-MSFF is trained by the Adam optimizer [57], which is computationally efficient, has little memory requirements, is invariant to diagonal rescaling of gradients, and is well suited for problems that are large in terms of data/parameters. For binary classification tasks, the loss function for training 6mAPred-MSFF is set as the binary cross-entropy for the measuring the difference between the target and the predicted output:

(19)$L\left(w\right)=-{\displaystyle \sum }_{i=1}^N{y}_i\log ({y^{\prime }}_i)+\left(1-y_i\right)\log (1-{y^{\prime }}_i)+\alpha \Vert w_2\Vert$

where $y_i$ is the true label, ${y^{\prime }}_i$ is the corresponding predicted value from 6mAPred-MSFF, and $\alpha \Vert w_2\Vert$ is a regularization term to avoid overfitting.

2.2.5. Evaluation Metrics

In our experiment, we used the following four indicators to evaluate the predictive performance of our proposed model, including Accuracy (ACC), Sensitivity (SN), Specificity (SP), and Mathew’s Correlation Coefficient (MCC). They are the four commonly used indicators for classifier performance evaluation in other bioinformatics fields [58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87]. They are formulated as follows:

$\mathrm{SN}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$

$\mathrm{SP}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}$

$\mathrm{Precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$

$\mathrm{ACC}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$

$\mathrm{MCC}=\frac{\mathrm{TP}\times \mathrm{TN}-\mathrm{FP}\times \mathrm{FN}}{\sqrt{\left(\mathrm{TP}+\mathrm{FP}\right)\times \left(\mathrm{TN}+\mathrm{FN}\right)\times \left(\mathrm{TP}+\mathrm{FN}\right)\times \left(\mathrm{TN}+\mathrm{FP}\right)}}$

where TP, TN, FP, and FN represent the numbers of true positives, true negatives, false positives, and false negatives, respectively. MCC and ACC are two metrics used to evaluate the overall prediction ability of a predictive model. The receiver operating characteristic curve (ROC), the area under the curve (AUC), and precision recall curves (PRC) are used to show the detailed performance of different methods [88]. The AUC ranges from 0.5 to 1. The higher the AUC score, the better the performance of the model [89,90,91].

3. Results and Discussion

3.1. Performance Comparison on Seven Benchmark Datasets

To demonstrate the effectiveness of the proposed method, we compared its performance with three other existing state-of-the-art methods, including MM-6mAPred [21], SNNRice6mA [33], and Deep6mA [34]. MM-6mAPred is traditional machine learning method training the model by hand-made features extracted from original DNA sequences. SNNRice6mA and Deep6mA are the deep learning methods. All the predictors are evaluated on the seven benchmark datasets by using 5-fold cross validation.

Table 2 lists the performances of the proposed method, namely 6mAPred-MSFF, and four existing predictors on the seven datasets. For the 6mA-rice-Lv dataset, when compared to the second-best method Deep6mA, our proposed method achieves an SN of 97.88%, an SP of 96.64%, an ACC of 96.26%, and an MCC of 0.926, yielding a relative improvement over Deep6mA of 3.15%, 0.92%, 2.03%, and 4.63%, respectively. Compared with the machine learning method MM-6mAPred, our proposed method improves the performance of SN, SP, ACC, and MCC by 4.61%, 5.8%, 5.12%, and 12.65%, respectively. For the 6mA-rice-Chen dataset, compared with the runner-up predictor Deep6mA, with the only exception that the value of SP is better than ours, our proposed method improves the performance of SN, ACC, and MCC by 9.89%, 2.79%, and 6.81%, respectively. For the A. thaliana and R. chinensis datasets, with the exception that the predictor Deep6mA achieves the best value of SP, and our proposed method improves the performance of SN, ACC, and MCC by 0.28–5.2%. Moreover, we can see that the performance of 6mAPred-MSFF is also better than all other methods on the F. vesca, H. sapiens, and D. melanogaster datasets.

For a more intuitive comparison, we further compared the ROC (receiver operating characteristic) curves and PR (precision-recall) curves of the four predictors, which are illustrated in Figure 2 and Figure 3, respectively. We can observe that our method achieves the best values for an auROC of 0.989 and an auPRC of 0.981 on the 6mA-rice-Lv dataset. We also achieve the best values for an auROC of 0.976 and auPRC of 0.973 on the 6mA-rice-Chen dataset. For the other species, our proposed method also achieves the best values of auROC and auPRC. Therefore, we can conclude that our proposed method can achieve the best predictive performance for detecting 6mA sites in rice species.

3.2. Validation on Other Species

To validate whether the model trained on rice data is applicable to predict the 6mA sites of other species. We take the datasets including A. thaliana, R. chinensis, F. vesca, H. sapiens, and D. melanogaster as the independent test datasets. We evaluate the performance of 6mAPred-MSFF, which was trained on the rice 6mA-rice-Lv dataset and performed the same evaluation for other methods including Deep6mA, SNNRice6mA-large, and MM-6mAPred.

The prediction results are listed in Table 3. For the two plant species (F. vesca and R. chinensis), with the only exception that the value of SP of our proposed method is lower than Deep6mA, and our proposed method outperforms all other methods in terms of SN, SP, ACC, and MCC. Specifically, for R. chinensis, our proposed method achieves an SN of 90.66%, an SP of 94.26%, an ACC of 92.46%, and an MCC of 0.850, where Deep6mA does have a higher SP of 95.29%. For F. vesca, our proposed method achieves an SN of 93.95%, an SP of 95.15%, an ACC of 94.55%, and an MCC of 0.891. However, Deep6mA does have a higher SP of 96.22%. Compared with the two plant species (F. vesca and R. chinensis), the performance of predicting 6mA sites of the two species (H. sapiens and D. melanogaster) is not very good. For the two animal species (H. sapiens and D. melanogaster), our proposed method outperforms all other methods in terms of SN, ACC, and MCC; however, the Deep6mA achieves the best value of SP. Surprisingly, all models trained on the 6mA-rice-Lv dataset fail in predicting the 6mA sites of A. thaliana, since the best value of SN is 58.66%; although there is a better value of SP, the performance of predicting 6mA sites of A. thaliana can not meet the requirement in practice. Therefore, we can conclude that the models trained with an 6mA-rice-Lv dataset can be used to predict the 6mA sites of other species and the features of 6mA sites also have certain similarities among different species.

3.3. Performance Impact by Different Encoding Schemes and Feature Descriptors

Our proposed model uses four different encoding schemes (1-gram, NAC, 2-grams, DNC) as the input of the embedding layer. The initial features will have an important impact on the training results of deep learning models. To explore the impact of different encoding methods, we will evaluate our proposed model with N-grams encoding scheme and different feature descriptors as the initial input.

N-grams encoding. To validate the impact of the N-grams encoding, we compare our proposed model with the encoding schemes (1-gram, 2-grams, 3-grams 1-gram and 2-grams, 1-gram and 3-grams, 2-grams and 3-grams) based on the 6mA-rice-Lv dataset. The prediction results are listed in Table 4. From Table 4, we can see that the 1-gram encoding scheme achieves an SN of 97.15% second only to our proposed method. The 1-gram and 2-gram encoding scheme achieves an SN of 94.47%, an SP of 94.67%, an ACC of 95.86, and an MCC of 0.917. Therefore, our proposed method uses the 1-gram and 2-gram encoding scheme integrating NAC and DNC to encode the input sequences.

Feature descriptors. There are several feature descriptors including Kmer, NCP, EIIP, and ENAC which are employed in machine learning methods. To validate the impact of the different feature descriptors, we take the feature descriptors as the initial input of our proposed model trained on 6mA-rice-Lv and perform the same evaluation on the Deep6mA, Random Forest (RF), and Artificial Neural Network (ANN). The prediction results are listed in Table 5. From Table 5, we can see that our proposed method trained by the NCP achieves the best values of an SN of 96.05%, an SP of 94.56%, an ACC of 95.31%, and an MCC of 0.906. The models trained by Kmer which concatenate 1mer, 2mer, and 3mer do not achieve satisfactory results. Meanwhile, we can observe that our proposed method outperforms all other methods in terms of SN, SP, ACC, and MCC under different feature descriptors.

For a more intuitive comparison, we further compared the ROC (receiver operating characteristic) curves and PR (precision-recall) curves of the four feature descriptors, which are illustrated in Figure 4 and Figure 5, respectively. We can observe that our proposed method achieves the best values of an auROC of 0.986 and an auPRC of 0.984 under the NCP. Our proposed method also achieves the best values of an auROC and auPRC under other feature descriptors. Therefore, we can conclude that our proposed method has a better ability to capture feature information.

4. Conclusions

In this study, we have proposed a novel predictor called 6mAPred-MSFF to predict the DNA 6mA sites. 6mAPred-MSFF is the first deep learning predictor, in which we integrate the global and local context by the inverted residual block and multi-scale channel attention module (MS-CAM). To fuse different feature vectors by fusion weights, 6mAPred-MSFF uses an attentional feature fusion (AFF) module based on the MS-CAM. As compared to existing predictors, our proposed method can automatically learn the global and local features and capture the characteristic specificity of 6mA sites. On the other hand, the experimental results demonstrate that our proposed method can effectively increase the accuracy and the generalization ability of the DNA 6mA sites prediction. Our proposed deep learning method shows better performance on the rice benchmark and the independent test of the other species.

Author Contributions

R.Z. and M.L.: conceptualization, investigation, methodology, formal analysis, supervision, project administration. R.Z.: software, validation, resources, data curation, writing—original draft preparation, writing—review and editing, visualization. M.L.: funding acquisition. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by “The National Natural Science Foundation of China, grant number 61702450”.

Institutional Review Board Statement

Not available.

Informed Consent Statement

Not available.

Data Availability Statement

Publicly available datasets were analyzed in this study. The data and source code can be found here: https://github.com/raozeng/6mAPred-MSFF (accessed on 20 August 2021).

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The flowchart of 6mAPred-MSFF. The sequence embedding module uses four different encoding schemes (1-gram, NAC, 2-grams, DNC) for the input of the embedding layer. Next, the feature matrices are fed into the feature extraction module to capture and combine the global and local features, respectively. Afterwards, the local and global contexts of the four different encoding schemes are aggregated into the feature fusion module. Finally, the output of the feature fusion module is fed into the prediction module to predict the 6mA site of a certain species.

Figure 2. The ROC curves (A–G) of four models (6mAPred-MSFF, Deep6mA, SNNRice6mA-large, and MM-6mAPred) based on the seven datasets (6mA-rice-Lv, 6mA-rice-Chen, A. thaliana, R. chinensis, F. vesca, H. sapiens, and D. melanogaster).

Figure 3. The PRC curves (A–G) of four models (6mAPred-MSFF, Deep6mA, SNNRice6mA-large, and MM-6mAPred) based on the seven datasets (6mA-rice-Lv, 6mA-rice-Chen, A. thaliana, R. chinensis, F. vesca, H. sapiens, and D. melanogaster).

Figure 4. The ROC curves (A–D) of four models (6mAPred-MSFF, Deep6mA, SNNRice6mA-large, and MM-6mAPred) based on four feature descriptors (Kmer, NCP, EIIP, ENAC).

Figure 5. The PR curves (A–D) of four models (6mAPred-MSFF, Deep6mA, SNNRice6mA-large, and MM-6mAPred) based on four feature descriptors (Kmer, NCP, EIIP, ENAC).

Table 1

Summary of seven benchmark datasets of six species.

Dataset	Positives	Negatives	Total
6mA-rice-Lv	154,000	154,000	300,800
6mA-rice-Chen	880	880	1760
A. thaliana	15,937	15,937	31,874
R. chinensis	11,815	11,815	23,630
F. vesca	22,700	22,700	45,400
H. sapiens	18,335	18,335	36,670
D. melanogaster	11,191	11,191	22,382

Table 2

Performance comparison of the proposed method and 6mA predictors based on seven datasets.

Benchmark Dataset	Method	SN (%)	SP (%)	ACC (%)	MCC	AUC
6mA-rice-Lv	6mAPred-MSFF	97.88	94.64	96.26	0.926	0.989
	Deep6mA	94.73	93.72	94.23	0.885	0.981
	SNNRice6mA-large	94.97	92.22	93.59	0.872	0.978
	MM-6mAPred	93.27	88.84	91.05	0.822	0.955
6mA-rice-Chen	6mAPred-MSFF	96.59	89.43	93.01	0.862	0.976
	Deep6mA	86.70	93.75	90.22	0.807	0.958
	SNNRice6mA-large	87.61	86.70	87.16	0.743	0.942
	MM-6mAPred	87.50	88.86	88.18	0.764	0.943
A. thaliana	6mAPred-MSFF	89.92	87.85	88.88	0.778	0.954
	Deep6mA	84.72	89.52	87.12	0.743	0.942
	SNNRice6mA-large	85.28	88.15	86.71	0.735	0.939
	MM-6mAPred	82.62	88.82	83.82	0.664	0.904
R. chinensis	6mAPred-MSFF	95.61	96.56	96.09	0.922	0.988
	Deep6mA	94.97	96.66	95.81	0.917	0.987
	SNNRice6mA-large	91.96	94.55	93.25	0.865	0.977
	MM-6mAPred	92.27	90.72	91.50	0.830	0.959
F. vesca	6mAPred-MSFF	98.79	98.44	98.62	0.972	0.998
	Deep6mA	98.35	97.10	98.48	0.970	0.998
	SNNRice6mA-large	97.72	97.65	97.72	0.955	0.996
	MM-6mAPred	95.13	94.72	94.93	0.899	0.980
H. sapiens	6mAPred-MSFF	94.65	94.91	94.78	0.896	0.985
	Deep6mA	90.16	95.18	92.67	0.854	0.985
	SNNRice6mA-large	88.62	93.56	91.09	0.823	0.972
	MM-6mAPred	88.08	87.22	87.65	0.753	0.937
D. melanogaster	6mAPred-MSFF	97.08	96.72	96.90	0.938	0.990
	Deep6mA	95.12	95.04	95.08	0.902	0.984
	SNNRice6mA-large	94.64	93.36	94.00	0.880	0.982
	MM-6mAPred	90.65	89.94	90.30	0.806	0.955

Table 3

Performance comparison of methods trained with the 6mA-rice-Lv dataset on cross species.

Species	Method	SN (%)	SP (%)	ACC (%)	MCC	AUC
R. chinensis	6mAPred-MSFF	90.66	94.26	92.46	0.850	0.966
	Deep6mA	88.87	95.29	92.07	0.843	0.966
	SNNRice6mA-large	89.79	94.13	91.96	0.840	0.964
	MM-6mAPred	88.51	91.35	89.93	0.799	0.949
F. vesca	6mAPred-MSFF	93.95	95.15	94.55	0.891	0.982
	Deep6mA	92.52	96.22	94.37	0.888	0.982
	SNNRice6mA-large	93.26	95.15	94.20	0.884	0.980
	MM-6mAPred	90.33	93.48	91.91	0.839	0.965
H. sapiens	6mAPred-MSFF	72.51	93.33	82.92	0.673	0.884
	Deep6mA	68.82	94.44	81.63	0.655	0.881
	SNNRice6mA-large	70.93	93.22	82.07	0.658	0.894
	MM-6mAPred	71.85	88.67	80.26	0.614	0.873
D. melanogaster	6mAPred-MSFF	74.30	92.62	83.46	0.680	0.903
	Deep6mA	68.81	93.78	81.30	0.646	0.901
	SNNRice6mA-large	70.73	92.27	81.50	0.645	0.906
	MM-6mAPred	64.96	88.84	76.90	0.554	0.868
A. thaliana	6mAPred-MSFF	58.66	94.28	76.47	0.567	0.837
	Deep6mA	54.56	95.26	74.91	0.545	0.834
	SNNRice6mA-large	57.92	94.11	76.01	0.558	0.853
	MM-6mAPred	57.14	91.44	74.30	0.517	0.845

Table 4

Performance comparison of the proposed method on different encoding schemes based on the 6mA-rice-Lv dataset.

Method	SN (%)	SP (%)	ACC (%)	MCC	AUC
6mAPred-MSFF	97.88	94.64	96.26	0.926	0.989
1-gram	97.15	93.68	95.42	0.909	0.987
2-grams	95.63	93.56	94.59	0.892	0.983
3-grams	96.10	93.60	94.85	0.897	0.984
1-gram and 2-grams	97.05	94.67	95.86	0.917	0.988
1-gram and 3-grams	95.97	93.49	94.73	0.895	0.983
2-grams and 3-grams	95.77	94.51	95.14	0.903	0.985

Table 5

Performance comparison of methods based on the 6mA-rice-Lv dataset under different feature descriptors.

Feature	Method	SN (%)	SP (%)	ACC (%)	MCC	AUC
Kmer(1mer + 2mer + 3mer)	6mAPred-MSFF	67.17	64.94	66.08	0.321	0.719
	Deep6mA	69.00	60.35	64.67	0.294	0.700
	RF	66.96	62.33	64.65	0.293	0.697
	ANN	66.87	62.35	64.61	0.293	0.697
NCP	6mAPred-MSFF	96.05	94.56	95.31	0.906	0.986
	Deep6mA	95.82	92.95	94.38	0.888	0.983
	RF	92.90	91.51	92.20	0.844	0.966
	ANN	92.92	91.53	92.22	0.845	0.966
EIIP	6mAPred-MSFF	95.35	91.10	93.22	0.865	0.975
	Deep6mA	94.20	88.77	91.48	0.831	0.966
	RF	92.56	91.21	91.89	0.838	0.964
	ANN	92.61	91.15	91.88	0.838	0.964
ENAC	6mAPred-MSFF	95.34	91.89	93.62	0.872	0.979
	Deep6mA	95.35	90.29	92.82	0.858	0.972
	RF	90.82	91.16	90.99	0.820	0.959
	ANN	90.80	91.17	90.99	0.820	0.959