Method for filling missing genotypes based on autoencoder sample matching

Abstract

The application provides a missing genotype filling method based on an automatic encoder sample matching, realizes low-cost and accurate filling of missing genotypes, and can provide more accurate genetic data support for various genetic analysis work. In the application, the genotype information value of each sample at each position in the target data file is converted and finally encoded into a one-hot encoding, a training set and a test set are divided, and then a convolutional denoising autoencoder model is constructed. The application uses an automatic preprocessing strategy to segment the sample data set participating in filling, reduces the device memory occupation, so that the user can successfully perform high-precision genotype filling using a low-cost device. The application has high filling precision, simple and reliable model structure, high training efficiency, and has a wide application prospect in the field of genetic sequence analysis, and can be used for subsequent biological whole genome association analysis and whole genome selection work.

Description

Technical Field

[0001] This invention belongs to the fields of computer and bioinformatics technology, specifically relating to a method for filling missing genotypes based on autoencoder sample matching. Background Technology

[0002] Whole-genome sequencing, utilizing high-throughput sequencing technology, can sequence the entire genome of an organism, revealing information such as sequence or structural variations. It is a crucial prerequisite for biological genetic analysis, providing data support for genome-wide association studies (GWAS) and genome-wide selection (GS) research. However, in modern genome sequencing, technological limitations and cost constraints lead to significant information loss in the detection of single nucleotide polymorphisms (SNPs), resulting in genotype deletions. In gene-related research, genotype deletions prevent researchers from obtaining accurate information about genes at that location, thus affecting the integrity and accuracy of downstream analyses (including genetic mapping, QTL mapping, and GWAS analysis). Therefore, genotype deletions pose an obstacle to GWAS, GS, and many other genetic analyses.

[0003] Genotype imputation is a crucial step in human and plant / animal genome sequence analysis, with results applicable to genetic research, including genome-wide association studies (GWAS) and genome-wide selection. Currently, solutions to genotype missing data fall into two main categories: one is imputation methods with a reference genome, which learn the distribution characteristics of data in the reference genome sequence and reconstruct the missing loci data from the actual gene sequencing process based on the dynamic linkage disequilibrium characteristics of genetics. Existing mainstream imputation methods require a reference panel and use the sample's haplotype information, such as the currently popular Beagle method, which estimates the missing genotype values ​​based on the reference template data of the genotypes to be imputed. However, because it is a linear imputation method based on Hidden Markov Models, it struggles to capture the non-linear relationships between genotypes at different loci, limiting its imputation accuracy. The other category is imputation methods without a reference genome, such as KNN, which fit parameters to a large number of gene sequences using a model to learn the overall genotype distribution characteristics under large sample sizes and impute missing genotypes. However, their imputation accuracy is not ideal. Deep learning-based methods can effectively learn the characteristics between data, thus solving the problem of missing data. Some researchers have used deep learning models for genotype imputation, but in practical applications, the accuracy of conventional deep learning models remains unsatisfactory for cross-population genotype data. Furthermore, due to the randomly scattered missing genotypes in real genotype files, ordinary deep learning methods struggle to capture the features of these random missing genotypes, resulting in low accuracy for imputation. Additionally, because genotype data files are typically large, conventional deep learning methods consume significant amounts of device memory, making it difficult for ordinary devices to handle the load, thus forcing the use of expensive equipment for training and prediction.

[0004] When studying the genes of biological populations, high-density sequencing of the population is often required. However, the cost of high-density sequencing remains high. Therefore, accurate genotyping of low-cost sequencing data obtained through simplified genome sequencing techniques is of significant practical importance. One effective technique is to use convolutional networks (CNNs), which can utilize convolutional kernels to learn the underlying structure and relationships in genotypic data. Convolutional kernels can learn various local patterns through filtering. Different regions of a biological genome sequence have different linkage characteristics between markers, and the closer the markers are, the stronger the linkage. If only a general denoising autoencoder is used to encode and model a large number of markers on chromosomes, problems such as the curse of dimensionality and overfitting will occur, making it impossible to effectively utilize the close linkage relationships between neighboring markers. Summary of the Invention

[0005] To address the aforementioned problems, this invention provides a missing genotype imputation method based on autoencoder sample matching, achieving low-cost and accurate imputation of missing genotypes and providing more accurate genetic data support for various genetic analyses. This invention converts and encodes the genotype information values ​​of samples at each location in the target data file into one-hot codes, divides the dataset into training and test sets, and then constructs a convolutional denoising autoencoder model. This invention uses an automatic preprocessing strategy to segment the sample dataset involved in imputation, reducing device memory usage and enabling users to perform high-precision genotype imputation smoothly using low-cost devices.

[0006] To achieve the above objectives, the technical solution of the present invention is as follows:

[0007] The missing genotype imputation method based on autoencoder sample matching includes the following steps:

[0008] Step 1: Select genomic sequence data of several biological samples from the gene bank, preprocess the selected genomic data, and obtain the original target genomic sequence data;

[0009] Step 2: Further process the raw target genome sequence data obtained in Step 1 to obtain a numerical matrix;

[0010] Step 3: Divide the numerical matrix obtained in Step 2 into a training set and a test set according to the samples to obtain the training set for model training and the test set for verifying accuracy evaluation and practical application effect.

[0011] Step 4: Hide the test set samples with scattered data of different proportions; if reconstruction is to be performed, the position index to be reconstructed needs to be provided, and genotype hidden values ​​will be added to all samples at the corresponding position index in the test set; after the test set samples are hidden, they are transformed into the final dataset to be predicted after the model is trained.

[0012] Step 5: Based on the test set, use the training and test sample self-matching algorithm to obtain a set of missing indices that are similar to the missing features of the test set;

[0013] Step 6: Construct the GCDA model, set the model hyperparameters, and select the initial data segment size and neuron size;

[0014] Step 7: The training set is segmented according to the segmentation intervals set in Step 6. The segmented data matrix is ​​then input into the input layer of the model. A convolutional autoencoder structure with encoder and decoder is used. The training set is masked based on the missing index obtained in Step 5 to learn the training set features. The mapping between the input space and the feature space is solved to minimize the reconstruction error of the input features, thus completing the feature learning of the input information and capturing the linkage relationship of the genotype data. Finally, the learned model is used to predict the test set data of the corresponding intervals, and the prediction results are output and saved to the collection.

[0015] Step 8: For the output data set obtained in Step 7, the filling result data of each group of predicted outputs are spliced ​​together in the segmented order to obtain the total genotype prediction result data matrix; From this genotype matrix, the predicted value and actual value of the hidden genotype corresponding to the test set are extracted, the prediction result is compared with the actual data, and the accuracy of filling or reconstruction is calculated.

[0016] Step 9: Repeat steps 7 and 8. Based on the calculated accuracy, optimize the hyperparameters and segmentation intervals of the autoencoder model to find the optimal model parameters.

[0017] Step 10: After preprocessing the dataset of missing real target genes to be predicted, input it into the model optimized in Step 9 for genotype prediction, output the predicted genotype data, complete the filling, and output the genotype file with the same format as the target file after filling.

[0018] Furthermore, step 2 supports input of multiple data file formats; for VCF format files, the data in the original VCF is extracted to obtain the GT genotype data in the corresponding sample class of the target chromosome, the genotype data is numerically encoded, and the data is integrated into a numerical matrix; for CSV format gene data files, the data in the original CSV file is extracted first using the CSV toolkit, the data is integrated into a matrix, and then the obtained genotype string matrix is ​​encoded and converted into a numerical matrix.

[0019] Furthermore, step 2 specifically includes the following sub-steps:

[0020] Step 2-1, Preprocessing: Automatically determine the specific chromosome information and locus information of the target sample based on the information of the target sample and training samples;

[0021] Step 2-2: Read the sample, and read the chromosomal gene information of each locus one by one in a loop;

[0022] Steps 2-3 involve numerically encoding alleles based on the similarities and differences between alleles and the reference genome, converting them into four numerical forms; and extracting genotype information (GT) data from the target sample as genotype sample data.

[0023] Steps 2-4: Organize the genotype sample data into a genotype matrix and convert it into tensor form;

[0024] Steps 2-5 involve using one-hot encoding to sparsify the tensor data obtained in Steps 2-4, and organizing all gene fragments into a numerical data set.

[0025] Furthermore, in steps 2-3, the specific encoding situations are as follows: if the locus genotype is a combination of two alleles 1, the corresponding encoding is "1"; if the locus genotype is a combination of allele 1 and allele 2, the corresponding encoding is "2"; if the locus genotype is a combination of two alleles 2, the corresponding encoding is "3"; if the locus genotype is empty, the corresponding encoding is "0".

[0026] Furthermore, step 4 specifically includes the following sub-steps:

[0027] Step 4-1: Traverse the target biological sample information to obtain the sample genotype information, and store the sample genotype locus information into a numerical matrix;

[0028] Step 4-2: Read the numerical matrix from Step 4-1 to obtain the total number of sites in the sample, and call the random() function to randomly select a certain number of sites from the sample data according to a certain proportion.

[0029] Step 4-3: Modify the randomly selected locus information to null values ​​to construct scattered and orderly genotype deletion loci. After the construction is completed, check the specific deletion ratio.

[0030] Step 4-4: Repeat the above steps several times to generate several sets of data matrices with different missing site locations.

[0031] Steps 4-5: Using Python's data processing techniques, attach VCF file information to different data matrix information to create VCF format missing genotype files.

[0032] Furthermore, step 5 specifically includes the following sub-steps:

[0033] Step 5-1: After preprocessing, the training sample set data matrix and the target sample data matrix are obtained;

[0034] Step 5-2: Calculate the average value of each position in the training set, and use the average value to initially fill the average value of each position in the target sample to obtain a complete pre-filled sample data matrix.

[0035] Step 5-3: Calculate the pairwise Euclidean distance between the test set samples and the training set samples after the mean pre-filling in Step 5-2, and save the calculated Euclidean distance values ​​as an Euclidean distance matrix.

[0036] Step 5-4: Calculate the row mean of the Euclidean distance matrix in Step 5-3 and store it in an array. This mean array is used to evaluate the overall similarity between each test set sample and the training set.

[0037] Step 5-5: Sort the mean array obtained in Step 5-4 in descending order, and retain the original data indices;

[0038] Steps 5-6: Based on the descending index of the mean array, find the test set sample corresponding to the index;

[0039] Steps 5-7: For the found test set samples, according to the Euclidean distance of the corresponding samples in Step 5-3, bind a training set sample with the smallest distance. The selected sample will not participate in the selection again.

[0040] Steps 5-8: Repeat steps 5-5 to 5-7 until all test sets are bound to a training set sample. If there are still training set samples that are not bound to any test set sample, start the binding process again according to the mean array index order in step 5-5 until there are no more training set samples that are not bound to test set samples. At this time, each training set sample is bound to one and only one test set sample.

[0041] Steps 5-9: Based on the binding results above, construct a genotype deletion that is consistent with the test set sample bound to each training set sample, and obtain a training set with specific deletions constructed according to the matching binding algorithm.

[0042] Steps 5-10: Obtain the index of the missing position based on the newly constructed training set with specific missing parts.

[0043] Furthermore, step 6 specifically includes the following sub-steps:

[0044] Step 6-1: The initial design of the GCDA model is a seven-layer convolutional architecture. The segment size for segmented training is set to 1000, and the filter size is 25.

[0045] Step 6-2: Determine the autoencoder mechanism, process the training set and test set into numerical matrices, and use their samples as feature vectors that can be input into the model for model training;

[0046] Step 6-3: Construct the encoder and decoder. The encoder is responsible for mapping the input sample data to the hidden layer representation. The input data is processed by the encoder to obtain a low-dimensional representation. The encoder parameters are learned by minimizing the reconstruction error, so that the hidden layer can capture the main features of the input data. The decoder consists of an upsampling layer and a convolutional layer. The upsampling layer is used to transform the feature data of the previous layer into the target size and map the hidden layer representation back to the original input space. The hidden layer representation is then mapped back to the original input space by the decoder and compared with the original input data. The convolutional layer obtains new feature data by fusing the feature data of the same layer and continues to process it in the next layer. After multiple layers of the same operation, the feature data compressed by the encoder can be successfully reconstructed into the target data. The reconstruction of the original data is completed through the above steps.

[0047] Step 6-4: Use a multi-class loss function to simultaneously evaluate the loss of the missing regions and the loss of the non-missing regions after reconstruction.

[0048] Furthermore, the formula and parameters of the autoencoder in step 6-2 specifically include:

[0049]

[0050] In the input vector I, M (i) is the output of the i-th tag; T is the convolution filter, and k represents the maximum limit of the convolution filter; the convolution operation in the equation is performed at each position of the input vector I, that is, this operation is performed for each genetic tag;

[0051] The encoding process of an encoder is the process of mapping raw data x from the input layer to the hidden layer to h through a mapping function:

[0052] h = f θ (x)=Φ1(W1x+b1)

[0053] Where θ = {W, b}, W is the weight matrix of the m×n encoder, b is the bias vector of the encoder, W1 is the weight matrix of the m×n encoder, b1 is the bias vector of the encoder, Φ1 is the activation function, and h is the latent variable obtained by the mapping function.

[0054] The decoder's decoding process involves restoring the latent variable h to a dimension similar to the input dimension using equations. The process:

[0055]

[0056] Where W2 is an n×m weight matrix, b2 is the bias vector of the decoder, and Φ2 is an activation function similar to Φ1. It is a matrix with the same format as the input x, obtained by the decoder function.

[0057] Furthermore, step 7 specifically includes the following sub-steps:

[0058] Step 7-1: Segment the training and test data, setting the initial segment size to 1000; set the length of each data set to 10000 points, that is, each data set will be divided into 10 segments; use the same process for each segment, train and fill segments separately, and each segment of data will be input into the model, into the encoder of the autoencoder, and extract features into the hidden layer.

[0059] Step 7-2: The data is upsampled by the decoder to transform the feature data of the previous layer into the target size, and the representation of the hidden layer is mapped back to the original input space. The representation of the hidden layer is mapped back to the original input space by the decoder and compared with the original input data.

[0060] Step 7-3: The convolutional layer obtains new feature data by fusing the feature data of the same layer and continues to process it in the next layer. After multiple layers of the same operation, the feature data compressed by the encoder is reconstructed into the target data, thus completing the reconstruction of the original data.

[0061] Furthermore, in step 8, the filling precision is the ratio of the number of correctly filled marks to the total number of filled marks.

[0062] The beneficial effects of this invention are as follows:

[0063] 1. This invention develops a missing genotype imputation method based on autoencoder sample self-matching, achieving effective learning of linkage features of each segment of biological gene markers. Compared with traditional methods, it can more accurately capture data features in gene data, significantly improving the imputation accuracy of gene markers, especially when the markers have a high missing rate. A new sample self-matching algorithm is invented to control missing masking. The GCDA model using the sample self-matching algorithm has significantly higher accuracy than the CDA model without this method. The missing genotype imputation method based on autoencoder sample self-matching not only has high imputation accuracy but also a simple and reliable model structure and high training efficiency. It has broad application prospects in the field of gene sequence analysis and can be used for subsequent genome-wide association analysis and genome-wide selection.

[0064] 2. This invention utilizes a training and testing sample self-matching algorithm based on Euclidean distance. By analyzing the correlation between cross-population training and testing samples and considering the missing structure of the testing samples, it masks the training samples, enabling the training process to more effectively capture data features between samples and improve the accuracy of model predictions. The Euclidean distance-based sample self-matching algorithm used in this invention solves the problem that ordinary deep learning algorithms cannot effectively obtain missing correlations between cross-population training and testing samples, significantly improving imputation accuracy. Furthermore, this invention converts the genotype information values ​​of samples at each position in the target data file into one-hot encodings, divides the training and testing sets, and designs an automatic preprocessing scheme to segment the sample datasets involved in imputation, reducing device memory usage and allowing users to perform high-precision genotype imputation smoothly using low-cost devices.

[0065] 3. The imputation validation was performed using maize and rice datasets. The results show that the segmentation strategy used in this invention can successfully complete the imputation task with low memory usage. It also achieves excellent accuracy for both random scattered missing and random neat missing genotype cases in the genotype file, which is higher than mainstream missing genotype imputation algorithms such as Beagle and KNN. At the same time, the model has fast prediction speed, high training efficiency and simple structure, and has broad application prospects in the field of gene sequence analysis. Attached Figure Description

[0066] Figure 1 This is a flowchart of the missing genotype imputation method based on autoencoder sample matching.

[0067] Figure 2 This is a schematic diagram of the autoencoder structure for a missing genotype imputation method based on autoencoder sample matching.

[0068] Figure 3 The flowchart shows the algorithm for self-matching training and test samples using Euclidean distance for genes.

[0069] Figure 4 This study compares the filling accuracy of three methods under different marker loss rates in rice.

[0070] Figure 5 This study compares the filling accuracy of three methods under different marker missing rates in maize. Detailed Implementation

[0071] The present invention will be further illustrated below with reference to the accompanying drawings and specific embodiments. It should be understood that the following specific embodiments are for illustrative purposes only and are not intended to limit the scope of the invention.

[0072] This invention provides a method for imputing missing genotypes based on autoencoder sample matching, the overall process of which is as follows: Figure 1 As shown, the specific steps include:

[0073] Step 1: Select genomic sequence data of several biological samples from the gene bank, preprocess the selected genomic data to obtain raw target genomic sequence data that can be further processed.

[0074] Step 2 further processes the raw data described in Step 1, supporting input of multiple data file formats. For VCF format files, data is extracted from the original VCF file to obtain the GT (genotype) genotype data of the corresponding sample class in the target chromosome. The genotype data is then numerically encoded and integrated into a numerical matrix. Additionally, CSV format gene data files are also supported. First, data is extracted from the original CSV file using a CSV toolkit and integrated into a matrix. Then, the obtained genotype string matrix is ​​encoded and converted into a numerical matrix (rows represent samples, columns represent loci). Specifically, this includes the following sub-steps:

[0075] Step 2-1, Preprocessing. The system automatically determines the specific chromosome and locus information for the selected target sample based on the information from the target sample and training samples.

[0076] Step 2-2: Read the sample, and read the chromosomal gene information of each site one by one in a loop.

[0077] Steps 2-3 convert alleles 1 and 2, defined based on their similarities and differences with the reference genome, from their original genotype representations into four numerical forms: 0 / 1 / 2 / 3, thus completing the numerical encoding of the alleles. Genotype information (GT) data is extracted from the target sample as the genotype sample data.

[0078] Specifically, the coding is as follows: A locus genotype consisting of two alleles 1 is encoded as "1". A locus genotype consisting of both allele 1 and allele 2 is encoded as "2". A locus genotype consisting of two alleles 2 is encoded as "3". A locus genotype that is empty is encoded as "0".

[0079] Steps 2-4: Organize the genotype sample data into a genotype matrix and convert it into tensor form.

[0080] Steps 2-5 involve using one-hot encoding to sparsify the tensor data obtained in Steps 2-4, and organizing all gene fragments into a numerical data set.

[0081] Step 3: Divide the numerical matrix obtained in Step 2 into a training set and a test set according to the samples, to obtain a training set for model training and a test set for verifying accuracy and evaluating the effect of practical application.

[0082] Step 4 involves hiding the test set samples with scattered data of varying proportions for later filling and validation. If reconstruction is required, the position index (POS) for reconstruction must be provided. Genotype hiding values ​​will be added to all training samples at the corresponding position index in the test set for subsequent reconstruction and validation. After hiding, the test set samples will be transformed into the dataset that the trained model will ultimately predict. This includes the following sub-steps:

[0083] Step 4-1: Traverse the target biological sample information to obtain the sample genotype information, and store the sample genotype locus information into a numerical matrix.

[0084] Step 4-2: Read the numerical matrix from Step 4-1 to obtain the total number of sites in the sample. Call the random() function to randomly select a certain number of sites from the sample data according to a certain proportion.

[0085] Step 4-3: Modify the randomly selected locus information to null values ​​to construct scattered and orderly genotype deletion loci. After the construction is completed, check the specific deletion ratio.

[0086] Step 4-4: Repeat the above steps 10 times to generate 10 sets of data matrices with different missing site locations.

[0087] Steps 4-5: Using Python's data processing techniques, attach VCF file information to different data matrix information to create VCF format missing genotype files.

[0088] Step 5: Based on the test set, use a training and test sample self-matching algorithm to obtain a set of missing indices similar to the missing features of the test set. These indices are used for data hiding of the training set in the autoencoder model, resulting in better feature learning performance. This includes the following sub-steps:

[0089] Step 5-1: After preprocessing, the training sample set data matrix and the target sample data matrix are obtained.

[0090] Step 5-2: Calculate the average value of each position in the training set, and use the average value to initially fill the average value of each position in the target sample to obtain a complete pre-filled sample data matrix.

[0091] Step 5-3: Calculate the pairwise Euclidean distance between the test set samples and the training set samples after the mean pre-filling in Step 5-2. Save the calculated Euclidean distance values ​​as an Euclidean distance matrix to evaluate the similarity between each training set sample and each test set sample.

[0092] Step 5-4: Calculate the row mean of the Euclidean distance matrix in Step 5-3 and store it in an array. This mean array can be used to evaluate the overall similarity between each test set sample and the training set.

[0093] Step 5-5: Sort the mean array obtained in step 5-4 in descending order and retain the original data index.

[0094] Steps 5-6: Find the test set sample with the corresponding index based on the index of the mean array in descending order.

[0095] Steps 5-7 involve binding the found test set samples to a training set sample with the smallest Euclidean distance according to the corresponding samples in Step 5-3. Once a sample is selected, it will not be selected again.

[0096] Steps 5-8 repeat steps 5-5 through 5-7 until all test set samples are bound to a training set sample. If there are still training set samples that are not bound to any test set sample, the binding process is restarted according to the mean array index order in step 5-5 until no training set sample is left unbound to a test set sample. At this point, each training set sample is bound to exactly one test set sample.

[0097] Steps 5-9: Based on the binding results described above, construct a genotype deletion for each training set sample that matches the bound test set sample. This yields a training set with specific deletions constructed using the matching binding algorithm.

[0098] Steps 5-10: Based on the newly constructed training set with missing parts, obtain the indices of the missing locations. These missing indices can be used to mask the corresponding locations in the training set during model training. Since the masked locations are closer to the actual masking situation in the test set, the training set can more effectively capture features during training, making model learning more efficient.

[0099] Step 6: Construct the GCDA model. Besides hyperparameters such as the dropout rate, which are adaptively set based on the features of the data to be filled, select initial model segment size (several consecutive genotype loci constitute one segment), neuron size, and other parameters. This includes the following sub-steps:

[0100] Step 6-1: The initial design of the GCDA model is a seven-layer convolutional architecture. The segment size for segmented training is set to 1000, and the filter size is 25.

[0101] Step 6-2: Determine the autoencoder mechanism by processing the training and test sets into numerical matrices, and using these samples as feature vectors that can be input into the model for training. The formula and parameters of the autoencoder are as follows:

[0102]

[0103] Where θ = {W, b}, W is the weight matrix of the m×n encoder, b is the bias vector of the encoder, and Φ is the encoder activation function, such as a linear or rectified linear unit (ReLU). The hidden representation h is also called the latent representation. The encoder uses the hidden representation h to map it into the reconstructed vector Rn; in the input vector I, M (i) This is the output of the i-th marker. T is the convolution filter, and k represents the maximum size of the convolution filter. The convolution operation in the equation is performed at each position of the input vector I, that is, this operation is performed for each genetic marker. Each convolution layer has a fixed number of convolution kernels, and the size of the convolution kernels determines the number of convolution filters.

[0104] The encoding process of an encoder is the process of mapping raw data x from the input layer to the hidden layer and then mapping it to h through a mapping function:

[0105] h = f θ (x)=Φ1(W1x+b1)

[0106] Where W1 is an m×n weight matrix, b1 is the encoder bias vector, Φ1 is the activation function, and h is the latent variable obtained from the mapping function.

[0107] The decoding process of the decoder involves restoring the latent variable h to a value similar to the input dimension using equations. The process:

[0108]

[0109] Where W2 is an n×m weight matrix, b2 is the bias vector of the decoder, and Φ2 is an activation function similar to Φ1. It is a matrix with the same format as the input x, obtained by the decoder function.

[0110] Step 6-3 involves constructing an encoder and decoder. The encoder maps the input sample data to the hidden layer representation, resulting in a low-dimensional representation of the input data. The encoder parameters are learned by minimizing the reconstruction error, enabling the hidden layer to capture the main features of the input data. The decoder consists of upsampling layers and convolutional layers. The upsampling layer transforms the feature data from the previous layer into the target size, mapping the hidden layer representation back to the original input space. The hidden layer representation is then mapped back to the original input space by the decoder and compared with the original input data. The convolutional layer fuses the feature data from the same layer to obtain new feature data, which is then processed by the next layer. After multiple layers of the same operation, the feature data compressed by the encoder can be successfully reconstructed into the target data. These steps complete the reconstruction of the original data.

[0111] Step 6-4: Use a multi-class loss function to simultaneously evaluate the loss of the missing regions and the loss of the non-missing regions after reconstruction.

[0112] Step 7: Segment the training set according to the intervals set in Step 6. Input the segmented data matrices sequentially into the model's input layer. Employ an encoder-decoder convolutional autoencoder structure. Mask the training set using the missing masking index obtained in Step 5 to learn training set features. Solve the mapping between the input space and the feature space to minimize the reconstruction error of the input features, thus completing feature learning of the input information, i.e., capturing the linkage relationships of genotype data. Finally, use the learned model to predict the test set data for the corresponding intervals and output the prediction results, saving them to a collection. Specifically, this includes the following sub-steps:

[0113] Step 7-1: Segment the training and test sets, setting the initial segment size to 1000. Each data set is set to a length of 10000 points, meaning each set will be divided into 10 segments. Each segment will undergo the same training and padding process to improve training efficiency and reduce memory usage. Each segment will be input into the model, specifically into the encoder of the autoencoder, to extract features for the hidden layers.

[0114] Step 7-2: The data is upsampled through the decoder to transform the feature data of the previous layer into the target size, and the representation of the hidden layer is mapped back to the original input space. The representation of the hidden layer is mapped back to the original input space through the decoder and compared with the original input data.

[0115] Step 7-3: The convolutional layer obtains new feature data by fusing the feature data of the same layer and then passes it to the next layer for processing. After multiple layers of the same operation, the feature data compressed by the encoder can be successfully reconstructed into the target data. The reconstruction of the original data is completed through the above steps.

[0116] Step 8: For the output dataset obtained in Step 7, concatenate the imputation results of each group of predicted outputs in segmented order to obtain the total genotype prediction result data matrix. From this genotype matrix, extract the predicted and actual values ​​of the genotypes corresponding to the hidden / masked genotypes in the test set, compare the prediction results with the actual data, and calculate the imputation or reconstruction accuracy. The imputation accuracy is the proportion of correctly imputed markers to the total number of imputed markers.

[0117] Step 9: Repeat steps 7 and 8. Based on the calculated accuracy, optimize the hyperparameters and segmentation intervals of the autoencoder model to find the optimal model parameters.

[0118] Step 10: After preprocessing the dataset of missing target genes to be predicted, input it into the optimized model completed in Step 9 for genotype prediction, output the predicted genotype data values, and complete the imputation. Output a genotype file with the same format as the target file after imputation.

[0119] The invention will be further verified and illustrated by the following example:

[0120] (1) First, select a chromosome of rice or corn, select gene segments of the selected chromosome, extract the first 10,000 loci, and identify the loci information to ensure that the obtained sample is our original sample set.

[0121] (2) Iteratively read the chromosome gene information for each locus. Extract the genotype information (GT) from the displayed information entries CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, sample_indexes, and samples as genotype sample data. Convert the genotype sample data into the original genotype matrix. For the genotype matrix, re-encode and define different genotype information into easily trainable coding data. For example, a locus genotype consisting of two alleles 1 is encoded as "1". A locus genotype consisting of both allele 1 and allele 2 is encoded as "2". A locus genotype consisting of two alleles 2 is encoded as "3". An empty locus genotype is encoded as "0".

[0122] (3) The numerical matrix data processed above is sparsified using one-hot encoding, and all gene fragments are organized into a numerical data set. The numerical data set is then transformed by modifying loci and individuals, and the encoding is re-customized for different gene data, converting it into a matrix format to obtain the organized sample set.

[0123] (4) Traverse the test set samples to obtain sample genotype information and store the sample genotype locus information as a numerical matrix. Obtain the total number of loci in the samples, call the random generation function random(), and randomly select a certain number of loci in the sample data at three ratios of 10%, 50%, and 90%, respectively, to mark the low, medium, and high deletion rates. Randomly generate hidden loci corresponding to the three different hiding ratios and store the hiding order of the hidden loci to ensure that the hidden positions of the missing loci processed by the three prediction models are consistent.

[0124] (5) Based on the hidden site sequence, the training set samples and test set samples are processed to hide gene sites. The randomly selected site information is modified to a missing value, corresponding to the code "0", in order to construct scattered and orderly genotype missing sites. After construction, the specific missing ratio is checked. The above operation is repeated ten times to obtain ten sets of test set data.

[0125] (6) Using a training and testing sample self-matching algorithm, calculate the average value of each position in the training set. Use the average value to initially fill in the missing values ​​at each position in the test sample, obtaining a complete pre-filled sample value matrix. Calculate the pairwise Euclidean distance between the pre-filled test set and training set samples, and save the Euclidean distance values ​​as an Euclidean distance matrix to evaluate the similarity between each training set sample and each test set sample. Next, calculate the row average of the Euclidean distance matrix and store it in an array. This average array can be used to evaluate the overall similarity between each test set sample and the training set. Sort the average array in descending order, retaining the original data indices. Based on the descending index of the average array, find the test set sample corresponding to the index. For the found test set sample, bind it to the training set sample with the smallest distance according to the Euclidean distance of the corresponding sample. The selected sample will not participate in the selection again. Repeat the above process until each test set sample is bound to at least one training set sample. If there are still training set samples that have not been bound to any test set samples, the binding process is restarted according to the mean array index order until each training set sample has been bound once. At this point, each training set sample is bound to exactly one test set sample. Based on the binding results, a genotype deletion consistent with the bound test set sample is constructed for each training set sample. A training set with specific deletions is obtained by constructing the matching binding algorithm. Based on the newly constructed training set with deletions, the index of the deletion position is obtained. The deletion index can be used to mask the corresponding positions in the training set during model training. Since the masked positions are closer to the actual masking situation in the test set, the training set can more effectively capture features during training, making model learning more effective.

[0126] (7) The GCDA model designed in this invention has a seven-layer convolutional architecture, with a segment size of 1000 and a filter size of 25. An autoencoder mechanism is determined, and the training and test sets are processed into numerical matrices, with their samples used as feature vectors input into the model for training. An encoder and decoder are constructed. The encoder maps the input sample data to the hidden layer representation, and the input data is represented in low dimension by the encoder. The encoder parameters are learned by minimizing the reconstruction error, enabling the hidden layer to capture the main features of the input data. The decoder consists of an upsampling layer and a convolutional layer. The upsampling layer converts the feature data from the previous layer into the target size, maps the hidden layer representation back to the original input space, and the hidden layer representation is mapped back to the original input space by the decoder and compared with the original input data. The convolutional layer obtains new feature data by fusing feature data at the same level and continues to process it in the next layer. After multiple layers of the same operation, the feature data compressed by the encoder can be successfully reconstructed into target data, completing the reconstruction of the original data through the above steps. A multi-class loss function is used to evaluate the loss in the missing region and the loss in the non-missing region after reconstruction.

[0127] The formula and parameters for GCDA specifically include:

[0128]

[0129] Where θ = {W, b}, W is the weight matrix of the m×n encoder, b is the bias vector of the encoder, and Φ is the encoder activation function, such as a linear or rectified linear unit (ReLU). The hidden representation h is also called the latent representation. The encoder uses the hidden representation h to map it into the reconstructed vector Rn; in the input vector I, M (i) This is the output of the i-th marker. T is the convolution filter, and k represents the maximum size of the convolution filter. The convolution operation in the equation is performed at each position of the input vector I, that is, this operation is performed for each genetic marker. Each convolution layer has a fixed number of convolution kernels, and the size of the convolution kernels determines the number of convolution filters.

[0130] The encoding process of an encoder is the process of mapping raw data x from the input layer to the hidden layer and then mapping it to h through a mapping function:

[0131] h = f θ (x)=Φ1(W1x+b1)

[0132] Where W1 is an m×n weight matrix, b1 is the encoder bias vector, Φ1 is the activation function, and h is the latent variable obtained from the mapping function.

[0133] The decoding process of the decoder involves restoring the latent variable h to a value similar to the input dimension using equation (3-3). The process:

[0134]

[0135] Where W2 is an n×m weight matrix, b2 is the bias vector of the decoder, and Φ2 is an activation function similar to Φ1. A matrix obtained by the decoder function with the same format as the input x.

[0136] (8) For each output dataset obtained in (7), insert the predicted results into the missing positions in the test set in order, and repeat the above experimental operation continuously to finally determine the optimal autoencoder model of the convolutional model. By comparing the accuracy of the convolutional denoising autoencoder training model with the accuracy of KNN and Beagle algorithms under the same conditions, the prediction effect of the three is finally obtained. At the same time, in the case of scattered missing data, the accuracy of RCDA using the sample self-matching algorithm in (6) of this invention is compared with CDA without using the algorithm. Figure 4 and Figure 5 As shown.

[0137] The above examples demonstrate that this invention features a reliable and highly accurate model, and can be trained under different loci and hidden ratios. By sparsifying and processing the sample dataset, the accuracy of label imputation is improved. Compared to traditional gene imputation methods such as KNN and Beagle, it exhibits a significant accuracy advantage regardless of the missing genotype rate. This invention employs a segmented processing approach, which can address local linkage relationships that are difficult for other algorithms to learn, enabling better feature extraction and thus achieving accurate imputation of missing genotypes.

[0138] It should be noted that the above content merely illustrates the technical concept of the present invention and should not be construed as limiting the scope of protection of the present invention. For those skilled in the art, various improvements and modifications can be made without departing from the principle of the present invention, and all such improvements and modifications fall within the scope of protection of the claims of the present invention.

Claims

1. A method for missing genotype imputation based on autoencoder sample matching, characterized in that, The method comprises the following steps: Step 1, selecting a plurality of biological sample genomic sequence data from a gene library, preprocessing the selected genomic data to obtain original target genomic sequence data; Step 2, further processing the original target genomic sequence data obtained in step 1 to obtain a numerical matrix; Step 3, dividing the numerical matrix obtained in step 2 into a training set and a test set according to samples, obtaining a training set for model training and a test set for verifying accuracy and evaluating actual application effect; Step 4, hiding the test set samples by different proportions of scattered data; if reconstruction is performed, the position index of the desired reconstruction is provided, and genotype hiding values will be added to all training samples at the corresponding position index of the test set; After the test set samples are hidden, they are converted into the data set to be finally predicted after the model training is completed; Step 5, using the training and test sample self-matching algorithm according to the test set to obtain a missing index set similar to the missing features of the test set; Step 6, constructing a GCDA model, setting model hyperparameters, selecting initial data segmentation size and neuron size; Step 7, segmenting the training set according to the segmentation interval set in step 6, sequentially inputting the segmented data matrix through the input layer of the model, using the convolutional autoencoder structure of the encoder-decoder, covering the training set according to the missing index obtained in step 5 to learn the features of the training set, solving the mapping of the input space and the feature space to minimize the reconstruction error of the input features, completing the feature learning of the input information, capturing the linkage relationship of the genotype data, and finally using the learned model to predict the test set data of the corresponding interval and output the prediction result to the set; Step 8, for the output data set obtained in step 7, the filling result data of each prediction output is spliced to obtain a total genotype prediction result data matrix according to the segmentation order; from the genotype matrix, the predicted value and the actual value of the genotype hidden and covered by the test set are extracted, the prediction result is compared with the actual data, and the accuracy of filling or reconstruction is calculated; Step 9, repeating steps 7 and 8, optimizing the autoencoder model hyperparameters and the size of the segmentation interval according to the calculated accuracy, and finding the best model parameters; Step 10, inputting the preprocessed missing real target gene data set into the model optimized in step 9 for genotype prediction, outputting the genotype data prediction value, completing the filling, and outputting the genotype file in the same format as the target file after filling.

2. The method of missing genotype imputation based on autoencoder sample matching according to claim 1, wherein, The step 2 supports input of multiple data file formats; for vcf format files, the data in the original vcf is extracted to obtain GT genotype data in the sample corresponding to the target chromosome, the genotype data is numerically encoded, and the data is integrated into a numerical matrix; for csv format gene data files, the data in the original csv file is extracted using a csv toolkit, and the data is integrated into a matrix, and then the genotype string matrix obtained is encoded and converted into a numerical matrix.

3. The method of missing genotype imputation based on autoencoder sample matching according to claim 1, wherein, The step 2 specifically comprises the following substeps: Step 2-1, preprocessing, automatically determine the specific chromosome information and site information of the selected target sample according to the target sample and training sample information; Step 2-2, reading the sample, cyclically reading the chromosome gene information of each site one by one; Step 2-3, according to the difference between the allele bases and the reference genome, the numerical encoding of the allele is completed, and it is converted into four numerical forms; the genotype information GT data in the target sample is extracted as the genotype sample data; Step 2-4, the genotype sample data is sorted into a genotype matrix and converted into a tensor form; Step 2-5, the tensor data obtained in step 2-4 is sparse using one-hot encoding, and all gene fragments are sorted as a numerical data set.

4. The method of missing genotype imputation based on autoencoder sample matching according to claim 3, wherein, In step 2-3, the encoding specific case is: the site genotype is composed of two alleles 1, which corresponds to the code "1"; the site genotype is composed of allele 1 and allele 2, which corresponds to the code "2"; the site genotype is composed of two alleles 2, which corresponds to the code "3"; the site genotype is empty, which corresponds to the code "0".

5. The method of missing genotype imputation based on autoencoder sample matching of claim 1, wherein, The step 4 specifically includes the following sub-steps: Step 4-1, traversing the target biological sample information to obtain sample genotype information, and storing the sample genotype site information in a numerical matrix; Step 4-2, reading the numerical matrix of step 4-1 to obtain the total number of sample sites, calling the random generation function random(), and randomly extracting a certain number of sites in the sample data at a certain proportion; Step 4-3, modifying the site information randomly extracted as null value, and constructing scattered and neat genotype missing sites, and checking the specific missing proportion after the construction is completed; Step 4-4, repeat the above steps several times to generate several groups of data matrices with different missing site positions; Step 4-5, using the data processing means of python language, appending the VCF file information to the different data matrix information, and making a missing genotype file in VCF format.

6. The method of missing genotype imputation based on autoencoder sample matching of claim 1, wherein, The step 5 specifically includes the following sub-steps: Step 5-1, after preprocessing, obtain the training sample set data matrix and the target sample data matrix; Step 5-2, calculate the average value of each position of the training set, use the average value to preliminarily fill the average value of each position of the target sample, and obtain the complete pre-filled sample data matrix; Step 5-3, calculate the Euclidean distance between each two of the test set samples and the training set samples after using the average value to pre-fill in step 5-2, and save the calculated Euclidean distance value as a Euclidean distance matrix; Step 5-4, calculate the row average of the Euclidean distance matrix in step 5-3 and store it in an array, which is used to evaluate the overall similarity between each test set sample and the training set; Step 5-5, according to the average array obtained in step 5-4, sort in descending order and keep the original data index; Step 5-6, according to the index of the average array after descending order, find the corresponding test set sample; Step 5-7, for the test set samples found, according to the Euclidean distance of the corresponding sample in step 5-3, bind a distance minimum training set sample to it, and the selected sample no longer participates in the selection; Step 5-8, loop step 5-5 to step 5-7 until all test sets are bound to a training set sample; if there are still remaining training set samples not bound to any test set sample at this time, then according to the mean array index order in step 5-5, start the binding process again, until there are no longer any training set samples not bound to test set samples, at this time, each training set sample has and only has one test set sample bound to it; Step 5-9, according to the above binding results, construct the genotype missing of each training set sample consistent with the bound test set sample, and obtain the training set with specific missing constructed according to the matching binding algorithm; Step 5-10, according to the newly constructed training set with specific missing, obtain the index of the missing position.

7. The method of missing genotype imputation based on autoencoder sample matching according to claim 1, wherein, The step 6 specifically comprises the following sub-steps: Step 6-1, initially design the GCDA model layer number as a seven-layer convolutional architecture, and set the segmentation size of the segmented training as 1000 and the filter size as 25; Step 6-2, determine the auto-encoder mechanism, process the training set and the test set into numerical matrices, and input the samples in the form of feature vectors into the model for training; Step 6-3, construct the encoder and the decoder, the encoder is responsible for mapping the input sample data to the hidden layer representation, the input data is obtained by the encoder to get a low-dimensional representation, the parameters of the encoder are learned by minimizing the reconstruction error, so that the hidden layer can capture the main features of the input data; The decoder is composed of an upsampling layer and a convolutional layer, the upsampling layer is used to convert the feature data of the previous layer to the target size, and the representation of the hidden layer is mapped back to the original input space, and the representation of the hidden layer is mapped back to the original input space through the decoder and compared with the original input data; the convolutional layer obtains new feature data by fusing the feature data of the same layer and continues to be processed by the next layer, and after the same operation of multiple layers, the compressed feature data of the encoder can be successfully reconstructed into target data, and the reconstruction of the original data is completed through the above steps; Step 6-4, a multi-classification loss function is used to evaluate the loss of the missing area after reconstruction and the loss of the non-missing area.

8. The method of missing genotype imputation based on autoencoder sample matching according to claim 7, wherein, The formula and parameters of the auto-encoder in step 6-2 specifically include: In the input vector I, M (i) is the output of the i marker; T is the convolution filter, k represents the maximum limit of the convolution filter; the convolution operation in the equation is performed for each position of the input vector I, that is, this operation is performed for each genetic marker; The encoding process of the encoder is the process of mapping the original data x from the input layer to the hidden layer through the mapping function to obtain h: h = f θ (x) = Φ1(W1x + b1) Where θ={W,b}, W is the weight matrix of the m*n encoder, b is the bias vector of the encoder, W1 is the weight matrix of m*n size, b1 is the bias vector of the encoder, Φ1 is the activation function, and h is the hidden variable obtained by the mapping function; The decoding process of the decoder is to reduce the latent variable h through equation into the same dimension as the input process: where W2 is an n x m weight matrix, b2 is a bias vector for the decoder, Φ2 is an activation function similar to Φ1, is a matrix of the same format as the input x resulting from the decoder function.

9. The method of missing genotype imputation based on autoencoder sample matching of claim 1, wherein, The step 7 specifically comprises the following sub-steps: Step 7-1, the training set and test set data are processed in segments, and the initial segment size is set to 1000; the length of each data set is set to 10000 sites, that is, each data set will be divided into 10 segments; the same process is used for each segment, and the training and padding are performed in segments; each segment of data is input into the model, input into the encoder of the autoencoder, and the features are extracted into the hidden layer; Step 7-2, the data is converted into the target size through the upsampling layer of the decoder, the hidden layer representation is mapped back to the original input space, and the hidden layer representation is mapped back to the original input space through the decoder and compared with the original input data; Step 7-3, the convolutional layer fuses the feature data of the same layer to obtain new feature data and continues to be processed by the next layer, and after multiple layers of the same operation, the compressed feature data of the encoder is reconstructed into target data, and the reconstruction of the original data is completed.

10. The method of missing genotype imputation based on autoencoder sample matching of claim 1, wherein, The filling accuracy in step 8 is the ratio of the number of correctly filled labels to the total number of filled labels.

Images