A method for automatically identifying centromeres of crops

Abstract

The application discloses a kind of crop centromere automatic identification method, belong to crop centromere identification field, comprising: according to the verified CENH3 binding area sequence and random non-centromere area sequence constructs training set;According to the multidimensional feature of DNA sequence extracted from training set;According to the multidimensional feature training CNN-RNN hybrid model and obtains the prediction model for identifying centromere;According to the prediction model, the target genomic sequence is identified, and the candidate centromere region is obtained.The application no longer depends on repeat sequence, but directly uses CENH3 binding feature, improves prediction accuracy.And can be applied across species, improve the universal applicability of centromere region identification.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of crop centromere recognition, and particularly relates to a crop centromere automatic recognition method. BACKGROUND

[0002] For crop centromere recognition, the existing technologies mainly include experimental verification means, sequence feature-based prediction methods and deep learning models developed in recent years. These technologies have promoted the development of centromere research to some extent, but still have many limitations, which restrict the process of crop molecular breeding.

[0003] For example, the CENH3 chromatin immunoprecipitation (ChIP-seq) technology as an experimental verification means can directly verify the binding of CENH3 on the centromere of the genome, but the experimental cost is high and the operation is complex.

[0004] The sequence feature-based prediction method indirectly judges by identifying the position of the tandem repeat sequence. However, taking corn as an example, existing studies have shown that the centromere region and the tandem repeat region do not always coincide, indicating that the prediction method based on the repeat sequence has the problems of low accuracy and poor generalization.

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

[0006] The existing technologies have significant limitations in accurately predicting the centromere position in multiple species, especially when the repeat sequence signal is weak or does not coincide with the functional region, the prediction accuracy is difficult to guarantee, so the repeat sequence method is unreliable, especially in plants such as corn; at the same time, there is a lack of scalable and general prediction model that combines sequence features and protein binding characteristics, which is difficult to adapt to the genomic characteristics of different crops; in addition, the existing methods cannot quickly locate the centromere region in newly assembled crop genomes (such as rice, wheat, corn, etc.), which directly restricts the process of fine mapping of important agronomic trait genes, chromosome engineering breeding and genome selection breeding, and becomes a key technical obstacle to break through the molecular breeding bottlenecks of crop yield improvement and stress resistance improvement. SUMMARY

[0007] To solve the above technical problems, the application provides a crop centromere automatic recognition method, which comprises:

[0008] constructing a training set according to the verified CENH3 binding region sequence and the random non-centromere region sequence;

[0009] extracting multi-dimensional features of the DNA sequences according to the training set;

[0010] training a CNN-RNN hybrid model according to the multi-dimensional features and obtaining a prediction model for identifying centromeres;

[0011] identifying a target genomic sequence according to the prediction model and obtaining a candidate centromere region.

[0012] Preferably, the process of constructing the training set according to the verified CENH3 binding region sequences and random non-centromere region sequences comprises:

[0013] cutting positive samples of a fixed length according to the CENH3 ChIP-seq peak sequences of published species;

[0014] randomly selecting non-centromere regions equal in length to the positive samples in the reference genome of the same species as negative samples;

[0015] merging the positive samples and the negative samples to obtain a labeled training set.

[0016] Preferably, the process of extracting multi-dimensional features of the DNA sequences according to the training set comprises:

[0017] converting the DNA sequences into numerical matrices according to One-hot encoding;

[0018] calculating DNA secondary structure features according to AT content, melting temperature, and curvature;

[0019] obtaining cross-species conservation scores according to multiple sequence alignment;

[0020] concatenating the numerical matrices, the DNA secondary structure features, and the cross-species conservation scores into a unified feature vector.

[0021] Preferably, the process of training a CNN-RNN hybrid model according to the multi-dimensional features comprises:

[0022] extracting local motif features according to a convolutional layer and obtaining a first feature map;

[0023] performing time-series modeling on the first feature map according to an LSTM layer and obtaining a second feature map;

[0024] weighting the second feature map according to an attention layer and obtaining key site representations;

[0025] mapping the key site representations to centromere probabilities according to a fully connected layer.

[0026] Preferably, the process of extracting local motif features according to a convolutional layer and obtaining a first feature map comprises:

[0027] The convolution kernel size is set to 10, and the number is 512.

[0028] One-dimensional convolution operation is performed on the input feature map and ReLU activation is performed, and the local motif feature map is output as the first feature map.

[0029] Preferably, the process of performing time series modeling on the first feature map according to the LSTM layer and obtaining the second feature map comprises:

[0030] The number of hidden units is set to 256, the first feature map is input into the LSTM layer in sequence order, and the hidden state of each time step is obtained, and the hidden state is spliced to obtain the second feature map.

[0031] Preferably, the process of weighting the second feature map according to the attention layer and obtaining the key site representation comprises:

[0032] The similarity score of each position of the second feature map and the query vector is calculated;

[0033] The similarity score is subjected to softmax normalization to obtain the attention weight;

[0034] The attention weight is weighted and summed with the second feature map to obtain the key site representation.

[0035] Preferably, the process of obtaining the prediction model for identifying the centromere comprises:

[0036] The training set is divided into a training subset and a validation subset by five-fold cross-validation;

[0037] The model parameters are iteratively updated according to the training subset;

[0038] The optimal model is selected as the prediction model according to the AUROC and AUPRC of the validation subset.

[0039] Preferably, the process of identifying the target genomic sequence according to the prediction model comprises:

[0040] The target genome is segmented into continuous fragments of the same length as the training set, and multi-dimensional features of each fragment are extracted;

[0041] The multi-dimensional features are input into the prediction model to obtain the centromere probability of each fragment;

[0042] The fragments with a centromere probability higher than a preset threshold are combined into candidate centromere regions.

[0043] Preferably, the method further comprises:

[0044] According to the prediction model trained on the first species, the genome of the second species is directly identified, and the candidate centromere region of the second species is output, realizing cross-species prediction.

[0045] Compared with the prior art, the present application has the following advantages and technical effects:

[0046] The present application no longer depends on repetitive sequences, but directly uses the CENH3 binding characteristics, for the first time realizing end-to-end automatic identification of crop centromeres at the computational level, outputting candidate regions without any experimental steps, and fundamentally breaking through the bottleneck of traditional experiments with high throughput, high cost and long cycle, thereby improving the prediction accuracy.

[0047] The method of the present application can be applied across species, and is still accurate in plants such as corn with weak repetitive signals. Compared with the existing prediction method which depends on repetitive sequences, the present application can improve the accuracy by 30-40% on a corn dataset, and has universality. And for a new species, only the genomic sequence needs to be input, and the candidate centromere region can be quickly obtained.

[0048] The present application uses three parallel sequences of One-hot, secondary structure and conservation, expands the original single information that only depends on tandem repeats to a three-dimensional description of "sequence-structure-evolution", and can provide sufficient discrimination signals even for complex genomes such as corn with weak or missing repetitive signals.

[0049] The hybrid model involved in the present application has a three-level progressive structure of capturing local motifs by CNN, capturing long-range dependencies by LSTM, and focusing on key sites by Attention, so that the model can not only identify short motif combinations, but also integrate remote context, thereby comprehensively improving the sensitivity and specificity of plant CENH3 binding sites. BRIEF DESCRIPTION OF DRAWINGS

[0050] The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application and are incorporated in and constitute a part of this application. The embodiments of this application and their description are used to explain the application and are not intended to limit the application. In the drawings:

[0051] Figure 1 The flowchart of the method of the embodiments of the present application is shown in the figure;

[0052] Figure 2 The system structure schematic diagram of the embodiments of the present application is shown in the figure. DETAILED DESCRIPTION

[0053] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0054] It should be noted that the steps shown in the flowchart of 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 herein can be executed in an order different from that shown.

[0055] As shown in the flowchart, the embodiment provides a method for automatically identifying centromeres of crops, comprising: Figure 1

[0056] According to the verified CENH3 binding region sequence and the random non-centromere region sequence, a training set is constructed;

[0057] According to the training set, multi-dimensional features of the DNA sequence are extracted;

[0058] According to the multi-dimensional features, a CNN-RNN hybrid model is trained to obtain a prediction model for identifying centromeres;

[0059] According to the prediction model, the target genomic sequence is identified to obtain a candidate centromere region.

[0060] Further, the process of constructing the training set according to the verified CENH3 binding region sequence and the random non-centromere region sequence comprises:

[0061] According to the CENH3 ChIP-seq peak sequence of the published species, a positive sample of a fixed length is intercepted;

[0062] In the same species reference genome, a non-centromere region equal in length to the positive sample is randomly selected as a negative sample;

[0063] The positive sample and the negative sample are combined to obtain a labeled training set.

[0064] Further, the embodiment is to construct a training set related to centromere prediction, starting from published species such as potato, sweet potato, rice, corn, soybean, and wheat. On the one hand, the centromere binding region sequence verified by CENH3 ChIP-seq technology is collected as a positive sample, and on the other hand, a non-centromere region is randomly selected in the genome of these species as a negative sample. At the same time, the reference genome sequences of the above-mentioned species are further collected, and finally integrated to form a training set for subsequent research.

[0065] Further, the process of extracting multi-dimensional features of the DNA sequence according to the training set comprises:

[0066] According to One-hot encoding, the DNA sequence is converted into a numerical matrix;

[0067] According to the AT content, the melting temperature, and the curvature, the DNA secondary structure features are calculated;

[0068] ​The cross-species conservation score is obtained according to the multiple sequence alignment;

[0069] The numerical matrix, the DNA secondary structure feature and the cross-species conservation score are spliced into a unified feature vector.

[0070] Further, the embodiment realizes the integration and extraction of the multi-dimensional features of the DNA sequence by One-hot encoding of the DNA sequence, introducing the DNA secondary structure features such as AT content and melting temperature, and fusing the information of the conserved regions identified by the multiple sequence alignment.

[0071] Specifically, the corresponding method and formula of the digital representation of the DNA sequence are as follows:

[0072] The DNA sequence is composed of four bases (A, C, G and T), and when One-hot encoding is used, each base is mapped to a 4-dimensional sparse vector:

[0073] ;

[0074] Given a DNA sequence with a length of L:

[0075] ;

[0076] In the formula, represents the vector representation of the whole DNA sequence with a length of L, represents the i-th base (i = 1, 2, …, L) in the sequence, represents the sequence length (i.e. the last base);

[0077] Its One-hot representation is a matrix:

[0078] ;

[0079] In the formula, represents the One-hot encoding matrix of the whole DNA sequence, represents the One-hot vector of the last base (the L-th base) in the sequence, represents the One-hot vector of the i-th base in the sequence;

[0080] The formula of the DNA secondary structure feature is as follows:

[0081] AT content:

[0082] ;

[0083] In the formula, Number of base A in the sequence, Number of base T in the sequence, Total length of the sequence;

[0084] Melting temperature T m Approximate formula:

[0085] ;

[0086] Where, Number of base G in the sequence, Number of base C in the sequence;

[0087] Curvature (DNA curvature):

[0088] Estimated by the triplex bending model, defined as:

[0089] ;

[0090] Where, Curvature value of the whole DNA sequence S, Total length of the DNA sequence S, Starting position of the current triplex base, Base at the i-th position of the sequence, Base at the i+1-th position of the sequence, Base at the i+2-th position of the sequence, and θ represents the bending angle corresponding to the specific triplex base combination.

[0091] Cross-species conservation region features:

[0092] Obtained by multiple sequence alignment (MSA, Multiple Sequence Alignment) to obtain conservation scores (such as PhastCons / GERP):

[0093] ;

[0094] Where, M is the number of reference species, Cross-species conservation score of the i-th base position of the DNA sequence S, Index of the reference species, Base of the j-th reference species at the i-th base position aligned with S, and δ represents the indicator function.

[0095] Multi-dimensional feature integration:

[0096] The final input features are spliced from the following parts:

[0097] ;

[0098] wherein, represents the final multi-dimensional integrated feature vector of DNA sequence S, represents the One-hot encoding matrix of DNA sequence, represents the melting temperature of DNA sequence, represents the curvature value of DNA sequence, represents the cross-species conservation score vector of DNA sequence;

[0099] Further, the process of training the CNN-RNN hybrid model according to the multi-dimensional features comprises:

[0100] extracting local motif features according to the convolutional layer and obtaining a first feature map;

[0101] performing time-series modeling on the first feature map according to the LSTM layer and obtaining a second feature map;

[0102] weighting the second feature map according to the attention layer and obtaining a key site representation;

[0103] mapping the key site representation to a centromere probability according to the fully connected layer.

[0104] Further, the process of extracting local motif features according to the convolutional layer and obtaining a first feature map comprises:

[0105] setting the convolution kernel size to 10 and the number to 512;

[0106] performing one-dimensional convolution operation on the input feature map and passing through ReLU activation, outputting the local motif feature map as the first feature map.

[0107] Further, the process of performing time-series modeling on the first feature map according to the LSTM layer and obtaining a second feature map comprises:

[0108] setting the number of hidden units to 256, inputting the first feature map into the LSTM layer in sequence order and obtaining the hidden state of each time step, and concatenating the hidden state to obtain the second feature map.

[0109] Further, the process of weighting the second feature map according to the attention layer and obtaining a key site representation comprises:

[0110] calculating the similarity score of each position of the second feature map with the query vector;

[0111] performing softmax normalization on the similarity score to obtain the attention weight;

[0112] performing weighted summation of the attention weight and the second feature map to obtain the key site representation.

[0113] 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.

[0114] The CNN-RNN hybrid architecture includes:

[0115] CNN layer: used to extract local motif features;

[0116] ;

[0117] 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.

[0118] RNN / LSTM layers: used to capture long-range dependencies;

[0119] LSTM cell recursive computation:

[0120] ;

[0121] ;

[0122] ;

[0123] 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 the cell state 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;

[0124] Attention mechanism: used to focus on key points;

[0125] ;

[0126] 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.

[0127] Output layer: Sigmoid classifier;

[0128] ;

[0129] 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.

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

[0131] Furthermore, the process of obtaining a predictive model for centromere identification includes:

[0132] The training set is divided into a training subset and a validation subset using five-fold cross-validation.

[0133] Iteratively update the model parameters based on a subset of the training data;

[0134] The optimal model is selected as the prediction model based on the AUROC and AUPRC of the validation subset.

[0135] Furthermore, in this embodiment

[0136] The specific methods and evaluation metrics for five-fold cross-validation are as follows:

[0137] 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:

[0138] ;

[0139] 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).

[0140] Evaluation indicators:

[0141] AUROC:

[0142] ;

[0143] 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;

[0144] AUPRC:

[0145] ;

[0146] in, This indicates the model's performance in scenarios with imbalanced positive and negative samples. Indicates accuracy. The inverse function of recall.

[0147] ;

[0148] 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).

[0149] Furthermore, the process of identifying target genome sequences based on prediction models includes:

[0150] 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.

[0151] By inputting multi-dimensional features into the prediction model, the centromere probability of each segment can be obtained;

[0152] Segments with a centromere probability higher than a preset threshold are merged into candidate centromere regions.

[0153] Furthermore, this embodiment includes cross-species predictive analysis.

[0154] Training is conducted using the source species, and testing is performed using the target species.

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

[0156] Evaluate generalization performance:

[0157] ;

[0158] 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.

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

[0160] Furthermore, the method also includes:

[0161] 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.

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

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

[0164] Build a CNN-LSTM model using PyTorch, with the following parameter settings:

[0165] The kernel size is 10, and the number of kernels is 512.

[0166] LSTM hidden layers 256;

[0167] Dropout rate 0.5;

[0168] The optimizer is Adam, with a learning rate of 1e-4.

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

[0170] 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 2 In 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.

[0171] 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; 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. 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. This indicates the total length of the sequence.

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 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.

4. The method according to claim 3, 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.

5. The method according to claim 3, 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.

6. The method according to claim 3, 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.

7. 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.

8. 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.

9. 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