Crop centromere automatic identification method

Abstract

The invention discloses a crop centromere automatic identification method, and belongs to the field of crop centromere identification, and the method comprises the steps: constructing a training set according to a verified CENH3 binding region sequence and a random non-centromere region sequence; extracting multi-dimensional features of the DNA sequence according to the training set; training a CNN-RNN hybrid model according to the multi-dimensional features and obtaining a prediction model used for identifying centromere; and identifying the target genome sequence according to the prediction model to obtain a candidate centromere region. The method does not depend on repetitive sequences any more, and directly utilizes the CENH3 binding characteristics, so that the prediction accuracy is improved. And cross-species application can be realized, so that the general applicability of centromere region identification is improved.

Description

Technical Field

[0001] This invention belongs to the field of crop centromere identification, and particularly relates to an automatic method for identifying crop centromeres. Background Technology

[0002] For centromere identification in crops, existing technologies mainly include experimental verification methods, sequence feature-based prediction methods, and the emerging deep learning models. These technologies have promoted the development of centromere research to some extent, but still have many limitations, which restrict the progress of molecular breeding of crops.

[0003] For example, the CENH3 chromatin immunoprecipitation (ChIP-seq) technique, used as an experimental verification method, can directly verify the binding of CENH3 to the centromere of the genome, but the experiment is costly and complex.

[0004] Prediction methods based on sequence features make indirect judgments by identifying the location of tandem repeat sequences. However, taking maize as an example, existing studies have shown that the centromere region and the tandem repeat region do not always coincide, indicating that prediction methods based on repeat sequences suffer from low accuracy and poor generalization.

[0005] In addition, deep learning, which has been developed in recent years, has included several models for predicting protein-DNA binding (such as DeepBind, DeepSEA, and BPNet). These models have been successfully applied to the prediction of transcription factor binding sites in mammals, but have not yet been systematically applied to the prediction of CENH3-DNA binding.

[0006] Existing technologies have significant limitations in accurately predicting centromere positions across multiple species, especially when repetitive sequence signals are weak or do not overlap with functional regions, making prediction accuracy difficult to guarantee. Therefore, repetitive sequence methods are unreliable, particularly performing poorly in plants such as maize. Furthermore, there is a lack of scalable, universal prediction models that incorporate fusion sequence features and protein-protein binding characteristics, making it difficult to adapt to the genomic characteristics of different crops. In addition, existing methods cannot rapidly locate centromere regions in newly assembled crop genomes (such as rice, wheat, and maize), directly hindering the fine mapping of genes for important agronomic traits, chromosome engineering breeding, and genome selection breeding. This has become a key technological obstacle to overcoming molecular breeding bottlenecks in improving crop yield and stress resistance. Summary of the Invention

[0007] To solve the above-mentioned technical problems, the present invention provides an automatic centromere identification method for crops, comprising: A training set was constructed based on the validated CENH3 binding region sequence and random noncentromere region sequences; Extract multi-dimensional features of DNA sequences from the training set; A CNN-RNN hybrid model is trained based on the multi-dimensional features to obtain a predictive model for centromere identification; The target genome sequence is identified based on the prediction model to obtain candidate centromere regions.

[0008] Preferably, the process of constructing a training set based on the verified CENH3 binding region sequence and random noncentromere region sequences includes: A fixed-length positive sample was extracted from the published CENH3 ChIP-seq peak sequences of the species. A noncentromere region of the same length as the positive sample was randomly selected from the reference genome of the same species as a negative sample. The positive and negative samples are merged to obtain a labeled training set.

[0009] Preferably, the process of extracting multi-dimensional features of DNA sequences based on the training set includes: The DNA sequence is converted into a numerical matrix according to One-hot encoding; DNA secondary structure characteristics were calculated based on AT content, melting temperature, and curvature. Cross-species conservation scores were obtained based on multiple sequence alignment. The numerical matrix, DNA secondary structure features, and cross-species conservation scores are concatenated into a unified feature vector.

[0010] Preferably, the process of training a CNN-RNN hybrid model based on the multi-dimensional features includes: Local motif features are extracted from the convolutional layer to obtain the first feature map; The first feature map is temporally modeled using the LSTM layer to obtain the second feature map; The key site representation is obtained by weighting the second feature map using the attention layer; The key sites are represented as centromere probabilities based on the fully connected layer.

[0011] Preferably, the process of extracting local motif features from the convolutional layer and obtaining the first feature map includes: Set the kernel size to 10 and the number of kernels to 512; The input feature map is subjected to a one-dimensional convolution operation and then activated by ReLU. The output local motif feature map is used as the first feature map.

[0012] Preferably, the process of performing temporal modeling on the first feature map based on the LSTM layer and obtaining the second feature map includes: The number of hidden units is set to 256. The first feature map is input into the LSTM layer in sequence and the hidden state at each time step is obtained. The hidden states are then concatenated to obtain the second feature map.

[0013] Preferably, the process of weighting the second feature map according to the attention layer and obtaining the key site representation includes: Calculate the similarity score between each position in the second feature map and the query vector; The similarity scores are then normalized using softmax to obtain the attention weights. The attention weights are weighted and summed with the second feature map to obtain the key site representation.

[0014] Preferably, the process of obtaining a predictive model for identifying centromeres includes: The training set is divided into a training subset and a validation subset using five-fold cross-validation. Iteratively update the model parameters based on the training subset; The optimal model is selected as the prediction model based on the AUROC and AUPRC of the validation subset.

[0015] Preferably, the process of identifying the target genome sequence based on the prediction model includes: The target genome is segmented into continuous segments of the same length as the training set, and multi-dimensional features of each segment are extracted. The multidimensional features are input into the prediction model to obtain the centromere probability of each segment; Segments with a centromere probability higher than a preset threshold are merged into candidate centromere regions.

[0016] Preferably, the method further includes: Based on the prediction model trained on the first species, the genome of the second species is directly identified, and candidate centromere regions of the second species are output, thus achieving cross-species prediction.

[0017] Compared with the prior art, the present invention has the following advantages and technical effects: This invention no longer relies on repetitive sequences, but directly utilizes CENH3 binding features to achieve end-to-end automatic identification of crop centromeres at the computational level for the first time. Candidate regions can be output without any experimental steps, fundamentally breaking through the bottlenecks of traditional experiments such as high throughput, low cost, and long cycle, and improving prediction accuracy.

[0018] The method of this invention can be applied across species, and remains accurate, especially in plants with weak repetitive signals such as maize. Compared with existing prediction methods that rely on repetitive sequences, this invention can improve the accuracy by 30–40% on maize datasets, demonstrating universality. Furthermore, for new species, only the genome sequence needs to be input to quickly obtain candidate centromere regions.

[0019] This invention expands the original single information relying solely on tandem repeats into a three-dimensional description of "sequence-structure-evolution" by using a three-pronged approach of one-hot sequence, secondary structure, and conservation. Even when faced with complex genomes such as maize with weak or missing repeat signals, it can provide sufficient discriminative signals.

[0020] The hybrid model involved in this invention uses a three-level progressive structure: CNN first captures local motifs, LSTM then captures long-range dependencies, and Attention finally focuses on key sites. This enables the model to recognize short motif combinations and integrate distant context, thus comprehensively improving the sensitivity and specificity of plant CENH3 binding sites. Attached Figure Description

[0021] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the system structure according to an embodiment of the present invention. Detailed Implementation

[0022] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0023] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0024] like Figure 1 As shown, this embodiment provides an automatic centromere identification method for crops, including: A training set was constructed based on the validated CENH3 binding region sequence and random noncentromere region sequences; Extract multi-dimensional features of DNA sequences from the training set; A CNN-RNN hybrid model is trained based on multi-dimensional features to obtain a predictive model for centromere identification; The target genome sequence is identified based on the prediction model to obtain candidate centromere regions.

[0025] Furthermore, the process of constructing a training set based on the validated CENH3 binding region sequence and random noncentromere region sequences includes: A fixed-length positive sample was extracted from the published CENH3 ChIP-seq peak sequences of the species. In the same species reference genome, a noncentromere region of the same length as the positive sample was randomly selected as a negative sample; Positive and negative samples are merged to obtain a labeled training set.

[0026] Furthermore, in this embodiment, to construct a training set related to centromere prediction, we start with published species such as potato, yam, rice, corn, soybean, and wheat. On the one hand, we collect centromere binding region sequences verified by CENH3 ChIP-seq technology as positive samples. On the other hand, we randomly select non-centromere regions from the genomes of these species as negative samples. At the same time, we further collect reference genome sequences of the above species and finally integrate them to form a training set for subsequent research.

[0027] Furthermore, the process of extracting multi-dimensional features from the DNA sequence based on the training set includes: Convert the DNA sequence into a numerical matrix based on One-hot encoding; DNA secondary structure characteristics were calculated based on AT content, melting temperature, and curvature. Cross-species conservation scores were obtained based on multiple sequence alignment. The numerical matrix, DNA secondary structure features, and cross-species conservation scores are concatenated into a unified feature vector.

[0028] Furthermore, this embodiment uses one-hot encoding to digitally represent the DNA sequence, while introducing DNA secondary structure features such as AT content, melting temperature, and curvature, and incorporating information on conserved regions identified through multi-species sequence alignment, thereby achieving the integrated extraction of multi-dimensional features of the DNA sequence.

[0029] Specifically, the corresponding methods and formulas for the digital representation of DNA sequences are as follows: DNA sequences are composed of four bases (A, C, G, T). When using one-hot encoding, each base is mapped to a 4-dimensional sparse vector: ; Given a DNA sequence of length L: ; In the formula, This represents a vector representation of a DNA sequence of length L. This represents the i-th base in the sequence (i = 1, 2, …, L). Indicates the sequence length (i.e., the last base); Its one-hot representation is a The matrix: ; In the formula, Represents the one-hot coding matrix of the entire DNA sequence. This represents the one-hot vector of the last base (the Lth base) in the sequence. Represents the one-hot vector of the i-th base in the sequence; The formula for the characteristics of DNA secondary structure is as follows: AT content: ; In the formula, This indicates the number of bases A in the sequence. This indicates the number of T bases in the sequence. Indicates the total length of the sequence; melting temperature T m Approximate formula: ; In the formula, This indicates the number of G bases in the sequence. Indicates the number of base C in the sequence; DNA curvature: Estimated by the triplet bending model, it is defined as: ; In the formula, This represents the curvature value of the entire DNA sequence S. This represents the total length of the DNA sequence S. This indicates the starting position of the current triplet base. This represents the base at the i-th position in the sequence. This represents the base at the (i+1)th position in the sequence. θ represents the base at position i+2 of the sequence, and θ represents the bending angle corresponding to a specific triplet base combination.

[0030] Characteristics of conserved regions across species: Conservation scores (such as PhaseCons / GERP) are obtained through multiple sequence alignment (MSA). ; Where M is the number of reference species, This represents the cross-species conservation score of the i-th base position in DNA sequence S. An index indicating the reference species. Let δ represent the base at the i-th base position aligned with S in the j-th reference species, and let δ represent the indicator function.

[0031] Multi-dimensional feature integration: The final input features are composed of the following parts: ; in, This represents the final multidimensional integrated feature vector of DNA sequence S. A one-hot coding matrix representing a DNA sequence. Indicates the melting temperature of the DNA sequence. The curvature value represents the DNA sequence. A cross-species conservation score vector representing a DNA sequence; Furthermore, the process of training a CNN-RNN hybrid model based on multi-dimensional features includes: Local motif features are extracted from the convolutional layer to obtain the first feature map; The first feature map is temporally modeled using the LSTM layer to obtain the second feature map; The key site representation is obtained by weighting the second feature map using the attention layer; The key sites are represented as centromere probabilities based on the fully connected layer.

[0032] Furthermore, the process of extracting local motif features from the convolutional layer and obtaining the first feature map includes: Set the kernel size to 10 and the number of kernels to 512; The input feature map is subjected to a one-dimensional convolution operation and then activated by ReLU. The output local motif feature map is used as the first feature map.

[0033] Furthermore, the process of performing temporal modeling on the first feature map based on the LSTM layer and obtaining the second feature map includes: The number of hidden units is set to 256. The first feature map is input into the LSTM layer in sequence and the hidden states at each time step are obtained. The hidden states are then concatenated to obtain the second feature map.

[0034] Furthermore, the process of weighting the second feature map according to the attention layer and obtaining the key site representation includes: Calculate the similarity score between each position in the second feature map and the query vector; Attention weights are obtained by applying softmax normalization to the similarity scores; The attention weights are summed with the second feature map to obtain the key site representation.

[0035] Furthermore, this embodiment constructs a hybrid architecture of convolutional neural networks (CNN) and recurrent neural networks (RNN / LSTM). This architecture learns DNA sequence features and contextual dependencies, while introducing an attention mechanism to focus on key binding sites.

[0036] The CNN-RNN hybrid architecture includes: CNN layer: used to extract local motif features; ; Where σ is ReLU, This represents the convolution output feature at position i. This represents the weight matrix of the convolution kernel. This represents an input feature segment of length k starting from the i-th position. σ represents the bias term of the CNN layer, σ represents the ReLU activation function, and * represents a one-dimensional convolution operation.

[0037] RNN / LSTM layers: used to capture long-range dependencies; LSTM cell recursive computation: ; ; ; in, This represents the output of the forget gate at time t (the value ranges from 0 to 1, where 1 means the cell state at the previous time step is completely preserved, and 0 means it is completely forgotten). The weight matrix represents the forget gate. This represents the hidden state of the LSTM at time t-1 (which stores the sequence information from the previous time step). This represents the input features at time t (from the output of the CNN layer). The bias term representing the forget gate. This represents the output of the input gate at time t. This represents the weight matrix of the input gate. This represents the bias term of the input gate. This represents the state of the candidate cells at time t. The weight matrix representing the candidate cell state. Bias terms representing the candidate cell state. This represents the state of the LSTM cell at time t. This represents the cell state at time t-1. This represents the output of the output gate at time t. This represents the weight matrix of the output gate. This represents the bias term of the output gate. Let represent the hidden state of the LSTM at time t. This represents a nonlinear transformation of the current cell state; Attention mechanism: used to focus on key points; ; in, This represents the attention weight at time step t. The hidden state h represents the state at time step t. t transpose, This represents the weight vector of the attention layer. The hidden state h represents the state at the k-th time step. k transpose, This represents the "key sequences" that the model focuses on; Indicates the sequence order, used for LSTM recursive computation; Indicates the total length of the input DNA sequence; This represents the hidden layer vector output by the LSTM layer, used to represent the feature encoding of the sequence at position t.

[0038] Output layer: Sigmoid classifier; ; in, This represents the predicted centromere probability value. This represents the weight matrix of the output layer (fully connected layer) (mapping the attention feature c to the probability space). σ represents the bias term of the output layer, and σ represents the Sigmoid activation function.

[0039] The model framework is similar to BPNet: CNN is responsible for motif detection, RNN captures dependencies, and Attention focuses on key regions.

[0040] Furthermore, the process of obtaining a predictive model for centromere identification includes: The training set is divided into a training subset and a validation subset using five-fold cross-validation. Iteratively update the model parameters based on a subset of the training data; The optimal model is selected as the prediction model based on the AUROC and AUPRC of the validation subset.

[0041] Furthermore, in this embodiment The specific methods and evaluation metrics for five-fold cross-validation are as follows: The dataset was divided into 5 parts, and 4 parts were used for training and 1 part for testing in a loop. The average result was: ; in, This represents the average evaluation score of the five-fold cross-validation. Let represent the evaluation score of the i-th fold cross-validation (i=1,2,3,4,5). Evaluation indicators: AUROC: ; in, This indicates the model's ability to distinguish between "centromere regions (positive samples)" and "non-centromere regions (negative samples)". Indicates the true positive rate (sensitivity). The inverse function representing the false positive rate; AUPRC: ; in, This indicates the model's performance in scenarios with imbalanced positive and negative samples. Indicates accuracy. The inverse function of recall. ; in, This represents the number of true cases (samples that are actually centromere regions and that the model predicts to be centromere regions). This represents the number of false positives (samples that are actually non-centromere regions but the model predicts to be centromere regions). This indicates the number of false negatives (samples that are actually in the centromere region but are predicted by the model to be in the non-centromere region). Furthermore, the process of identifying target genome sequences based on prediction models includes: The target genome is segmented into continuous segments of the same length as the training set, and multi-dimensional features of each segment are extracted. By inputting multi-dimensional features into the prediction model, the centromere probability of each segment can be obtained; Segments with a centromere probability higher than a preset threshold are merged into candidate centromere regions.

[0042] Furthermore, this embodiment includes cross-species predictive analysis. Training is conducted using the source species, and testing is performed using the target species.

[0043] The feature layer remains consistent, and the input target species sequence features Ftarget are used.

[0044] Evaluate generalization performance: ; in, A multi-dimensional feature vector representing the genome sequence of the target species. θ represents the evaluation score for cross-species prediction. source This represents the model parameters obtained from training the source species. Indicated by the source species parameter θ source Defined prediction function.

[0045] Use the same AUROC / AUPRC to measure cross-species performance.

[0046] Furthermore, the method also includes: Based on the prediction model trained on the first species, the genome of the second species is directly identified, and candidate centromere regions of the second species are output, thus achieving cross-species prediction.

[0047] As a preferred implementation, this embodiment collects CENH3ChIP-seq peak sequences (200bp in length) from rice (MSU7), Arabidopsis thaliana (TAIR10), and maize (B73RefGen_v5).

[0048] Data split: 70% training, 15% validation, 15% testing.

[0049] Build a CNN-LSTM model using PyTorch, with the following parameter settings: The kernel size is 10, and the number of kernels is 512. LSTM hidden layers 256; Dropout rate 0.5; The optimizer is Adam, with a learning rate of 1e-4.

[0050] The training process is repeated 50 times, and the model that performs best on the validation set is selected.

[0051] As another preferred implementation method, such as Figure 2 As shown, this embodiment focuses on cross-species prediction, using a model trained on Arabidopsis thaliana and rice to predict the centromere region of the soybean genome; the prediction results are compared with existing soybean CENH3ChIP-seq data to calculate AUROC. Here, CENH3-Xfer is the name of the program designed in this embodiment, and quarTeT is the name of a similar program; Figure 2In this context, 'a' represents accuracy, which refers to the number of centromeres of the species that are identified; 'b' represents precision, which is the proportion of correctly identified positive samples out of all predicted positive samples; and 'c' represents recall, which is the proportion of correctly identified positive samples out of all actual positive samples. By comparing the program CENH3-Xfer designed in this embodiment with the similar program quarTeT, the results show that the algorithm innovation of this embodiment significantly improves the centromere prediction accuracy compared to quarTeT, demonstrating the effectiveness of the innovative algorithm designed in this embodiment.

[0052] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for automatic identification of centromeres in crops, characterized in that, include: A training set was constructed based on the validated CENH3 binding region sequence and random noncentromere region sequences; Extract multi-dimensional features of DNA sequences from the training set; A CNN-RNN hybrid model is trained based on the multi-dimensional features to obtain a predictive model for centromere identification; The target genome sequence is identified based on the prediction model to obtain candidate centromere regions.

2. The method according to claim 1, characterized in that, The process of constructing a training set based on verified CENH3 binding region sequences and random noncentromere region sequences includes: A fixed-length positive sample was extracted from the published CENH3 ChIP-seq peak sequences of the species. A noncentromere region of the same length as the positive sample was randomly selected from the reference genome of the same species as a negative sample. The positive and negative samples are merged to obtain a labeled training set.

3. The method according to claim 1, characterized in that, The process of extracting multi-dimensional features of DNA sequences from the training set includes: The DNA sequence is converted into a numerical matrix according to One-hot encoding; DNA secondary structure characteristics were calculated based on AT content, melting temperature, and curvature. Cross-species conservation scores were obtained based on multiple sequence alignment. The numerical matrix, DNA secondary structure features, and cross-species conservation scores are concatenated into a unified feature vector.

4. The method according to claim 1, characterized in that, The process of training a CNN-RNN hybrid model based on the multi-dimensional features includes: Local motif features are extracted from the convolutional layer to obtain the first feature map; The first feature map is temporally modeled using the LSTM layer to obtain the second feature map; The key site representation is obtained by weighting the second feature map using the attention layer; The key sites are represented as centromere probabilities based on the fully connected layer.

5. The method according to claim 4, characterized in that, The process of extracting local motif features from convolutional layers and obtaining the first feature map includes: Set the kernel size to 10 and the number of kernels to 512; The input feature map is subjected to a one-dimensional convolution operation and then activated by ReLU. The output local motif feature map is used as the first feature map.

6. The method according to claim 4, characterized in that, The process of performing temporal modeling on the first feature map based on the LSTM layer and obtaining the second feature map includes: The number of hidden units is set to 256. The first feature map is input into the LSTM layer in sequence and the hidden state at each time step is obtained. The hidden states are then concatenated to obtain the second feature map.

7. The method according to claim 4, characterized in that, The process of weighting the second feature map using the attention layer to obtain the key site representation includes: Calculate the similarity score between each position in the second feature map and the query vector; The similarity scores are then normalized using softmax to obtain the attention weights. The attention weights are weighted and summed with the second feature map to obtain the key site representation.

8. The method according to claim 1, characterized in that, The process of obtaining a predictive model for centromere identification includes: The training set is divided into a training subset and a validation subset using five-fold cross-validation. Iteratively update the model parameters based on the training subset; The optimal model is selected as the prediction model based on the AUROC and AUPRC of the validation subset.

9. The method according to claim 1, characterized in that, The process of identifying the target genome sequence based on the prediction model includes: The target genome is segmented into continuous segments of the same length as the training set, and multi-dimensional features of each segment are extracted. The multidimensional features are input into the prediction model to obtain the centromere probability of each segment; Segments with a centromere probability higher than a preset threshold are merged into candidate centromere regions.

10. The method according to claim 1, characterized in that, The method further includes: Based on the prediction model trained on the first species, the genome of the second species is directly identified, and candidate centromere regions of the second species are output, thus achieving cross-species prediction.

Images