Prediction method based on combined feature coding and DNA binding sites

Abstract

The invention discloses a prediction method based on combined feature coding and DNA binding sites, and belongs to the field of site prediction. The method comprises the following steps of: eliminating low variance features of an initial data set and setting a recombination rate threshold marking sample; determining a task type according to the input label and using a classification or regression algorithm; optimizing the parameters by using Optuna and distributing the parameters to corresponding parameters; randomly distributing training and testing data sets; respectively inputting the data set into a plurality of algorithms including a random forest, a deep forest and the like for training; selecting a model with the best performance as a final prediction model according to the evaluation indexes; and finally, predicting the test set through the optimal model. According to the method, excellent high precision and low error are achieved, meanwhile, better portability is achieved, and a new thought is provided for research of an efficient recombination system in drug design.

Description

Technical Field

[0001] The invention belongs to the field of site prediction, and mainly relates to a prediction method for DNA binding sites, and specifically relates to a prediction method based on combined feature coding and DNA binding sites. Background Art In bioinformatics, the interaction between DNA and proteins is an important component of cell function, especially in cell growth, differentiation and gene regulation. Gene regulation mechanisms rely on the binding of transcription factors to specific DNA sequences, which directly affects the transcription level of genes and thus regulates gene expression. Identifying DNA binding sites, which usually consist of 5 to 20 base pairs, can help us gain a deeper understanding of these regulatory mechanisms and play an important role in drug development and biomedical research. However, traditional experimental methods, such as chromatin immunoprecipitation sequencing (ChIP-seq), although they can provide high-resolution binding site information, are limited in their widespread application due to their high resource consumption and long time requirements.

[0002] Existing computational prediction methods usually use position weight matrices (PWMs) or motif recognition techniques, but these methods fail to fully capture the complexity of DNA binding. The specificity of binding sites is also affected by multiple factors such as neighboring sequences, GC content, and DNA shape. With the rapid development of deep learning technology in fields such as image processing and natural language processing, its application in bioinformatics is also increasing. By using deep learning architectures such as convolutional neural networks (CNNs), researchers can efficiently extract important features from massive DNA sequences and significantly improve prediction accuracy. In addition, machine learning has also shown good prospects in the prediction of DNA binding sites, being able to handle complex data structures and effectively extract key information, especially in small sample cases. Although these algorithms have improved in accuracy, they still face some challenges, especially when dealing with complex biological features, they often rely on a single feature or model, which limits their scope of application and accuracy. Therefore, researchers are actively exploring new methods such as ensemble learning and feature optimization to enhance the reliability and accuracy of predictions and promote further development in this field. Summary of the invention

[0003] In view of the shortcomings of existing prediction methods for DNA binding sites, the present invention provides a prediction method OptimDase based on combined feature coding and DNA binding sites. OptimDase uses a combined feature representation method to better represent the potential features of DNA sequences, and also combines the three algorithms of deep forest, random forest and XGBoost to improve the prediction accuracy. Its purpose is to solve the problems of unsatisfactory results and time-consuming training process in the current prediction methods for DNA binding sites.

[0004] The present invention is achieved by the following measures: A prediction method based on combined feature coding and DNA binding sites comprises the following steps: (1) Initial dataset of DNA binding sites , 1≤j≤k, k represents the total number of samples in the initial data set, where X j represents a single DNA fragment in the dataset, Y j Indicates whether this DNA sequence fragment can become a binding site. Its value is 1 or 0, representing a binding site and a non-binding site, respectively. On this basis, the data set is expanded. In order to obtain the sequence variant of each data, the sequence variants of each data sequence are obtained by performing reverse, complementary and complementary reverse operations through basic binary operations, and then they are spliced ​​with the original data set to finally obtain the data set. ; The k is a natural number greater than or equal to 1; (2) Labeling the dataset by defining a recombination rate threshold and adding a feature column Whether the sample in is a positive sample or a negative sample, where the feature column is named "Class"; attC sites with a recombination rate greater than the threshold are marked as 1, which is a positive sample, and those less than the threshold are marked as 0, which is a negative sample; (3) Through the one-hot encoding table For the dataset The sequence information in is converted into feature vector C1, and then feature vector C2 is obtained by combining feature encoding representation. The feature vectors C1 and C2 are concatenated to obtain the combined feature vector C. Encode the categorized results again; (4) The dataset after feature representation in step (3) 60~95% of the samples in the training set D x , the remaining data is used as the test set D y ; (5) Determine the task type based on the input label Flag: If the task is a classification task, use the classification algorithm, otherwise use the regression prediction algorithm; then, use Optuna to optimize the parameters of random forest and XGBoost, and reallocate the optimized parameters to the corresponding model; (6) Decision tree analysis D x Get the weight vector W={ }, 1≤i≤c, where c represents the total number of features, and the calculation formula is as follows: Where c is the total number of features, is the value of the feature in the i-th column of the weight vector W, and the calculation formula is as follows: In the formula is the value of the t-th tree node, T is the total number of decision trees, and the calculation formula is as follows: In the formula and is the Gini coefficient of all child nodes marked as 0 and child nodes marked as 1; is the Gini coefficient, and the calculation formula is as follows: Where N is D x Sample size, Respectively represent the number of samples marked as “0” and “1” in the node branch; (7) Using the weighted scanning method to x The corresponding feature vector C of the sample in is processed, and the execution steps are as follows: the feature vectors are processed respectively through a sliding window of length m and weight vector Slide, and the window sliding step size is set to p during the sliding process. After sliding, the feature matrix f is obtained respectively (1,n-m+1,m) and the weight matrix w (1,n-m+1,m) ; The feature matrix f (1,n-m+1,m) With the weight matrix w (1,n-m+1,m) Multiply to get the pre-input matrix tf (1,n-m+1,m) , where n represents the number of base codes 0 or 1, and m is the length of the sliding window, ranging from 10 to 1000; (8) tf (1,n-m+1,m) Input into the full random forest and random forest, through the above operations, the feature matrix fa is obtained respectively (n-m+1,2) and FB (n-m+1,2) ; Finally, the feature matrix fa captured by two random forests (n-m+1,2) and FB (n -m+1,2) Concatenate into a one-dimensional feature vector F (2*2(n-m+1)) , the feature vector is input into the cascade forest for training, and a prediction model for DNA binding site classification is obtained; (9) The training set D x and the test set D y The models are input into XGBoost and Random Forest algorithms for training respectively, and the best performing model is selected as the final prediction model M according to the evaluation indicators; Finally, the test set is input into model M for prediction.

[0005] In the above prediction method, preferably, k in step (1) is a natural number between 1000 and 50000.

[0006] In the above prediction method, preferably, the recombination rate threshold of step (2) is 0.46, the number of positive samples screened is 1762, and the number of negative samples is 11117.

[0007] In the above prediction method, preferably, in the combined feature encoding method in step (3), the feature length Length is given by the formula Obtain, where b represents the length of the combined bases and takes a value of 3. The bases A, T, C, and G can be randomly selected to form a sequence set by permuting three repeatable bases, each of which has a single length of 3 bp. Each element in the set C is set as a feature column, for a total of 64 feature columns, and its elements are the feature names of the feature columns; The feature vector C2 is calculated as follows: first, 64 feature columns are constructed with 3 bases as a group, and then a window with a length of 3 is used to slide the entire DNA sequence with a step size of 1, and the feature columns corresponding to the 3 base groups in the window are marked as 1. The feature vector C2 obtained by combining feature encoding has a length of 64.

[0008] In the above prediction method, preferably, the number of samples in the training set of step (4) accounts for 80% of the total number of samples in the data set.

[0009] In the above prediction method, preferably, in step (6), the total number of decision trees is 462, and the maximum depth of each decision tree is 11; and c represents the total number of features.

[0010] In the above prediction method, preferably, in step (7), the value of m is 50 and the value of p is 1.

[0011] Compared with the prior art methods, the present invention has the following beneficial effects: Since the identification of DNA binding sites is a meaningful and difficult task, we proposed an algorithm OptimDase to identify DNA binding sites and predict the recombination frequency of DNA recombination sites. The algorithm improves the accuracy of classification results by extracting the combined features of bases and their adjacent bases. At the same time, the OptimDase algorithm is based on the optimal decision and is implemented by using XGBoost, Deep Forest and Random Forest algorithms, which effectively improves the prediction ability and avoids the shortcomings of a single model. Among them, the Deep Forest algorithm is optimized by the weighted scanning method, and its prediction results are improved on all data sets. Studies have shown that the prediction performance of OptimDase is better than other existing algorithms. More importantly, it overcomes the problem that existing algorithms only focus on low prediction accuracy. In short, OptimDase has significant anti-interference ability and good prediction accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 Flow chart of the prediction method based on combined feature coding and DNA binding sites; Figure 2 Expand the example for the DNA sequence dataset; Figure 3 is the one-hot encoding correspondence table of DNA sequences; Figure 4 A feature representation for converting DNA sequence data into a combined feature encoding; Figure 5 Comparison of the prediction performance of mainstream algorithms on TFBSs site datasets; Figure 6 Comparison of the prediction performance results of mainstream algorithms on the attC site dataset. DETAILED DESCRIPTION

[0013] In order to clearly illustrate the technical solution of the present invention, the present invention is described below in conjunction with the accompanying drawings (1-6) and examples. The examples herein are only used to explain the present invention and do not limit the present invention.

[0014] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present invention belongs.

[0015] Figure 1 The model workflow for identifying DNA binding sites using the OptimDase method is demonstrated. First, the DNA sequence or related feature data is input. Next, the embedding layer converts the input data into a format suitable for further analysis. The system then determines whether the task is a classification or regression problem, and selects the corresponding algorithm based on this decision: if it is a classification task, a classification algorithm is used; if it is a regression task, a regression algorithm is used. After determining the type of task, the model performs parameter optimization to improve performance. The optimized models include random forests, deep forests, and XGBoost, and all models are evaluated to determine their effectiveness. Finally, the model that performs best in the evaluation is selected for prediction. At the end of the entire process, the prediction results generated by the selected model are output, emphasizing the integration of algorithms and optimization techniques to improve prediction accuracy. The data sets selected in this embodiment are the SP1 binding site data set and the attC site data set.

[0016] The input file type of the OptimDase method is CSV type. The Origin.csv file contains two mutant attCr0 sequences, a total of 12,879 attCr0 mutants and 2,400 samples of SP1 binding sites, recorded as the initial data set D. Each attCr0 mutant contains 292 structural features, of which 9 are global features and 283 are basic features. Each SP1 binding site data contains a DNA sequence and its label (to determine whether it is a DNA binding site). Initial data preprocessing is performed based on this data set; the output file of the OptimDase method contains an output type file. The Result.output file is the prediction result of each site in the test set output by the OptimDase method.

[0017] The OptimDase prediction method can be divided into the following steps: 1. Data preprocessing In this example, for the attCr0 mutant database, we first filtered out useless features and features with small variance. Here, 14 features with zero or small variance were deleted, where features with small variance were defined as values ​​with more than 80% of the values ​​being the same in individual features. The number of features remaining after processing was 278. Next, the feature values ​​were linearly normalized to [0,1]. Finally, we defined a recombination rate threshold (0.46) and added a "Class" feature column to mark whether the sample is a positive sample. attC sites with a recombination rate greater than the threshold are marked as 1 (positive sample), and those less than the threshold are marked as 0 (negative sample). Finally, we screened out 1762 positive samples and 11117 negative samples, and constructed a balanced dataset using oversampling operations. Specifically, we obtained 11117 positive samples on 1762 positive samples by replacement sampling, and formed a balanced dataset containing 22234 samples with 11117 negative samples.

[0018] This example pre-targets the SP1 binding site dataset Process the data. Since the number of DNA sequence samples is limited, we first need to expand the dataset by adding some additional sequences. Specifically, we calculated the inversion sequence, complementary sequence, and reverse complementary sequence for each DNA sequence. These three newly calculated DNA sequences will be used as a new data. The label of the new data is the same as the original DNA sequence. After the above operations, the number of positive and negative samples in the TFBS dataset is expanded from 1200 to four times the number. ( Figure 2 ). Then, through the one-hot encoding table Transform the DNA sequence in the data set D* to obtain the feature vector C1( Figure 3). We then encode each base in association with the context base to achieve multi-base feature encoding. This encoding method first constructs 64 feature columns in groups of 3 bases, then uses a window of length 3 to slide the entire DNA sequence with a step size of 1 (different from the window of weighted multi-granularity scanning), and marks the feature columns corresponding to the 3 base groups in the window as 1. This constitutes the feature vector C2( Figure 4 ). Then concatenate C1 and C2 to get the feature vector C.

[0019] In this embodiment, the data is preprocessed. It contains a total of 9600 samples, of which half are positive and half are negative. Each sample data consists of 120 feature fields and 1 result category. Finally, the positive and negative samples are rearranged in order.

[0020] 2. Dataset Division The dataset after feature representation in step (3) 80% of the data is used as the training set Dx, and the remaining 20% ​​of the data is used as the test set Dy.

[0021] 3. Task classification The task type is determined based on the input label Flag. If the task is a classification task, the classification algorithm is used, otherwise the regression prediction algorithm is used. Then, Optuna is used to optimize the parameters of random forest and XGBoost, and the optimized parameters are reassigned to the corresponding models.

[0022] 4. Calculate weights Through the decision tree, the processed training set D x Get the weight vector W. According to the formula Get the Gini coefficient of the node , where N is D x The total number of samples, N postive and N negative Respectively represent the number of samples marked as "0" and "1" in the node branch. According to the formula Computing the value of nodes , where and is the Gini coefficient of all child nodes marked as 0 and child nodes marked as 1 under this node, is the Gini coefficient of the node itself. (t) Evaluate the value of the feature in column i , where T is the total number of decision trees. According to the formula Calculate the weight W of each feature i , where is the feature value of the i-th column, and c is the total number of features.

[0023] In this example, the top 10 site features and their weights in attC sites and TFBSs sites are: 5. Weighted scanning process The training set D is scanned by weighted x The corresponding feature vector C of the sample in is processed, and the execution steps are as follows: the feature vectors are processed respectively through a sliding window of length m and weight vector Slide, the sliding step is 1. After sliding, the feature matrix f is obtained respectively (1,n-m+1,m) and the weight matrix w (1,n-m+1,m) . The feature matrix f (1,n-m+1,m) With the weight matrix w (1,n-m+1,m) Multiply to get the pre-input matrix tf (1,n-m+1,m) , n represents the number of base codes 0 or 1. m is the length of the sliding window, which is 50.

[0024] 6. Predict DNA binding sites tf (1,n-m+1,m) Input into the full random forest and random forest. Through the above operations, the feature matrix fa is obtained respectively (n-m+1,2) and FB (n-m+1,2) Finally, the feature matrix fa captured by the two random forests (n-m+1,2) and FB (n-m+1,2) Concatenate into a one-dimensional feature vector F (2*2(n-m+1)) , input the feature vector into the cascade forest for training, and obtain a prediction model for DNA binding site classification. Then the training set D x and the test set D y They are respectively input into the XGBoost and random forest algorithms for training, and based on the evaluation indicators, the best performing model is selected as the final prediction model M. Finally, the test set is input into model M for prediction.

[0025] 7. Performance Evaluation We use multiple evaluation methods to represent the performance of the algorithm in the classification task: including accuracy, F1 value (F1 Measure) and AUC, which are expressed by the formula ,formula The prediction accuracy and F1 value are obtained. These indicators are calculated by comparing the confusion matrix of the actual category and the model prediction category. In the confusion matrix, TP, TN, FP and FN represent the number of true positive samples, true negative samples, false positive samples and false negative samples respectively. P and N represent the number of samples whose prediction results are positive samples and negative samples respectively. The precision value and recall value are calculated by the formula and formula Calculated; the accuracy ranges from 0 to 1. The closer the value is to 1, the more samples are correctly classified, and vice versa. The larger the F1 value, the closer the algorithm is to the expected value. The ROC curve is composed of the probability of true positive samples on the Y axis and the probability of false positive samples on the X axis. The area called AUC value is the area between the ROC curve and the coordinate axis, which provides an assessment of the model's ability. The larger the AUC value, the stronger the model's performance. The above formula can be used to obtain the test set D x The accuracy, F1 value and AUC are 0.8943, 0.8920 and 0.9219 respectively.

[0026] We compare it with some existing machine learning algorithms, such as Random Forest and LightGBM. The hyperparameters of all algorithms are optimized by the Optuna framework based on 5-fold cross validation. Experimental results ( Figure 5 ) shows the prediction performance results of all algorithms on the TFBSs dataset.

[0027] In the regression prediction experiment, we also selected four evaluation indicators to evaluate the performance of the algorithm: MAE (mean absolute error), RMSE (root mean square error), PCC (Pearson correlation coefficient) and VarScore (variance score). The calculation formula of MAE is , which represents the mean absolute difference between the predicted value and the true value. The calculation formula of RMSE is , which is used to indicate the distance between the predicted value and the true value. PCC indicates the correlation between the predicted result and the true result, and its calculation method is The calculation formula of VarScore is This indicator measures the impact of data set fluctuations on the algorithm. V(·) represents the variance calculation function. Experimental results ( Figure 6 ) shows the prediction performance results of all algorithms on the attC dataset.

[0028] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned implementation cases, those skilled in the art can still modify the technical solutions recorded in the aforementioned implementation cases or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A prediction method based on combined feature coding and DNA binding sites, characterized in that: The steps include: (1) Basic data set of DNA binding sites , 1≤j≤k, k represents the total number of samples in the initial data set, where X j is a single DNA fragment in the dataset, Y j Whether the DNA sequence fragment can be used as a binding site, its value is 1 or 0, representing the binding site and non-binding site respectively. On this basis, the data set is expanded. In order to obtain the sequence variant of each data, the sequence variants of each data sequence are obtained by performing reverse, complementary and complementary reverse operations through binary basic operations, and then spliced ​​with the original data set to finally obtain the data set. ; The k is a natural number greater than or equal to 1; (2) Labeling the dataset by defining a recombination rate threshold and adding a feature column Whether the sample in is a positive sample or a negative sample, where the feature column is named "Class"; attC sites with a recombination rate greater than the threshold are marked as 1, which is a positive sample, and those less than the threshold are marked as 0, which is a negative sample; (3) Through the one-hot encoding table For the dataset The sequence information in is converted to obtain the feature vector C1, and then the feature vector C2 is obtained by combining the feature encoding representation, and the combined feature vector C is obtained by concatenating C1 and C2. Encode the categorized results again; (4) The dataset after feature representation in step (3) 60~95% of the samples in the training set D x , the remaining data is used as the test set D y ; (5) Determine the task type based on the input label Flag: If the task is a classification task, use the classification algorithm, otherwise use the regression prediction algorithm; then, use Optuna to optimize the parameters of random forest and XGBoost, and reallocate the optimized parameters to the corresponding model; (6) Decision tree analysis D x Get the weight vector W={ }, 1≤i≤c, where c represents the total number of features, and the calculation formula is as follows: , where c is the total number of features, is the value of the feature in the i-th column of the weight vector W, and the calculation formula is as follows: ， In the formula is the value of the node of the tth tree, T is the total number of decision trees, and the calculation formula is as follows: In the formula and is the Gini coefficient of all child nodes marked as 0 and child nodes marked as 1; It represents the Gini coefficient, and the calculation formula is as follows: Where N represents D x The total number of Respectively represent the number of samples marked as "0" and "1" in the node branch; (7) Using the weighted scanning method to x The corresponding feature vector C of the sample in is processed, and the execution steps are as follows: the feature vectors are processed respectively through a sliding window of length m and weight vector Slide, and the window sliding step size is set to p during the sliding process. After sliding, the feature matrix f is obtained respectively (1,n-m+1,m) and the weight matrix w (1,n-m+1,m) ; The feature matrix f (1 ,n-m+1,m) With the weight matrix w (1,n-m+1,m) Multiply to get the pre-input matrix tf (1,n-m+1,m) , where n represents the number of base codes 0 or 1, and m is the length of the sliding window, ranging from 10 to 1000; (8) tf (1,n-m+1,m) Input into the full random forest and random forest, through the above operations, the feature matrix fa is obtained respectively (n-m+1,2) and FB (n-m+1,2) ; Finally, the feature matrix fa captured by two random forests (n-m+1,2) and FB (n-m+1,2) Concatenate into a one-dimensional feature vector F (2*2(n-m+1)) , the feature vector is input into the cascade forest for training, and a prediction model for DNA binding site classification is obtained; (9) The training set D x and the test set D y The models are input into XGBoost and Random Forest algorithms for training respectively, and the best performing model is selected as the final prediction model M according to the evaluation indicators; Finally, the test set is input into model M for prediction.

2. The prediction method based on combined feature coding and DNA binding sites according to claim 1, characterized in that: In step (1), k is a natural number between 1000 and 50000.

3. The prediction method based on combined feature coding and DNA binding sites according to claim 1, characterized in that: The recombination rate threshold of step (2) is 0.46, the number of positive samples screened out is 1762, and the number of negative samples is 11117.

4. The prediction method based on combined feature coding and DNA binding sites according to claim 1, characterized in that: In the combined feature encoding method in step (3), the feature length Length is given by the formula The result is: b represents the length of the combined bases, and its value is 3. The bases A, T, C, and G can be randomly selected to form a sequence set with a single length of 3 bp. Each element in the set C is set as a feature column, for a total of 64 feature columns, and its elements are the feature names of the feature columns. The feature vector C2 is calculated as follows: first, 64 feature columns are constructed with 3 bases as a group, and then a window with a length of 3 is used to slide the entire DNA sequence with a step size of 1, and the feature columns corresponding to the 3 base groups in the window are marked as 1. The feature vector C2 obtained by combining feature encoding has a length of 64.

5. The prediction method based on combined feature coding and DNA binding sites according to claim 1, characterized in that: The number of samples in the training set of step (4) accounts for 80% of the total number of samples in the data set.

6. The prediction method based on combined feature coding and DNA binding sites according to claim 1, characterized in that: In step (6), the total number of decision trees is 462, and the maximum depth of each decision tree is 11; c represents the total number of features.

7. The prediction method based on combined feature coding and DNA binding sites according to claim 1, characterized in that: In step (7), the value of m is 50 and the value of p is 1.

Images