RNA binding site prediction method based on Mangbar and graph neural network

Abstract

The invention belongs to the field of bioinformatics, and particularly relates to an RNA binding site prediction method based on a Mangbar and graph neural network. The method comprises the following steps: firstly, extracting sequence and structure features from an RNA secondary structure file, and constructing RNA graph structure data; secondly, on the basis of node features in the graph, a bidirectional Mangban state space model is adopted for sequence coding; then, based on an edge connection relationship, constructing a multi-scale convolutional pyramid network to extract structural features; then, carrying out adaptive fusion on the sequence features and the multi-scale structure features; and finally, based on the fusion features, generating a prediction result by adopting a dynamic routing mechanism. According to the method, RNA sequence and structural characteristics can be automatically analyzed, and high-precision binding site prediction is realized through a dynamic multi-modal fusion mechanism. The invention aims to solve the problems of low sequence modeling efficiency, single feature fusion mechanism and insufficient structural information utilization in the prior art, and realizes efficient and accurate RNA binding site prediction.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics, specifically relating to an RNA binding site prediction method based on Mamba and graph neural networks. This method can automatically analyze RNA sequence and structural features and achieve high-precision binding site prediction through a dynamic multimodal fusion mechanism. Background Technology

[0002] RNA-protein interactions play a crucial role in bioinformatics processes such as gene expression regulation and post-transcriptional modifications. Accurate prediction of RNA binding sites is essential for understanding cellular function and disease mechanisms. Existing prediction methods are mainly based on traditional machine learning or deep learning models, such as graph neural networks like MVRBind. However, these methods still have the following problems: ① Limited sequence modeling capability: Traditional Transformer models have high computational complexity when processing long RNA sequences and struggle to capture long-term dependencies; ② Simple feature fusion mechanism: Multimodal feature fusion typically uses simple splicing or fixed weights, failing to adaptively adjust the contributions of different features; ③ Insufficient structural information: Multi-scale feature extraction of RNA secondary structures is insufficient, lacking hierarchical modeling. Therefore, there is an urgent need in this field for a prediction method that can efficiently process RNA sequences, fully mine structural information, and provide adaptive feature fusion. Summary of the Invention

[0003] This invention proposes an RNA binding site prediction method based on Mamba and graph neural networks, aiming to solve the problems of low sequence modeling efficiency, single feature fusion mechanism, and insufficient utilization of structural information in existing technologies, and to achieve efficient and accurate RNA binding site prediction. The specific technical solution includes five steps: Step 1, RNA Feature Extraction and Graph Structure Construction: Extract sequence information and structural features from the RNA secondary structure file to construct RNA graph structure data; Step 2, Bidirectional State-Space Sequence Encoding: Based on the node features in the RNA graph structure constructed in Step 1, a bidirectional Mamba state-space model is used for sequence feature encoding; Step 3, Multi-scale graph convolution feature extraction: Based on the edge connection relationships of the RNA graph structure in Step 1, a multi-scale convolutional pyramid network is constructed to extract structural features; Step 4, Cross-modal attention fusion: Adaptively fuse the sequence obtained in Step 2 with the multi-scale structural features obtained in Step 3; Step 5: Dynamic route prediction output: Based on the fused features obtained in Step 4, the final prediction result is generated using a dynamic routing mechanism.

[0004] A method for predicting RNA binding sites based on Mamba and graph neural networks, the implementation process of step 1 is as follows: The system reads sequence and base pairing information from the RNA secondary structure file and calculates a six-dimensional feature vector for each nucleotide. For the first nucleotide in the sequence... Each nucleotide, its node feature vector The structure is as follows: ,in, For one-hot coding functions of nucleotide types, This represents vector concatenation. An indicator function representing the pairing status. The normalized pairing strength is then used; subsequently, the edges of the graph are constructed based on sequence adjacency, base pairing, and spatial proximity; finally, the node features and edge indices are encapsulated into a graph data object.

[0005] A method for predicting RNA binding sites based on Mamba and graph neural networks, the implementation process of step 2 is as follows: First, the node features are projected into a high-dimensional space, and then input into the forward and backward Mamba state space modules for sequence encoding. The Mamba module captures the long-term dependencies of the sequence through the following discretized state space equation: ,in, for The hidden state at all times For input, For output, For continuous parameters Discretized using the zero-order ZOH rule, we obtain: ,in The time step parameter is used; the outputs from both directions are concatenated and fused to obtain the final sequence feature representation. 。

[0006] A method for predicting RNA binding sites based on Mamba and graph neural networks, the implementation process of step 3 is as follows: First, structural features are extracted using a multi-scale graph convolutional pyramid network. Then, multiple cascaded graph attention convolutional layers (TransformerConv) are used, and the original node features are fused into the input of each layer. Features of layer graph convolution The updated formula is: ,in, It is the ReLU activation function. This indicates the splicing or averaging of multiple attention streams. This is the normalized adjacency matrix calculated based on edge indices. For the first Each attention head has a learnable weight matrix. Finally, the outputs from each scale are stacked to form a multi-scale structural feature. .

[0007] A method for predicting RNA binding sites based on the Mamba model and graph neural networks, the implementation process of step 4 is as follows: First, multi-scale structural features are pooled to obtain... Then, using sequence features As a query Structural features after pooling As a key Sum Multi-head attention calculation is performed. The formula for calculating attention weights is: Finally, a gating fusion mechanism is used to adaptively weight and fuse the original sequence features with the attention output. (Gating weights) The calculation formula is: To obtain the aligned features ,in This indicates element-wise multiplication.

[0008] A method for predicting RNA binding sites based on Mamba and graph neural networks, the implementation process of step 5 is as follows: First, a dynamic routing mechanism is adopted to calculate the routing weights between sequence and structural features through multiple capsule networks; second, for the... Each capsule, its routing score The cosine similarity matrix between the transformed sequence features and the graph features is calculated and its mean is taken. The routing weight α is obtained by softmax normalization: Features are weighted and fused based on their weights. After global pooling and an output network consisting of fully connected layers, normalization layers, and activation functions, the predicted probability of the binding site is finally output through the sigmoid function. . Attached Figure Description

[0009] Figure 1 This is a diagram illustrating the overall architecture of an RNA binding site prediction method based on Mamba and graph neural networks.

[0010] Figure 2 This is a flowchart of the RNA feature extraction and graph structure construction module.

[0011] Figure 3 This is a flowchart of the bidirectional state-space sequence encoding module.

[0012] Figure 4 This is a flowchart of the multi-scale graph convolution feature extraction module.

[0013] Figure 5 This is a flowchart of the cross-modal attention fusion module.

[0014] Figure 6 This is a flowchart of the dynamic route prediction output module. Detailed Implementation

[0015] The present invention will now be described in detail with reference to the accompanying drawings and examples.

[0016] RNA feature extraction and graph structure construction: such as Figure 1 and Figure 2As shown, the specific implementation process of this step is as follows: The system reads multiple structure files from the RNA secondary structure directory. First, it reads the bpseq format base pairing file, which is in text format. Each line contains three fields: nucleotide number, nucleotide type, and paired nucleotide number. The system uses Python's open function to open the file in read-only mode, reads the content line by line, removes leading and trailing whitespace characters from each line using the strip method, and splits it by whitespace characters using the split method. If the length of the split list is greater than or equal to 3 and it is not a comment line starting with #, the first element is extracted as the nucleotide number; the second element is the nucleotide type (A, U, G, C); and the third element is the paired nucleotide number. All nucleotide types are stored in the list sequence, and pairing information is stored in the list pairings. Next, it reads the ct format connection table file, which describes in detail the connection relationship of each nucleotide in the RNA sequence. The system skips the first header row and starts reading from the second row. Each row is split using the `split` method, extracting the first element as the current nucleotide sequence number and the fifth element as the paired nucleotide sequence number. When the fifth element is not 0, (current sequence number, paired sequence number) is stored as a key-value pair in the dictionary `ct_pairs`. Simultaneously, multiple CSV-formatted interaction files are read, including `2JUKA-canonical.csv` (canonical base pair), `2JUKA-stacking.csv` (base stacking), `2JUKA-base-phosphate.csv` (base-phosphate interaction), and `2JUKA-base-ribose.csv` (base-ribose interaction). These files are read using the `read_csv` function from the pandas library, checking for the existence of columns `position1` and `position2` (for canonical base pair), `base1` and `base2` (for base stacking), `base` and `phosphate` (for base-phosphate), or `base` and `ribose` (for base-ribose). The corresponding column data is extracted and stored as a list of tuples. After data reading is complete, the system performs data cleaning. First, filter out self-loop edges: traverse the pairings list, and if the pairing number at a certain position is equal to the current number plus 1 (since bpseq starts counting from 1, while Python lists start from 0), then set the pairing number at that position to 0; then verify all numerical characteristics: check whether the pairing number is a valid number, ensuring that it is greater than or equal to 0 and less than or equal to the sequence length.Finally, sequence consistency is checked: ensure that the sequence length in the bpseq file is consistent with the sequence length in the ct file; during the feature calculation stage, the system constructs a six-dimensional feature vector for each nucleotide; for each position i in the sequence (from 0 to the sequence length minus 1): create a list one_hot of length 4, initialized to [0,0,0,0]; if sequence[i] equals 'A', set one_hot[0] to 1; if sequence[i] equals 'U' or 'T', set one_hot[1] to 1; if sequence[i] equals 'G', set one_hot[2] to 1; if sequence[i] equals 'C', set one_hot[3] to 1. Pairing state feature: if pairings[i] is greater than 0, then paired is 1, otherwise it is 0. Normalized pairing strength feature: normalized_strength equals pairs[i] divided by sequence length; concatenate these three features into a six-dimensional vector: node_feature equals one_hot plus [paired, normalized_strength]; during the graph structure construction phase, the system creates an empty edge list edge_list; add sequence adjacency edges: for i from 0 to sequence length minus 2, add [i, i+1] and [i+1, i] to edge_list; add base pairing edges: traverse the pairs list, for each position i, if pairs[i] is greater than 0 and i+1 is less than pairs[i], then j equals pairs[i] minus 1 (converted to zero-based index), add [i, j] and [j, i] to edge_list; add spatial proximity edges: for each position i, j from max(0, i-3) to min(sequence length, i+4), if i is not equal to j, add [i, j] to edge_list. Convert `edge_list` to a PyTorch tensor: `edge_index` equals `torch.tensor(edge_list, dtype=torch.long).t().contiguous()`, with shape [2, number of edges]. Convert `node_features` to a PyTorch tensor: `node_features` equals `torch.tensor(node_features, dtype=torch.float32)`, with shape [sequence length, 6]. Wrap the graph data using the PyTorch Geometric Data class: `graph_data` equals `Data(x=node_features, edge_index=edge_index)`.

[0017] Bidirectional state-space sequence encoding: such as Figure 3 As shown, the specific implementation process of this step is as follows: The system first creates an input projection layer: using nn.Linear(6, 128) to project the six-dimensional node features onto a 128-dimensional hidden space; the node features are then passed through the projection layer: x_proj equals self.input_proj(graph_data.x), and the shape changes from [sequence length, 6] to [sequence length, 128]. The MambaBlock class is defined, containing the following components: layer normalization uses nn.LayerNorm(128); input projection uses nn.Linear(128, 256); one-dimensional convolution uses nn.Conv1d(256, 256, kernel_size=4,groups=256, padding=3); state space parameter A_log uses nn.Parameter(torch.randn(256,16)); output projection uses nn.Linear(256, 128). Forward processing: x_forward equals x_proj, passing through a forward Mamba block; inside the Mamba block: the input is saved as a residual; the input is normalized; the 128-dimensional array is mapped to 256-dimensional arrays through an input projection layer; one-dimensional convolution is performed, first transposing the dimension to [1, 256, sequence length], the convolution output is clipped to the original sequence length, then transposed back to [sequence length, 256], and the SiLU activation function is used. State space model calculation: A equals negative exponential A_log, the hidden state h is initialized to zeros(1, 256, 16); for each time step t of the sequence: the discretization parameter dA is calculated to be equal to the exponential A, the state update h is equal to dA multiplied by h plus the current input, the output y_t is equal to the sum of h plus D multiplied by the current input; the outputs of all time steps are stacked and passed through an output projection layer; finally, the output is connected with residuals. Reverse processing: x_reversed is equivalent to flipping x_proj in the sequence dimension; x_backward_temp is obtained through the same Mamba block; x_backward is equivalent to flipping x_backward_temp again and aligning it with the original sequence; Bidirectional fusion: combined is equivalent to concatenating x_forward and x_backward in the feature dimension with a shape of [sequence length, 256]; a fusion layer nn.Linear(256, 128) is created; seq_features is equivalent to the fusion layer (combined) with a shape of [sequence length, 128].

[0018] Multi-scale graph convolution feature extraction: such as Figure 4As shown, the specific implementation process of this step is as follows: The system constructs three Transformer graph convolutional layers; the first layer uses TransformerConv(128, 128, heads=4, concat=False), the second layer uses TransformerConv(256, 128, heads=4, concat=False), and the third layer uses TransformerConv(256, 128, heads=4, concat=False). First-scale graph convolution: Input x_proj has shape [sequence length, 128]; passes through the first Transformer convolutional layer, using graph_data.edge_index as the edge index; output is activated by ReLU, resulting in x_scale1 with shape [sequence length, 128]. Second-scale graph convolution: Input x_input2 is equal to concatenating x_scale1 and x_proj along the feature dimension, with shape [sequence length, 256]; passes through the second Transformer convolutional layer; output is activated by ReLU, resulting in x_scale2 with shape [sequence length, 128]. Third-scale graph convolution: Input x_input3 is equal to concatenating x_scale2 and x_proj along the feature dimension, with shape [sequence length, 256]; passes through the third Transformer convolutional layer; output is activated by ReLU, resulting in x_scale3 with shape [sequence length, 256]. 128]; Multi-scale feature stacking: graph_features is equivalent to stacking x_scale1, x_scale2, and x_scale3 in dimension 1 using torch.stack([x_scale1, x_scale2, x_scale3], dim=1), with the shape [sequence length, 3, 128].

[0019] Cross-modal attention fusion: such as Figure 5As shown, the specific implementation process of this step is as follows: Graph feature pooling: graph_features_pooled is equal to taking the mean of graph_features in dimension 1, using torch.mean(graph_features, dim=1), changing the shape from [sequence length, 3, 128] to [sequence length, 128]; Creating attention projection layers: query projection layer nn.Linear(128, 128), key projection layer nn.Linear(128, 128), value projection layer nn.Linear(128, 128), output projection layer nn.Linear(128, 128); Multi-head attention calculation: Q is equal to reshaping seq_features into [sequence length, 4, 32] after passing through the query projection layer and transposing it to [4, sequence length, 32]; K is equal to reshaping graph_features_pooled into [sequence length, 4, 32] after passing through the key projection layer and transposing it to [4, sequence length, 32] 32]; V equals reshaping graph_features_pooled through the value projection layer to [sequence length, 4, 32] and transposing it to [4, sequence length, 32]; For each attention head i from 0 to 3: attention_weights[i] equals matmul(Q[i], K[i] transpose) divided by the square root of 32, shape [sequence length, sequence length]; attention_weights[i] equals softmax(attention_weights[i], dim=-1); head_output[i] equals matmul(attention_weights[i], V[i]), shape [sequence length, 32]; Transpose all head outputs to [sequence length, 4, 32] and reshape them to [sequence length, 128]; attention_output equals passing through the output projection layer, shape [sequence length, 128]; Create a gated fusion network: nn.Sequential(nn.Linear(256, 128), nn.Sigmoid()); combined_features is equal to concatenating seq_features and attention_output along the feature dimension, with shape [sequence length, 256]; gate_weights is equal to fusion through a gated network, with shape [sequence length, 128]; aligned_features is equal to gate_weights multiplied by seq_features plus (1-gate_weights) multiplied by attention_output, with shape [sequence length, 128].

[0020] Dynamic route prediction output: such as Figure 6As shown, the specific implementation process of this step is as follows: Define dynamic routing parameters: num_capsules equals 4; routing_weights equals nn.Parameter(torch.randn(4, 128, 128)); capsule_biases equals nn.Parameter(torch.randn(4, 128)); Graph feature pooling: graph_pooled equals the mean of graph_features in dimension 1, shape [sequence length, 128]; Routing score calculation: Initialize an empty list routing_scores; For i from 0 to 3: transformed_seq equals matmul(aligned_features, routing_weights[i]), shape [sequence length, 128]; transformed_graph equals matmul(graph_pooled, routing_weights[i]), shape [sequence length, 128]; similarity equals cosine_similarity(transformed_seq increases dimension 1, transformed_graph increases dimension 0, dim=-1), shape[sequence length, sequence length]; routing_scores.append(mean(similarity)); routing_scores equals stack(routing_scores), shape[4]; routing_weights equals softmax(routing_scores, dim=0), shape[4]; weighted feature fusion: fused_features initialized to zeros_like(aligned_features), shape[sequence length, 128]; for i from 0 to 3: weight equals routing_weights[i]; transformed_seq equals matmul(aligned_features, routing_weights[i]), shape[sequence length, 128]; capsule_output equals transformed_seq plus capsule_biases[i], shape[sequence length,

[128] ; fused_features equals fused_features plus weight multiplied by capsule_output; global pooling: global_feature equals the mean of fused_features at dimension 0, using torch.`mean(fused_features, dim=0, keepdim=True)`, shape [1, 128]; Create the output network: `nn.Sequential(nn.Linear(128, 256), nn.LayerNorm(256), nn.ReLU(), nn.Dropout(0.3), nn.Linear(256, 128), nn.LayerNorm(128), nn.ReLU(), nn.Dropout(0.2), nn.Linear(128, 1), nn.Sigmoid())`. Final prediction: `prediction` equals the output network (global_feature), then `squeeze()` to obtain the scalar prediction probability value.

[0021] To verify the effectiveness of this invention, tests were conducted on the 2JUK benchmark dataset. Compared with the MVPBind method, this invention improved the ACC value from 0.86 to 0.90, an improvement of 4%, demonstrating the effectiveness and superiority of this invention.

[0022] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A method for predicting RNA binding sites based on Manba and graph neural networks, characterized in that, The specific technical solutions include five steps: Step 1, RNA feature extraction and graph structure construction: extract sequence information and structure features from RNA secondary structure files, and construct RNA graph structure data; Step 2, bidirectional state space sequence encoding: based on the node features in the RNA graph structure constructed in step 1, use the bidirectional Mamba state space model for sequence feature encoding; Step 3, multi-scale graph convolution feature extraction: based on the edge connection relationship of the RNA graph structure in step 1, construct a multi-scale convolution pyramid network to extract structure features; Step 4, cross-modal attention fusion: adaptively fuse the sequence obtained in step 2 and the multi-scale structure features obtained in step 3; Step 5, dynamic routing prediction output: based on the fusion features obtained in step 4, use the dynamic routing mechanism to generate the final prediction result.

2. The RNA binding site prediction method based on the Manba and graph neural network according to claim 1, wherein, The implementation process of step 1 is as follows: the system reads the sequence and base pairing information from the RNA secondary structure file, calculates a six-dimensional feature vector for each nucleotide; for the first nucleotide in the sequence, its node feature vector is constructed as follows: is the normalized pairing strength; then, the edges of the graph are constructed based on sequence adjacency, base pairing and spatial proximity; finally, the node features and edge indices are packaged into a graph data object.​​​​​​ 3. The method of claim 1, wherein, The implementation process of step 2 is as follows: first, project the node features into a high-dimensional space, and then input the forward and reverse Mamba state space modules respectively for sequence encoding. The Mamba module captures the long-term dependencies of the sequence through the following discrete state space equation: wherein, is the hidden state at time t, is the input, is the output, is the output, is the output, is discretized by the zero-order hold (ZOH) rule: wherein is the time step parameter; the outputs of the two directions are spliced and fused to obtain the final sequence feature representation: 。 4. The method of claim 1, wherein, The implementation process of step 3 is as follows: first, structural features are extracted using a multi-scale graph convolutional pyramid network, then a plurality of serial graph attention convolution layers TransformerConv are passed through, and the original node features are fused in the input of each layer, the first Layer graph convolution features The update formula is: wherein, is a ReLU activation function, denotes a concatenation or averaging operation of multi-head attention, is a normalized adjacency matrix calculated based on edge index, is a learnable weight matrix of the first head; finally, the outputs of each scale are stacked to form a multi-scale structural feature .

5. The method of claim 1, wherein, The implementation process of step 4 is as follows: first, the multi-scale structure features are pooled to obtain Then, the sequence features are taken as the query , the pooled structure features are taken as the key and the value , and multi-head attention calculation is performed; the calculation formula of the attention weight is: Finally, through a gating fusion mechanism, the original sequence features and the attention output are adaptively weighted and fused; the calculation formula of the gating weight is: , and the aligned features are obtained, wherein represents element-wise multiplication.

6. The method of claim 1, wherein, The implementation process of step 5 is as follows: first, a dynamic routing mechanism is used to calculate the routing weights between the sequence and structure features through multiple capsule networks; second, for the first capsule, the routing score is calculated as follows: The cosine similarity matrix between the transformed sequence features and the graph features is calculated and averaged to obtain: The routing weight α is obtained by softmax normalization: The features are weighted and fused according to the weight: After global pooling and an output network composed of a fully connected layer, a normalization layer and an activation function, the prediction probability of the binding site is finally output by the Sigmoid function .​

Images